ProblemStat.inc.hpp 14.7 KB
Newer Older
1
2
#pragma once

3
4
5
6
#include <map>
#include <string>
#include <utility>

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
12
#include <dune/typetree/childextraction.hh>

13
#include <amdis/AdaptInfo.hpp>
14
#include <amdis/BackupRestore.hpp>
15
#include <amdis/GridFunctionOperator.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
16
#include <amdis/io/FileWriterCreator.hpp>
17
#include <amdis/linearalgebra/SymmetryStructure.hpp>
18

19
20
namespace AMDiS {

21
22
template <class Traits>
void ProblemStat<Traits>::initialize(
23
24
25
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
26
{
27
  // create grids
28
  if (!grid_) {
29
30
31
    if (initFlag.isSet(CREATE_MESH) ||
        (!adoptFlag.isSet(INIT_MESH) &&
        (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
32
      createGrid();
33
    }
34

35
36
37
38
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
39
      adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
40
    }
41
  }
42

43
  if (!grid_)
44
    warning("no grid created");
45

46
  // create fespace
47
  if (!globalBasis_) {
48
49
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
50
      createGlobalBasis();
51
    }
52

53
54
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Praetorius, Simon's avatar
Praetorius, Simon committed
55
      adoptGlobalBasis(adoptProblem->globalBasis_);
56
    }
57
  }
58

59
  if (!globalBasis_)
60
    warning("no globalBasis created\n");
61

62
63
64
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
65

66
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
67
    systemMatrix_ = adoptProblem->systemMatrix_;
68
69
    solution_ = adoptProblem->solution_;
    rhs_ = adoptProblem->rhs_;
Praetorius, Simon's avatar
Praetorius, Simon committed
70
    estimates_ = adoptProblem->estimates_;
71
  }
72

73

74
  // create solver
75
  if (!linearSolver_) {
76
77
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
78

79
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
80
81
      test_exit(!linearSolver_, "solver already created\n");
      linearSolver_ = adoptProblem->linearSolver_;
82
    }
83
84
  }

85
  if (!linearSolver_) {
86
    warning("no solver created\n");
87
88
  }

89
90
91
92
93
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
94
      marker_ = adoptProblem->marker_;
95

96

97
98
99
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
100

101
  solution_->resizeZero();
102
}
103

104

Praetorius, Simon's avatar
Praetorius, Simon committed
105
106
107
108
109
110
111
112
113
template <class Traits>
void ProblemStat<Traits>::
restore(Flag initFlag)
{
  std::string grid_filename = Parameters::get<std::string>(name_ + "->restore->grid").value();
  std::string solution_filename = Parameters::get<std::string>(name_ + "->restore->solution").value();
  test_exit(filesystem::exists(grid_filename), "Restore file '{}' not found.", grid_filename);
  test_exit(filesystem::exists(solution_filename), "Restore file '{}' not found.", solution_filename);

Praetorius, Simon's avatar
Praetorius, Simon committed
114
115
116
  // TODO(SP): implement BAckupRestore independent of wrapped grid
  using HostGrid = typename Grid::HostGrid;

Praetorius, Simon's avatar
Praetorius, Simon committed
117
  // restore grid from file
118
  if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
Praetorius, Simon's avatar
Praetorius, Simon committed
119
    adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
Praetorius, Simon's avatar
Praetorius, Simon committed
120
  else
Praetorius, Simon's avatar
Praetorius, Simon committed
121
    adoptGrid(std::shared_ptr<HostGrid>(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename)));
Praetorius, Simon's avatar
Praetorius, Simon committed
122
123
124
125
126
127
128
129
130
131

  // create fespace
  if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))
    createGlobalBasis();

  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();

  // create solver
132
  if (!linearSolver_ && initFlag.isSet(INIT_SOLVER))
Praetorius, Simon's avatar
Praetorius, Simon committed
133
134
135
136
137
138
139
140
141
142
    createSolver();

  // create marker
  if (initFlag.isSet(INIT_MARKER))
    createMarker();

  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();

143
  solution_->resize(sizeInfo(*globalBasis_));
Praetorius, Simon's avatar
Praetorius, Simon committed
144
145
146
147
  solution_->restore(solution_filename);
}


Praetorius, Simon's avatar
Praetorius, Simon committed
148
149
150
151
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
Praetorius, Simon's avatar
Praetorius, Simon committed
152
153

  MeshCreator<Grid> creator(gridName_);
154
  grid_ = creator.create();
155

156
  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
157
158
  if (!creator.boundaryIds().empty())
    boundaryManager_->setBoundaryIds(creator.boundaryIds());
Praetorius, Simon's avatar
Praetorius, Simon committed
159

160
161
162
163
164
165
166
  info(3,"Create grid:");
  info(3,"#elements = {}"   , grid_->size(0));
  info(3,"#faces/edges = {}", grid_->size(1));
  info(3,"#vertices = {}"   , grid_->size(dim));
  info(3,"overlap-size = {}", grid_->leafGridView().overlapSize(0));
  info(3,"ghost-size = {}"  , grid_->leafGridView().ghostSize(0));
  info(3,"");
Praetorius, Simon's avatar
Praetorius, Simon committed
167
168
169
170
171
172
173
174
175
176
177
}


template <class T, class GV>
using HasCreate = decltype(T::create(std::declval<GV>()));


template <class Traits>
void ProblemStat<Traits>::createGlobalBasis()
{
  createGlobalBasisImpl(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
178
  initGlobalBasis();
Praetorius, Simon's avatar
Praetorius, Simon committed
179
180
181
182
183
184
185
}


template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type)
{
  assert( bool(grid_) );
186
  static_assert(std::is_same_v<GridView, typename Grid::LeafGridView>, "");
187
188
  auto basis = Traits::create(name_, grid_->leafGridView());
  globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis));
Praetorius, Simon's avatar
Praetorius, Simon committed
189
190
191
192
193
194
195
196
197
198
199
}


template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::false_type)
{
  error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
}


template <class Traits>
200
void ProblemStat<Traits>::initGlobalBasis() {}
Praetorius, Simon's avatar
Praetorius, Simon committed
201
202
203
204
205


template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
206
  systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
207
208
209
210
  std::string symmetryStr = "unknown";
  Parameters::get(name_ + "->symmetry", symmetryStr);
  systemMatrix_->setSymmetryStructure(symmetryStr);

Praetorius, Simon's avatar
Praetorius, Simon committed
211
  solution_ = std::make_shared<SolutionVector>(globalBasis_);
212
  rhs_ = std::make_shared<SystemVector>(globalBasis_);
213

Praetorius, Simon's avatar
Praetorius, Simon committed
214
  auto localView = globalBasis_->localView();
215
  Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void
Praetorius, Simon's avatar
Praetorius, Simon committed
216
217
218
219
220
221
222
223
224
225
226
227
228
  {
    std::string i = to_string(treePath);
    estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
    for (std::size_t j = 0; j < estimates_[i].size(); j++)
      estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
  });
}


template <class Traits>
void ProblemStat<Traits>::createSolver()
{
  std::string solverName = "default";
229
  Parameters::get(name_ + "->solver", solverName);
Praetorius, Simon's avatar
Praetorius, Simon committed
230

231
  linearSolver_ = std::make_shared<LinearSolver>(solverName, name_ + "->solver");
Praetorius, Simon's avatar
Praetorius, Simon committed
232
233
234
}


235
236
237
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
Praetorius, Simon's avatar
Praetorius, Simon committed
238
  marker_.clear();
239
  auto localView = globalBasis_->localView();
240
  Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void
241
  {
242
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
243
244
245
246

    if (!Parameters::get<std::string>(componentName + "->strategy"))
      return;

247
    std::string tp = to_string(treePath);
Praetorius, Simon's avatar
Praetorius, Simon committed
248
249
    auto newMarker
      = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
250
    assert(bool(newMarker));
251
    this->addMarker(std::move(newMarker));
252
253
254
255
  });
}


256
257
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
258
{
Praetorius, Simon's avatar
Praetorius, Simon committed
259
  FileWriterCreator<SolutionVector> creator(solution_, boundaryManager_);
Praetorius, Simon's avatar
Praetorius, Simon committed
260

Praetorius, Simon's avatar
Praetorius, Simon committed
261
  filewriter_.clear();
262
  auto localView = globalBasis_->localView();
263
  Traversal::forEachNode(localView.tree(), [&](auto const& /*node*/, auto treePath) -> void
264
  {
265
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
Praetorius, Simon's avatar
Praetorius, Simon committed
266
267
268
269
270
271
272
    auto format = Parameters::get<std::vector<std::string>>(componentName + "->format");

    if (!format && to_string(treePath).empty()) {
      // alternative for root treepath
      componentName = name_ + "->output";
      format = Parameters::get<std::vector<std::string>>(componentName + "->format");
    }
273

Praetorius, Simon's avatar
Praetorius, Simon committed
274
    if (!format)
275
276
      return;

Praetorius, Simon's avatar
Praetorius, Simon committed
277
278
279
280
281
    for (std::string const& type : format.value()) {
      auto writer = creator.create(type, componentName, treePath);
      if (writer)
        filewriter_.push_back(std::move(writer));
    }
282
283
284
285
  });
}


286
// Adds a Dirichlet boundary condition
287
template <class Traits>
288
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
289
void ProblemStat<Traits>::
290
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
291
{
292
  static constexpr bool isValidPredicate = Concepts::Functor<Predicate, bool(WorldVector)>;
293
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
294
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
295

296
297
298
299
300
301
  static constexpr bool isValidTreePath =
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
  static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addDirichletBC!");

  if constexpr (isValidPredicate && isValidTreePath) {
Praetorius, Simon's avatar
Praetorius, Simon committed
302
    auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col));
303
    auto valueGridFct = makeGridFunction(values, this->gridView());
304

Praetorius, Simon's avatar
Praetorius, Simon committed
305
306
    constraints_.push_back(DirichletBC{
      std::move(basis), {predicate}, valueGridFct});
307
  }
308
309
310
311
312
313
314
315
316
}


// Adds a Dirichlet boundary condition
template <class Traits>
  template <class RowTreePath, class ColTreePath, class Values>
void ProblemStat<Traits>::
addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& values)
{
317
318
319
320
321
322
  static constexpr bool isValidTreePath =
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
  static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addDirichletBC!");

  if constexpr (isValidTreePath) {
Praetorius, Simon's avatar
Praetorius, Simon committed
323
    auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath(col));
324
325
    auto valueGridFct = makeGridFunction(values, this->gridView());

Praetorius, Simon's avatar
Praetorius, Simon committed
326
327
    constraints_.push_back(DirichletBC{
      std::move(basis), {*boundaryManager_, id}, valueGridFct});
328
  }
329
330
331
332
333
334
335
336
}


template <class Traits>
void ProblemStat<Traits>::
addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
{
  auto localView = globalBasis_->localView();
337
  auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath());
Praetorius, Simon's avatar
Praetorius, Simon committed
338
339
  constraints_.push_back(makePeriodicBC(
    std::move(basis), {*boundaryManager_, id}, {matrix, vector}));
340
}
341

342

Praetorius, Simon's avatar
Praetorius, Simon committed
343

344
345
template <class Traits>
void ProblemStat<Traits>::
346
solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData)
347
{
348
  Dune::InverseOperatorResult stat;
349

350
  solution_->resize();
351

352
353
  if (createMatrixData)
    linearSolver_->init(systemMatrix_->impl());
354

355
356
357
358
359
360
361
362
363
  // solve the linear system
  linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);

  if (!storeMatrixData)
    linearSolver_->finish();

  info(2, "solution of discrete system needed {} seconds", stat.elapsed);
  if (stat.reduction >= 0.0)
    info(1, "Residual norm: ||b-Ax||/||b-Ax^0|| = {}", stat.reduction);
364

365
366
367
  bool ignoreConverged = false;
  Parameters::get(name_ + "->solver->ignore converged", ignoreConverged);
  test_exit(stat.converged || ignoreConverged, "Could not solver the linear system!");
368
}
369

370

371
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
372
373
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
374
375
376
377
{
  Dune::Timer t;

  Flag markFlag = 0;
378
  for (auto& currentMarker : marker_)
379
    markFlag |= currentMarker.second->markGrid(adaptInfo);
380

381
  msg("markElements needed {} seconds", t.elapsed());
382
383
384
385
386

  return markFlag;
}


387
388
389
390
391
392
template <class Traits>
Flag ProblemStat<Traits>::
globalCoarsen(int n)
{
  Dune::Timer t;
  bool adapted = false;
393
  // TODO(FM): Find a less expensive alternative to the loop adaption
394
  for (int i = 0; i < n; ++i) {
395
    // mark all entities for coarsening
396
397
398
    for (const auto& element : elements(grid_->leafGridView()))
      grid_->mark(-1, element);

399
400
401
402
403
404
405
    bool adaptedInLoop = grid_->preAdapt();
    adaptedInLoop |= grid_->adapt();
    grid_->postAdapt();
    if (!adaptedInLoop)
      break;
    else
      adapted = true;
406
407
408
409
410
411
412
413
414
  }

  msg("globalCoarsen needed {} seconds", t.elapsed());
  return adapted ? MESH_ADAPTED : Flag(0);
}


// grid has globalRefine(int, AdaptDataHandleInterface&)
template <class G>
415
416
using HasGlobalRefineADHI = decltype(
  std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
417
418
419

template <class Traits>
Flag ProblemStat<Traits>::
420
globalRefine(int n)
421
422
{
  Dune::Timer t;
423
424
425
426
  if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
    grid_->globalRefine(n, globalBasis_->globalRefineCallback());
  else
    grid_->globalRefine(n);
427
428

  msg("globalRefine needed {} seconds", t.elapsed());
429
  return n > 0 ? MESH_ADAPTED : Flag(0);
430
431
432
}


433
434
template <class Traits>
Flag ProblemStat<Traits>::
435
adaptGrid(AdaptInfo& /*adaptInfo*/)
436
437
438
{
  Dune::Timer t;

439
440
441
  bool adapted = grid_->preAdapt();
  adapted |= grid_->adapt();
  grid_->postAdapt();
442

443
  msg("adaptGrid needed {} seconds", t.elapsed());
444
  return adapted ? MESH_ADAPTED : Flag(0);
445
446
447
}


448
449
template <class Traits>
void ProblemStat<Traits>::
450
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
451
{
452
  Dune::Timer t;
453
  Dune::Timer t2;
454

Praetorius, Simon's avatar
Praetorius, Simon committed
455
456
457
  // 0. initialize boundary condition and other constraints
  for (auto& bc : constraints_)
    bc.init();
458
459
460
461

  t2.reset();

  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
462
  systemMatrix_->init();
463
  rhs_->init(sizeInfo(*globalBasis_), asmVector);
464
465
466
467
468
469

  // statistic about system size
  if (Environment::mpiSize() > 1)
    msg("{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
  else
    msg("{} local DOFs", rhs_->localSize());
470

471
472
473
  auto localMatOperators = localOperators(systemMatrix_->operators());
  auto localVecOperators = localOperators(rhs_->operators());

474
  // 2. traverse grid and assemble operators on the elements
Praetorius, Simon's avatar
Praetorius, Simon committed
475
  auto localView = globalBasis_->localView();
Praetorius, Simon's avatar
Praetorius, Simon committed
476
  for (auto const& element : elements(gridView(), PartitionSet{})) {
477
    localView.bind(element);
478
479

    if (asmMatrix)
480
      systemMatrix_->assemble(localView, localView, localMatOperators);
481
    if (asmVector)
482
      rhs_->assemble(localView, localVecOperators);
483

484
    localView.unbind();
485
486
487
  }

  // 3. finish matrix insertion and apply dirichlet boundary conditions
488
489
  systemMatrix_->finish();
  rhs_->finish();
490

491
492
493
  info(2,"  assemble operators needed {} seconds", t2.elapsed());
  t2.reset();

494
  solution_->resize(sizeInfo(*globalBasis_));
Praetorius, Simon's avatar
Praetorius, Simon committed
495
496
497
498

  // 4. apply boundary condition and constraints to matrix, solution, and rhs
  for (auto& bc : constraints_)
    bc.apply(*systemMatrix_, *solution_, *rhs_);
499

500
501
  info(2,"  assemble boundary conditions needed {} seconds", t2.elapsed());

502
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
503
  msg("assemble needed {} seconds", t.elapsed());
504
}
505

506

507
508
template <class Traits>
void ProblemStat<Traits>::
509
writeFiles(AdaptInfo& adaptInfo, bool force)
510
{
511
  Dune::Timer t;
512
  for (auto writer : filewriter_)
Praetorius, Simon's avatar
Praetorius, Simon committed
513
    writer->write(adaptInfo, force);
514
  msg("writeFiles needed {} seconds", t.elapsed());
515
}
516

517
} // end namespace AMDiS