ProblemStat.inc.hpp 8.45 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
namespace AMDiS {

9
10
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::initialize(
11
12
13
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
14
{
15
16
17
  // create grides
  if (grid) {
    warning("grid already created");
18
19
20
21
22
23
  }
  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
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
30
      grid = adoptProblem->getMesh();
31
    }
32

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

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

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


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

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

62
63
  if (!globalBasis)
    warning("no globalBasis 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 GlobalBasis>
void ProblemStat<GlobalBasis>::
108
addMatrixOperator(ElementOperator const& op, int i, int j, double* factor, double* estFactor)
109
{
110
111
112
113
  matrixOperators[i][j].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
}

114
115
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
116
117
118
119
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
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
124
125
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 GlobalBasis>
void ProblemStat<GlobalBasis>::
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 GlobalBasis>
void ProblemStat<GlobalBasis>::
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 GlobalBasis>
void ProblemStat<GlobalBasis>::
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
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
165
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
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
173
174
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 GlobalBasis>
void ProblemStat<GlobalBasis>::
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 GlobalBasis>
void ProblemStat<GlobalBasis>::
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
// Adds a Dirichlet boundary condition
201
template <class GlobalBasis>
202
  template <class Predicate, class Values>
203
void ProblemStat<GlobalBasis>::
204
205
206
207
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
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
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
264
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
265
266
267
{
  Timer t;

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

  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
template <class GlobalBasis>
void ProblemStat<GlobalBasis>::
writeFiles(AdaptInfo& adaptInfo, bool /*force*/)
282
283
284
285
286
{
  Timer t;
  filewriter->write(adaptInfo.getTime(), *solution);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
287

288
} // end namespace AMDiS