ProblemStat.inc.hpp 9.64 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
8
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
9
10
#include <dune/amdis/linear_algebra/HierarchicWrapper.hpp>
#include <dune/amdis/utility/TreePath.hpp>
11

12
13
namespace AMDiS {

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

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

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

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

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


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

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

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


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

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

81

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

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

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

100

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

105
106
  solution->compress();
}
107

108

109
// add matrix/vector operator terms
110

111
template <class GlobalBasis>
112
113
114
115
116
  template <class RowTreePath, class ColTreePath>
void ProblemStat<GlobalBasis>::addMatrixOperator(
    ElementOperator const& op,
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
117
{
118
119
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
120
121
122
123
  matrixOperators[i][j].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
}

124

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

138

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

153
#if 0
154
155
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
156
addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
157
158
159
160
{
  for (auto op : ops){
    int r = op.first.first;
    int c = op.first.second;
161
162
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second)});
    matrixOperators[r][c].changing = true;
163
  }
164
}
165

166
167
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
168
addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
169
170
171
172
{
  for (auto op : ops) {
    int r = op.first.first;
    int c = op.first.second;
173
174
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    matrixOperators[r][c].changing = matrixOperators[r][c].changing || op.second.second;
175
  }
176
}
177
#endif
178
179


180
template <class GlobalBasis>
181
182
183
184
185
  template <class TreePath>
void ProblemStat<GlobalBasis>::addVectorOperator(
    ElementOperator const& op,
    TreePath path,
    double* factor, double* estFactor)
186
{
187
  auto i = child(globalBasis->localView().tree(), makeTreePath(path));
188
189
  rhsOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
190
191
}

192

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

205

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

219
#if 0
220
221
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
222
addVectorOperator(std::map< int, ElementOperator> ops)
223
224
{
  for (auto op : ops) {
225
226
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second)});
    rhsOperators[op.first].changing = true;
227
228
  }
}
229

230
231
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
232
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
233
234
{
  for (auto op : ops) {
235
236
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    rhsOperators[op.first].changing = rhsOperators[op.first].changing || op.second.second;
237
  }
238
}
239
#endif
240
241


242
// Adds a Dirichlet boundary condition
243
template <class GlobalBasis>
244
  template <class Predicate, class Values>
245
void ProblemStat<GlobalBasis>::
246
247
248
249
addDirichletBC(Predicate const& predicate, int row, int col, Values const& values)
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
    "Function passed to addDirichletBC for predicate does not model the Functor<bool(WorldVector)> concept");
250

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

254
#if 0
255
256
  test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
  test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);
257

258
//     boundaryConditionSet = true;
259

260
  using BcType = DirichletBC<WorldVector>;
261
  matrixOperators[row][col].dirichlet.emplace_back( new BcType{predicate, values} );
262
#endif
263
}
264

265

266
267
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
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());
    }
  }
292

293
294
295
296
297
298
299
300
301
302
303
  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!");
  }
}
304

305

306
307
308
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
309
310
311
{
  Timer t;

312
  Assembler<GlobalBasis> assembler(*globalBasis, matrixOperators, rhsOperators);
313

314
315
316
  auto gv = leafGridView();
  assembler.update(gv);
  assembler.assemble(*systemMatrix, *solution, *rhs, asmMatrix, asmVector);
317

318
319
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
320

321

322
323
324
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
writeFiles(AdaptInfo& adaptInfo, bool /*force*/)
325
326
327
328
329
{
  Timer t;
  filewriter->write(adaptInfo.getTime(), *solution);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
330

331
} // end namespace AMDiS