ProblemStat.inc.hpp 12.4 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
135
136
137
138
139
140
141
142
143
144
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
  auto localView = globalBasis->localView();
  forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
  {
    std::string componentName = name + "->marker[" + to_string(treePath) + "]";

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

    int i = std::stoi(to_string(treePath)); // TODO: To be removed
145
    // replace i with treePath once supported
146
147
148
149
150
151
    auto newMarker = Marker<Traits>::createMarker(
      componentName, std::to_string(i), estimates[std::to_string(i)], grid);
    if (newMarker)
    {
      marker.push_back(std::move(newMarker));
    }
152
153
154
155
156
157
158
159
160
    // If there is more than one marker, and all components are defined
    // on the same grid, the maximum marking has to be enabled.
    // TODO: What about two markers each for two grids?
    if (marker.size() > 1 && nGrids == 1)
 	    marker.back()->setMaximumMarking(true);
  });
}


161
162
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
163
{
164
  AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
165
  {
166
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
167
168
169
170

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

171
    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
172
    filewriter_.push_back(std::move(writer));
173
174
175
176
  });
}


177
// add matrix/vector operator terms
178

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

191
192
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
193

194
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
195
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
196

197
198
  matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
  matrixOperators_[i][j].changing = true;
199
}
200

201

202
template <class Traits>
203
  template <class Operator, class RowTreePath, class ColTreePath>
204
void ProblemStat<Traits>::addMatrixOperator(
205
    BoundaryType b,
206
    Operator const& preOp,
207
208
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
209
{
210
  static_assert( Concepts::PreTreePath<RowTreePath>,
211
      "row must be a valid treepath, or an integer/index-constant");
212
  static_assert( Concepts::PreTreePath<ColTreePath>,
213
214
      "col must be a valid treepath, or an integer/index-constant");

215
216
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
217
  using Intersection = typename GridView::Intersection;
218

219
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
220
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
221

222
223
  matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
  matrixOperators_[i][j].changing = true;
224
}
225
226


227
template <class Traits>
228
  template <class Operator, class TreePath>
229
void ProblemStat<Traits>::addVectorOperator(
230
    Operator const& preOp,
231
232
    TreePath path,
    double* factor, double* estFactor)
233
{
234
  static_assert( Concepts::PreTreePath<TreePath>,
235
236
      "path must be a valid treepath, or an integer/index-constant");

237
  auto i = child(localView_->tree(), makeTreePath(path));
238

239
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
240
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
241

242
243
  rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
  rhsOperators_[i].changing = true;
244
}
245

246

247
template <class Traits>
248
  template <class Operator, class TreePath>
249
void ProblemStat<Traits>::addVectorOperator(
250
    BoundaryType b,
251
    Operator const& preOp,
252
253
    TreePath path,
    double* factor, double* estFactor)
254
{
255
  static_assert( Concepts::PreTreePath<TreePath>,
256
257
      "path must be a valid treepath, or an integer/index-constant");

258
  auto i = child(localView_->tree(), makeTreePath(path));
259
  using Intersection = typename GridView::Intersection;
260

261
  auto op = makeGridOperator(preOp, globalBasis_->gridView());
262
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
263

264
265
  rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
  rhsOperators_[i].changing = true;
266
}
267
268


269
// Adds a Dirichlet boundary condition
270
template <class Traits>
271
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
272
void ProblemStat<Traits>::
273
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
274
275
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
276
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
277

278
279
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
280

281
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
282

283
  using Range = RangeType_t<decltype(i)>;
284
  using BcType = DirichletBC<WorldVector,Range>;
285
  auto bc = std::make_shared<BcType>(predicate, valueGridFct);
286
  constraints_[i][j].push_back(bc);
287
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
288
}
289

290

291
292
template <class Traits>
void ProblemStat<Traits>::
293
294
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
295
  Dune::Timer t;
296

297
  SolverInfo solverInfo(name_ + "->solver");
298
299
300
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

301
  solution_->compress();
302

303
  linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
304
305
306
307
308
309
310
311
312
313
314
315
316
                      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());
    }
  }
317

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

  Flag markFlag = 0;
367
368
  for (auto currentMarker : marker)
    markFlag |= currentMarker->markGrid(adaptInfo);
369
370
371
372
373
374
375
376
377
378

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

  return markFlag;
}


template <class Traits>
Flag ProblemStat<Traits>::refineMesh(AdaptInfo& adaptInfo)
{
379
  // TODO: actual element data, other components
380
381

  grid->preAdapt();
382
383
384

  std::map<typename Grid::LocalIdSet::IdType, double> vertexData;
//  std::map<typename Grid::LocalIdSet::IdType, double> elementData;
385
  const auto& gridView = globalBasis->gridView();
386
  const auto& indexSet = gridView.indexSet();
387
  const auto& idSet = grid->localIdSet();
388
389
390

  // Save data in container during adaptation
  for (const auto& v : vertices(gridView)) {
391
392
393
394
395
    vertexData[idSet.id(v)] = (*solution)[indexSet.index(v)];
  }/*
  for (const auto& e : elements(gridView)) {
    elementData[idSet.id(e)] = something[indexSet.index(e)];
  }*/
396
397
398
399

  Flag flag = grid->adapt();
/*
  // Unpack data
400
  solution->compress();
401
  for (const auto& v : vertices(gridView)) {
402
403
404
405
    (*solution)[indexSet.index(v)] = vertexData[idSet.id(v)];
  }/*
  for (const auto& e : elements(gridView)) {
    something[indexSet.index(e)] = elementData[idSet.id(e)];
406
407
408
409
410
411
412
413
  }*/

  grid->postAdapt();

  return flag;
}


414
415
template <class Traits>
void ProblemStat<Traits>::
416
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
417
{
418
  Dune::Timer t;
419

420
421
422
423
  Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
  auto gv = leafGridView();
  globalBasis_->update(gv);
  assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);
424

425
426
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
427

428

429
430
template <class Traits>
void ProblemStat<Traits>::
431
writeFiles(AdaptInfo& adaptInfo, bool force)
432
{
433
  Dune::Timer t;
434
  for (auto writer : filewriter_)
435
    writer->writeFiles(adaptInfo, force);
436
437
  msg("writeFiles needed ", t.elapsed(), " seconds");
}
438

439

440
} // end namespace AMDiS