ProblemStat.inc.hpp 15.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/Assembler.hpp>
16
#include <amdis/GridFunctionOperator.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
17
#include <amdis/io/FileWriterCreator.hpp>
18
#include <amdis/linearalgebra/SymmetryStructure.hpp>
19

20
21
namespace AMDiS {

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

36
37
38
39
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        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
50
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
51
      createGlobalBasis();
52
    }
53

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

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

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

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

74

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

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

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

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

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

97

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

102
  solution_->resizeZero();
103
}
104

105

Praetorius, Simon's avatar
Praetorius, Simon committed
106
107
108
109
110
111
112
113
114
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
115
116
117
  // TODO(SP): implement BAckupRestore independent of wrapped grid
  using HostGrid = typename Grid::HostGrid;

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

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

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

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

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

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

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


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

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

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

161
162
163
164
165
166
167
  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
168
169
170
171
172
173
174
175
176
177
178
}


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>{});
179
  initGlobalBasis();
Praetorius, Simon's avatar
Praetorius, Simon committed
180
181
182
183
184
185
186
}


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


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


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


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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
215
  auto localView = globalBasis_->localView();
216
  for_each_node(localView.tree(), [&,this](auto&&, auto treePath) -> void
Praetorius, Simon's avatar
Praetorius, Simon committed
217
218
219
220
221
222
223
224
225
226
227
228
229
  {
    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";
230
  Parameters::get(name_ + "->solver", solverName);
Praetorius, Simon's avatar
Praetorius, Simon committed
231
232

  auto solverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
233
    = named(CreatorMap<LinearSolver>::getCreator(solverName, name_ + "->solver"));
Praetorius, Simon's avatar
Praetorius, Simon committed
234

235
  linearSolver_ = solverCreator->createWithString(name_ + "->solver");
Praetorius, Simon's avatar
Praetorius, Simon committed
236
237
238
}


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

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

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


260
261
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
262
{
Praetorius, Simon's avatar
Praetorius, Simon committed
263
  FileWriterCreator<SolutionVector> creator(solution_, boundaryManager_);
Praetorius, Simon's avatar
Praetorius, Simon committed
264

Praetorius, Simon's avatar
Praetorius, Simon committed
265
  filewriter_.clear();
266
  auto localView = globalBasis_->localView();
Praetorius, Simon's avatar
Praetorius, Simon committed
267
  for_each_node(localView.tree(), [&](auto const& /*node*/, auto treePath) -> void
268
  {
269
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
Praetorius, Simon's avatar
Praetorius, Simon committed
270
271
272
273
274
275
276
    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");
    }
277

Praetorius, Simon's avatar
Praetorius, Simon committed
278
    if (!format)
279
280
      return;

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


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

300
301
302
303
304
305
306
307
308
309
310
  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) {
    auto localView = globalBasis_->localView();
    auto i = makeTreePath(row);
    auto j = makeTreePath(col);
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, i);
    auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, j);
311

312
    auto valueGridFct = makeGridFunction(values, this->gridView());
313

314
315
316
317
    auto bc = makeDirichletBC<SystemMatrix, SolutionVector, SystemVector>(
                std::move(rowBasis), std::move(colBasis), {predicate}, valueGridFct);
    boundaryConditions_[i][j].push_back(makeUniquePtr(std::move(bc)));
  }
318
319
320
321
322
323
324
325
326
}


// 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)
{
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
  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) {
    auto localView = globalBasis_->localView();
    auto i = makeTreePath(row);
    auto j = makeTreePath(col);
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, i);
    auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, j);

    auto valueGridFct = makeGridFunction(values, this->gridView());

    auto bc = makeDirichletBC<SystemMatrix, SolutionVector, SystemVector>(
                std::move(rowBasis), std::move(colBasis), {*boundaryManager_, id}, valueGridFct);
    boundaryConditions_[i][j].push_back(makeUniquePtr(std::move(bc)));
  }
345
346
347
348
349
350
351
352
}


template <class Traits>
void ProblemStat<Traits>::
addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
{
  auto localView = globalBasis_->localView();
353
  auto basis = Dune::Functions::subspaceBasis(*globalBasis_, makeTreePath());
354
355
  auto bc = makePeriodicBC<SystemMatrix, SolutionVector, SystemVector>(
              std::move(basis), {*boundaryManager_, id}, {matrix, vector});
356
  boundaryConditions_[makeTreePath()][makeTreePath()].push_back(makeUniquePtr(std::move(bc)));
357
}
358

359

360
361
template <class Traits>
void ProblemStat<Traits>::
362
solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData)
363
{
364
  Dune::Timer t;
365

366
  SolverInfo solverInfo(name_ + "->solver");
367
368
369
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

370
  solution_->resize();
Praetorius, Simon's avatar
Praetorius, Simon committed
371
  linearSolver_->solve(*systemMatrix_, *solution_, *rhs_, solverInfo);
372

Praetorius, Simon's avatar
Praetorius, Simon committed
373
  if (solverInfo.info() > 0) {
374
    msg("solution of discrete system needed {} seconds", t.elapsed());
375

Praetorius, Simon's avatar
Praetorius, Simon committed
376
377
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
378
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
379
          solverInfo.absResidual(), solverInfo.relResidual());
380
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
381
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
382
383
    }
  }
384

385
  test_exit(!solverInfo.doBreak() && !solverInfo.error(), "Could not solver the linear system!");
386
}
387

388

389
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
390
391
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
392
393
394
395
{
  Dune::Timer t;

  Flag markFlag = 0;
396
  for (auto& currentMarker : marker_)
397
    markFlag |= currentMarker.second->markGrid(adaptInfo);
398

399
  msg("markElements needed {} seconds", t.elapsed());
400
401
402
403
404

  return markFlag;
}


405
406
407
408
409
410
template <class Traits>
Flag ProblemStat<Traits>::
globalCoarsen(int n)
{
  Dune::Timer t;
  bool adapted = false;
411
  // TODO(FM): Find a less expensive alternative to the loop adaption
412
  for (int i = 0; i < n; ++i) {
413
    // mark all entities for coarsening
414
415
416
    for (const auto& element : elements(grid_->leafGridView()))
      grid_->mark(-1, element);

417
418
419
420
421
422
423
    bool adaptedInLoop = grid_->preAdapt();
    adaptedInLoop |= grid_->adapt();
    grid_->postAdapt();
    if (!adaptedInLoop)
      break;
    else
      adapted = true;
424
425
426
427
428
429
430
431
432
  }

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


// grid has globalRefine(int, AdaptDataHandleInterface&)
template <class G>
433
434
using HasGlobalRefineADHI = decltype(
  std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
435
436
437

template <class Traits>
Flag ProblemStat<Traits>::
438
globalRefine(int n)
439
440
{
  Dune::Timer t;
441
442
443
444
  if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
    grid_->globalRefine(n, globalBasis_->globalRefineCallback());
  else
    grid_->globalRefine(n);
445
446

  msg("globalRefine needed {} seconds", t.elapsed());
447
  return n > 0 ? MESH_ADAPTED : Flag(0);
448
449
450
}


451
452
template <class Traits>
Flag ProblemStat<Traits>::
453
adaptGrid(AdaptInfo& /*adaptInfo*/)
454
455
456
{
  Dune::Timer t;

457
458
459
  bool adapted = grid_->preAdapt();
  adapted |= grid_->adapt();
  grid_->postAdapt();
460

461
  msg("adaptGrid needed {} seconds", t.elapsed());
462
  return adapted ? MESH_ADAPTED : Flag(0);
463
464
465
}


466
467
template <class Traits>
void ProblemStat<Traits>::
468
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
469
{
470
  Dune::Timer t;
471
  Dune::Timer t2;
472

473
  auto localView = globalBasis_->localView();
474
475
476
  for_each_node(localView.tree(), [&](auto&&, auto rowTp) -> void {
    for_each_node(localView.tree(), [&](auto&&, auto colTp) -> void {
      for (auto bc : boundaryConditions_[rowTp][colTp])
477
        bc->init();
478
479
    });
  });
480
481
482
483

  t2.reset();

  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
484
  systemMatrix_->init();
485
  rhs_->init(sizeInfo(*globalBasis_), asmVector);
486
487
488
489
490
491

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

  // 2. traverse grid and assemble operators on the elements
Praetorius, Simon's avatar
Praetorius, Simon committed
494
  for (auto const& element : elements(gridView(), PartitionSet{})) {
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
507
  systemMatrix_->finish();
  rhs_->finish();
508

509
510
511
  info(2,"  assemble operators needed {} seconds", t2.elapsed());
  t2.reset();

512
  solution_->resize(sizeInfo(*globalBasis_));
513
514
  for_each_node(localView.tree(), [&](auto&&, auto rowTp) -> void {
    for_each_node(localView.tree(), [&](auto&&, auto colTp) -> void {
515
      // finish boundary condition
516
      for (auto bc : boundaryConditions_[rowTp][colTp])
517
        bc->apply(*systemMatrix_, *solution_, *rhs_);
518
519
    });
  });
520

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

523
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
524
  msg("assemble needed {} seconds", t.elapsed());
525
}
526

527

528
529
template <class Traits>
void ProblemStat<Traits>::
530
writeFiles(AdaptInfo& adaptInfo, bool force)
531
{
532
  Dune::Timer t;
533
  for (auto writer : filewriter_)
Praetorius, Simon's avatar
Praetorius, Simon committed
534
    writer->write(adaptInfo, force);
535
  msg("writeFiles needed {} seconds", t.elapsed());
536
}
537

538
} // end namespace AMDiS