ProblemStat.inc.hpp 8.84 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()
{
123
  AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
124
  {
125
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
126
127
128
129

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

130
131
    std::string tp = to_string(treePath);
    auto newMarker = Marker<Traits>::createMarker(componentName, tp, estimates_[tp], grid_);
132
    if (newMarker)
133
134
      marker_.push_back(std::move(newMarker));

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


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

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

154
    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
155
    filewriter_.push_back(std::move(writer));
156
157
158
159
  });
}


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

169
170
  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
171

172
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
173

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

181

182
183
template <class Traits>
void ProblemStat<Traits>::
184
185
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
186
  Dune::Timer t;
187

188
  SolverInfo solverInfo(name_ + "->solver");
189
190
191
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

192
  solution_->compress();
193

194
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
195
196
197
                      solverInfo);

  if (solverInfo.getInfo() > 0) {
198
    msg("solution of discrete system needed {} seconds", t.elapsed());
199
200
201

    if (solverInfo.getAbsResidual() >= 0.0) {
      if (solverInfo.getRelResidual() >= 0.0)
202
203
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
          solverInfo.getAbsResidual(), solverInfo.getRelResidual());
204
      else
205
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.getAbsResidual());
206
207
    }
  }
208

209
210
211
212
213
214
215
216
  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() << " ";
217
    error_exit("Tolerance {} could not be reached!", tol_str.str());
218
219
  }
}
220

221

222
223
224
225
226
227
template <class Traits>
Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  Flag markFlag = 0;
228
  for (auto& currentMarker : marker_)
229
    markFlag |= currentMarker->markGrid(adaptInfo);
230

231
  msg("markElements needed {} seconds", t.elapsed());
232
233
234
235
236

  return markFlag;
}


237
238
239
240
241
242
243
244
245
246
247
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  grid_->preAdapt();
  grid_->adapt();
  globalBasis_->update(gridView());
  grid_->postAdapt();

248
  msg("adaptGrid needed {} seconds", t.elapsed());
249
250
251
252
  return 0;
}


253
254
template <class Traits>
void ProblemStat<Traits>::
255
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
256
{
257
  Dune::Timer t;
258

259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);
  solution_->resize(*globalBasis_);

  forEachNode_(localView_->tree(), [&,this](auto const& rowNode, auto rowTreePath) {
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTreePath);
    forEachNode_(localView_->tree(), [&,this](auto const& colNode, auto colTreePath) {
      auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTreePath);
      for (auto bc : constraints_[rowNode][colNode])
        bc->init(*systemMatrix_, *solution_, *rhs_, rowBasis, colBasis);
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

  // 2. traverse grid and assemble operators on the elements
  for (auto const& element : elements(gridView())) {
    localView_->bind(element);
    auto geometry = element.geometry();

    if (asmMatrix)
      systemMatrix_->assemble(element, geometry, *localView_, *localView_);
    if (asmVector)
      rhs_->assemble(element, geometry, *localView_);

    localView_->unbind();
  }

  // 3. finish matrix insertion and apply dirichlet boundary conditions
  systemMatrix_->finish(asmMatrix);
  rhs_->finish(asmVector);

  forEachNode_(localView_->tree(), [&,this](auto const& rowNode, auto rowTreePath) {
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTreePath);
    forEachNode_(localView_->tree(), [&,this](auto const& colNode, auto colTreePath) {
      auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTreePath);
      // finish boundary condition
      for (auto bc : constraints_[rowNode][colNode])
        bc->finish(*systemMatrix_, *solution_, *rhs_, rowBasis, colBasis);
    });
  });
300

301
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
302
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
303
}
304

305

306
307
template <class Traits>
void ProblemStat<Traits>::
308
writeFiles(AdaptInfo& adaptInfo, bool force)
309
{
310
  Dune::Timer t;
311
  for (auto writer : filewriter_)
312
    writer->writeFiles(adaptInfo, force);
313
  msg("writeFiles needed {} seconds", t.elapsed());
314
}
315

316

317
} // end namespace AMDiS