ProblemStat.inc.hpp 16.5 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
10
#include <dune/typetree/childextraction.hh>

11
#include <amdis/AdaptInfo.hpp>
12
#include <amdis/BackupRestore.hpp>
13
#include <amdis/FileWriter.hpp>
14
#include <amdis/Assembler.hpp>
15
#include <amdis/GridFunctionOperator.hpp>
16
#include <amdis/GridTransferManager.hpp>
17

18
19
namespace AMDiS {

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

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

45
  if (!grid_)
46
    warning("no grid created");
47

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

58
  // create fespace
59
  if (globalBasis_) {
60
    warning("globalBasis already created");
61
62
63
64
  }
  else {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
65
      createGlobalBasis();
66
    }
67

68
69
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
Praetorius, Simon's avatar
Praetorius, Simon committed
70
      adoptGlobalBasis(adoptProblem->globalBasis_);
71
    }
72
  }
73

74
  if (!globalBasis_)
75
    warning("no globalBasis created\n");
76
77


78
79
80
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
81

82
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
83
    systemMatrix_ = adoptProblem->systemMatrix_;
84
85
    solution_ = adoptProblem->solution_;
    rhs_ = adoptProblem->rhs_;
Praetorius, Simon's avatar
Praetorius, Simon committed
86
    estimates_ = adoptProblem->estimates_;
87
  }
88

89

90
  // create solver
91
  if (linearSolver_) {
92
93
94
95
96
    warning("solver already created\n");
  }
  else {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
97

98
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
99
100
      test_exit(!linearSolver_, "solver already created\n");
      linearSolver_ = adoptProblem->linearSolver_;
101
    }
102
103
  }

104
  if (!linearSolver_) {
105
    warning("no solver created\n");
106
107
  }

108
109
110
111
112
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
113
      marker_ = adoptProblem->marker_;
114

115

116
117
118
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
119

120
  solution_->compress();
121
}
122

123

Praetorius, Simon's avatar
Praetorius, Simon committed
124
125
126
127
128
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
  grid_ = MeshCreator<Grid>::create(gridName_);
129
  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
Praetorius, Simon's avatar
Praetorius, Simon committed
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169

  msg("Create grid:");
  msg("#elements = {}"   , grid_->size(0));
  msg("#faces/edges = {}", grid_->size(1));
  msg("#vertices = {}"   , grid_->size(dim));
  msg("");
}


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>{});
  initGlobalBasis(*globalBasis_);
}


template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type)
{
  assert( bool(grid_) );
  static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
  globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
}


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


template <class Traits>
void ProblemStat<Traits>::initGlobalBasis(GlobalBasis const& globalBasis)
{
170
171
  dirichletBCs_.init(*globalBasis_, *globalBasis_);
  periodicBCs_.init(*globalBasis_, *globalBasis_);
Praetorius, Simon's avatar
Praetorius, Simon committed
172
173
174
175
176
177
178
}


template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
  systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
179
180
  solution_ = std::make_shared<SystemVector>(*globalBasis_, DataTransferOperation::INTERPOLATE);
  rhs_ = std::make_shared<SystemVector>(*globalBasis_, DataTransferOperation::NO_OPERATION);
181

Praetorius, Simon's avatar
Praetorius, Simon committed
182
  auto localView = globalBasis_->localView();
183
  for_each_node(localView.tree(), [&,this](auto const& node, auto treePath)
Praetorius, Simon's avatar
Praetorius, Simon committed
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
  {
    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";
  Parameters::get(name_ + "->solver->name", solverName);

  auto solverCreator
    = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));

202
  linearSolver_ = solverCreator->createWithString(name_ + "->solver");
Praetorius, Simon's avatar
Praetorius, Simon committed
203
204
205
}


206
207
208
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
Praetorius, Simon's avatar
Praetorius, Simon committed
209
  marker_.clear();
210
  auto localView = globalBasis_->localView();
211
  for_each_node(localView.tree(), [&,this](auto const& node, auto treePath)
212
  {
213
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
214
215
216
217

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

218
    std::string tp = to_string(treePath);
Praetorius, Simon's avatar
Praetorius, Simon committed
219
220
    auto newMarker
      = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
221
    if (newMarker)
222
223
      marker_.push_back(std::move(newMarker));

224
225
    // If there is more than one marker, and all components are defined
    // on the same grid, the maximum marking has to be enabled.
226
227
228
    // TODO: needs to be checked!
    if (marker_.size() > 1)
      marker_.back()->setMaximumMarking(true);
229
230
231
232
  });
}


233
234
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
235
{
Praetorius, Simon's avatar
Praetorius, Simon committed
236
  filewriter_.clear();
237
  auto localView = globalBasis_->localView();
238
  for_each_node(localView.tree(), [&,this](auto const& node, auto treePath)
239
  {
240
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
241
242
243
244

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

Praetorius, Simon's avatar
Praetorius, Simon committed
245
    auto writer = makeFileWriterPtr(componentName, this->solution(treePath));
246
    filewriter_.push_back(std::move(writer));
247
248
249
250
  });
}


251
// Adds a Dirichlet boundary condition
252
template <class Traits>
253
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
254
void ProblemStat<Traits>::
255
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
256
257
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
258
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
259

260
261
262
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));
263

264
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
265

266
  using Range = RangeType_t<decltype(i)>;
267
  using BcType = DirichletBC<WorldVector,Range>;
268
  auto bc = std::make_unique<BcType>(predicate, valueGridFct);
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
  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));

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

  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));
301
}
302

303

304
305
template <class Traits>
void ProblemStat<Traits>::
306
307
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
308
  Dune::Timer t;
309

310
  SolverInfo solverInfo(name_ + "->solver");
311
312
313
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

314
  solution_->compress();
315

316
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
317
318
                      solverInfo);

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

Praetorius, Simon's avatar
Praetorius, Simon committed
322
323
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
324
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
325
          solverInfo.absResidual(), solverInfo.relResidual());
326
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
327
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
328
329
    }
  }
330

331
332
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
Praetorius, Simon's avatar
Praetorius, Simon committed
333
334
335
336
337
338
    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() << " ";
339
    error_exit("Tolerance {} could not be reached!", tol_str.str());
340
341
  }
}
342

343

344
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
345
346
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
347
348
349
350
{
  Dune::Timer t;

  Flag markFlag = 0;
351
  for (auto& currentMarker : marker_)
352
    markFlag |= currentMarker->markGrid(adaptInfo);
353

354
  msg("markElements needed {} seconds", t.elapsed());
355
356
357
358
359

  return markFlag;
}


360
361
362
363
364
365
366
367
368
369
370
template <class Traits>
Flag ProblemStat<Traits>::
globalCoarsen(int n)
{
  Dune::Timer t;
  bool adapted = false;
  for (int i = 0; i < n; ++i) {
    // mark all entities for grid refinement
    for (const auto& element : elements(grid_->leafGridView()))
      grid_->mark(-1, element);

371
    adapted |= GridTransferManager::adapt(*grid_, *globalBasis_);
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
  }

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


// grid has globalRefine(int, AdaptDataHandleInterface&)
template <class G>
using HasGlobalRefineADHI = decltype(std::declval<G>().globalRefine(1,std::declval<GridTransfer<G>&>()));

template <class Traits>
Flag ProblemStat<Traits>::
globalRefine(int refCount)
{
  Dune::Timer t;
  bool adapted = false;
  Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalRefineADHI, Grid>{},
  /*then*/ [&](auto id) {
391
    // TODO(FM): Add a way to pass a GridTransfer as ADH with *globalBasis_ argument
392
393
394
395
396
397
398
399
400
    id(grid_)->globalRefine(refCount, GridTransferManager::gridTransfer(*grid_));
    globalBasis_->update(this->gridView());
  },
  /*else*/ [&](auto id) {
    for (int i = 0; i < refCount; ++i) {
      // mark all entities for grid refinement
      for (const auto& element : elements(grid_->leafGridView()))
        grid_->mark(1, element);

401
      adapted |= GridTransferManager::adapt(*id(grid_), *id(globalBasis_));
402
403
404
405
406
407
408
409
    }
  });

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


410
411
412
413
414
415
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

416
  bool adapted = GridTransferManager::adapt(*grid_, *globalBasis_);
417

418
  msg("adaptGrid needed {} seconds", t.elapsed());
419
  return adapted ? MESH_ADAPTED : Flag(0);
420
421
422
}


423
424
template <class Traits>
void ProblemStat<Traits>::
425
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
426
{
427
  Dune::Timer t;
428

429
430
431
432
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

433
  auto localView = globalBasis_->localView();
434
  for_each_node(localView.tree(), [&,this](auto const& rowNode, auto rowTp) {
435
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTp);
436
    for_each_node(localView.tree(), [&,this](auto const& colNode, auto colTp) {
437
438
439
440
441
      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);
442
443
444
445
446
447
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

  // 2. traverse grid and assemble operators on the elements
  for (auto const& element : elements(gridView())) {
448
    localView.bind(element);
449
450

    if (asmMatrix)
451
      systemMatrix_->assemble(localView, localView);
452
    if (asmVector)
453
      rhs_->assemble(localView);
454

455
    localView.unbind();
456
457
458
459
460
461
  }

  // 3. finish matrix insertion and apply dirichlet boundary conditions
  systemMatrix_->finish(asmMatrix);
  rhs_->finish(asmVector);

462
463
  for_each_node(localView.tree(), [&,this](auto const& rowNode, auto row_tp) {
    for_each_node(localView.tree(), [&,this](auto const& colNode, auto col_tp) {
464
      // finish boundary condition
465
      for (auto bc : dirichletBCs_[rowNode][colNode])
466
        bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, row_tp, colNode, col_tp);
467
      for (auto bc : periodicBCs_[rowNode][colNode])
468
        bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, row_tp, colNode, col_tp);
469
470
    });
  });
471

472
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
473
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
474
}
475

476

477
478
template <class Traits>
void ProblemStat<Traits>::
479
writeFiles(AdaptInfo& adaptInfo, bool force)
480
{
481
  Dune::Timer t;
482
  for (auto writer : filewriter_)
483
    writer->writeFiles(adaptInfo, force);
484
  msg("writeFiles needed {} seconds", t.elapsed());
485
}
486

487

488
489
template <class Traits>
void ProblemStat<Traits>::
490
backup(std::string const& filename) const
491
{
492
493
494
  try {
    Dune::BackupRestoreFacility<Grid>::backup(*grid_, filename);
  } catch (Dune::NotImplemented const&) {
495
    warning("Falling back to backup of gridview.");
496
497
    BackupRestoreByGridFactory<Grid>::backup(gridView(), filename);
  }
498
499
500
501
}

template <class Traits>
void ProblemStat<Traits>::
502
backup(AdaptInfo& adaptInfo)
503
{
504
505
  auto param = Parameters::get<std::string>(name_ + "->backup->grid");
  std::string grid_filename = param ? *param : name_ + "_" + std::to_string(adaptInfo.timestepNumber()) + ".grid";
506

507
508
509
510
511
512
513
514
515
516
517
518
519
  auto param2 = Parameters::get<std::string>(name_ + "->backup->solution");
  std::string solution_filename = param2 ? *param2 : name_ + "_" + std::to_string(adaptInfo.timestepNumber()) + ".solution";

  backup(grid_filename);
  solution_->backup(solution_filename);
  msg("Problem backuped to files '{}' and '{}'.", grid_filename, solution_filename);
}


template <class Traits>
auto ProblemStat<Traits>::
restore(std::string const& filename)
{
520
  // restore grid from file
521
  std::unique_ptr<Grid> grid;
522
523
524
525
526
  try {
    grid.reset(Dune::BackupRestoreFacility<Grid>::restore(filename));
  } catch (Dune::NotImplemented const&) {
    grid.reset(BackupRestoreByGridFactory<Grid>::restore(filename));
  }
527
528
529
530
531
532
533
534
535
536
537
538
539
540
  return std::move(grid);
}

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);

  // restore grid from file
  adoptGrid(restore(grid_filename));
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564

  // 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_->compress();
565
  solution_->restore(solution_filename);
566
567
}

568
} // end namespace AMDiS