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

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

7
#include <dune/common/timer.hh>
8
9
#include <dune/typetree/childextraction.hh>

10
11
#include <amdis/AdaptInfo.hpp>
#include <amdis/FileWriter.hpp>
12
#include <amdis/Assembler.hpp>
13
#include <amdis/GridFunctionOperator.hpp>
14
#include <amdis/GridTransferManager.hpp>
15
#include <amdis/common/Loops.hpp>
16

17
18
namespace AMDiS {

19
20
template <class Traits>
void ProblemStat<Traits>::initialize(
21
22
23
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
24
{
25
  // create grids
26
  if (grid_) {
27
    warning("grid already created");
28
29
30
31
32
  }
  else {
    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
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

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

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

71
  if (!globalBasis_)
72
    warning("no globalBasis created\n");
73
74


75
76
77
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
78

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

86

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

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

101
  if (!linearSolver_) {
102
    warning("no solver created\n");
103
104
  }

105
106
107
108
109
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
110
      marker_ = adoptProblem->marker_;
111

112

113
114
115
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
116

117
  solution_->compress();
118
}
119

120

Praetorius, Simon's avatar
Praetorius, Simon committed
121
122
123
124
125
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
  grid_ = MeshCreator<Grid>::create(gridName_);
126
  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
Praetorius, Simon's avatar
Praetorius, Simon committed
127
128
129
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

  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)
{
167
168
  dirichletBCs_.init(*globalBasis_, *globalBasis_);
  periodicBCs_.init(*globalBasis_, *globalBasis_);
Praetorius, Simon's avatar
Praetorius, Simon committed
169
170
171
172
173
174
175
}


template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
  systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
176
177
178
  solution_ = std::make_shared<SystemVector>(*globalBasis_, INTERPOLATE);
  rhs_ = std::make_shared<SystemVector>(*globalBasis_, NO_OPERATION);

Praetorius, Simon's avatar
Praetorius, Simon committed
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
  auto localView = globalBasis_->localView();
  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
  {
    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"));

199
  linearSolver_ = solverCreator->createWithString(name_ + "->solver");
Praetorius, Simon's avatar
Praetorius, Simon committed
200
201
202
}


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

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

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

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


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

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

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


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

257
258
259
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));
260

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

263
  using Range = RangeType_t<decltype(i)>;
264
  using BcType = DirichletBC<WorldVector,Range>;
265
  auto bc = std::make_unique<BcType>(predicate, valueGridFct);
266
267
268
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
  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));
298
}
299

300

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

307
  SolverInfo solverInfo(name_ + "->solver");
308
309
310
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

311
  solution_->compress();
312

313
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
314
315
                      solverInfo);

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

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

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

340

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

  Flag markFlag = 0;
348
  for (auto& currentMarker : marker_)
349
    markFlag |= currentMarker->markGrid(adaptInfo);
350

351
  msg("markElements needed {} seconds", t.elapsed());
352
353
354
355
356

  return markFlag;
}


357
358
359
360
361
362
363
364
365
366
367
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);

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

  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) {
388
    // TODO(FM): Add a way to pass a GridTransfer as ADH with *globalBasis_ argument
389
390
391
392
393
394
395
396
397
    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);

398
      adapted |= GridTransferManager::adapt(*id(grid_), *id(globalBasis_));
399
400
401
402
403
404
405
406
    }
  });

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


407
408
409
410
411
412
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

413
  bool adapted = GridTransferManager::adapt(*grid_, *globalBasis_);
414

415
  msg("adaptGrid needed {} seconds", t.elapsed());
416
  return adapted ? MESH_ADAPTED : Flag(0);
417
418
419
}


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

426
427
428
429
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

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

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

    if (asmMatrix)
448
      systemMatrix_->assemble(localView, localView);
449
    if (asmVector)
450
      rhs_->assemble(localView);
451

452
    localView.unbind();
453
454
455
456
457
458
  }

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

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

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

473

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

484

485
} // end namespace AMDiS