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/io/FileWriterCreator.hpp> 17 #include <amdis/linearalgebra/SymmetryStructure.hpp> 21 template <
class Traits>
29 if (initFlag.
isSet(CREATE_MESH) ||
30 (!adoptFlag.
isSet(INIT_MESH) &&
31 (initFlag.
isSet(INIT_SYSTEM) || initFlag.
isSet(INIT_FE_SPACE)))) {
36 (adoptFlag.
isSet(INIT_MESH) ||
37 adoptFlag.
isSet(INIT_SYSTEM) ||
38 adoptFlag.
isSet(INIT_FE_SPACE))) {
44 warning(
"no grid created");
48 if (initFlag.
isSet(INIT_FE_SPACE) ||
49 (initFlag.
isSet(INIT_SYSTEM) && !adoptFlag.
isSet(INIT_FE_SPACE))) {
54 (adoptFlag.
isSet(INIT_FE_SPACE) || adoptFlag.
isSet(INIT_SYSTEM))) {
60 warning(
"no globalBasis created\n");
63 if (initFlag.
isSet(INIT_SYSTEM))
64 createMatricesAndVectors();
66 if (adoptProblem && adoptFlag.
isSet(INIT_SYSTEM)) {
69 rhs_ = adoptProblem->
rhs_;
76 if (initFlag.
isSet(INIT_SOLVER))
79 if (adoptProblem && adoptFlag.
isSet(INIT_SOLVER)) {
80 test_exit(!linearSolver_,
"solver already created\n");
86 warning(
"no solver created\n");
90 if (initFlag.
isSet(INIT_MARKER))
93 if (adoptProblem && adoptFlag.
isSet(INIT_MARKER))
94 marker_ = adoptProblem->
marker_;
98 if (initFlag.
isSet(INIT_FILEWRITER))
101 solution_->resizeZero();
105 template <
class Traits>
109 std::string grid_filename = Parameters::get<std::string>(name_ +
"->restore->grid").value();
110 std::string solution_filename = Parameters::get<std::string>(name_ +
"->restore->solution").value();
111 test_exit(filesystem::exists(grid_filename),
"Restore file '{}' not found.", grid_filename);
112 test_exit(filesystem::exists(solution_filename),
"Restore file '{}' not found.", solution_filename);
115 using HostGrid =
typename Grid::HostGrid;
118 if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
119 adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
124 if (initFlag.
isSet(INIT_FE_SPACE) || initFlag.
isSet(INIT_SYSTEM))
128 if (initFlag.
isSet(INIT_SYSTEM))
129 createMatricesAndVectors();
132 if (!linearSolver_ && initFlag.
isSet(INIT_SOLVER))
136 if (initFlag.
isSet(INIT_MARKER))
140 if (initFlag.
isSet(INIT_FILEWRITER))
143 solution_->resize(sizeInfo(*globalBasis_));
144 solution_->restore(solution_filename);
148 template <
class Traits>
156 boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
158 boundaryManager_->setBoundaryIds(creator.
boundaryIds());
160 info(3,
"Create grid:");
161 info(3,
"#elements = {}" , grid_->size(0));
162 info(3,
"#faces/edges = {}", grid_->size(1));
163 info(3,
"#vertices = {}" , grid_->size(dim));
164 info(3,
"overlap-size = {}", grid_->leafGridView().overlapSize(0));
165 info(3,
"ghost-size = {}" , grid_->leafGridView().ghostSize(0));
170 template <
class T,
class GV>
171 using HasCreate = decltype(T::create(std::declval<GV>()));
174 template <
class Traits>
177 createGlobalBasisImpl(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
182 template <
class Traits>
185 assert(
bool(grid_) );
186 static_assert(std::is_same_v<GridView, typename Grid::LeafGridView>,
"");
187 auto basis = Traits::create(name_, grid_->leafGridView());
188 globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis));
192 template <
class Traits>
195 error_exit(
"Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
199 template <
class Traits>
203 template <
class Traits>
206 systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
207 std::string symmetryStr =
"unknown";
209 systemMatrix_->setSymmetryStructure(symmetryStr);
211 solution_ = std::make_shared<SolutionVector>(globalBasis_);
212 rhs_ = std::make_shared<SystemVector>(globalBasis_);
214 auto localView = globalBasis_->localView();
215 Traversal::forEachNode(localView.tree(), [&,
this](
auto&&,
auto treePath) ->
void 217 std::string i = to_string(treePath);
218 estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
219 for (std::size_t j = 0; j < estimates_[i].size(); j++)
220 estimates_[i][j] = 0.0;
225 template <
class Traits>
228 std::string solverName =
"default";
231 linearSolver_ = std::make_shared<LinearSolver>(solverName, name_ +
"->solver");
235 template <
class Traits>
239 auto localView = globalBasis_->localView();
240 Traversal::forEachNode(localView.tree(), [&,
this](
auto&&,
auto treePath) ->
void 242 std::string componentName = name_ +
"->marker[" + to_string(treePath) +
"]";
244 if (!Parameters::get<std::string>(componentName +
"->strategy"))
247 std::string tp = to_string(treePath);
250 assert(
bool(newMarker));
251 this->addMarker(std::move(newMarker));
256 template <
class Traits>
262 auto localView = globalBasis_->localView();
263 Traversal::forEachNode(localView.tree(), [&](
auto const& ,
auto treePath) ->
void 265 std::string componentName = name_ +
"->output[" + to_string(treePath) +
"]";
266 auto format = Parameters::get<std::vector<std::string>>(componentName +
"->format");
268 if (!format && to_string(treePath).empty()) {
270 componentName = name_ +
"->output";
271 format = Parameters::get<std::vector<std::string>>(componentName +
"->format");
277 for (std::string
const& type : format.value()) {
278 auto writer = creator.
create(type, componentName, treePath);
280 filewriter_.push_back(std::move(writer));
287 template <
class Traits>
288 template <
class Predicate,
class RowTreePath,
class ColTreePath,
class Values>
292 static constexpr
bool isValidPredicate = Concepts::Functor<Predicate, bool(WorldVector)>;
294 "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
296 static constexpr
bool isValidTreePath =
297 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
298 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
299 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
301 if constexpr (isValidPredicate && isValidTreePath) {
302 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col));
306 std::move(basis), {predicate}, valueGridFct});
312 template <
class Traits>
313 template <
class RowTreePath,
class ColTreePath,
class Values>
317 static constexpr
bool isValidTreePath =
318 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
319 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
320 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
322 if constexpr (isValidTreePath) {
323 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col));
327 std::move(basis), {*boundaryManager_,
id}, valueGridFct});
332 template <
class Traits>
336 auto localView = globalBasis_->localView();
337 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath());
338 constraints_.push_back(makePeriodicBC(
339 std::move(basis), {*boundaryManager_,
id}, {matrix, vector}));
344 template <
class Traits>
348 Dune::InverseOperatorResult stat;
352 if (createMatrixData)
353 linearSolver_->init(systemMatrix_->impl());
356 linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);
358 if (!storeMatrixData)
359 linearSolver_->finish();
361 info(2,
"solution of discrete system needed {} seconds", stat.elapsed);
362 info(1,
"Residual norm: ||b-Ax|| = {}",
363 residuum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
365 if (stat.reduction >= 0.0)
366 info(2,
"Relative residual norm: ||b-Ax||/||b|| = {}", stat.reduction);
368 info(2,
"Relative residual norm: ||b-Ax||/||b|| = {}",
369 relResiduum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
371 bool ignoreConverged =
false;
372 Parameters::get(name_ +
"->solver->ignore converged", ignoreConverged);
373 test_exit(stat.converged || ignoreConverged,
"Could not solver the linear system!");
377 template <
class Traits>
384 for (
auto& currentMarker : marker_)
385 markFlag |= currentMarker.second->markGrid(adaptInfo);
388 int markFlagValue =
int(markFlag);
389 Dune::MPIHelper::getCollectiveCommunication().max(&markFlagValue, 1);
390 markFlag =
Flag(markFlagValue);
392 msg(
"markElements needed {} seconds", t.elapsed());
398 template <
class Traits>
403 bool adapted =
false;
405 for (
int i = 0; i < n; ++i) {
407 for (
const auto& element : elements(grid_->leafGridView()))
408 grid_->mark(-1, element);
410 bool adaptedInLoop = grid_->preAdapt();
411 adaptedInLoop |= grid_->adapt();
419 msg(
"globalCoarsen needed {} seconds", t.elapsed());
420 return adapted ? MESH_ADAPTED :
Flag(0);
426 using HasGlobalRefineADHI = decltype(
427 std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
429 template <
class Traits>
434 if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
435 grid_->globalRefine(n, globalBasis_->globalRefineCallback());
437 grid_->globalRefine(n);
439 msg(
"globalRefine needed {} seconds", t.elapsed());
440 return n > 0 ? MESH_ADAPTED :
Flag(0);
444 template <
class Traits>
450 bool adapted = grid_->preAdapt();
451 adapted |= grid_->adapt();
454 msg(
"adaptGrid needed {} seconds", t.elapsed());
455 return adapted ? MESH_ADAPTED :
Flag(0);
459 template <
class Traits>
467 for (
auto& bc : constraints_)
473 systemMatrix_->init();
474 rhs_->init(sizeInfo(*globalBasis_), asmVector);
478 msg(
"{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
480 msg(
"{} local DOFs", rhs_->localSize());
482 auto localMatOperators =
localOperators(systemMatrix_->operators());
486 auto localView = globalBasis_->localView();
487 for (
auto const& element : elements(gridView(), PartitionSet{})) {
488 localView.bind(element);
491 systemMatrix_->assemble(localView, localView, localMatOperators);
493 rhs_->assemble(localView, localVecOperators);
499 systemMatrix_->finish();
502 info(2,
" assemble operators needed {} seconds", t2.elapsed());
505 solution_->resize(sizeInfo(*globalBasis_));
508 for (
auto& bc : constraints_)
509 bc.apply(*systemMatrix_, *solution_, *rhs_);
511 info(2,
" assemble boundary conditions needed {} seconds", t2.elapsed());
513 msg(
"fill-in of assembled matrix: {}", systemMatrix_->nnz());
514 msg(
"assemble needed {} seconds", t.elapsed());
518 template <
class Traits>
523 for (
auto writer : filewriter_)
524 writer->write(adaptInfo, force);
525 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:334
The Flag class encapsulates flags which represents simple information. Used e.g. while mesh traversal...
Definition: Flag.hpp:13
auto localOperators(Container const &c)
Definition: Operator.hpp:125
std::map< std::string, std::vector< double > > estimates_
Definition: ProblemStat.hpp:555
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:290
std::shared_ptr< SystemMatrix > systemMatrix_
Matrix that is filled during assembling.
Definition: ProblemStat.hpp:544
constexpr bool Predicate
A predicate is a function that returns a boolean.
Definition: Concepts.hpp:142
Implements a boundary condition of Dirichlet-type.
Definition: DirichletBC.hpp:39
std::shared_ptr< SystemVector > rhs_
Definition: ProblemStat.hpp:551
std::shared_ptr< SolutionVector > solution_
Vector with the solution components.
Definition: ProblemStat.hpp:547
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:52
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:22
std::shared_ptr< GlobalBasis > globalBasis_
FE space of this problem.
Definition: ProblemStat.hpp:529
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:107
std::map< std::string, std::shared_ptr< Marker< Grid > > > marker_
Pointer to the adaptation markers.
Definition: ProblemStat.hpp:535
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:379
Flag adaptGrid(AdaptInfo &adaptInfo) override
Implementation of ProblemStatBase::refineMesh.
Definition: ProblemStat.inc.hpp:446
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:520
std::shared_ptr< Grid > grid_
Grid of this problem.
Definition: ProblemStat.hpp:520
void solve(AdaptInfo &adaptInfo, bool createMatrixData=true, bool storeMatrixData=false) override
Implementation of ProblemStatBase::solve.
Definition: ProblemStat.inc.hpp:346
Flag globalRefine(int n) override
Uniform global refinement by n level.
Definition: ProblemStat.inc.hpp:431
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:526
void buildAfterAdapt(AdaptInfo &adaptInfo, Flag flag, bool asmMatrix=true, bool asmVector=true) override
Implementation of ProblemStatBase::buildAfterCoarse.
Definition: ProblemStat.inc.hpp:461
std::shared_ptr< LinearSolverInterface > linearSolver_
Pointer to the estimators for this problem.
Definition: ProblemStat.hpp:541
Flag globalCoarsen(int n) override
Uniform global grid coarsening by up to n level.
Definition: ProblemStat.inc.hpp:400
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