ProblemStat.inc.hpp 8.74 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
#include <amdis/AdaptInfo.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/LocalAssembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/common/Loops.hpp>
15

16
17
namespace AMDiS {

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

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

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

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


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

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

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


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

78
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
79
    solution_ = adoptProblem->solution_;
80
    estimates_ = adoptProblem->estimates_;
81
82
    rhs_ = adoptProblem->rhs_;
    systemMatrix_ = adoptProblem->systemMatrix_;
83
  }
84

85

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

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

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

104
105
106
107
108
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();

    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
109
      marker_ = adoptProblem->marker_;
110

111

112
113
114
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
115

116
  solution_->compress();
117
}
118

119

120
121
122
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
Praetorius, Simon's avatar
Praetorius, Simon committed
123
  marker_.clear();
124
125
  auto localView = globalBasis_->localView();
  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
126
  {
127
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
128
129
130
131

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

132
    std::string tp = to_string(treePath);
Praetorius, Simon's avatar
Praetorius, Simon committed
133
    auto newMarker = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
134
    if (newMarker)
135
136
      marker_.push_back(std::move(newMarker));

137
138
    // If there is more than one marker, and all components are defined
    // on the same grid, the maximum marking has to be enabled.
139
140
141
    // TODO: needs to be checked!
    if (marker_.size() > 1)
      marker_.back()->setMaximumMarking(true);
142
143
144
145
  });
}


146
147
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
148
{
Praetorius, Simon's avatar
Praetorius, Simon committed
149
  filewriter_.clear();
150
151
  auto localView = globalBasis_->localView();
  forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
152
  {
153
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
154
155
156
157

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

158
    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
159
    filewriter_.push_back(std::move(writer));
160
161
162
163
  });
}


164
// Adds a Dirichlet boundary condition
165
template <class Traits>
166
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
167
void ProblemStat<Traits>::
168
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
169
170
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
171
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
172

173
174
175
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));
176

177
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
178

179
  using Range = RangeType_t<decltype(i)>;
180
  using BcType = DirichletBC<WorldVector,Range>;
181
  auto bc = std::make_shared<BcType>(predicate, valueGridFct);
182
  constraints_[i][j].push_back(bc);
183
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
184
}
185

186

187
188
template <class Traits>
void ProblemStat<Traits>::
189
190
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
191
  Dune::Timer t;
192

193
  SolverInfo solverInfo(name_ + "->solver");
194
195
196
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

197
  solution_->compress();
198

199
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
200
201
202
                      solverInfo);

  if (solverInfo.getInfo() > 0) {
203
    msg("solution of discrete system needed {} seconds", t.elapsed());
204
205
206

    if (solverInfo.getAbsResidual() >= 0.0) {
      if (solverInfo.getRelResidual() >= 0.0)
207
208
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
          solverInfo.getAbsResidual(), solverInfo.getRelResidual());
209
      else
210
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.getAbsResidual());
211
212
    }
  }
213

214
215
216
217
218
219
220
221
  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() << " ";
222
    error_exit("Tolerance {} could not be reached!", tol_str.str());
223
224
  }
}
225

226

227
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
228
229
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
230
231
232
233
{
  Dune::Timer t;

  Flag markFlag = 0;
234
  for (auto& currentMarker : marker_)
235
    markFlag |= currentMarker->markGrid(adaptInfo);
236

237
  msg("markElements needed {} seconds", t.elapsed());
238
239
240
241
242

  return markFlag;
}


243
244
245
246
247
248
249
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  grid_->preAdapt();
Praetorius, Simon's avatar
Praetorius, Simon committed
250
251
252
253
254
255
256
257
  bool changed = grid_->adapt();

  // update data
  if (changed) {
    globalBasis_->update(gridView());
    solution_->resize(*globalBasis_);
  }

258
259
  grid_->postAdapt();

260
  msg("adaptGrid needed {} seconds", t.elapsed());
Praetorius, Simon's avatar
Praetorius, Simon committed
261
  return changed ? MESH_ADAPTED : Flag(0);
262
263
264
}


265
266
template <class Traits>
void ProblemStat<Traits>::
267
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
268
{
269
  Dune::Timer t;
270

271
272
273
274
275
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);
  solution_->resize(*globalBasis_);

276
277
278
  auto localView = globalBasis_->localView();
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
279
      for (auto bc : constraints_[rowNode][colNode])
280
        bc->init(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
281
282
283
284
285
286
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

  // 2. traverse grid and assemble operators on the elements
  for (auto const& element : elements(gridView())) {
287
    localView.bind(element);
288
289

    if (asmMatrix)
290
      systemMatrix_->assemble(localView, localView);
291
    if (asmVector)
292
      rhs_->assemble(localView);
293

294
    localView.unbind();
295
296
297
298
299
  }

  // 3. finish matrix insertion and apply dirichlet boundary conditions
  systemMatrix_->finish(asmMatrix);
  rhs_->finish(asmVector);
300
  (*solution_) = 0;
301

302
303
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
304
305
      // finish boundary condition
      for (auto bc : constraints_[rowNode][colNode])
306
        bc->finish(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
307
308
    });
  });
309

310
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
311
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
312
}
313

314

315
316
template <class Traits>
void ProblemStat<Traits>::
317
writeFiles(AdaptInfo& adaptInfo, bool force)
318
{
319
  Dune::Timer t;
320
  for (auto writer : filewriter_)
321
    writer->writeFiles(adaptInfo, force);
322
  msg("writeFiles needed {} seconds", t.elapsed());
323
}
324

325

326
} // end namespace AMDiS