7 #include <dune/common/hybridutilities.hh> 8 #include <dune/common/timer.hh> 9 #include <dune/functions/functionspacebases/subspacebasis.hh> 10 #include <dune/grid/common/capabilities.hh> 11 #include <dune/typetree/childextraction.hh> 13 #include <amdis/AdaptInfo.hpp> 14 #include <amdis/BackupRestore.hpp> 15 #include <amdis/Assembler.hpp> 16 #include <amdis/GridFunctionOperator.hpp> 17 #include <amdis/io/FileWriterCreator.hpp> 18 #include <amdis/linearalgebra/SymmetryStructure.hpp> 22 template <
class Traits>
30 if (initFlag.
isSet(CREATE_MESH) ||
31 (!adoptFlag.
isSet(INIT_MESH) &&
32 (initFlag.
isSet(INIT_SYSTEM) || initFlag.
isSet(INIT_FE_SPACE)))) {
37 (adoptFlag.
isSet(INIT_MESH) ||
38 adoptFlag.
isSet(INIT_SYSTEM) ||
39 adoptFlag.
isSet(INIT_FE_SPACE))) {
40 adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
45 warning(
"no grid created");
49 if (initFlag.
isSet(INIT_FE_SPACE) ||
50 (initFlag.
isSet(INIT_SYSTEM) && !adoptFlag.
isSet(INIT_FE_SPACE))) {
55 (adoptFlag.
isSet(INIT_FE_SPACE) || adoptFlag.
isSet(INIT_SYSTEM))) {
56 adoptGlobalBasis(adoptProblem->globalBasis_);
61 warning(
"no globalBasis created\n");
64 if (initFlag.
isSet(INIT_SYSTEM))
65 createMatricesAndVectors();
67 if (adoptProblem && adoptFlag.
isSet(INIT_SYSTEM)) {
68 systemMatrix_ = adoptProblem->systemMatrix_;
69 solution_ = adoptProblem->solution_;
70 rhs_ = adoptProblem->rhs_;
71 estimates_ = adoptProblem->estimates_;
77 if (initFlag.
isSet(INIT_SOLVER))
80 if (adoptProblem && adoptFlag.
isSet(INIT_SOLVER)) {
81 test_exit(!linearSolver_,
"solver already created\n");
82 linearSolver_ = adoptProblem->linearSolver_;
87 warning(
"no solver created\n");
91 if (initFlag.
isSet(INIT_MARKER))
94 if (adoptProblem && adoptFlag.
isSet(INIT_MARKER))
95 marker_ = adoptProblem->marker_;
99 if (initFlag.
isSet(INIT_FILEWRITER))
102 solution_->resizeZero();
106 template <
class Traits>
110 std::string grid_filename = Parameters::get<std::string>(name_ +
"->restore->grid").value();
111 std::string solution_filename = Parameters::get<std::string>(name_ +
"->restore->solution").value();
112 test_exit(filesystem::exists(grid_filename),
"Restore file '{}' not found.", grid_filename);
113 test_exit(filesystem::exists(solution_filename),
"Restore file '{}' not found.", solution_filename);
116 using HostGrid =
typename Grid::HostGrid;
119 if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
120 adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
125 if (initFlag.
isSet(INIT_FE_SPACE) || initFlag.
isSet(INIT_SYSTEM))
129 if (initFlag.
isSet(INIT_SYSTEM))
130 createMatricesAndVectors();
133 if (!linearSolver_ && initFlag.
isSet(INIT_SOLVER))
137 if (initFlag.
isSet(INIT_MARKER))
141 if (initFlag.
isSet(INIT_FILEWRITER))
144 solution_->resize(sizeInfo(*globalBasis_));
145 solution_->restore(solution_filename);
149 template <
class Traits>
157 boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
159 boundaryManager_->setBoundaryIds(creator.
boundaryIds());
161 info(3,
"Create grid:");
162 info(3,
"#elements = {}" , grid_->size(0));
163 info(3,
"#faces/edges = {}", grid_->size(1));
164 info(3,
"#vertices = {}" , grid_->size(dim));
165 info(3,
"overlap-size = {}", grid_->leafGridView().overlapSize(0));
166 info(3,
"ghost-size = {}" , grid_->leafGridView().ghostSize(0));
171 template <
class T,
class GV>
172 using HasCreate = decltype(T::create(std::declval<GV>()));
175 template <
class Traits>
178 createGlobalBasisImpl(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
183 template <
class Traits>
186 assert(
bool(grid_) );
187 static_assert(std::is_same_v<GridView, typename Grid::LeafGridView>,
"");
188 auto basis = Traits::create(name_, grid_->leafGridView());
189 globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis));
193 template <
class Traits>
196 error_exit(
"Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
200 template <
class Traits>
204 template <
class Traits>
207 systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
208 std::string symmetryStr =
"unknown";
210 systemMatrix_->setSymmetryStructure(symmetryStr);
212 solution_ = std::make_shared<SolutionVector>(globalBasis_);
213 rhs_ = std::make_shared<SystemVector>(globalBasis_);
215 auto localView = globalBasis_->localView();
216 for_each_node(localView.tree(), [&,
this](
auto&&,
auto treePath) ->
void 218 std::string i = to_string(treePath);
219 estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
220 for (std::size_t j = 0; j < estimates_[i].size(); j++)
221 estimates_[i][j] = 0.0;
226 template <
class Traits>
229 std::string solverName =
"default";
235 linearSolver_ = solverCreator->createWithString(name_ +
"->solver");
239 template <
class Traits>
243 auto localView = globalBasis_->localView();
244 for_each_node(localView.tree(), [&,
this](
auto&&,
auto treePath) ->
void 246 std::string componentName = name_ +
"->marker[" + to_string(treePath) +
"]";
248 if (!Parameters::get<std::string>(componentName +
"->strategy"))
251 std::string tp = to_string(treePath);
254 assert(
bool(newMarker));
255 this->addMarker(std::move(newMarker));
260 template <
class Traits>
266 auto localView = globalBasis_->localView();
267 for_each_node(localView.tree(), [&](
auto const& ,
auto treePath) ->
void 269 std::string componentName = name_ +
"->output[" + to_string(treePath) +
"]";
270 auto format = Parameters::get<std::vector<std::string>>(componentName +
"->format");
272 if (!format && to_string(treePath).empty()) {
274 componentName = name_ +
"->output";
275 format = Parameters::get<std::vector<std::string>>(componentName +
"->format");
281 for (std::string
const& type : format.value()) {
282 auto writer = creator.
create(type, componentName, treePath);
284 filewriter_.push_back(std::move(writer));
291 template <
class Traits>
292 template <
class Predicate,
class RowTreePath,
class ColTreePath,
class Values>
296 static constexpr
bool isValidPredicate = Concepts::Functor<Predicate, bool(WorldVector)>;
298 "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
300 static constexpr
bool isValidTreePath =
301 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
302 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
303 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
305 if constexpr (isValidPredicate && isValidTreePath) {
306 auto localView = globalBasis_->localView();
307 auto i = makeTreePath(row);
308 auto j = makeTreePath(col);
309 auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, i);
310 auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, j);
314 auto bc = makeDirichletBC<SystemMatrix, SolutionVector, SystemVector>(
315 std::move(rowBasis), std::move(colBasis), {predicate}, valueGridFct);
316 boundaryConditions_[i][j].push_back(
makeUniquePtr(std::move(bc)));
322 template <
class Traits>
323 template <
class RowTreePath,
class ColTreePath,
class Values>
327 static constexpr
bool isValidTreePath =
328 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
329 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
330 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
332 if constexpr (isValidTreePath) {
333 auto localView = globalBasis_->localView();
334 auto i = makeTreePath(row);
335 auto j = makeTreePath(col);
336 auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, i);
337 auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, j);
341 auto bc = makeDirichletBC<SystemMatrix, SolutionVector, SystemVector>(
342 std::move(rowBasis), std::move(colBasis), {*boundaryManager_,
id}, valueGridFct);
343 boundaryConditions_[i][j].push_back(
makeUniquePtr(std::move(bc)));
348 template <
class Traits>
352 auto localView = globalBasis_->localView();
353 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath());
354 auto bc = makePeriodicBC<SystemMatrix, SolutionVector, SystemVector>(
355 std::move(basis), {*boundaryManager_,
id}, {matrix, vector});
356 boundaryConditions_[makeTreePath()][makeTreePath()].push_back(
makeUniquePtr(std::move(bc)));
360 template <
class Traits>
371 linearSolver_->solve(*systemMatrix_, *solution_, *rhs_, solverInfo);
373 if (solverInfo.
info() > 0) {
374 msg(
"solution of discrete system needed {} seconds", t.elapsed());
378 msg(
"Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
385 test_exit(!solverInfo.doBreak() && !solverInfo.
error(),
"Could not solver the linear system!");
389 template <
class Traits>
396 for (
auto& currentMarker : marker_)
397 markFlag |= currentMarker.second->markGrid(adaptInfo);
399 msg(
"markElements needed {} seconds", t.elapsed());
405 template <
class Traits>
410 bool adapted =
false;
412 for (
int i = 0; i < n; ++i) {
414 for (
const auto& element : elements(grid_->leafGridView()))
415 grid_->mark(-1, element);
417 bool adaptedInLoop = grid_->preAdapt();
418 adaptedInLoop |= grid_->adapt();
426 msg(
"globalCoarsen needed {} seconds", t.elapsed());
427 return adapted ? MESH_ADAPTED :
Flag(0);
433 using HasGlobalRefineADHI = decltype(
434 std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
436 template <
class Traits>
441 if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
442 grid_->globalRefine(n, globalBasis_->globalRefineCallback());
444 grid_->globalRefine(n);
446 msg(
"globalRefine needed {} seconds", t.elapsed());
447 return n > 0 ? MESH_ADAPTED :
Flag(0);
451 template <
class Traits>
457 bool adapted = grid_->preAdapt();
458 adapted |= grid_->adapt();
461 msg(
"adaptGrid needed {} seconds", t.elapsed());
462 return adapted ? MESH_ADAPTED :
Flag(0);
466 template <
class Traits>
473 auto localView = globalBasis_->localView();
474 for_each_node(localView.tree(), [&](
auto&&,
auto rowTp) ->
void {
475 for_each_node(localView.tree(), [&](
auto&&,
auto colTp) ->
void {
476 for (
auto bc : boundaryConditions_[rowTp][colTp])
484 systemMatrix_->init();
485 rhs_->init(sizeInfo(*globalBasis_), asmVector);
489 msg(
"{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
491 msg(
"{} local DOFs", rhs_->localSize());
494 for (
auto const& element : elements(gridView(), PartitionSet{})) {
495 localView.bind(element);
498 systemMatrix_->assemble(localView, localView);
500 rhs_->assemble(localView);
506 systemMatrix_->finish();
509 info(2,
" assemble operators needed {} seconds", t2.elapsed());
512 solution_->resize(sizeInfo(*globalBasis_));
513 for_each_node(localView.tree(), [&](
auto&&,
auto rowTp) ->
void {
514 for_each_node(localView.tree(), [&](
auto&&,
auto colTp) ->
void {
516 for (
auto bc : boundaryConditions_[rowTp][colTp])
517 bc->apply(*systemMatrix_, *solution_, *rhs_);
521 info(2,
" assemble boundary conditions needed {} seconds", t2.elapsed());
523 msg(
"fill-in of assembled matrix: {}", systemMatrix_->nnz());
524 msg(
"assemble needed {} seconds", t.elapsed());
528 template <
class Traits>
533 for (
auto writer : filewriter_)
534 writer->write(adaptInfo, force);
535 msg(
"writeFiles needed {} seconds", t.elapsed());
constexpr bool isSet(Flag const &f) const
Checks whether all set bits of f.flags_ are set in flags_ too.
Definition: Flag.hpp:156
double relResidual() const
Returns relResidual_.
Definition: SolverInfo.hpp:50
void addPeriodicBC(BoundaryType id, WorldMatrix const &A, WorldVector const &b)
Definition: ProblemStat.inc.hpp:350
The Flag class encapsulates flags which represents simple information. Used e.g. while mesh traversal...
Definition: Flag.hpp:13
constexpr bool Functor
A Functor is a function F with signature Signature.
Definition: Concepts.hpp:134
void addDirichletBC(Predicate const &predicate, RowTreePath row, ColTreePath col, Values const &values)
Add boundary conditions to the system.
Definition: ProblemStat.inc.hpp:294
void error_exit(std::string const &str, Args &&... args)
print a message and exit
Definition: Output.hpp:142
constexpr bool Predicate
A predicate is a function that returns a boolean.
Definition: Concepts.hpp:142
decltype(auto) makeGridFunction(PreGridFct const &preGridFct, GridView const &gridView)
Generator for Gridfunctions from Expressions (PreGridfunctions)
Definition: GridFunction.hpp:154
A CreatorMap is used to construct objects, which types depends on key words determined at run time...
Definition: CreatorMap.hpp:29
Contains all classes needed for solving linear and non linear equation systems.
Definition: AdaptBase.hpp:6
Definition: BackupRestore.hpp:15
void setCreateMatrixData(bool b)
Sets createMatrixData_.
Definition: SolverInfo.hpp:104
A creator class for dune grids.
Definition: MeshCreator.hpp:52
Definition: ProblemStat.hpp:52
int info() const
Returns info.
Definition: SolverInfo.hpp:38
Creator class for filewriters depending on a given type name.
Definition: FileWriterCreator.hpp:25
void initialize(Flag initFlag, Self *adoptProblem=nullptr, Flag adoptFlag=INIT_NOTHING)
Initialisation of the problem.
Definition: ProblemStat.inc.hpp:23
int error() const
Returns error code in last run of an iterative solver.
Definition: SolverInfo.hpp:32
static int mpiSize()
Return the MPI_Size of the group created by Dune::MPIHelper.
Definition: Environment.hpp:74
void restore(Flag initFlag)
Read the grid and solution from backup files and initialize the problem.
Definition: ProblemStat.inc.hpp:108
void msg(std::string const &str, Args &&... args)
print a message
Definition: Output.hpp:98
static std::shared_ptr< Grid > create(std::string name)
Static create mthod. See create()
Definition: MeshCreator.hpp:71
static std::optional< T > get(std::string const &key)
Get parameter-values from parameter-tree.
Definition: Initfile.hpp:25
Flag markElements(AdaptInfo &adaptInfo) override
Implementation of ProblemStatBase::markElements.
Definition: ProblemStat.inc.hpp:391
CreatorInterfaceName< BaseClass > * named(CreatorInterface< BaseClass > *ptr)
cast a ptr of CreatorInterface to CreatorInterfaceName
Definition: CreatorInterface.hpp:64
Flag adaptGrid(AdaptInfo &adaptInfo) override
Implementation of ProblemStatBase::refineMesh.
Definition: ProblemStat.inc.hpp:453
void writeFiles(AdaptInfo &adaptInfo, bool force=false)
Writes output files. If force=true write even if timestep out of write rhythm.
Definition: ProblemStat.inc.hpp:530
void solve(AdaptInfo &adaptInfo, bool createMatrixData=true, bool storeMatrixData=false) override
Implementation of ProblemStatBase::solve.
Definition: ProblemStat.inc.hpp:362
Flag globalRefine(int n) override
Uniform global refinement by n level.
Definition: ProblemStat.inc.hpp:438
Holds adapt parameters and infos about the problem.
Definition: AdaptInfo.hpp:25
void info(int noInfoLevel, std::string const &str, Args &&... args)
prints a message, if Environment::infoLevel() >= noInfoLevel
Definition: Output.hpp:105
void for_each_node(Tree &&tree, PreFunc &&preFunc, LeafFunc &&leafFunc, PostFunc &&postFunc)
Traverse tree and visit each node.
Definition: Traversal.hpp:81
Definition: SolverInfo.hpp:11
void buildAfterAdapt(AdaptInfo &adaptInfo, Flag flag, bool asmMatrix=true, bool asmVector=true) override
Implementation of ProblemStatBase::buildAfterCoarse.
Definition: ProblemStat.inc.hpp:468
Flag globalCoarsen(int n) override
Uniform global grid coarsening by up to n level.
Definition: ProblemStat.inc.hpp:407
auto makeUniquePtr(Obj &&obj)
Create a unique_ptr by copy/move construction.
Definition: TypeTraits.hpp:107
void test_exit(bool condition, std::string const &str, Args &&... args)
test for condition and in case of failure print message and exit
Definition: Output.hpp:163
void setStoreMatrixData(bool b)
Sets storeMatrixData_.
Definition: SolverInfo.hpp:110
static std::unique_ptr< EstimatorMarker< Grid > > createMarker(std::string const &name, std::string const &component, Estimates const &est, std::shared_ptr< Grid > const &grid)
Creates a scalar marker depending on the strategy set in parameters.
Definition: Marker.inc.hpp:102
std::vector< int > const & boundaryIds() const
Return the filled vector of boundary ids. NOTE: not all creators support reading this.
Definition: MeshCreator.hpp:140
std::unique_ptr< FileWriterInterface > create(std::string type, std::string prefix, Indices... ii) const
Create a new FileWriter of type type
Definition: FileWriterCreator.hpp:47
double absResidual() const
Returns absResidual_.
Definition: SolverInfo.hpp:44