AMDiS  2.10
The Adaptive Multi-Dimensional Simulation Toolbox
ProblemStat.inc.hpp
1 #pragma once
2 
3 #include <map>
4 #include <string>
5 #include <utility>
6 
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>
12 
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>
19 
20 namespace AMDiS {
21 
22 template <class Traits>
24  Flag initFlag,
25  Self* adoptProblem,
26  Flag adoptFlag)
27 {
28  // create grids
29  if (!grid_) {
30  if (initFlag.isSet(CREATE_MESH) ||
31  (!adoptFlag.isSet(INIT_MESH) &&
32  (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
33  createGrid();
34  }
35 
36  if (adoptProblem &&
37  (adoptFlag.isSet(INIT_MESH) ||
38  adoptFlag.isSet(INIT_SYSTEM) ||
39  adoptFlag.isSet(INIT_FE_SPACE))) {
40  adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
41  }
42  }
43 
44  if (!grid_)
45  warning("no grid created");
46 
47  // create fespace
48  if (!globalBasis_) {
49  if (initFlag.isSet(INIT_FE_SPACE) ||
50  (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
51  createGlobalBasis();
52  }
53 
54  if (adoptProblem &&
55  (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
56  adoptGlobalBasis(adoptProblem->globalBasis_);
57  }
58  }
59 
60  if (!globalBasis_)
61  warning("no globalBasis created\n");
62 
63  // create system
64  if (initFlag.isSet(INIT_SYSTEM))
65  createMatricesAndVectors();
66 
67  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
68  systemMatrix_ = adoptProblem->systemMatrix_;
69  solution_ = adoptProblem->solution_;
70  rhs_ = adoptProblem->rhs_;
71  estimates_ = adoptProblem->estimates_;
72  }
73 
74 
75  // create solver
76  if (!linearSolver_) {
77  if (initFlag.isSet(INIT_SOLVER))
78  createSolver();
79 
80  if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
81  test_exit(!linearSolver_, "solver already created\n");
82  linearSolver_ = adoptProblem->linearSolver_;
83  }
84  }
85 
86  if (!linearSolver_) {
87  warning("no solver created\n");
88  }
89 
90  // create marker
91  if (initFlag.isSet(INIT_MARKER))
92  createMarker();
93 
94  if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
95  marker_ = adoptProblem->marker_;
96 
97 
98  // create file writer
99  if (initFlag.isSet(INIT_FILEWRITER))
100  createFileWriter();
101 
102  solution_->resizeZero();
103 }
104 
105 
106 template <class Traits>
108 restore(Flag initFlag)
109 {
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);
114 
115  // TODO(SP): implement BAckupRestore independent of wrapped grid
116  using HostGrid = typename Grid::HostGrid;
117 
118  // restore grid from file
119  if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
120  adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
121  else
122  adoptGrid(std::shared_ptr<HostGrid>(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename)));
123 
124  // create fespace
125  if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))
126  createGlobalBasis();
127 
128  // create system
129  if (initFlag.isSet(INIT_SYSTEM))
130  createMatricesAndVectors();
131 
132  // create solver
133  if (!linearSolver_ && initFlag.isSet(INIT_SOLVER))
134  createSolver();
135 
136  // create marker
137  if (initFlag.isSet(INIT_MARKER))
138  createMarker();
139 
140  // create file writer
141  if (initFlag.isSet(INIT_FILEWRITER))
142  createFileWriter();
143 
144  solution_->resize(*globalBasis_);
145  solution_->restore(solution_filename);
146 }
147 
148 
149 template <class Traits>
151 {
152  Parameters::get(name_ + "->mesh", gridName_);
153 
154  MeshCreator<Grid> creator(gridName_);
155  grid_ = creator.create();
156 
157  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
158  if (!creator.boundaryIds().empty())
159  boundaryManager_->setBoundaryIds(creator.boundaryIds());
160 
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));
167  info(3,"");
168 }
169 
170 
171 template <class T, class GV>
172 using HasCreate = decltype(T::create(std::declval<GV>()));
173 
174 
175 template <class Traits>
177 {
178  createGlobalBasisImpl(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
179  initGlobalBasis();
180 }
181 
182 
183 template <class Traits>
185 {
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));
190 }
191 
192 
193 template <class Traits>
194 void ProblemStat<Traits>::createGlobalBasisImpl(std::false_type)
195 {
196  error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
197 }
198 
199 
200 template <class Traits>
202 
203 
204 template <class Traits>
206 {
207  systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
208  std::string symmetryStr = "unknown";
209  Parameters::get(name_ + "->symmetry", symmetryStr);
210  systemMatrix_->setSymmetryStructure(symmetryStr);
211 
212  solution_ = std::make_shared<SolutionVector>(globalBasis_);
213  rhs_ = std::make_shared<SystemVector>(globalBasis_);
214 
215  auto localView = globalBasis_->localView();
216  Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void
217  {
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; // TODO: Remove when estimate() is implemented
222  });
223 }
224 
225 
226 template <class Traits>
228 {
229  std::string solverName = "default";
230  Parameters::get(name_ + "->solver", solverName);
231 
232  linearSolver_ = std::make_shared<LinearSolver>(solverName, name_ + "->solver");
233 }
234 
235 
236 template <class Traits>
238 {
239  marker_.clear();
240  auto localView = globalBasis_->localView();
241  Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void
242  {
243  std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
244  auto strategy = Parameters::get<std::string>(componentName + "->strategy");
245 
246  if (!strategy && to_string(treePath).empty()) {
247  // alternative for root treepath
248  componentName = name_ + "->strategy";
249  strategy = Parameters::get<std::string>(componentName + "->strategy");
250  }
251 
252  if (!strategy)
253  return;
254 
255  std::string tp = to_string(treePath);
256  auto newMarker
257  = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
258  assert(bool(newMarker));
259  this->addMarker(std::move(newMarker));
260  });
261 }
262 
263 
264 template <class Traits>
266 {
267  FileWriterCreator<SolutionVector> creator(solution_, boundaryManager_);
268 
269  filewriter_.clear();
270  auto localView = globalBasis_->localView();
271  Traversal::forEachNode(localView.tree(), [&](auto const& /*node*/, auto treePath) -> void
272  {
273  std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
274  auto format = Parameters::get<std::vector<std::string>>(componentName + "->format");
275 
276  if (!format && to_string(treePath).empty()) {
277  // alternative for root treepath
278  componentName = name_ + "->output";
279  format = Parameters::get<std::vector<std::string>>(componentName + "->format");
280  }
281 
282  if (!format)
283  return;
284 
285  for (std::string const& type : format.value()) {
286  auto writer = creator.create(type, componentName, treePath);
287  if (writer)
288  filewriter_.push_back(std::move(writer));
289  }
290  });
291 }
292 
293 
294 // Adds a Dirichlet boundary condition
295 template <class Traits>
296  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
298 addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
299 {
300  static constexpr bool isValidPredicate = Concepts::Functor<Predicate, bool(WorldVector)>;
301  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
302  "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
303 
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!");
308 
309  if constexpr (isValidPredicate && isValidTreePath) {
310  auto valueGridFct = makeGridFunction(values, this->gridView());
311 
312  constraints_.push_back(DirichletBC{
313  globalBasis_, makeTreePath(row), makeTreePath(col), {predicate}, valueGridFct});
314  }
315 }
316 
317 
318 // Adds a Dirichlet boundary condition
319 template <class Traits>
320  template <class RowTreePath, class ColTreePath, class Values>
322 addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& values)
323 {
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!");
328 
329  if constexpr (isValidTreePath) {
330  auto valueGridFct = makeGridFunction(values, this->gridView());
331 
332  constraints_.push_back(DirichletBC{
333  globalBasis_, makeTreePath(row), makeTreePath(col), {boundaryManager_, id}, valueGridFct});
334  }
335 }
336 
337 
338 template <class Traits>
340 addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
341 {
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}));
346  });
347 }
348 
349 
350 
351 template <class Traits>
353 solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData)
354 {
355  Dune::Timer t;
356  Dune::InverseOperatorResult stat;
357 
358  solution_->resize();
359 
360  if (createMatrixData)
361  linearSolver_->init(systemMatrix_->impl());
362 
363  // solve the linear system
364  linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);
365 
366  if (!storeMatrixData)
367  linearSolver_->finish();
368 
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()));
372 
373  if (stat.reduction >= 0.0)
374  info(2, "Relative residual norm: ||b-Ax||/||b|| = {}", stat.reduction);
375  else
376  info(2, "Relative residual norm: ||b-Ax||/||b|| = {}",
377  relResiduum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
378 
379  bool ignoreConverged = false;
380  Parameters::get(name_ + "->solver->ignore converged", ignoreConverged);
381  test_exit(stat.converged || ignoreConverged, "Could not solver the linear system!");
382 }
383 
384 
385 template <class Traits>
388 {
389  Dune::Timer t;
390 
391  Flag markFlag = 0;
392  for (auto& currentMarker : marker_)
393  markFlag |= currentMarker.second->markGrid(adaptInfo);
394 
395  // synchronize mark flag over processors
396  int markFlagValue = int(markFlag);
397  grid_->comm().max(&markFlagValue, 1);
398  markFlag = Flag(markFlagValue);
399 
400  info(2, "markElements needed {} seconds", t.elapsed());
401 
402  return markFlag;
403 }
404 
405 
406 template <class Traits>
409 {
410  Dune::Timer t;
411  bool adapted = false;
412  // TODO(FM): Find a less expensive alternative to the loop adaption
413  for (int i = 0; i < n; ++i) {
414  // mark all entities for coarsening
415  for (const auto& element : elements(grid_->leafGridView()))
416  grid_->mark(-1, element);
417 
418  bool adaptedInLoop = grid_->preAdapt();
419  adaptedInLoop |= grid_->adapt();
420  grid_->postAdapt();
421  if (!adaptedInLoop)
422  break;
423  else
424  adapted = true;
425  }
426 
427  info(2, "globalCoarsen needed {} seconds", t.elapsed());
428  return adapted ? MESH_ADAPTED : Flag(0);
429 }
430 
431 
432 // grid has globalRefine(int, AdaptDataHandleInterface&)
433 template <class G>
434 using HasGlobalRefineADHI = decltype(
435  std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
436 
437 template <class Traits>
440 {
441  Dune::Timer t;
442  if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
443  grid_->globalRefine(n, globalBasis_->globalRefineCallback());
444  else
445  grid_->globalRefine(n);
446 
447  info(2, "globalRefine needed {} seconds", t.elapsed());
448  return n > 0 ? MESH_ADAPTED : Flag(0);
449 }
450 
451 
452 template <class Traits>
454 adaptGrid(AdaptInfo& /*adaptInfo*/)
455 {
456  Dune::Timer t;
457 
458  bool adapted = grid_->preAdapt();
459  adapted |= grid_->adapt();
460  grid_->postAdapt();
461 
462  info(2, "adaptGrid needed {} seconds", t.elapsed());
463  return adapted ? MESH_ADAPTED : Flag(0);
464 }
465 
466 
467 template <class Traits>
469 buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
470 {
471  Dune::Timer t;
472  Dune::Timer t2;
473 
474  // 0. initialize boundary condition and other constraints
475  for (auto& bc : constraints_)
476  bc.init();
477 
478  t2.reset();
479 
480  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
481  if (asmMatrix)
482  systemMatrix_->init();
483  if (asmVector)
484  rhs_->init(*globalBasis_, true);
485 
486  // statistic about system size
487  if (gridView().comm().size() > 1)
488  msg("{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
489  else
490  msg("{} local DOFs", rhs_->localSize());
491 
492  // 2. traverse grid and assemble operators on the elements
493  auto localView = globalBasis_->localView();
494  for (auto const& element : entitySet(*globalBasis_)) {
495  localView.bind(element);
496 
497  if (asmMatrix)
498  systemMatrix_->assemble(localView, localView);
499  if (asmVector)
500  rhs_->assemble(localView);
501 
502  localView.unbind();
503  }
504 
505  // 3. finish matrix insertion and apply dirichlet boundary conditions
506  if (asmMatrix)
507  systemMatrix_->finish();
508  if (asmVector)
509  rhs_->finish();
510 
511  info(2," assemble operators needed {} seconds", t2.elapsed());
512  t2.reset();
513 
514  solution_->resize(*globalBasis_);
515 
516  // 4. apply boundary condition and constraints to matrix, solution, and rhs
517  for (auto& bc : constraints_)
518  bc.apply(*systemMatrix_, *solution_, *rhs_);
519 
520  rhs_->finish();
521  solution_->finish();
522  info(2," assemble boundary conditions needed {} seconds", t2.elapsed());
523 
524  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
525  msg("assemble needed {} seconds", t.elapsed());
526 }
527 
528 
529 template <class Traits>
531 writeFiles(AdaptInfo& adaptInfo, bool force)
532 {
533  Dune::Timer t;
534  for (auto writer : filewriter_)
535  writer->write(adaptInfo, force);
536  msg("writeFiles needed {} seconds", t.elapsed());
537 }
538 
539 } // end namespace AMDiS
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