ProblemStat.inc.hpp 13.9 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
      adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
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
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
  grid_ = MeshCreator<Grid>::create(gridName_);
126
  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
Praetorius, Simon's avatar
Praetorius, Simon committed
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

  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)
{
167 168
  dirichletBCs_.init(*globalBasis_, *globalBasis_);
  periodicBCs_.init(*globalBasis_, *globalBasis_);
Praetorius, Simon's avatar
Praetorius, Simon committed
169 170 171 172 173 174 175
}


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

Praetorius, Simon's avatar
Praetorius, Simon committed
179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202
  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");
}


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

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

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

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


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

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

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


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

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

261
  auto valueGridFct = makeGridFunction(values, globalBasis_->gridView());
262

263
  using Range = RangeType_t<decltype(i)>;
264
  using BcType = DirichletBC<WorldVector,Range>;
265
  auto bc = std::make_unique<BcType>(predicate, valueGridFct);
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
  dirichletBCs_[i][j].push_back(std::move(bc));
}


// Adds a Dirichlet boundary condition
template <class Traits>
  template <class RowTreePath, class ColTreePath, class Values>
void ProblemStat<Traits>::
addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& values)
{
  auto localView = globalBasis_->localView();
  auto i = child(localView.tree(), makeTreePath(row));
  auto j = child(localView.tree(), makeTreePath(col));

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

  using Range = RangeType_t<decltype(i)>;
  using BcType = DirichletBC<WorldVector,Range>;
  auto bc = std::make_unique<BcType>(boundaryManager_, id, valueGridFct);
  dirichletBCs_[i][j].push_back(std::move(bc));
}


template <class Traits>
void ProblemStat<Traits>::
addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
{
  auto localView = globalBasis_->localView();
  using BcType = PeriodicBC<WorldVector, typename GlobalBasis::MultiIndex>;
  using FaceTrafo = FaceTransformation<typename WorldVector::field_type, WorldVector::dimension>;
  auto bc = std::make_unique<BcType>(boundaryManager_, id, FaceTrafo{matrix, vector});
  periodicBCs_[localView.tree()][localView.tree()].push_back(std::move(bc));
298
}
299

300

301 302
template <class Traits>
void ProblemStat<Traits>::
303 304
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
305
  Dune::Timer t;
306

307
  SolverInfo solverInfo(name_ + "->solver");
308 309 310
  solverInfo.setCreateMatrixData(createMatrixData);
  solverInfo.setStoreMatrixData(storeMatrixData);

311
  solution_->compress();
312

313
  linearSolver_->solve(systemMatrix_->matrix(), solution_->vector(), rhs_->vector(),
314 315
                      solverInfo);

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

Praetorius, Simon's avatar
Praetorius, Simon committed
319 320
    if (solverInfo.absResidual() >= 0.0) {
      if (solverInfo.relResidual() >= 0.0)
321
        msg("Residual norm: ||b-Ax|| = {}, ||b-Ax||/||b|| = {}",
Praetorius, Simon's avatar
Praetorius, Simon committed
322
          solverInfo.absResidual(), solverInfo.relResidual());
323
      else
Praetorius, Simon's avatar
Praetorius, Simon committed
324
        msg("Residual norm: ||b-Ax|| = {}", solverInfo.absResidual());
325 326
    }
  }
327

328 329
  if (solverInfo.doBreak()) {
    std::stringstream tol_str;
Praetorius, Simon's avatar
Praetorius, Simon committed
330 331 332 333 334 335
    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() << " ";
336
    error_exit("Tolerance {} could not be reached!", tol_str.str());
337 338
  }
}
339

340

341
template <class Traits>
Praetorius, Simon's avatar
Praetorius, Simon committed
342 343
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
344 345 346 347
{
  Dune::Timer t;

  Flag markFlag = 0;
348
  for (auto& currentMarker : marker_)
349
    markFlag |= currentMarker->markGrid(adaptInfo);
350

351
  msg("markElements needed {} seconds", t.elapsed());
352 353 354 355 356

  return markFlag;
}


357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407
template <class Traits>
Flag ProblemStat<Traits>::
globalCoarsen(int n)
{
  Dune::Timer t;
  bool adapted = false;
  for (int i = 0; i < n; ++i) {
    // mark all entities for grid refinement
    for (const auto& element : elements(grid_->leafGridView()))
      grid_->mark(-1, element);

    adapted |= GridTransferManager::adapt(*grid_);
    globalBasis_->update(gridView());
  }

  msg("globalCoarsen needed {} seconds", t.elapsed());
  return adapted ? MESH_ADAPTED : Flag(0);
}


// grid has globalRefine(int, AdaptDataHandleInterface&)
template <class G>
using HasGlobalRefineADHI = decltype(std::declval<G>().globalRefine(1,std::declval<GridTransfer<G>&>()));

template <class Traits>
Flag ProblemStat<Traits>::
globalRefine(int refCount)
{
  Dune::Timer t;
  bool adapted = false;
  Dune::Hybrid::ifElse(Dune::Std::is_detected<HasGlobalRefineADHI, Grid>{},
  /*then*/ [&](auto id) {
    id(grid_)->globalRefine(refCount, GridTransferManager::gridTransfer(*grid_));
    globalBasis_->update(this->gridView());
  },
  /*else*/ [&](auto id) {
    for (int i = 0; i < refCount; ++i) {
      // mark all entities for grid refinement
      for (const auto& element : elements(grid_->leafGridView()))
        grid_->mark(1, element);

      adapted |= GridTransferManager::adapt(*id(grid_));
      globalBasis_->update(this->gridView());
    }
  });

  msg("globalRefine needed {} seconds", t.elapsed());
  return adapted ? MESH_ADAPTED : Flag(0);
}


408 409 410 411 412 413
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

414 415
  bool adapted = GridTransferManager::adapt(*grid_);
  globalBasis_->update(gridView());
416

417
  msg("adaptGrid needed {} seconds", t.elapsed());
418
  return adapted ? MESH_ADAPTED : Flag(0);
419 420 421
}


422 423
template <class Traits>
void ProblemStat<Traits>::
424
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
425
{
426
  Dune::Timer t;
427

428 429 430 431
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  systemMatrix_->init(asmMatrix);
  rhs_->init(asmVector);

432
  auto localView = globalBasis_->localView();
433 434 435 436 437 438 439 440
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto rowTp) {
    auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTp);
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto colTp) {
      auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTp);
      for (auto bc : dirichletBCs_[rowNode][colNode])
        bc->init(rowBasis, colBasis);
      for (auto bc : periodicBCs_[rowNode][colNode])
        bc->init(rowBasis, colBasis);
441 442 443 444 445 446
    });
  });
  msg("{} total DOFs", globalBasis_->dimension());

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

    if (asmMatrix)
450
      systemMatrix_->assemble(localView, localView);
451
    if (asmVector)
452
      rhs_->assemble(localView);
453

454
    localView.unbind();
455 456 457 458 459 460
  }

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

461 462
  forEachNode_(localView.tree(), [&,this](auto const& rowNode, auto) {
    forEachNode_(localView.tree(), [&,this](auto const& colNode, auto) {
463
      // finish boundary condition
464 465 466 467
      for (auto bc : dirichletBCs_[rowNode][colNode])
        bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
      for (auto bc : periodicBCs_[rowNode][colNode])
        bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, colNode);
468 469
    });
  });
470

471
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
472
  msg("buildAfterAdapt needed {} seconds", t.elapsed());
473
}
474

475

476 477
template <class Traits>
void ProblemStat<Traits>::
478
writeFiles(AdaptInfo& adaptInfo, bool force)
479
{
480
  Dune::Timer t;
481
  for (auto writer : filewriter_)
482
    writer->writeFiles(adaptInfo, force);
483
  msg("writeFiles needed {} seconds", t.elapsed());
484
}
485

486

487
} // end namespace AMDiS