ProblemStat.inc.hpp 10.3 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

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
  int globalRefinements = 0;
48
  Parameters::get(gridName_ + "->global refinements", globalRefinements);
49
  if (globalRefinements > 0) {
50
    grid_->globalRefine(globalRefinements);
51
  }
52
53


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

64
65
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
66
67
      globalBasis_ = adoptProblem->globalBasis_;
      initGlobalBasis(*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)) {
80
    solution_ = adoptProblem->solution_;
81
    estimates_ = adoptProblem->estimates_;
82
83
    rhs_ = adoptProblem->rhs_;
    systemMatrix_ = adoptProblem->systemMatrix_;
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

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

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

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

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


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

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

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


161
// add matrix/vector operator terms
162

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

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

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

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

184

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

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

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

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


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

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

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

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

227

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

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

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

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


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

258
259
  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_shared<BcType>(predicate, valueGridFct);
266
  constraints_[i][j].push_back(bc);
267
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
268
}
269

270

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

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

281
  solution_->compress();
282

283
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
284
285
286
                      solverInfo);

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

    if (solverInfo.getAbsResidual() >= 0.0) {
      if (solverInfo.getRelResidual() >= 0.0)
291
292
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
          solverInfo.getAbsResidual(), solverInfo.getRelResidual());
293
      else
294
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.getAbsResidual());
295
296
    }
  }
297

298
299
300
301
302
303
304
305
  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() << " ";
306
    error_exit("Tolerance {} could not be reached!", tol_str.str());
307
308
  }
}
309

310

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

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

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

  return markFlag;
}


326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  grid_->preAdapt();
  grid_->adapt();
  globalBasis_->update(gridView());
  grid_->postAdapt();

  msg("Grid adaption needed {} seconds", t.elapsed());
  return 0;
}


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

348
349
  Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
  assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
350

351
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
352
}
353

354

355
356
template <class Traits>
void ProblemStat<Traits>::
357
writeFiles(AdaptInfo& adaptInfo, bool force)
358
{
359
  Dune::Timer t;
360
  for (auto writer : filewriter_)
361
    writer->writeFiles(adaptInfo, force);
362
  msg("writeFiles needed {} seconds", t.elapsed());
363
}
364

365

366
} // end namespace AMDiS