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

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

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

10
11
12
13
14
15
#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/LocalAssembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/common/Loops.hpp>
16
// #include <amdis/Estimator.hpp>
17

18
19
namespace AMDiS {

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

37
38
39
40
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
41
      grid_ = adoptProblem->grid_;
42
    }
43
  }
44

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

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


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

65
66
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
67
68
      globalBasis_ = adoptProblem->globalBasis_;
      initGlobalBasis(*globalBasis_);
69
    }
70
  }
71

72
  if (!globalBasis_)
73
    warning("no globalBasis created\n");
74
75


76
77
78
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
79

80
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
81
    solution_ = adoptProblem->solution_;
82
    estimates = adoptProblem->estimates;
83
84
    rhs_ = adoptProblem->rhs_;
    systemMatrix_ = adoptProblem->systemMatrix_;
85
  }
86

87

88
  // create solver
89
  if (linearSolver_) {
90
91
92
93
94
    warning("solver already created\n");
  }
  else {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
95

96
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
97
98
      test_exit(!linearSolver_, "solver already created\n");
      linearSolver_ = adoptProblem->linearSolver_;
99
    }
100
101
  }

102
  if (!linearSolver_) {
103
    warning("no solver created\n");
104
105
  }

106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
/*
  // create estimator
  estimator.resize(nComponents, NULL);

  if (initFlag.isSet(INIT_ESTIMATOR))
    createEstimator();

  if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
    estimator = adoptProblem->estimator;
*/

  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
      marker = adoptProblem->marker;

124

125
126
127
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
128

129
  solution_->compress();
130
}
131

132

133
134
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
135
{
136
  AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
137
  {
138
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
139
140
141
142

    if (!Parameters::get<std::string>(componentName + "->filename"))
      return;

143
    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
144
    filewriter_.push_back(std::move(writer));
145
146
147
148
  });
}


149
// add matrix/vector operator terms
150

151
template <class Traits>
152
  template <class Operator, class RowTreePath, class ColTreePath>
153
void ProblemStat<Traits>::addMatrixOperator(
154
    Operator const& preOp,
155
156
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
157
{
158
  static_assert( Concepts::PreTreePath<RowTreePath>,
159
      "row must be a valid treepath, or an integer/index-constant");
160
  static_assert( Concepts::PreTreePath<ColTreePath>,
161
162
      "col must be a valid treepath, or an integer/index-constant");

163
164
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
165

166
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
167
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
168

169
170
  matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
  matrixOperators_[i][j].changing = true;
171
}
172

173

174
template <class Traits>
175
  template <class Operator, class RowTreePath, class ColTreePath>
176
void ProblemStat<Traits>::addMatrixOperator(
177
    BoundaryType b,
178
    Operator const& preOp,
179
180
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
181
{
182
  static_assert( Concepts::PreTreePath<RowTreePath>,
183
      "row must be a valid treepath, or an integer/index-constant");
184
  static_assert( Concepts::PreTreePath<ColTreePath>,
185
186
      "col must be a valid treepath, or an integer/index-constant");

187
188
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
189
  using Intersection = typename GridView::Intersection;
190

191
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
192
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
193

194
195
  matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
  matrixOperators_[i][j].changing = true;
196
}
197
198


199
template <class Traits>
200
  template <class Operator, class TreePath>
201
void ProblemStat<Traits>::addVectorOperator(
202
    Operator const& preOp,
203
204
    TreePath path,
    double* factor, double* estFactor)
205
{
206
  static_assert( Concepts::PreTreePath<TreePath>,
207
208
      "path must be a valid treepath, or an integer/index-constant");

209
  auto i = child(localView_->tree(), makeTreePath(path));
210

211
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
212
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
213

214
215
  rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
  rhsOperators_[i].changing = true;
216
}
217

218

219
template <class Traits>
220
  template <class Operator, class TreePath>
221
void ProblemStat<Traits>::addVectorOperator(
222
    BoundaryType b,
223
    Operator const& preOp,
224
225
    TreePath path,
    double* factor, double* estFactor)
226
{
227
  static_assert( Concepts::PreTreePath<TreePath>,
228
229
      "path must be a valid treepath, or an integer/index-constant");

230
  auto i = child(localView_->tree(), makeTreePath(path));
231
  using Intersection = typename GridView::Intersection;
232

233
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
234
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
235

236
237
  rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
  rhsOperators_[i].changing = true;
238
}
239
240


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

250
251
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
252

253
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
254

255
  using Range = RangeType_t<decltype(i)>;
256
  using BcType = DirichletBC<WorldVector,Range>;
257
  auto bc = std::make_shared<BcType>(predicate, valueGridFct);
258
  constraints_[i][j].push_back(bc);
259
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
260
}
261

262

263
264
template <class Traits>
void ProblemStat<Traits>::
265
266
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
267
  Dune::Timer t;
268

269
  SolverInfo solverInfo(name_ + "->solver");
270
271
272
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

273
  solution_->compress();
274

275
  linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
276
277
278
279
280
281
282
283
284
285
286
287
288
                      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());
    }
  }
289

290
291
292
293
294
295
296
297
298
299
300
  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!");
  }
}
301

302

303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
template <class Traits>
void ProblemStat<Traits>::
estimate(AdaptInfo& adaptInfo)
{
  Dune::Timer t;
/*
  for (std::size_t i = 0; i < nComponents; i++) {
    Estimator *scalEstimator = estimator[i];

    if (scalEstimator) {
      scalEstimator->estimate(adaptInfo->getTimestep());

      adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
      adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);

      if (adaptInfo->getRosenbrockMode() == false) {
        adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
        adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
      }
	  }
  }

//#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
//    MPI::COMM_WORLD.Barrier();
//#endif
*/
  msg("estimation of the error needed ", t.elapsed(), "seconds\n");
}


template <class Traits>
Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  test_exit(static_cast<unsigned int>(nComponents) == marker.size(), "Wrong number of markers!\n");

  Flag markFlag = 0;
  for (std::size_t i = 0; i < nComponents; i++)
    if (marker[i])
      markFlag |= marker[i]->markGrid(adaptInfo);

  msg("markElements needed ", t.elapsed(), " seconds");

  return markFlag;
}


template <class Traits>
Flag ProblemStat<Traits>::refineMesh(AdaptInfo& adaptInfo)
{
  // TODO: data transfer

  grid->preAdapt();
/*
  const auto& gridView = grid->leafGridView();
  const auto& indexSet = gridView.indexSet();
  const auto& idSet = grid−>localIdSet();

  // Save data in container during adaptation
  for (const auto& v : vertices(gridView)) {
    persistentContainer[idSet.id(v)] = data[indexSet.index(v)];
  } */

  Flag flag = grid->adapt();
/*
  // Unpack data
  data.resize(gridView.size(dim));
  for (const auto& v : vertices(gridView)) {
    data[indexSet.index(v)] = persistentContainer[idSet.id(v)];
  }*/

  grid->postAdapt();

  return flag;
}


381
382
template <class Traits>
void ProblemStat<Traits>::
383
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
384
{
385
  Dune::Timer t;
386

387
388
389
390
  Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
  auto gv = leafGridView();
  globalBasis_->update(gv);
  assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
391

392
393
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
394

395

396
397
template <class Traits>
void ProblemStat<Traits>::
398
writeFiles(AdaptInfo& adaptInfo, bool force)
399
{
400
  Dune::Timer t;
401
  for (auto writer : filewriter_)
402
    writer->writeFiles(adaptInfo, force);
403
404
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
405

406

407
} // end namespace AMDiS