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

3
4
5
#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
6

7
8
9
namespace AMDiS {

template <class Traits>
10
11
12
13
void ProblemStat<Traits>::initialize(
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
14
{
15
16
17
18
19
20
21
22
23
  // create meshes
  if (mesh) {
    warning("mesh already created");
  }
  else {
    if (initFlag.isSet(CREATE_MESH) ||
        (!adoptFlag.isSet(INIT_MESH) &&
        (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
      createMesh();
24
    }
25

26
27
28
29
30
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
      mesh = adoptProblem->getMesh();
31
    }
32

33
34
    componentMeshes.resize(nComponents, mesh.get());
  }
35

36
37
  if (!mesh)
    warning("no mesh created");
38

39
40
41
42
43
  int globalRefinements = 0;
  Parameters::get(meshName + "->global refinements", globalRefinements);
  if (globalRefinements > 0) {
    mesh->globalRefine(globalRefinements);
  }
44
45


46
47
48
49
50
51
52
53
  // create fespace
  if (feSpaces) {
    warning("feSpaces already created");
  }
  else {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
      createFeSpaces();
54
    }
55

56
57
58
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
      feSpaces = adoptProblem->getFeSpaces();
59
    }
60
  }
61

62
63
  if (!feSpaces)
    warning("no feSpaces created\n");
64
65


66
67
68
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
69

70
71
72
73
74
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    solution = adoptProblem->getSolution();
    rhs = adoptProblem->getRhs();
    systemMatrix = adoptProblem->getSystemMatrix();
  }
75

76

77
78
79
80
81
82
83
  // create solver
  if (linearSolver) {
    warning("solver already created\n");
  }
  else {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
84

85
86
87
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
      test_exit(!linearSolver, "solver already created\n");
      linearSolver = adoptProblem->getSolver();
88
    }
89
90
  }

91
92
  if (!linearSolver) {
    warning("no solver created\n");
93
94
  }

95

96
97
98
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
99

100
101
  solution->compress();
}
102

103

104
// add matrix/vector operator terms
105

106
107
template <class Traits>
void ProblemStat<Traits>::
108
addMatrixOperator(ElementOperator const& op, int i, int j, double* factor, double* estFactor)
109
{
110
111
112
113
114
115
116
117
118
119
  matrixOperators[i][j].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
{
  matrixOperators[i][j].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
120
}
121

122
123
124
125
template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
{
126
127
  matrixOperators[i][j].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  matrixOperators[i][j].changing = true;
128
}
129

130
131
template <class Traits>
void ProblemStat<Traits>::
132
addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
133
134
135
136
{
  for (auto op : ops){
    int r = op.first.first;
    int c = op.first.second;
137
138
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second)});
    matrixOperators[r][c].changing = true;
139
  }
140
}
141

142
143
template <class Traits>
void ProblemStat<Traits>::
144
addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
145
146
147
148
{
  for (auto op : ops) {
    int r = op.first.first;
    int c = op.first.second;
149
150
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    matrixOperators[r][c].changing = matrixOperators[r][c].changing || op.second.second;
151
  }
152
}
153
154


155
156
template <class Traits>
void ProblemStat<Traits>::
157
158
addVectorOperator(ElementOperator const& op, int i, double* factor, double* estFactor)
{
159
160
  rhsOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
161
162
163
164
165
}

template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(IntersectionOperator const& op, int i, double* factor, double* estFactor)
166
{
167
168
  rhsOperators[i].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
169
}
170

171
172
173
174
template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
{
175
176
  rhsOperators[i].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  rhsOperators[i].changing = true;
177
}
178

179
180
template <class Traits>
void ProblemStat<Traits>::
181
addVectorOperator(std::map< int, ElementOperator> ops)
182
183
{
  for (auto op : ops) {
184
185
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second)});
    rhsOperators[op.first].changing = true;
186
187
  }
}
188

189
190
template <class Traits>
void ProblemStat<Traits>::
191
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
192
193
{
  for (auto op : ops) {
194
195
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    rhsOperators[op.first].changing = rhsOperators[op.first].changing || op.second.second;
196
  }
197
}
198
199


200
201
202
203
204
205
206
207
// Adds a Dirichlet boundary condition
template <class Traits>
  template <class Predicate, class Values>
void ProblemStat<Traits>::
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");
208

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

212
213
  test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
  test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);
214

215
//     boundaryConditionSet = true;
216

217
  using BcType = DirichletBC<WorldVector>;
218
  matrixOperators[row][col].dirichlet.emplace_back( new BcType{predicate, values} );
219
}
220

221

222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
template <class Traits>
void ProblemStat<Traits>::
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());
    }
  }
248

249
250
251
252
253
254
255
256
257
258
259
  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!");
  }
}
260

261

262
263
template <class Traits>
void ProblemStat<Traits>::
264
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix, bool asmVector)
265
266
267
{
  Timer t;

268
269
270
271
272
273
  Assembler<FeSpaces> assembler(getFeSpaces(), asmMatrix, asmVector);

  auto gridView = leafGridView();
  assembler.assemble(gridView,
                     *systemMatrix, *solution, *rhs,
                     matrixOperators, rhsOperators);
274

275
276
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
277

278

279
280
281
282
283
284
285
286
template <class Traits>
void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
  Timer t;
  filewriter->write(adaptInfo.getTime(), *solution);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
287

288
} // end namespace AMDiS