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

3
4
5
6
#include <map>
#include <string>
#include <utility>

7
8
#include <dune/typetree/childextraction.hh>

9
#include <dune/amdis/AdaptInfo.hpp>
10
#include <dune/amdis/Assembler.hpp>
11
#include <dune/amdis/FileWriter.hpp>
12
13
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>
14

15
16
namespace AMDiS {

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

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

41
    componentGrids.resize(nComponents, grid.get());
42
  }
43

44
45
  if (!grid)
    warning("no grid created");
46

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


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

64
65
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
66
      globalBasis = adoptProblem->getGlobalBasis();
67
    }
68
  }
69

70
71
  if (!globalBasis)
    warning("no globalBasis created\n");
72
73


74
75
76
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
77

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

84

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

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

99
100
  if (!linearSolver) {
    warning("no solver created\n");
101
102
  }

103

104
105
106
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
107

108
109
  solution->compress();
}
110

111

112
113
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
114
115
116
117
118
119
120
121
122
{
  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;

123
    auto writer = makeFileWriterPtr<Traits>(componentName, globalBasis, treePath, solution);
124
125
126
127
128
    filewriter.push_back(writer);
  });
}


129
// add matrix/vector operator terms
130

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

144

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

158

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

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

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


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

212

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

225

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

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

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


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

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

274
275
  auto i = child(globalBasis->localView().tree(), makeTreePath(row));
  auto j = child(globalBasis->localView().tree(), makeTreePath(col));
276

277
278
  using Range = decltype(values(std::declval<WorldVector>()));
  using BcType = DirichletBC<WorldVector,Range>;
279

280
281
  auto bc = std::make_shared<BcType>(predicate, values);
  constraints[i][i].push_back(bc);
282
}
283

284

285
286
template <class Traits>
void ProblemStat<Traits>::
287
288
289
290
291
292
293
294
295
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
  Timer t;

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

  solution->compress();
296
297
298
299
300

  msg("matrix: ",num_rows(systemMatrix->getMatrix()), " x ", num_cols(systemMatrix->getMatrix()));
  msg("x: ",num_rows(solution->getVector()));
  msg("rhs: ",num_rows(rhs->getVector()));

301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
  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());
    }
  }
316

317
318
319
320
321
322
323
324
325
326
327
  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!");
  }
}
328

329

330
331
template <class Traits>
void ProblemStat<Traits>::
332
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
333
334
335
{
  Timer t;

336
  Assembler<Traits> assembler(*globalBasis, matrixOperators, rhsOperators, constraints);
337

338
339
340
  auto gv = leafGridView();
  assembler.update(gv);
  assembler.assemble(*systemMatrix, *solution, *rhs, asmMatrix, asmVector);
341

342
343
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
344

345

346
347
template <class Traits>
void ProblemStat<Traits>::
348
writeFiles(AdaptInfo& adaptInfo, bool force)
349
350
{
  Timer t;
351
352
  for (auto writer : filewriter)
    writer->writeFiles(adaptInfo, force);
353
354
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
355

356
} // end namespace AMDiS