ProblemStat.inc.hpp 16.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
    warning("grid already created");
31
32
33
34
35
  }
  else {
    if (initFlag.isSet(CREATE_MESH) ||
        (!adoptFlag.isSet(INIT_MESH) &&
        (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
36
      createGrid();
37
    }
38

39
40
41
42
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
43
      adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
44
    }
45
  }
46

47
  if (!grid_)
48
    warning("no grid created");
49

50
51
52
  if (initFlag.isSet(INIT_MESH)) {
    int globalRefinements = 0;
    Parameters::get(gridName_ + "->global refinements", globalRefinements);
53
    if (globalRefinements > 0)
54
      grid_->globalRefine(globalRefinements);
55
56
57
58

    bool loadBalance = false;
    Parameters::get(gridName_ + "->load balance", loadBalance);
    if (loadBalance)
59
      loadBalance = grid_->loadBalance();
60

61
    if (globalBasis_ && (globalRefinements > 0 || loadBalance))
62
      globalBasis_->update(globalBasis_->gridView());
63
  }
64

65
  // create fespace
66
  if (globalBasis_) {
67
    warning("globalBasis already created");
68
69
70
71
  }
  else {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
72
      createGlobalBasis();
73
    }
74

75
76
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Praetorius, Simon's avatar
Praetorius, Simon committed
77
      adoptGlobalBasis(adoptProblem->globalBasis_);
78
    }
79
  }
80

81
  if (!globalBasis_)
82
    warning("no globalBasis created\n");
83

84
85
86
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
87

88
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
89
    systemMatrix_ = adoptProblem->systemMatrix_;
90
91
    solution_ = adoptProblem->solution_;
    rhs_ = adoptProblem->rhs_;
Praetorius, Simon's avatar
Praetorius, Simon committed
92
    estimates_ = adoptProblem->estimates_;
93
  }
94

95

96
  // create solver
97
  if (linearSolver_) {
98
99
100
101
102
    warning("solver already created\n");
  }
  else {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
103

104
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
105
106
      test_exit(!linearSolver_, "solver already created\n");
      linearSolver_ = adoptProblem->linearSolver_;
107
    }
108
109
  }

110
  if (!linearSolver_) {
111
    warning("no solver created\n");
112
113
  }

114
115
116
117
118
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
119
      marker_ = adoptProblem->marker_;
120

121

122
123
124
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
125

126
  solution_->resizeZero();
127
}
128

129

Praetorius, Simon's avatar
Praetorius, Simon committed
130
131
132
133
134
135
136
137
138
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
139
140
141
  // TODO(SP): implement BAckupRestore independent of wrapped grid
  using HostGrid = typename Grid::HostGrid;

Praetorius, Simon's avatar
Praetorius, Simon committed
142
  // restore grid from file
143
  if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
Praetorius, Simon's avatar
Praetorius, Simon committed
144
    adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
Praetorius, Simon's avatar
Praetorius, Simon committed
145
  else
Praetorius, Simon's avatar
Praetorius, Simon committed
146
    adoptGrid(std::shared_ptr<HostGrid>(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename)));
Praetorius, Simon's avatar
Praetorius, Simon committed
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174

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

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

  // create solver
  if (linearSolver_)
    warning("solver already created\n");
  else if (initFlag.isSet(INIT_SOLVER))
    createSolver();

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

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

  solution_->resize();
  solution_->restore(solution_filename);
}


Praetorius, Simon's avatar
Praetorius, Simon committed
175
176
177
178
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
Praetorius, Simon's avatar
Praetorius, Simon committed
179
180

  MeshCreator<Grid> creator(gridName_);
181
  grid_ = creator.create();
182
183

  Dune::Timer t;
184
185
186
  bool loadBalance = grid_->loadBalance();
  if (loadBalance)
    info(2,"load balance needed {} seconds", t.elapsed());
187

188
  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
189
190
  if (!creator.boundaryIds().empty())
    boundaryManager_->setBoundaryIds(creator.boundaryIds());
Praetorius, Simon's avatar
Praetorius, Simon committed
191

192
193
194
195
196
197
198
  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
199
200
201
202
203
204
205
206
207
208
209
}


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>{});
210
  initGlobalBasis();
Praetorius, Simon's avatar
Praetorius, Simon committed
211
212
213
214
215
216
217
}


template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type)
{
  assert( bool(grid_) );
Praetorius, Simon's avatar
Praetorius, Simon committed
218
  static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
219
220
  auto basis = Traits::create(name_, grid_->leafGridView());
  globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis));
Praetorius, Simon's avatar
Praetorius, Simon committed
221
222
223
224
225
226
227
228
229
230
231
}


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


template <class Traits>
232
void ProblemStat<Traits>::initGlobalBasis()
Praetorius, Simon's avatar
Praetorius, Simon committed
233
{
234
235
  dirichletBCs_.init(*globalBasis_, *globalBasis_);
  periodicBCs_.init(*globalBasis_, *globalBasis_);
Praetorius, Simon's avatar
Praetorius, Simon committed
236
237
238
239
240
241
}


template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
242
243
244
  systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
  solution_ = std::make_shared<CoefficientVector>(globalBasis_);
  rhs_ = std::make_shared<SystemVector>(globalBasis_);
245

Praetorius, Simon's avatar
Praetorius, Simon committed
246
  auto localView = globalBasis_->localView();
247
  for_each_node(localView.tree(), [&,this](auto const& node, auto treePath) -> void
Praetorius, Simon's avatar
Praetorius, Simon committed
248
249
250
251
252
253
254
255
256
257
258
259
260
  {
    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";
261
  Parameters::get(name_ + "->solver", solverName);
Praetorius, Simon's avatar
Praetorius, Simon committed
262
263

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

266
  linearSolver_ = solverCreator->createWithString(name_ + "->solver");
Praetorius, Simon's avatar
Praetorius, Simon committed
267
268
269
}


270
271
272
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
Praetorius, Simon's avatar
Praetorius, Simon committed
273
  marker_.clear();
274
  auto localView = globalBasis_->localView();
275
  for_each_node(localView.tree(), [&,this](auto const& node, auto treePath) -> void
276
  {
277
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
278
279
280
281

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

282
    std::string tp = to_string(treePath);
Praetorius, Simon's avatar
Praetorius, Simon committed
283
284
    auto newMarker
      = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
285
    assert(bool(newMarker));
286
    this->addMarker(std::move(newMarker));
287
288
289
290
  });
}


291
292
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
293
{
Praetorius, Simon's avatar
Praetorius, Simon committed
294
295
  FileWriterCreator<CoefficientVector> creator(solution_, boundaryManager_);

Praetorius, Simon's avatar
Praetorius, Simon committed
296
  filewriter_.clear();
297
  auto localView = globalBasis_->localView();
Praetorius, Simon's avatar
Praetorius, Simon committed
298
  for_each_node(localView.tree(), [&](auto const& /*node*/, auto treePath) -> void
299
  {
300
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
Praetorius, Simon's avatar
Praetorius, Simon committed
301
302
303
304
305
306
307
    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");
    }
308

Praetorius, Simon's avatar
Praetorius, Simon committed
309
    if (!format)
310
311
      return;

Praetorius, Simon's avatar
Praetorius, Simon committed
312
313
314
315
316
    for (std::string const& type : format.value()) {
      auto writer = creator.create(type, componentName, treePath);
      if (writer)
        filewriter_.push_back(std::move(writer));
    }
317
318
319
320
  });
}


321
// Adds a Dirichlet boundary condition
322
template <class Traits>
323
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
324
void ProblemStat<Traits>::
325
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
326
327
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
328
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
329

330
331
332
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));
333

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

336
  using Range = RangeType_t<decltype(i)>;
337
  using BcType = DirichletBC<WorldVector,Range>;
338
  auto bc = std::make_unique<BcType>(predicate, valueGridFct);
339
340
341
342
343
344
345
346
347
348
349
350
351
352
  dirichletBCs_[i][j].push_back(std::move(bc));
}


// 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)
{
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));

353
  auto valueGridFct = makeGridFunction(values, this->gridView());
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370

  using Range = RangeType_t<decltype(i)>;
  using BcType = DirichletBC<WorldVector,Range>;
  auto bc = std::make_unique<BcType>(boundaryManager_, id, valueGridFct);
  dirichletBCs_[i][j].push_back(std::move(bc));
}


template <class Traits>
void ProblemStat<Traits>::
addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
{
  auto localView = globalBasis_->localView();
  using BcType = PeriodicBC<WorldVector, typename GlobalBasis::MultiIndex>;
  using FaceTrafo = FaceTransformation<typename WorldVector::field_type, WorldVector::dimension>;
  auto bc = std::make_unique<BcType>(boundaryManager_, id, FaceTrafo{matrix, vector});
  periodicBCs_[localView.tree()][localView.tree()].push_back(std::move(bc));
371
}
372

373

374
375
template <class Traits>
void ProblemStat<Traits>::
376
377
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
378
  Dune::Timer t;
379

380
  SolverInfo solverInfo(name_ + "->solver");
381
382
383
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

384
  solution_->resize();
Praetorius, Simon's avatar
Praetorius, Simon committed
385
386
  linearSolver_->solve(systemMatrix_->impl().matrix(), solution_->impl().vector(),
                       rhs_->impl().vector(), globalBasis_->communicator(), solverInfo);
387

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

Praetorius, Simon's avatar
Praetorius, Simon committed
391
392
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
393
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
394
          solverInfo.absResidual(), solverInfo.relResidual());
395
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
396
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
397
398
    }
  }
399

400
401
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
Praetorius, Simon's avatar
Praetorius, Simon committed
402
403
404
405
406
407
    if (solverInfo.absTolerance() > 0
        && solverInfo.absResidual() > solverInfo.absTolerance())
      tol_str << "absTol = " << solverInfo.absTolerance() << " ";
    if (solverInfo.relTolerance() > 0
        && solverInfo.relResidual() > solverInfo.relTolerance())
      tol_str << "relTol = " << solverInfo.relTolerance() << " ";
408
    error_exit("Tolerance {} could not be reached!", tol_str.str());
409
410
  }
}
411

412

413
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
414
415
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
416
417
418
419
{
  Dune::Timer t;

  Flag markFlag = 0;
420
  for (auto& currentMarker : marker_)
421
    markFlag |= currentMarker.second->markGrid(adaptInfo);
422

423
  msg("markElements needed {} seconds", t.elapsed());
424
425
426
427
428

  return markFlag;
}


429
430
431
432
433
434
template <class Traits>
Flag ProblemStat<Traits>::
globalCoarsen(int n)
{
  Dune::Timer t;
  bool adapted = false;
435
  // TODO(FM): Find a less expensive alternative to the loop adaption
436
  for (int i = 0; i < n; ++i) {
437
    // mark all entities for coarsening
438
439
440
    for (const auto& element : elements(grid_->leafGridView()))
      grid_->mark(-1, element);

441
442
443
444
445
446
447
    bool adaptedInLoop = grid_->preAdapt();
    adaptedInLoop |= grid_->adapt();
    grid_->postAdapt();
    if (!adaptedInLoop)
      break;
    else
      adapted = true;
448
449
450
451
452
453
454
455
456
  }

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


// grid has globalRefine(int, AdaptDataHandleInterface&)
template <class G>
457
458
using HasGlobalRefineADHI = decltype(
  std::declval<G>().globalRefine(1,std::declval<typename G::ADHI&>()));
459
460
461

template <class Traits>
Flag ProblemStat<Traits>::
462
globalRefine(int n)
463
464
465
{
  Dune::Timer t;
  Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalRefineADHI, Grid>{},
466
  /*then*/ [&](auto id) -> void {
467
    id(grid_)->globalRefine(n, globalBasis_->globalRefineCallback());
468
  },
469
  /*else*/ [&](auto id) -> void {
470
    id(grid_)->globalRefine(n);
471
472
473
  });

  msg("globalRefine needed {} seconds", t.elapsed());
474
  return n > 0 ? MESH_ADAPTED : Flag(0);
475
476
477
}


478
479
480
481
482
483
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

484
485
486
  bool adapted = grid_->preAdapt();
  adapted |= grid_->adapt();
  grid_->postAdapt();
487

488
  msg("adaptGrid needed {} seconds", t.elapsed());
489
  return adapted ? MESH_ADAPTED : Flag(0);
490
491
492
}


493
494
template <class Traits>
void ProblemStat<Traits>::
495
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
496
{
497
  Dune::Timer t;
498
  Dune::Timer t2;
499

500
  auto localView = globalBasis_->localView();
501
  for_each_node(localView.tree(), [&,this](auto const& rowNode, auto rowTp) -> void {
502
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTp);
503
    for_each_node(localView.tree(), [&,this](auto const& colNode, auto colTp) -> void {
504
505
506
507
508
      auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTp);
      for (auto bc : dirichletBCs_[rowNode][colNode])
        bc->init(rowBasis, colBasis);
      for (auto bc : periodicBCs_[rowNode][colNode])
        bc->init(rowBasis, colBasis);
509
510
    });
  });
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525

  t2.reset();

  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  std::string symmetryStr = "unknown";
  Parameters::get(name_ + "->symmetry", symmetryStr);

  systemMatrix_->init(symmetryStructure(symmetryStr));
  rhs_->init(asmVector);

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

  // 2. traverse grid and assemble operators on the elements
528
  for (auto const& element : elements(gridView(), typename LinearSolverTraits::PartitionSet{})) {
529
    localView.bind(element);
530
531

    if (asmMatrix)
532
      systemMatrix_->assemble(localView, localView);
533
    if (asmVector)
534
      rhs_->assemble(localView);
535

536
    localView.unbind();
537
538
539
  }

  // 3. finish matrix insertion and apply dirichlet boundary conditions
540
541
  systemMatrix_->finish();
  rhs_->finish();
542

543
544
545
546
  info(2,"  assemble operators needed {} seconds", t2.elapsed());
  t2.reset();

  solution_->resize();
547
548
  for_each_node(localView.tree(), [&,this](auto const& rowNode, auto row_tp) -> void {
    for_each_node(localView.tree(), [&,this](auto const& colNode, auto col_tp) -> void {
549
      // finish boundary condition
550
      for (auto bc : dirichletBCs_[rowNode][colNode])
551
        bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, row_tp, colNode, col_tp);
552
      for (auto bc : periodicBCs_[rowNode][colNode])
553
        bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, row_tp, colNode, col_tp);
554
555
    });
  });
556

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

559
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
560
  msg("assemble needed {} seconds", t.elapsed());
561
}
562

563

564
565
template <class Traits>
void ProblemStat<Traits>::
566
writeFiles(AdaptInfo& adaptInfo, bool force)
567
{
568
  Dune::Timer t;
569
  for (auto writer : filewriter_)
Praetorius, Simon's avatar
Praetorius, Simon committed
570
    writer->write(adaptInfo, force);
571
  msg("writeFiles needed {} seconds", t.elapsed());
572
}
573

574
} // end namespace AMDiS