ProblemStat.inc.hpp 8.33 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
10
11
12
namespace AMDiS {

template <class Traits>
void ProblemStat<Traits>::initialize(Flag initFlag,
                                      Self* adoptProblem,
                                      Flag adoptFlag)
13
{
14
15
16
17
18
19
20
21
22
  // 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();
23
    }
24

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

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

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

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


45
46
47
48
49
50
51
52
  // 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();
53
    }
54

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

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


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

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

75

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

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

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

94

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

99
100
  solution->compress();
}
101

102

103
// add matrix/vector operator terms
104

105
106
template <class Traits>
void ProblemStat<Traits>::
107
addMatrixOperator(ElementOperator const& op, int i, int j, double* factor, double* estFactor)
108
{
109
110
111
112
113
114
115
116
117
118
  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;
119
}
120

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

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

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


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

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

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

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

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


199
200
201
202
203
204
205
206
// 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");
207

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

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

214
//     boundaryConditionSet = true;
215

216
  using BcType = DirichletBC<WorldVector>;
217
  matrixOperators[row][col].dirichlet.emplace_back( new BcType{predicate, values} );
218
}
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
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());
    }
  }
247

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

260

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

267
268
269
270
  Assembler<FeSpaces> assembler(getFeSpaces());
  assembler.assemble(*systemMatrix, *solution, *rhs,
                     matrixOperators, rhsOperators,
                     asmMatrix_, asmVector_);
271

272
273
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
274

275

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

285
} // end namespace AMDiS