ProblemStat.inc.hpp 8.99 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 grides
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
81
82
    solution_ = adoptProblem->solution_;
    rhs_ = adoptProblem->rhs_;
    systemMatrix_ = adoptProblem->systemMatrix_;
83
  }
84

85

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

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

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

104

105
106
107
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
108

109
  solution_->compress();
110
}
111

112

113
114
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
115
{
116
  AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
117
  {
118
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
119
120
121
122

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

123
    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
124
    filewriter_.push_back(std::move(writer));
125
126
127
128
  });
}


129
// add matrix/vector operator terms
130

131
template <class Traits>
132
  template <class Operator, class RowTreePath, class ColTreePath>
133
void ProblemStat<Traits>::addMatrixOperator(
134
    Operator const& preOp,
135
136
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
137
{
138
  static_assert( Concepts::PreTreePath<RowTreePath>,
139
      "row must be a valid treepath, or an integer/index-constant");
140
  static_assert( Concepts::PreTreePath<ColTreePath>,
141
142
      "col must be a valid treepath, or an integer/index-constant");

143
144
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
145

146
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
147
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
148

149
150
  matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
  matrixOperators_[i][j].changing = true;
151
}
152

153

154
template <class Traits>
155
  template <class Operator, class RowTreePath, class ColTreePath>
156
void ProblemStat<Traits>::addMatrixOperator(
157
    BoundaryType b,
158
    Operator const& preOp,
159
160
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
161
{
162
  static_assert( Concepts::PreTreePath<RowTreePath>,
163
      "row must be a valid treepath, or an integer/index-constant");
164
  static_assert( Concepts::PreTreePath<ColTreePath>,
165
166
      "col must be a valid treepath, or an integer/index-constant");

167
168
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
169
  using Intersection = typename GridView::Intersection;
170

171
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
172
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
173

174
175
  matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
  matrixOperators_[i][j].changing = true;
176
}
177
178


179
template <class Traits>
180
  template <class Operator, class TreePath>
181
void ProblemStat<Traits>::addVectorOperator(
182
    Operator const& preOp,
183
184
    TreePath path,
    double* factor, double* estFactor)
185
{
186
  static_assert( Concepts::PreTreePath<TreePath>,
187
188
      "path must be a valid treepath, or an integer/index-constant");

189
  auto i = child(localView_->tree(), makeTreePath(path));
190

191
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
192
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
193

194
195
  rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
  rhsOperators_[i].changing = true;
196
}
197

198

199
template <class Traits>
200
  template <class Operator, class TreePath>
201
void ProblemStat<Traits>::addVectorOperator(
202
    BoundaryType b,
203
    Operator const& preOp,
204
205
    TreePath path,
    double* factor, double* estFactor)
206
{
207
  static_assert( Concepts::PreTreePath<TreePath>,
208
209
      "path must be a valid treepath, or an integer/index-constant");

210
  auto i = child(localView_->tree(), makeTreePath(path));
211
  using Intersection = typename GridView::Intersection;
212

213
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
214
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
215

216
217
  rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
  rhsOperators_[i].changing = true;
218
}
219
220


221
// Adds a Dirichlet boundary condition
222
template <class Traits>
223
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
224
void ProblemStat<Traits>::
225
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
226
227
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
228
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
229

230
231
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
232

233
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
234

235
  using Range = RangeType_t<decltype(i)>;
236
  using BcType = DirichletBC<WorldVector,Range>;
237
  auto bc = std::make_shared<BcType>(predicate, valueGridFct);
238
  constraints_[i][j].push_back(bc);
239
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
240
}
241

242

243
244
template <class Traits>
void ProblemStat<Traits>::
245
246
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
247
  Dune::Timer t;
248

249
  SolverInfo solverInfo(name_ + "->solver");
250
251
252
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

253
  solution_->compress();
254

255
  linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
256
257
258
259
260
261
262
263
264
265
266
267
268
                      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());
    }
  }
269

270
271
272
273
274
275
276
277
278
279
280
  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!");
  }
}
281

282

283
284
template <class Traits>
void ProblemStat<Traits>::
285
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
286
{
287
  Dune::Timer t;
288

289
290
291
292
  Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
  auto gv = leafGridView();
  globalBasis_->update(gv);
  assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
293

294
295
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
296

297

298
299
template <class Traits>
void ProblemStat<Traits>::
300
writeFiles(AdaptInfo& adaptInfo, bool force)
301
{
302
  Dune::Timer t;
303
  for (auto writer : filewriter_)
304
    writer->writeFiles(adaptInfo, force);
305
306
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
307

308
} // end namespace AMDiS