Liebe Gitlab-Nutzer, lieber Gitlab-Nutzer, es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Ein Anmelden über dieses erzeugt ein neues Konto. Das alte Konto ist über den Reiter "Standard" erreichbar. Die Administratoren

Dear Gitlab user, it is now possible to log in to our service using the ZIH login/LDAP. Logging in via this will create a new account. The old account can be accessed via the "Standard" tab. The administrators

ProblemStat.inc.hpp 10.8 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
#include <amdis/AdaptInfo.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/GridFunctionOperator.hpp>
13
#include <amdis/GridTransferManager.hpp>
14
#include <amdis/LocalAssembler.hpp>
15
#include <amdis/common/Loops.hpp>
16

17 18
namespace AMDiS {

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

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

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

47 48 49 50 51 52
  if (initFlag.isSet(INIT_MESH)) {
    int globalRefinements = 0;
    Parameters::get(gridName_ + "->global refinements", globalRefinements);
    if (globalRefinements > 0) {
      grid_->globalRefine(globalRefinements);
    }
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))) {
Praetorius, Simon's avatar
Praetorius, Simon committed
67
      adoptGlobalBasis(adoptProblem->globalBasis_);
68
    }
69
  }
70

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


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

79
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
Praetorius, Simon's avatar
Praetorius, Simon committed
80
    systemMatrix_ = adoptProblem->systemMatrix_;
81 82
    solution_ = adoptProblem->solution_;
    rhs_ = adoptProblem->rhs_;
Praetorius, Simon's avatar
Praetorius, Simon committed
83
    estimates_ = adoptProblem->estimates_;
84
  }
85

86

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

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

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

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

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

112

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

117
  solution_->compress();
118
}
119

120

Praetorius, Simon's avatar
Praetorius, Simon committed
121 122 123 124 125 126 127 128 129 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 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
  grid_ = MeshCreator<Grid>::create(gridName_);

  msg("Create grid:");
  msg("#elements = {}"   , grid_->size(0));
  msg("#faces/edges = {}", grid_->size(1));
  msg("#vertices = {}"   , grid_->size(dim));
  msg("");
}


template <class T, class GV>
using HasCreate = decltype(T::create(std::declval<GV>()));


template <class Traits>
void ProblemStat<Traits>::createGlobalBasis()
{
  createGlobalBasisImpl(Dune::Std::is_detected<HasCreate,Traits,GridView>{});
  initGlobalBasis(*globalBasis_);
}


template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type)
{
  assert( bool(grid_) );
  static_assert(std::is_same<GridView, typename Grid::LeafGridView>::value, "");
  globalBasis_ = std::make_shared<GlobalBasis>(Traits::create(grid_->leafGridView()));
}


template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::false_type)
{
  error_exit("Cannot create GlobalBasis from type. Pass a BasisCreator instead!");
}


template <class Traits>
void ProblemStat<Traits>::initGlobalBasis(GlobalBasis const& globalBasis)
{
  constraints_.init(*globalBasis_, *globalBasis_);
}


template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
  systemMatrix_ = std::make_shared<SystemMatrix>(*globalBasis_, *globalBasis_);
174 175 176
  solution_ = std::make_shared<SystemVector>(*globalBasis_, INTERPOLATE);
  rhs_ = std::make_shared<SystemVector>(*globalBasis_, NO_OPERATION);

Praetorius, Simon's avatar
Praetorius, Simon committed
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
  auto localView = globalBasis_->localView();
  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
  {
    std::string i = to_string(treePath);
    estimates_[i].resize(globalBasis_->gridView().indexSet().size(0));
    for (std::size_t j = 0; j < estimates_[i].size(); j++)
      estimates_[i][j] = 0.0; // TODO: Remove when estimate() is implemented
  });
}


template <class Traits>
void ProblemStat<Traits>::createSolver()
{
  std::string solverName = "default";
  Parameters::get(name_ + "->solver->name", solverName);

  auto solverCreator
    = named(CreatorMap<LinearSolverType>::getCreator(solverName, name_ + "->solver->name"));

  linearSolver_ = solverCreator->create(name_ + "->solver");
}


201 202 203
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
Praetorius, Simon's avatar
Praetorius, Simon committed
204
  marker_.clear();
205 206
  auto localView = globalBasis_->localView();
  AMDiS::forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
207
  {
208
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
209 210 211 212

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

213
    std::string tp = to_string(treePath);
Praetorius, Simon's avatar
Praetorius, Simon committed
214 215
    auto newMarker
      = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
216
    if (newMarker)
217 218
      marker_.push_back(std::move(newMarker));

219 220
    // If there is more than one marker, and all components are defined
    // on the same grid, the maximum marking has to be enabled.
221 222 223
    // TODO: needs to be checked!
    if (marker_.size() > 1)
      marker_.back()->setMaximumMarking(true);
224 225 226 227
  });
}


228 229
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
230
{
Praetorius, Simon's avatar
Praetorius, Simon committed
231
  filewriter_.clear();
232 233
  auto localView = globalBasis_->localView();
  forEachNode_(localView.tree(), [&,this](auto const& node, auto treePath)
234
  {
235
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
236 237 238 239

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

Praetorius, Simon's avatar
Praetorius, Simon committed
240
    auto writer = makeFileWriterPtr(componentName, this->solution(treePath));
241
    filewriter_.push_back(std::move(writer));
242 243 244 245
  });
}


246
// Adds a Dirichlet boundary condition
247
template <class Traits>
248
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
249
void ProblemStat<Traits>::
250
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
251 252
{
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
253
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
254

255 256 257
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));
258

259
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
260

261
  using Range = RangeType_t<decltype(i)>;
262
  using BcType = DirichletBC<WorldVector,Range>;
263 264
  auto bc = std::make_unique<BcType>(predicate, valueGridFct);
  constraints_[i][j].push_back(std::move(bc));
265
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
266
}
267

268

269 270
template <class Traits>
void ProblemStat<Traits>::
271 272
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
273
  Dune::Timer t;
274

275
  SolverInfo solverInfo(name_ + "->solver");
276 277 278
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

279
  solution_->compress();
280

281
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
282 283
                      solverInfo);

Praetorius, Simon's avatar
Praetorius, Simon committed
284
  if (solverInfo.info() > 0) {
285
    msg("solution of discrete system needed {} seconds", t.elapsed());
286

Praetorius, Simon's avatar
Praetorius, Simon committed
287 288
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
289
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
290
          solverInfo.absResidual(), solverInfo.relResidual());
291
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
292
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
293 294
    }
  }
295

296 297
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
Praetorius, Simon's avatar
Praetorius, Simon committed
298 299 300 301 302 303
    if (solverInfo.absTolerance() > 0
        && solverInfo.absResidual() > solverInfo.absTolerance())
      tol_str << "absTol = " << solverInfo.absTolerance() << " ";
    if (solverInfo.relTolerance() > 0
        && solverInfo.relResidual() > solverInfo.relTolerance())
      tol_str << "relTol = " << solverInfo.relTolerance() << " ";
304
    error_exit("Tolerance {} could not be reached!", tol_str.str());
305 306
  }
}
307

308

309
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
310 311
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
312 313 314 315
{
  Dune::Timer t;

  Flag markFlag = 0;
316
  for (auto& currentMarker : marker_)
317
    markFlag |= currentMarker->markGrid(adaptInfo);
318

319
  msg("markElements needed {} seconds", t.elapsed());
320 321 322 323 324

  return markFlag;
}


325 326 327 328 329 330
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

331 332
  bool adapted = GridTransferManager::adapt(*grid_);
  globalBasis_->update(gridView());
333

334
  msg("adaptGrid needed {} seconds", t.elapsed());
335
  return adapted ? MESH_ADAPTED : Flag(0);
336 337 338
}


339 340
template <class Traits>
void ProblemStat<Traits>::
341
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
342
{
343
  Dune::Timer t;
344

345 346 347 348
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

349 350 351
  auto localView = globalBasis_->localView();
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
352
      for (auto bc : constraints_[rowNode][colNode])
353
        bc->init(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
354 355 356 357 358 359
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

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

    if (asmMatrix)
363
      systemMatrix_->assemble(localView, localView);
364
    if (asmVector)
365
      rhs_->assemble(localView);
366

367
    localView.unbind();
368 369 370 371 372 373
  }

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

374 375
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
376 377
      // finish boundary condition
      for (auto bc : constraints_[rowNode][colNode])
378
        bc->finish(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
379 380
    });
  });
381

382
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
383
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
384
}
385

386

387 388
template <class Traits>
void ProblemStat<Traits>::
389
writeFiles(AdaptInfo& adaptInfo, bool force)
390
{
391
  Dune::Timer t;
392
  for (auto writer : filewriter_)
393
    writer->writeFiles(adaptInfo, force);
394
  msg("writeFiles needed {} seconds", t.elapsed());
395
}
396

397

398
} // end namespace AMDiS