ProblemStat.inc.hpp 8.88 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
14
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
15

16
17
namespace AMDiS {

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

35
36
37
38
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
39
      grid = adoptProblem->getGrid();
40
    }
41

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

45
46
  if (!grid)
    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
57
  if (globalBasis) {
    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
      globalBasis = adoptProblem->getGlobalBasis();
68
    }
69
  }
70

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


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

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

85

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

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

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

104

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

109
110
  solution->compress();
}
111

112

113
114
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
115
116
117
118
119
120
121
122
123
{
  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;

124
    auto writer = makeFileWriterPtr<Traits>(componentName, globalBasis, treePath, solution);
125
126
127
128
129
    filewriter.push_back(writer);
  });
}


130
// add matrix/vector operator terms
131

132
template <class Traits>
133
  template <class Operator, class RowTreePath, class ColTreePath>
134
void ProblemStat<Traits>::addMatrixOperator(
135
    Operator const& op,
136
137
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
138
{
139
140
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
141

142
  auto localAssembler = makeLocalAssemblerPtr<Element>(op, i, j);
143

144
  matrixOperators[i][j].element.push_back({localAssembler, factor, estFactor});
145
  matrixOperators[i][j].changing = true;
146
}
147

148

149
template <class Traits>
150
  template <class Operator, class RowTreePath, class ColTreePath>
151
void ProblemStat<Traits>::addMatrixOperator(
152
    BoundaryType b,
153
    Operator const& op,
154
155
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
156
{
157
158
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
159
  using Intersection = typename GridView::Intersection;
160

161
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i, j);
162

163
164
  matrixOperators[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
  matrixOperators[i][j].changing = true;
165
}
166
167


168
template <class Traits>
169
  template <class Operator, class TreePath>
170
void ProblemStat<Traits>::addVectorOperator(
171
    Operator const& op,
172
173
    TreePath path,
    double* factor, double* estFactor)
174
{
175
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
176

177
  auto localAssembler = makeLocalAssemblerPtr<Element>(op, i);
178

179
  rhsOperators[i].element.push_back({localAssembler, factor, estFactor});
180
  rhsOperators[i].changing = true;
181
}
182

183

184
template <class Traits>
185
  template <class Operator, class TreePath>
186
void ProblemStat<Traits>::addVectorOperator(
187
    BoundaryType b,
188
    Operator const& op,
189
190
    TreePath path,
    double* factor, double* estFactor)
191
{
192
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
193
  using Intersection = typename GridView::Intersection;
194

195
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i);
196

197
198
  rhsOperators[i].boundary.push_back({localAssembler, factor, estFactor, b});
  rhsOperators[i].changing = true;
199
}
200
201


202
// Adds a Dirichlet boundary condition
203
template <class Traits>
204
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
205
void ProblemStat<Traits>::
206
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
207
208
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
209
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
210

211
212
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
213

214
215
216
217
218
219
220
221
222
  addDirichletBCImpl(predicate, i, j, values);
}

template <class Traits>
  template <class Predicate, class RowNode, class ColNode>
void ProblemStat<Traits>::
addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, RangeType<RowNode> const& value)
{
  using Range = RangeType<RowNode>;
223
  using BcType = DirichletBC<WorldVector,Range>;
224
225
226
227
228
229
230
231
232
233
234
235
236
  auto bc = std::make_shared<BcType>(predicate, [value](auto const&) { return value; });
  constraints[row][col].push_back(bc);
}

template <class Traits>
  template <class Predicate, class RowNode, class ColNode, class Values>
    std::enable_if_t<not std::is_convertible<Values, RangeType<RowNode>>::value>
ProblemStat<Traits>::
addDirichletBCImpl(Predicate const& predicate, RowNode const& row, ColNode const& col, Values const& values)
{
  using Range = RangeType<RowNode>;
  static_assert( Concepts::Functor<Values, Range(WorldVector)>,
    "Function passed to addDirichletBC for `values` does not model the Functor<Range(WorldVector)> concept");
237

238
  using BcType = DirichletBC<WorldVector,Range>;
239
  auto bc = std::make_shared<BcType>(predicate, values);
240
  constraints[row][col].push_back(bc);
241
}
242

243

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

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

  solution->compress();
255

256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
  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());
    }
  }
271

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

284

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

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

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

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

300

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

311
} // end namespace AMDiS