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))) {
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(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 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";
235 linearSolver_ = solverCreator->createWithString(name_ +
"->solver");
239 template <
class Traits>
243 auto localView = globalBasis_->localView();
244 Traversal::forEachNode(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 Traversal::forEachNode(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 basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col));
310 std::move(basis), {predicate}, valueGridFct});
316 template <
class Traits>
317 template <
class RowTreePath,
class ColTreePath,
class Values>
321 static constexpr
bool isValidTreePath =
322 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
323 Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
324 static_assert(isValidTreePath,
"Invalid row and/or col treepath passed to addDirichletBC!");
326 if constexpr (isValidTreePath) {
327 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col));
331 std::move(basis), {*boundaryManager_,
id}, valueGridFct});
336 template <
class Traits>
340 auto localView = globalBasis_->localView();
341 auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath());
343 std::move(basis), {*boundaryManager_,
id}, {matrix, vector}));
348 template <
class Traits>
359 linearSolver_->solve(*systemMatrix_, *solution_, *rhs_, solverInfo);
361 if (solverInfo.
info() > 0) {
362 msg(
"solution of discrete system needed {} seconds", t.elapsed());
366 msg(
"Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
373 test_exit(!solverInfo.doBreak() || !solverInfo.
error(),
"Could not solver the linear system!");
377 template <
class Traits>
384 for (
auto& currentMarker : marker_)
385 markFlag |= currentMarker.second->markGrid(adaptInfo);
387 msg(
"markElements needed {} seconds", t.elapsed());
393 template <
class Traits>
398 bool adapted =
false;
400 for (
int i = 0; i < n; ++i) {
402 for (
const auto& element : elements(grid_->leafGridView()))
403 grid_->mark(-1, element);
405 bool adaptedInLoop = grid_->preAdapt();
406 adaptedInLoop |= grid_->adapt();
414 msg(
"globalCoarsen needed {} seconds", t.elapsed());
415 return adapted ? MESH_ADAPTED :
Flag(0);
421 using HasGlobalRefineADHI = decltype(
422 std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
424 template <
class Traits>
429 if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
430 grid_->globalRefine(n, globalBasis_->globalRefineCallback());
432 grid_->globalRefine(n);
434 msg(
"globalRefine needed {} seconds", t.elapsed());
435 return n > 0 ? MESH_ADAPTED :
Flag(0);
439 template <
class Traits>
445 bool adapted = grid_->preAdapt();
446 adapted |= grid_->adapt();
449 msg(
"adaptGrid needed {} seconds", t.elapsed());
450 return adapted ? MESH_ADAPTED :
Flag(0);
454 template <
class Traits>
462 for (
auto& bc : constraints_)
468 systemMatrix_->init();
469 rhs_->init(sizeInfo(*globalBasis_), asmVector);
473 msg(
"{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
475 msg(
"{} local DOFs", rhs_->localSize());
478 auto localView = globalBasis_->localView();
479 for (
auto const& element : elements(gridView(), PartitionSet{})) {
480 localView.bind(element);
483 systemMatrix_->assemble(localView, localView);
485 rhs_->assemble(localView);
491 systemMatrix_->finish();
494 info(2,
" assemble operators needed {} seconds", t2.elapsed());
497 solution_->resize(sizeInfo(*globalBasis_));
500 for (
auto& bc : constraints_)
501 bc.apply(*systemMatrix_, *solution_, *rhs_);
503 info(2,
" assemble boundary conditions needed {} seconds", t2.elapsed());
505 msg(
"fill-in of assembled matrix: {}", systemMatrix_->nnz());
506 msg(
"assemble needed {} seconds", t.elapsed());
510 template <
class Traits>
515 for (
auto writer : filewriter_)
516 writer->write(adaptInfo, force);
517 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:338
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:553
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
std::shared_ptr< SystemMatrix > systemMatrix_
Matrix that is filled during assembling.
Definition: ProblemStat.hpp:542
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:549
std::shared_ptr< SolutionVector > solution_
Vector with the solution components.
Definition: ProblemStat.hpp:545
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
std::shared_ptr< GlobalBasis > globalBasis_
FE space of this problem.
Definition: ProblemStat.hpp:527
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
std::map< std::string, std::shared_ptr< Marker< Grid > > > marker_
Pointer to the adaptation markers.
Definition: ProblemStat.hpp:533
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:379
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:441
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:512
std::shared_ptr< Grid > grid_
Grid of this problem.
Definition: ProblemStat.hpp:518
void solve(AdaptInfo &adaptInfo, bool createMatrixData=true, bool storeMatrixData=false) override
Implementation of ProblemStatBase::solve.
Definition: ProblemStat.inc.hpp:350
Flag globalRefine(int n) override
Uniform global refinement by n level.
Definition: ProblemStat.inc.hpp:426
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
Definition: SolverInfo.hpp:11
std::shared_ptr< BoundaryManager< Grid > > boundaryManager_
Management of boundary conditions.
Definition: ProblemStat.hpp:524
void buildAfterAdapt(AdaptInfo &adaptInfo, Flag flag, bool asmMatrix=true, bool asmVector=true) override
Implementation of ProblemStatBase::buildAfterCoarse.
Definition: ProblemStat.inc.hpp:456
Flag globalCoarsen(int n) override
Uniform global grid coarsening by up to n level.
Definition: ProblemStat.inc.hpp:395
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::shared_ptr< LinearSolver > linearSolver_
Pointer to the estimators for this problem.
Definition: ProblemStat.hpp:539
auto makePeriodicBC(Dune::Functions::SubspaceBasis< B, TP > basis, BoundarySubset< typename B::GridView::Intersection > boundarySubset, Impl::FaceTrafo< B > trafo)
Make a PeriodicBC from a subspacebasis.
Definition: PeriodicBC.hpp:137
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