ProblemStat.inc.hpp 10.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
12
13
14
#include <amdis/AdaptInfo.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/LocalAssembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
#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))) {
Praetorius, Simon's avatar
Praetorius, Simon committed
39
      adoptGrid(adoptProblem->grid_);
40
    }
41
  }
42

43
  if (!grid_)
44
    warning("no grid created");
45

46
  int globalRefinements = 0;
47
  Parameters::get(gridName_ + "->global refinements", globalRefinements);
48
  if (globalRefinements > 0) {
49
    grid_->globalRefine(globalRefinements);
50
  }
51
52


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

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

69
  if (!globalBasis_)
70
    warning("no globalBasis created\n");
71
72


73
74
75
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
76

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

84

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

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

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

103
104
105
106
107
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
108
      marker_ = adoptProblem->marker_;
109

110

111
112
113
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
114

115
  solution_->compress();
116
}
117

118

Praetorius, Simon's avatar
Praetorius, Simon committed
119
120
121
122
123
124
125
126
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
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
  grid_ = MeshCreator<Grid>::create(gridName_);

  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_);
  solution_ = std::make_shared<SystemVector>(*globalBasis_);
  rhs_ = std::make_shared<SystemVector>(*globalBasis_);

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


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

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

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

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


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

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

Praetorius, Simon's avatar
Praetorius, Simon committed
238
    auto writer = makeFileWriterPtr(componentName, this->solution(treePath));
239
    filewriter_.push_back(std::move(writer));
240
241
242
243
  });
}


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

253
254
255
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));
256

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

259
  using Range = RangeType_t<decltype(i)>;
260
  using BcType = DirichletBC<WorldVector,Range>;
261
  auto bc = std::make_shared<BcType>(predicate, valueGridFct);
262
  constraints_[i][j].push_back(bc);
263
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
264
}
265

266

267
268
template <class Traits>
void ProblemStat<Traits>::
269
270
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
271
  Dune::Timer t;
272

273
  SolverInfo solverInfo(name_ + "->solver");
274
275
276
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

277
  solution_->compress();
278

279
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
280
281
                      solverInfo);

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

Praetorius, Simon's avatar
Praetorius, Simon committed
285
286
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
287
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
288
          solverInfo.absResidual(), solverInfo.relResidual());
289
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
290
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
291
292
    }
  }
293

294
295
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
Praetorius, Simon's avatar
Praetorius, Simon committed
296
297
298
299
300
301
    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() << " ";
302
    error_exit("Tolerance {} could not be reached!", tol_str.str());
303
304
  }
}
305

306

307
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
308
309
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
310
311
312
313
{
  Dune::Timer t;

  Flag markFlag = 0;
314
  for (auto& currentMarker : marker_)
315
    markFlag |= currentMarker->markGrid(adaptInfo);
316

317
  msg("markElements needed {} seconds", t.elapsed());
318
319
320
321
322

  return markFlag;
}


323
324
325
326
327
328
329
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  grid_->preAdapt();
Praetorius, Simon's avatar
Praetorius, Simon committed
330
331
332
333
334
335
336
337
  bool changed = grid_->adapt();

  // update data
  if (changed) {
    globalBasis_->update(gridView());
    solution_->resize(*globalBasis_);
  }

338
339
  grid_->postAdapt();

340
  msg("adaptGrid needed {} seconds", t.elapsed());
Praetorius, Simon's avatar
Praetorius, Simon committed
341
  return changed ? MESH_ADAPTED : Flag(0);
342
343
344
}


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

351
352
353
354
355
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);
  solution_->resize(*globalBasis_);

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

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

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

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

  // 3. finish matrix insertion and apply dirichlet boundary conditions
  systemMatrix_->finish(asmMatrix);
  rhs_->finish(asmVector);
380
  (*solution_) = 0;
381

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