ProblemStat.inc.hpp 12.4 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
12
#include <amdis/AdaptInfo.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/GridFunctionOperator.hpp>
13
#include <amdis/GridTransferManager.hpp>
14
#include <amdis/LocalAssembler.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
199
200
201
202
  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"));

  linearSolver_ = solverCreator->create(name_ + "->solver");
}


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
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

363
364
  bool adapted = GridTransferManager::adapt(*grid_);
  globalBasis_->update(gridView());
365

366
  msg("adaptGrid needed {} seconds", t.elapsed());
367
  return adapted ? MESH_ADAPTED : Flag(0);
368
369
370
}


371
372
template <class Traits>
void ProblemStat<Traits>::
373
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
374
{
375
  Dune::Timer t;
376

377
378
379
380
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

381
  auto localView = globalBasis_->localView();
382
383
384
385
386
387
388
389
  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);
390
391
392
393
394
395
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

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

    if (asmMatrix)
399
      systemMatrix_->assemble(localView, localView);
400
    if (asmVector)
401
      rhs_->assemble(localView);
402

403
    localView.unbind();
404
405
406
407
408
409
  }

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

410
411
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
412
      // finish boundary condition
413
414
415
416
      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);
417
418
    });
  });
419

420
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
421
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
422
}
423

424

425
426
template <class Traits>
void ProblemStat<Traits>::
427
writeFiles(AdaptInfo& adaptInfo, bool force)
428
{
429
  Dune::Timer t;
430
  for (auto writer : filewriter_)
431
    writer->writeFiles(adaptInfo, force);
432
  msg("writeFiles needed {} seconds", t.elapsed());
433
}
434

435

436
} // end namespace AMDiS