ProblemStat.inc.hpp 10.8 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
      grid_ = adoptProblem->grid_;
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
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
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_);
174
175
176
  solution_ = std::make_shared<SystemVector>(*globalBasis_, INTERPOLATE);
  rhs_ = std::make_shared<SystemVector>(*globalBasis_, NO_OPERATION);

Praetorius, Simon's avatar
Praetorius, Simon committed
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
  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");
}


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

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

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

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


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

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

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


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

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

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

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

268

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

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

279
  solution_->compress();
280

281
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
282
283
                      solverInfo);

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

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

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

308

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

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

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

  return markFlag;
}


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

331
332
  bool adapted = GridTransferManager::adapt(*grid_);
  globalBasis_->update(gridView());
333

334
  msg("adaptGrid needed {} seconds", t.elapsed());
335
  return adapted ? MESH_ADAPTED : Flag(0);
336
337
338
}


339
340
template <class Traits>
void ProblemStat<Traits>::
341
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
342
{
343
  Dune::Timer t;
344

345
346
347
348
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

349
350
351
  auto localView = globalBasis_->localView();
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
352
      for (auto bc : constraints_[rowNode][colNode])
353
        bc->init(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
354
355
356
357
358
359
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

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

    if (asmMatrix)
363
      systemMatrix_->assemble(localView, localView);
364
    if (asmVector)
365
      rhs_->assemble(localView);
366

367
    localView.unbind();
368
369
370
371
372
373
  }

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

374
375
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
376
377
      // finish boundary condition
      for (auto bc : constraints_[rowNode][colNode])
378
        bc->finish(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
379
380
    });
  });
381

382
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
383
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
384
}
385

386

387
388
template <class Traits>
void ProblemStat<Traits>::
389
writeFiles(AdaptInfo& adaptInfo, bool force)
390
{
391
  Dune::Timer t;
392
  for (auto writer : filewriter_)
393
    writer->writeFiles(adaptInfo, force);
394
  msg("writeFiles needed {} seconds", t.elapsed());
395
}
396

397

398
} // end namespace AMDiS