ProblemStat.inc.hpp 12.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
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
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
    // replace with treePath once supported
    marker.push_back(Marker<Traits>::createMarker(
      componentName, i, estimates[i], componentGrids[i]));

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


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

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

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


174
// add matrix/vector operator terms
175

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

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

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

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

198

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

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

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

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


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

234
  auto i = child(localView_->tree(), makeTreePath(path));
235

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

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

243

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

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

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

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


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

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

278
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
279

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

287

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

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

298
  solution_->compress();
299

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

315
316
317
318
319
320
321
322
323
324
325
  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!");
  }
}
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
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;
364
365
  for (auto marker_ : marker)
    markFlag |= marker_->markGrid(adaptInfo);
366
367
368
369
370
371
372
373
374
375

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

  return markFlag;
}


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

  grid->preAdapt();
379
380
381

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

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

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

  grid->postAdapt();

  return flag;
}


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

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

422
423
  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}
424

425

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

436

437
} // end namespace AMDiS