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
146
    // replace i with treePath once supported
    auto marker_ = Marker<Traits>::createMarker(
147
      componentName, std::to_string(i), estimates[std::to_string(i)], componentGrids[i]);
148
149
    if (marker_)
      marker.push_back(marker_);
150
151
152
153
154
155
156
157
158
159

    // 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);
  });
}


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

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

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


176
// add matrix/vector operator terms
177

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

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

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

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

200

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

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

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

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


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

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

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

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

245

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

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

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

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


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

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

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

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

289

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

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

300
  solution_->compress();
301

302
  linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
303
304
305
306
307
308
309
310
311
312
313
314
315
                      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
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
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;
366
367
  for (auto marker_ : marker)
    markFlag |= marker_->markGrid(adaptInfo);
368
369
370
371
372
373
374
375
376
377

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

  return markFlag;
}


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

  grid->preAdapt();
381
382
383

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

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

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

  grid->postAdapt();

  return flag;
}


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

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

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

427

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

438

439
} // end namespace AMDiS