ProblemStat.inc.hpp 9.13 KB
Newer Older
1
2
#pragma once

3
4
5
6
#include <map>
#include <string>
#include <utility>

7
8
#include <dune/typetree/childextraction.hh>

9
#include <dune/amdis/AdaptInfo.hpp>
10
#include <dune/amdis/Assembler.hpp>
11
#include <dune/amdis/FileWriter.hpp>
12
#include <dune/amdis/LocalAssembler.hpp>
13
#include <dune/amdis/GridFunctionOperator.hpp>
14
15
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.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
26
27
  // create grides
  if (grid) {
    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->getGrid();
41
    }
42

43
    componentGrids.resize(nComponents, grid.get());
44
  }
45

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

49
  int globalRefinements = 0;
50
  Parameters::get(gridName + "->global refinements", globalRefinements);
51
  if (globalRefinements > 0) {
52
    grid->globalRefine(globalRefinements);
53
  }
54
55


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

66
67
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
68
      globalBasis = adoptProblem->getGlobalBasis();
69
    }
70
  }
71

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


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

80
81
82
83
84
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    solution = adoptProblem->getSolution();
    rhs = adoptProblem->getRhs();
    systemMatrix = adoptProblem->getSystemMatrix();
  }
85

86

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

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

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

105

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

110
111
  solution->compress();
}
112

113

114
115
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
116
117
118
119
120
121
122
123
124
{
  auto localView = globalBasis->localView();
  forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
  {
    std::string componentName = name + "->output[" + to_string(treePath) + "]";

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

Praetorius, Simon's avatar
Praetorius, Simon committed
125
126
    auto writer = makeFileWriterPtr<Traits>(componentName, this->getSolution(treePath));
    filewriter.push_back(std::move(writer));
127
128
129
130
  });
}


131
// add matrix/vector operator terms
132

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

145
146
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
147

148
  auto op = makeGridOperator(preOp, globalBasis->gridView());
149
  auto localAssembler = makeLocalAssemblerPtr<Element>(op, i, j);
150

151
  matrixOperators[i][j].element.push_back({localAssembler, factor, estFactor});
152
  matrixOperators[i][j].changing = true;
153
}
154

155

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

169
170
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
171
  using Intersection = typename GridView::Intersection;
172

173
  auto op = makeGridOperator(preOp, globalBasis->gridView());
174
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i, j);
175

176
177
  matrixOperators[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
  matrixOperators[i][j].changing = true;
178
}
179
180


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

191
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
192

193
  auto op = makeGridOperator(preOp, globalBasis->gridView());
194
  auto localAssembler = makeLocalAssemblerPtr<Element>(op, i);
195

196
  rhsOperators[i].element.push_back({localAssembler, factor, estFactor});
197
  rhsOperators[i].changing = true;
198
}
199

200

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

212
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
213
  using Intersection = typename GridView::Intersection;
214

215
  auto op = makeGridOperator(preOp, globalBasis->gridView());
216
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i);
217

218
219
  rhsOperators[i].boundary.push_back({localAssembler, factor, estFactor, b});
  rhsOperators[i].changing = true;
220
}
221
222


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

232
233
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
234

235
  auto valueGridFct = makeGridFunction(values, globalBasis->gridView());
236

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

244

245
246
template <class Traits>
void ProblemStat<Traits>::
247
248
249
250
251
252
253
254
255
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
  Timer t;

  SolverInfo solverInfo(name + "->solver");
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

  solution->compress();
256

257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
  linearSolver->solve(systemMatrix->getMatrix(),
                      solution->getVector(), rhs->getVector(),
                      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());
    }
  }
272

273
274
275
276
277
278
279
280
281
282
283
  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!");
  }
}
284

285

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

292
  Assembler<Traits> assembler(*globalBasis, matrixOperators, rhsOperators, constraints);
293

294
295
296
  auto gv = leafGridView();
  assembler.update(gv);
  assembler.assemble(*systemMatrix, *solution, *rhs, asmMatrix, asmVector);
297

298
299
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
300

301

302
303
template <class Traits>
void ProblemStat<Traits>::
304
writeFiles(AdaptInfo& adaptInfo, bool force)
305
306
{
  Timer t;
307
308
  for (auto writer : filewriter)
    writer->writeFiles(adaptInfo, force);
309
310
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
311

312
} // end namespace AMDiS