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

3
4
#include <dune/typetree/childextraction.hh>

5
#include <dune/amdis/AdaptInfo.hpp>
6
#include <dune/amdis/Assembler.hpp>
7
#include <dune/amdis/FileWriter.hpp>
8
9
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
10

11
12
namespace AMDiS {

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

30
31
32
33
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
34
      grid = adoptProblem->getGrid();
35
    }
36

37
    componentGrids.resize(nComponents, grid.get());
38
  }
39

40
41
  if (!grid)
    warning("no grid created");
42

43
  int globalRefinements = 0;
44
  Parameters::get(gridName + "->global refinements", globalRefinements);
45
  if (globalRefinements > 0) {
46
    grid->globalRefine(globalRefinements);
47
  }
48
49


50
  // create fespace
51
52
  if (globalBasis) {
    warning("globalBasis already created");
53
54
55
56
  }
  else {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
57
      createGlobalBasis();
58
    }
59

60
61
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
62
      globalBasis = adoptProblem->getGlobalBasis();
63
    }
64
  }
65

66
67
  if (!globalBasis)
    warning("no globalBasis created\n");
68
69


70
71
72
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
73

74
75
76
77
78
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    solution = adoptProblem->getSolution();
    rhs = adoptProblem->getRhs();
    systemMatrix = adoptProblem->getSystemMatrix();
  }
79

80

81
82
83
84
85
86
87
  // create solver
  if (linearSolver) {
    warning("solver already created\n");
  }
  else {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
88

89
90
91
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
      test_exit(!linearSolver, "solver already created\n");
      linearSolver = adoptProblem->getSolver();
92
    }
93
94
  }

95
96
  if (!linearSolver) {
    warning("no solver created\n");
97
98
  }

99

100
101
102
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
103

104
105
  solution->compress();
}
106

107

108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::createFileWriter()
{
  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;

    using TP = decltype(treePath);
    auto writer = std::make_shared<FileWriter<GlobalBasis,TP>>(componentName, globalBasis, treePath, solution);
    filewriter.push_back(writer);
  });
}


126
// add matrix/vector operator terms
127

128
template <class GlobalBasis>
129
130
131
132
133
  template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator(
    ElementOperator const& op,
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
134
{
135
136
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
137
138
139
140
  matrixOperators[i][j].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
}

141

142
template <class GlobalBasis>
143
144
145
146
147
  template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator(
    IntersectionOperator const& op,
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
148
{
149
150
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
151
152
  matrixOperators[i][j].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
153
}
154

155

156
template <class GlobalBasis>
157
158
159
160
161
162
  template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator(
    BoundaryType b,
    IntersectionOperator const& op,
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
163
{
164
165
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
166
167
  matrixOperators[i][j].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  matrixOperators[i][j].changing = true;
168
}
169

170
#if 0
171
172
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
173
addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
174
175
176
177
{
  for (auto op : ops){
    int r = op.first.first;
    int c = op.first.second;
178
179
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second)});
    matrixOperators[r][c].changing = true;
180
  }
181
}
182

183
184
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
185
addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
186
187
188
189
{
  for (auto op : ops) {
    int r = op.first.first;
    int c = op.first.second;
190
191
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    matrixOperators[r][c].changing = matrixOperators[r][c].changing || op.second.second;
192
  }
193
}
194
#endif
195
196


197
template <class GlobalBasis>
198
199
200
201
202
  template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator(
    ElementOperator const& op,
    TreePath path,
    double* factor, double* estFactor)
203
{
204
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
205
206
  rhsOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
207
208
}

209

210
template <class GlobalBasis>
211
212
213
214
215
  template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator(
    IntersectionOperator const& op,
    TreePath path,
    double* factor, double* estFactor)
216
{
217
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
218
219
  rhsOperators[i].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
220
}
221

222

223
template <class GlobalBasis>
224
225
226
227
228
229
  template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator(
    BoundaryType b,
    IntersectionOperator const& op,
    TreePath path,
    double* factor, double* estFactor)
230
{
231
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
232
233
  rhsOperators[i].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  rhsOperators[i].changing = true;
234
}
235

236
#if 0
237
238
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
239
addVectorOperator(std::map< int, ElementOperator> ops)
240
241
{
  for (auto op : ops) {
242
243
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second)});
    rhsOperators[op.first].changing = true;
244
245
  }
}
246

247
248
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
249
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
250
251
{
  for (auto op : ops) {
252
253
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    rhsOperators[op.first].changing = rhsOperators[op.first].changing || op.second.second;
254
  }
255
}
256
#endif
257
258


259
// Adds a Dirichlet boundary condition
260
template <class GlobalBasis>
261
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
262
void ProblemStat<GlobalBasis>::
263
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
264
265
266
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
    "Function passed to addDirichletBC for predicate does not model the Functor<bool(WorldVector)> concept");
267

268
269
  static_assert( Concepts::Callable<Values, WorldVector>,
    "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");
270

271
272
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
273

274
275
  using Range = decltype(values(std::declval<WorldVector>()));
  using BcType = DirichletBC<WorldVector,Range>;
276

277
278
  auto bc = std::make_shared<BcType>(predicate, values);
  constraints[i][i].push_back(bc);
279
}
280

281

282
283
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
  Timer t;

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

  solution->compress();
  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());
    }
  }
308

309
310
311
312
313
314
315
316
317
318
319
  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!");
  }
}
320

321

322
323
324
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
325
326
327
{
  Timer t;

328
  Assembler<GlobalBasis> assembler(*globalBasis, matrixOperators, rhsOperators, constraints);
329

330
331
332
  auto gv = leafGridView();
  assembler.update(gv);
  assembler.assemble(*systemMatrix, *solution, *rhs, asmMatrix, asmVector);
333

334
335
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
336

337

338
339
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
340
writeFiles(AdaptInfo& adaptInfo, bool force)
341
342
{
  Timer t;
343
344
  for (auto writer : filewriter)
    writer->writeFiles(adaptInfo, force);
345
346
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
347

348
} // end namespace AMDiS