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/GridFunctionOperator.hpp> 16 #include <amdis/functions/EntitySet.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))) {
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))) {
61 warning(
"no globalBasis created\n");
64 if (initFlag.
isSet(INIT_SYSTEM))
65 createMatricesAndVectors();
67 if (adoptProblem && adoptFlag.
isSet(INIT_SYSTEM)) {
70 rhs_ = adoptProblem->
rhs_;
77 if (initFlag.
isSet(INIT_SOLVER))
80 if (adoptProblem && adoptFlag.
isSet(INIT_SOLVER)) {
81 test_exit(!linearSolver_,
"solver already created\n");
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(*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 Traversal::forEachNode(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";
232 linearSolver_ = std::make_shared<LinearSolver>(solverName, name_ +
"->solver");
236 template <
class Traits>
240 auto localView = globalBasis_->localView();
241 Traversal::forEachNode(localView.tree(), [&,
this](
auto&&,
auto treePath) ->
void 243 std::string componentName = name_ +
"->marker[" + to_string(treePath) +
"]";
244 auto strategy = Parameters::get<std::string>(componentName +
"->strategy");
246 if (!strategy && to_string(treePath).empty()) {
248 componentName = name_ +
"->strategy";
249 strategy = Parameters::get<std::string>(componentName +
"->strategy");
255 std::string tp = to_string(treePath);
258 assert(
bool(newMarker));
259 this->addMarker(std::move(newMarker));
264 template <
class Traits>
270 auto localView = globalBasis_->localView();
271 Traversal::forEachNode(localView.tree(), [&](
auto const& ,
auto treePath) ->
void 273 std::string componentName = name_ +
"->output[" + to_string(treePath) +
"]";
274 auto format = Parameters::get<std::vector<std::string>>(componentName +
"->format");
276 if (!format && to_string(treePath).empty()) {
278 componentName = name_ +
"->output";
279 format = Parameters::get<std::vector<std::string>>(componentName +
"->format");
285 for (std::string
const& type : format.value()) {
286 auto writer = creator.
create(type, componentName, treePath);
288 filewriter_.push_back(std::move(writer));
295 template <
class Traits>
296 template <
class Predicate,
class RowTreePath,
class ColTreePath,
class Values>
300 static constexpr
bool isValidPredicate = Concepts::Functor<Predicate, bool(WorldVector)>;
302 "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
304 static constexpr
bool isValidTreePath =
305 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
306 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
307 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
309 if constexpr (isValidPredicate && isValidTreePath) {
313 globalBasis_, makeTreePath(row), makeTreePath(col), {predicate}, valueGridFct});
319 template <
class Traits>
320 template <
class RowTreePath,
class ColTreePath,
class Values>
324 static constexpr
bool isValidTreePath =
325 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
326 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
327 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
329 if constexpr (isValidTreePath) {
333 globalBasis_, makeTreePath(row), makeTreePath(col), {boundaryManager_,
id}, valueGridFct});
338 template <
class Traits>
342 Traversal::forEachLeafNode(globalBasis_->localView().tree(), [&](
auto&&,
auto tp) {
343 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, tp);
344 constraints_.push_back(makePeriodicBC(
345 std::move(basis), {boundaryManager_,
id}, {matrix, vector}));
351 template <
class Traits>
356 Dune::InverseOperatorResult stat;
360 if (createMatrixData)
361 linearSolver_->init(systemMatrix_->impl());
364 linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);
366 if (!storeMatrixData)
367 linearSolver_->finish();
369 info(2,
"solution of discrete system needed {} seconds", t.elapsed());
370 info(1,
"Residual norm: ||b-Ax|| = {}",
371 residuum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
373 if (stat.reduction >= 0.0)
374 info(2,
"Relative residual norm: ||b-Ax||/||b|| = {}", stat.reduction);
376 info(2,
"Relative residual norm: ||b-Ax||/||b|| = {}",
377 relResiduum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
379 bool ignoreConverged =
false;
380 Parameters::get(name_ +
"->solver->ignore converged", ignoreConverged);
381 test_exit(stat.converged || ignoreConverged,
"Could not solver the linear system!");
385 template <
class Traits>
392 for (
auto& currentMarker : marker_)
393 markFlag |= currentMarker.second->markGrid(adaptInfo);
396 int markFlagValue =
int(markFlag);
397 grid_->comm().max(&markFlagValue, 1);
398 markFlag =
Flag(markFlagValue);
400 info(2,
"markElements needed {} seconds", t.elapsed());
406 template <
class Traits>
411 bool adapted =
false;
413 for (
int i = 0; i < n; ++i) {
415 for (
const auto& element : elements(grid_->leafGridView()))
416 grid_->mark(-1, element);
418 bool adaptedInLoop = grid_->preAdapt();
419 adaptedInLoop |= grid_->adapt();
427 info(2,
"globalCoarsen needed {} seconds", t.elapsed());
428 return adapted ? MESH_ADAPTED :
Flag(0);
434 using HasGlobalRefineADHI = decltype(
435 std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
437 template <
class Traits>
442 if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
443 grid_->globalRefine(n, globalBasis_->globalRefineCallback());
445 grid_->globalRefine(n);
447 info(2,
"globalRefine needed {} seconds", t.elapsed());
448 return n > 0 ? MESH_ADAPTED :
Flag(0);
452 template <
class Traits>
458 bool adapted = grid_->preAdapt();
459 adapted |= grid_->adapt();
462 info(2,
"adaptGrid needed {} seconds", t.elapsed());
463 return adapted ? MESH_ADAPTED :
Flag(0);
467 template <
class Traits>
475 for (
auto& bc : constraints_)
482 systemMatrix_->init();
484 rhs_->init(*globalBasis_,
true);
487 if (gridView().comm().size() > 1)
488 msg(
"{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
490 msg(
"{} local DOFs", rhs_->localSize());
493 auto localView = globalBasis_->localView();
494 for (
auto const& element : entitySet(*globalBasis_)) {
495 localView.bind(element);
498 systemMatrix_->assemble(localView, localView);
500 rhs_->assemble(localView);
507 systemMatrix_->finish();
511 info(2,
" assemble operators needed {} seconds", t2.elapsed());
514 solution_->resize(*globalBasis_);
517 for (
auto& bc : constraints_)
518 bc.apply(*systemMatrix_, *solution_, *rhs_);
522 info(2,
" assemble boundary conditions needed {} seconds", t2.elapsed());
524 msg(
"fill-in of assembled matrix: {}", systemMatrix_->nnz());
525 msg(
"assemble needed {} seconds", t.elapsed());
529 template <
class Traits>
534 for (
auto writer : filewriter_)
535 writer->write(adaptInfo, force);
536 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
void addPeriodicBC(BoundaryType id, WorldMatrix const &A, WorldVector const &b)
Definition: ProblemStat.inc.hpp:340
The Flag class encapsulates flags which represents simple information. Used e.g. while mesh traversal...
Definition: Flag.hpp:13
std::map< std::string, std::vector< double > > estimates_
Definition: ProblemStat.hpp:565
constexpr bool Functor
A Functor is a function F with signature Signature.
Definition: Concepts.hpp:133
void addDirichletBC(Predicate const &predicate, RowTreePath row, ColTreePath col, Values const &values)
Add boundary conditions to the system.
Definition: ProblemStat.inc.hpp:298
std::shared_ptr< SystemMatrix > systemMatrix_
Matrix that is filled during assembling.
Definition: ProblemStat.hpp:554
constexpr bool Predicate
A predicate is a function that returns a boolean.
Definition: Concepts.hpp:141
Implements a boundary condition of Dirichlet-type.
Definition: DirichletBC.hpp:37
std::shared_ptr< SystemVector > rhs_
Definition: ProblemStat.hpp:561
std::shared_ptr< SolutionVector > solution_
Vector with the solution components.
Definition: ProblemStat.hpp:557
decltype(auto) makeGridFunction(PreGridFct const &preGridFct, GridView const &gridView)
Generator for Gridfunctions from Expressions (PreGridfunctions)
Definition: GridFunction.hpp:168
Definition: AdaptBase.hpp:6
Definition: BackupRestore.hpp:15
A creator class for dune grids.
Definition: MeshCreator.hpp:51
Definition: ProblemStat.hpp:52
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
std::shared_ptr< GlobalBasis > globalBasis_
FE space of this problem.
Definition: ProblemStat.hpp:539
void restore(Flag initFlag)
Read the grid and solution from backup files and initialize the problem.
Definition: ProblemStat.inc.hpp:108
std::map< std::string, std::shared_ptr< Marker< Grid > > > marker_
Pointer to the adaptation markers.
Definition: ProblemStat.hpp:545
static std::shared_ptr< Grid > create(std::string name)
Static create mthod. See create()
Definition: MeshCreator.hpp:70
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:387
Flag adaptGrid(AdaptInfo &adaptInfo) override
Implementation of ProblemStatBase::refineMesh.
Definition: ProblemStat.inc.hpp:454
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:531
std::shared_ptr< Grid > grid_
Grid of this problem.
Definition: ProblemStat.hpp:530
void solve(AdaptInfo &adaptInfo, bool createMatrixData=true, bool storeMatrixData=false) override
Implementation of ProblemStatBase::solve.
Definition: ProblemStat.inc.hpp:353
Flag globalRefine(int n) override
Uniform global refinement by n level.
Definition: ProblemStat.inc.hpp:439
Holds adapt parameters and infos about the problem.
Definition: AdaptInfo.hpp:25
std::shared_ptr< BoundaryManager< Grid > > boundaryManager_
Management of boundary conditions.
Definition: ProblemStat.hpp:536
void buildAfterAdapt(AdaptInfo &adaptInfo, Flag flag, bool asmMatrix=true, bool asmVector=true) override
Implementation of ProblemStatBase::buildAfterCoarse.
Definition: ProblemStat.inc.hpp:469
std::shared_ptr< LinearSolverInterface > linearSolver_
Pointer to the estimators for this problem.
Definition: ProblemStat.hpp:551
Flag globalCoarsen(int n) override
Uniform global grid coarsening by up to n level.
Definition: ProblemStat.inc.hpp:408
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:139
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