ProblemStat.inc.hpp 10.2 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
15
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/LocalAssembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/common/Loops.hpp>
16
// #include <amdis/Estimator.hpp>
17

18
19
namespace AMDiS {

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

37
38
39
40
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
41
      grid_ = adoptProblem->grid_;
42
    }
43
  }
44

45
  if (!grid_)
46
    warning("no grid created");
47

48
  int globalRefinements = 0;
49
  Parameters::get(gridName_ + "->global refinements", globalRefinements);
50
  if (globalRefinements > 0) {
51
    grid_->globalRefine(globalRefinements);
52
  }
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))) {
67
68
      globalBasis_ = adoptProblem->globalBasis_;
      initGlobalBasis(*globalBasis_);
69
    }
70
  }
71

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


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

80
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
81
    solution_ = adoptProblem->solution_;
82
    estimates_ = adoptProblem->estimates_;
83
84
    rhs_ = adoptProblem->rhs_;
    systemMatrix_ = adoptProblem->systemMatrix_;
85
  }
86

87

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

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

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

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

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

113

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

118
  solution_->compress();
119
}
120

121

122
123
124
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
125
  AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
126
  {
127
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
128
129
130
131

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

132
133
    std::string tp = to_string(treePath);
    auto newMarker = Marker<Traits>::createMarker(componentName, tp, estimates_[tp], grid_);
134
    if (newMarker)
135
136
      marker_.push_back(std::move(newMarker));

137
138
    // If there is more than one marker, and all components are defined
    // on the same grid, the maximum marking has to be enabled.
139
140
141
    // TODO: needs to be checked!
    if (marker_.size() > 1)
      marker_.back()->setMaximumMarking(true);
142
143
144
145
  });
}


146
147
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
148
{
Praetorius, Simon's avatar
Praetorius, Simon committed
149
  forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
150
  {
151
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
152
153
154
155

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

156
    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
157
    filewriter_.push_back(std::move(writer));
158
159
160
161
  });
}


162
// add matrix/vector operator terms
163

164
template <class Traits>
165
  template <class Operator, class RowTreePath, class ColTreePath>
166
void ProblemStat<Traits>::addMatrixOperator(
167
    Operator const& preOp,
Praetorius, Simon's avatar
Praetorius, Simon committed
168
    RowTreePath row, ColTreePath col)
169
{
170
  static_assert( Concepts::PreTreePath<RowTreePath>,
171
      "row must be a valid treepath, or an integer/index-constant");
172
  static_assert( Concepts::PreTreePath<ColTreePath>,
173
174
      "col must be a valid treepath, or an integer/index-constant");

175
176
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
177

178
  auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
179
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
180

Praetorius, Simon's avatar
Praetorius, Simon committed
181
  matrixOperators_[i][j].element.push_back({localAssembler, nullptr, nullptr});
182
  matrixOperators_[i][j].changing = true;
183
}
184

185

186
template <class Traits>
187
  template <class Operator, class RowTreePath, class ColTreePath>
188
void ProblemStat<Traits>::addMatrixOperator(
189
    BoundaryType b,
190
    Operator const& preOp,
Praetorius, Simon's avatar
Praetorius, Simon committed
191
    RowTreePath row, ColTreePath col)
192
{
193
  static_assert( Concepts::PreTreePath<RowTreePath>,
194
      "row must be a valid treepath, or an integer/index-constant");
195
  static_assert( Concepts::PreTreePath<ColTreePath>,
196
197
      "col must be a valid treepath, or an integer/index-constant");

198
199
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
200
  using Intersection = typename GridView::Intersection;
201

202
  auto op = makeLocalOperator<Intersection>(preOp, globalBasis_->gridView());
203
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
204

Praetorius, Simon's avatar
Praetorius, Simon committed
205
  matrixOperators_[i][j].boundary.push_back({localAssembler, nullptr, nullptr, b});
206
  matrixOperators_[i][j].changing = true;
207
}
208
209


210
template <class Traits>
211
  template <class Operator, class TreePath>
212
void ProblemStat<Traits>::addVectorOperator(
213
    Operator const& preOp,
Praetorius, Simon's avatar
Praetorius, Simon committed
214
    TreePath path)
215
{
216
  static_assert( Concepts::PreTreePath<TreePath>,
217
218
      "path must be a valid treepath, or an integer/index-constant");

219
  auto i = child(localView_->tree(), makeTreePath(path));
220

221
  auto op = makeLocalOperator<Element>(preOp, globalBasis_->gridView());
222
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
223

Praetorius, Simon's avatar
Praetorius, Simon committed
224
  rhsOperators_[i].element.push_back({localAssembler, nullptr, nullptr});
225
  rhsOperators_[i].changing = true;
226
}
227

228

229
template <class Traits>
230
  template <class Operator, class TreePath>
231
void ProblemStat<Traits>::addVectorOperator(
232
    BoundaryType b,
233
    Operator const& preOp,
Praetorius, Simon's avatar
Praetorius, Simon committed
234
    TreePath path)
235
{
236
  static_assert( Concepts::PreTreePath<TreePath>,
237
238
      "path must be a valid treepath, or an integer/index-constant");

239
  auto i = child(localView_->tree(), makeTreePath(path));
240
  using Intersection = typename GridView::Intersection;
241

242
  auto op = makeLocalOperator<Intersection>(preOp, globalBasis_->gridView());
243
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
244

Praetorius, Simon's avatar
Praetorius, Simon committed
245
  rhsOperators_[i].boundary.push_back({localAssembler, nullptr, nullptr, b});
246
  rhsOperators_[i].changing = true;
247
}
248
249


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

259
260
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
261

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

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

271

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

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

282
  solution_->compress();
283

284
  linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
285
286
287
288
289
290
291
292
293
294
295
296
297
                      solverInfo);

  if (solverInfo.getInfo() > 0) {
    msg("solution of discrete system needed ", t.elapsed(), " seconds");

    if (solverInfo.getAbsResidual() >= 0.0) {
      if (solverInfo.getRelResidual() >= 0.0)
        msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual(),
                          ", ||b-Ax||/||b|| = ", solverInfo.getRelResidual());
      else
        msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual());
    }
  }
298

299
300
301
302
303
304
305
306
307
308
309
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
    if (solverInfo.getAbsTolerance() > 0
        && solverInfo.getAbsResidual() > solverInfo.getAbsTolerance())
      tol_str << "absTol = " << solverInfo.getAbsTolerance() << " ";
    if (solverInfo.getRelTolerance() > 0
        && solverInfo.getRelResidual() > solverInfo.getRelTolerance())
      tol_str << "relTol = " << solverInfo.getRelTolerance() << " ";
    error_exit("Tolerance ", tol_str.str(), " could not be reached!");
  }
}
310

311

312
313
314
315
316
317
template <class Traits>
Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  Flag markFlag = 0;
318
  for (auto& currentMarker : marker_)
319
    markFlag |= currentMarker->markGrid(adaptInfo);
320
321
322
323
324
325
326

  msg("markElements needed ", t.elapsed(), " seconds");

  return markFlag;
}


327
328
template <class Traits>
void ProblemStat<Traits>::
329
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
330
{
331
  Dune::Timer t;
332

333
  Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
334
  auto gv = grid_->leafGridView();
335
336
  globalBasis_->update(gv);
  assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
337

338
339
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
340

341

342
343
template <class Traits>
void ProblemStat<Traits>::
344
writeFiles(AdaptInfo& adaptInfo, bool force)
345
{
346
  Dune::Timer t;
347
  for (auto writer : filewriter_)
348
    writer->writeFiles(adaptInfo, force);
349
350
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
351

352

353
} // end namespace AMDiS