ProblemStat.inc.hpp 11.1 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/LocalAssembler.hpp>
14
#include <amdis/common/Loops.hpp>
15

16
17
namespace AMDiS {

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

35
36
37
38
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
39
40
      grid_ = adoptProblem->grid_;
      gridTransfer_ = adoptProblem->gridTransfer_;
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
  gridTransfer_ = std::make_shared<GridTransfer<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
167
168
169
170
171
172
173
174

  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)
{
  constraints_.init(*globalBasis_, *globalBasis_);
}


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

  assert(bool(gridTransfer_));
  gridTransfer_->attach(solution_.get());
  gridTransfer_->attach(rhs_.get());
Praetorius, Simon's avatar
Praetorius, Simon committed
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205

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


206
207
208
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
Praetorius, Simon's avatar
Praetorius, Simon committed
209
  marker_.clear();
210
211
  auto localView = globalBasis_->localView();
  AMDiS::forEachNode_(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
238
  auto localView = globalBasis_->localView();
  forEachNode_(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
269
  auto bc = std::make_unique<BcType>(predicate, valueGridFct);
  constraints_[i][j].push_back(std::move(bc));
270
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
271
}
272

273

274
275
template <class Traits>
void ProblemStat<Traits>::
276
277
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
278
  Dune::Timer t;
279

280
  SolverInfo solverInfo(name_ + "->solver");
281
282
283
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

284
  solution_->compress();
285

286
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
287
288
                      solverInfo);

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

Praetorius, Simon's avatar
Praetorius, Simon committed
292
293
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
294
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
295
          solverInfo.absResidual(), solverInfo.relResidual());
296
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
297
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
298
299
    }
  }
300

301
302
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
Praetorius, Simon's avatar
Praetorius, Simon committed
303
304
305
306
307
308
    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() << " ";
309
    error_exit("Tolerance {} could not be reached!", tol_str.str());
310
311
  }
}
312

313

314
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
315
316
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
317
318
319
320
{
  Dune::Timer t;

  Flag markFlag = 0;
321
  for (auto& currentMarker : marker_)
322
    markFlag |= currentMarker->markGrid(adaptInfo);
323

324
  msg("markElements needed {} seconds", t.elapsed());
325
326
327
328
329

  return markFlag;
}


330
331
332
333
334
335
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

336
337
338
  bool adapted = gridTransfer_->preAdapt();
  adapted |= gridTransfer_->adapt(); // NOTE: |= does not short-circuit and works with bools as ||=
  if (adapted)
Praetorius, Simon's avatar
Praetorius, Simon committed
339
    globalBasis_->update(gridView());
340
  gridTransfer_->postAdapt();
341

342
  msg("adaptGrid needed {} seconds", t.elapsed());
343
  return adapted ? MESH_ADAPTED : Flag(0);
344
345
346
}


347
348
template <class Traits>
void ProblemStat<Traits>::
349
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
350
{
351
  Dune::Timer t;
352

353
354
355
356
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

357
358
359
  auto localView = globalBasis_->localView();
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
360
      for (auto bc : constraints_[rowNode][colNode])
361
        bc->init(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
362
363
364
365
366
367
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

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

    if (asmMatrix)
371
      systemMatrix_->assemble(localView, localView);
372
    if (asmVector)
373
      rhs_->assemble(localView);
374

375
    localView.unbind();
376
377
378
379
380
381
  }

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

382
383
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
384
385
      // finish boundary condition
      for (auto bc : constraints_[rowNode][colNode])
386
        bc->finish(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
387
388
    });
  });
389

390
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
391
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
392
}
393

394

395
396
template <class Traits>
void ProblemStat<Traits>::
397
writeFiles(AdaptInfo& adaptInfo, bool force)
398
{
399
  Dune::Timer t;
400
  for (auto writer : filewriter_)
401
    writer->writeFiles(adaptInfo, force);
402
  msg("writeFiles needed {} seconds", t.elapsed());
403
}
404

405

406
} // end namespace AMDiS