ProblemStat.inc.hpp

#pragma once
 
#include <map>
#include <string>
#include <utility>
 
#include <dune/common/hybridutilities.hh>
#include <dune/common/timer.hh>
#include <dune/functions/functionspacebases/subspacebasis.hh>
#include <dune/grid/common/capabilities.hh>
#include <dune/typetree/childextraction.hh>
 
#include <amdis/AdaptInfo.hpp>
#include <amdis/BackupRestore.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/functions/EntitySet.hpp>
#include <amdis/io/FileWriterCreator.hpp>
#include <amdis/linearalgebra/SymmetryStructure.hpp>
 
namespace AMDiS {
 
template <class Traits>
void ProblemStat<Traits>::initialize(
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
{
  // create grids
  if (!grid_) {
    if (initFlag.isSet(CREATE_MESH) ||
        (!adoptFlag.isSet(INIT_MESH) &&
        (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
      createGrid();
    }
 
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
      adoptGrid(adoptProblem->grid_, adoptProblem->boundaryManager_);
    }
  }
 
  if (!grid_)
    warning("no grid created");
 
  // create fespace
  if (!globalBasis_) {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
      createGlobalBasis();
    }
 
    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
      adoptGlobalBasis(adoptProblem->globalBasis_);
    }
  }
 
  if (!globalBasis_)
    warning("no globalBasis created\n");
 
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
 
  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    systemMatrix_ = adoptProblem->systemMatrix_;
    solution_ = adoptProblem->solution_;
    rhs_ = adoptProblem->rhs_;
    estimates_ = adoptProblem->estimates_;
  }
 
 
  // create solver
  if (!linearSolver_) {
    if (initFlag.isSet(INIT_SOLVER))
      createSolver();
 
    if (adoptProblem && adoptFlag.isSet(INIT_SOLVER)) {
      test_exit(!linearSolver_, "solver already created\n");
      linearSolver_ = adoptProblem->linearSolver_;
    }
  }
 
  if (!linearSolver_) {
    warning("no solver created\n");
  }
 
  // create marker
    if (initFlag.isSet(INIT_MARKER))
      createMarker();
 
    if (adoptProblem && adoptFlag.isSet(INIT_MARKER))
      marker_ = adoptProblem->marker_;
 
 
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
 
  solution_->resizeZero();
}
 
 
template <class Traits>
void ProblemStat<Traits>::
restore(Flag initFlag)
{
  std::string grid_filename = Parameters::get<std::string>(name_ + "->restore->grid").value();
  std::string solution_filename = Parameters::get<std::string>(name_ + "->restore->solution").value();
  test_exit(filesystem::exists(grid_filename), "Restore file '{}' not found.", grid_filename);
  test_exit(filesystem::exists(solution_filename), "Restore file '{}' not found.", solution_filename);
 
  // TODO(SP): implement BAckupRestore independent of wrapped grid
  using HostGrid = typename Grid::HostGrid;
 
  // restore grid from file
  if (Dune::Capabilities::hasBackupRestoreFacilities<HostGrid>::v)
    adoptGrid(std::shared_ptr<HostGrid>(Dune::BackupRestoreFacility<HostGrid>::restore(grid_filename)));
  else
    adoptGrid(std::shared_ptr<HostGrid>(BackupRestoreByGridFactory<HostGrid>::restore(grid_filename)));
 
  // create fespace
  if (initFlag.isSet(INIT_FE_SPACE) || initFlag.isSet(INIT_SYSTEM))
    createGlobalBasis();
 
  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();
 
  // create solver
  if (!linearSolver_ && initFlag.isSet(INIT_SOLVER))
    createSolver();
 
  // create marker
  if (initFlag.isSet(INIT_MARKER))
    createMarker();
 
  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();
 
  solution_->resize(*globalBasis_);
  solution_->restore(solution_filename);
}
 
 
template <class Traits>
void ProblemStat<Traits>::createGrid()
{
  Parameters::get(name_ + "->mesh", gridName_);
 
  MeshCreator<Grid> creator(gridName_);
  grid_ = creator.create();
 
  boundaryManager_ = std::make_shared<BoundaryManager<Grid>>(grid_);
  if (!creator.boundaryIds().empty())
    boundaryManager_->setBoundaryIds(creator.boundaryIds());
 
  info(3,"Create grid:");
  info(3,"#elements = {}"   , grid_->size(0));
  info(3,"#faces/edges = {}", grid_->size(1));
  info(3,"#vertices = {}"   , grid_->size(dim));
  info(3,"overlap-size = {}", grid_->leafGridView().overlapSize(0));
  info(3,"ghost-size = {}"  , grid_->leafGridView().ghostSize(0));
  info(3,"");
}
 
 
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();
}
 
 
template <class Traits>
void ProblemStat<Traits>::createGlobalBasisImpl(std::true_type)
{
  assert( bool(grid_) );
  static_assert(std::is_same_v<GridView, typename Grid::LeafGridView>, "");
  auto basis = Traits::create(name_, grid_->leafGridView());
  globalBasis_ = std::make_shared<GlobalBasis>(std::move(basis));
}
 
 
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() {}
 
 
template <class Traits>
void ProblemStat<Traits>::createMatricesAndVectors()
{
  systemMatrix_ = std::make_shared<SystemMatrix>(globalBasis_, globalBasis_);
  std::string symmetryStr = "unknown";
  Parameters::get(name_ + "->symmetry", symmetryStr);
  systemMatrix_->setSymmetryStructure(symmetryStr);
 
  solution_ = std::make_shared<SolutionVector>(globalBasis_);
  rhs_ = std::make_shared<SystemVector>(globalBasis_);
 
  auto localView = globalBasis_->localView();
  Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void
  {
    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", solverName);
 
  linearSolver_ = std::make_shared<LinearSolver>(solverName, name_ + "->solver");
}
 
 
template <class Traits>
void ProblemStat<Traits>::createMarker()
{
  marker_.clear();
  auto localView = globalBasis_->localView();
  Traversal::forEachNode(localView.tree(), [&,this](auto&&, auto treePath) -> void
  {
    std::string componentName = name_ + "->marker[" + to_string(treePath) + "]";
    auto strategy = Parameters::get<std::string>(componentName + "->strategy");
 
    if (!strategy && to_string(treePath).empty()) {
      // alternative for root treepath
      componentName = name_ + "->strategy";
      strategy = Parameters::get<std::string>(componentName + "->strategy");
    }
 
    if (!strategy)
      return;
 
    std::string tp = to_string(treePath);
    auto newMarker
      = EstimatorMarker<Grid>::createMarker(componentName, tp, estimates_[tp], grid_);
    assert(bool(newMarker));
    this->addMarker(std::move(newMarker));
  });
}
 
 
template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
  FileWriterCreator<SolutionVector> creator(solution_, boundaryManager_);
 
  filewriter_.clear();
  auto localView = globalBasis_->localView();
  Traversal::forEachNode(localView.tree(), [&](auto const& /*node*/, auto treePath) -> void
  {
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";
    auto format = Parameters::get<std::vector<std::string>>(componentName + "->format");
 
    if (!format && to_string(treePath).empty()) {
      // alternative for root treepath
      componentName = name_ + "->output";
      format = Parameters::get<std::vector<std::string>>(componentName + "->format");
    }
 
    if (!format)
      return;
 
    for (std::string const& type : format.value()) {
      auto writer = creator.create(type, componentName, treePath);
      if (writer)
        filewriter_.push_back(std::move(writer));
    }
  });
}
 
 
// Adds a Dirichlet boundary condition
template <class Traits>
  template <class Predicate, class RowTreePath, class ColTreePath, class Values>
void ProblemStat<Traits>::
addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Values const& values)
{
  static constexpr bool isValidPredicate = Concepts::Functor<Predicate, bool(WorldVector)>;
  static_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");
 
  static constexpr bool isValidTreePath =
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
  static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addDirichletBC!");
 
  if constexpr (isValidPredicate && isValidTreePath) {
    auto valueGridFct = makeGridFunction(values, this->gridView());
 
    constraints_.push_back(DirichletBC{
      globalBasis_, makeTreePath(row), makeTreePath(col), {predicate}, valueGridFct});
  }
}
 
 
// 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)
{
  static constexpr bool isValidTreePath =
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, RowTreePath> &&
    Concepts::ValidTreePath<typename GlobalBasis::LocalView::Tree, ColTreePath>;
  static_assert(isValidTreePath, "Invalid row and/or col treepath passed to addDirichletBC!");
 
  if constexpr (isValidTreePath) {
    auto valueGridFct = makeGridFunction(values, this->gridView());
 
    constraints_.push_back(DirichletBC{
      globalBasis_, makeTreePath(row), makeTreePath(col), {boundaryManager_, id}, valueGridFct});
  }
}
 
 
template <class Traits>
void ProblemStat<Traits>::
addPeriodicBC(BoundaryType id, WorldMatrix const& matrix, WorldVector const& vector)
{
  Traversal::forEachLeafNode(globalBasis_->localView().tree(), [&](auto&&, auto tp) {
    auto basis = Dune::Functions::subspaceBasis(*globalBasis_, tp);
    constraints_.push_back(makePeriodicBC(
      std::move(basis), {boundaryManager_, id}, {matrix, vector}));
  });
}
 
 
 
template <class Traits>
void ProblemStat<Traits>::
solve(AdaptInfo& /*adaptInfo*/, bool createMatrixData, bool storeMatrixData)
{
  Dune::Timer t;
  Dune::InverseOperatorResult stat;
 
  solution_->resize();
 
  if (createMatrixData)
    linearSolver_->init(systemMatrix_->impl());
 
  // solve the linear system
  linearSolver_->apply(solution_->impl(), rhs_->impl(), stat);
 
  if (!storeMatrixData)
    linearSolver_->finish();
 
  info(2, "solution of discrete system needed {} seconds", t.elapsed());
  info(1, "Residual norm: ||b-Ax|| = {}",
    residuum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
 
  if (stat.reduction >= 0.0)
    info(2, "Relative residual norm: ||b-Ax||/||b|| = {}", stat.reduction);
  else
    info(2, "Relative residual norm: ||b-Ax||/||b|| = {}",
      relResiduum(systemMatrix_->impl(),solution_->impl(), rhs_->impl()));
 
  bool ignoreConverged = false;
  Parameters::get(name_ + "->solver->ignore converged", ignoreConverged);
  test_exit(stat.converged || ignoreConverged, "Could not solver the linear system!");
}
 
 
template <class Traits>
Flag ProblemStat<Traits>::
markElements(AdaptInfo& adaptInfo)
{
  Dune::Timer t;
 
  Flag markFlag = 0;
  for (auto& currentMarker : marker_)
    markFlag |= currentMarker.second->markGrid(adaptInfo);
 
  // synchronize mark flag over processors
  int markFlagValue = int(markFlag);
  grid_->comm().max(&markFlagValue, 1);
  markFlag = Flag(markFlagValue);
 
  info(2, "markElements needed {} seconds", t.elapsed());
 
  return markFlag;
}
 
 
template <class Traits>
Flag ProblemStat<Traits>::
globalCoarsen(int n)
{
  Dune::Timer t;
  bool adapted = false;
  // TODO(FM): Find a less expensive alternative to the loop adaption
  for (int i = 0; i < n; ++i) {
    // mark all entities for coarsening
    for (const auto& element : elements(grid_->leafGridView()))
      grid_->mark(-1, element);
 
    bool adaptedInLoop = grid_->preAdapt();
    adaptedInLoop |= grid_->adapt();
    grid_->postAdapt();
    if (!adaptedInLoop)
      break;
    else
      adapted = true;
  }
 
  info(2, "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<typename G::ADHI&>()));
 
template <class Traits>
Flag ProblemStat<Traits>::
globalRefine(int n)
{
  Dune::Timer t;
  if constexpr (Dune::Std::is_detected<HasGlobalRefineADHI, Grid>::value)
    grid_->globalRefine(n, globalBasis_->globalRefineCallback());
  else
    grid_->globalRefine(n);
 
  info(2, "globalRefine needed {} seconds", t.elapsed());
  return n > 0 ? MESH_ADAPTED : Flag(0);
}
 
 
template <class Traits>
Flag ProblemStat<Traits>::
adaptGrid(AdaptInfo& /*adaptInfo*/)
{
  Dune::Timer t;
 
  bool adapted = grid_->preAdapt();
  adapted |= grid_->adapt();
  grid_->postAdapt();
 
  info(2, "adaptGrid needed {} seconds", t.elapsed());
  return adapted ? MESH_ADAPTED : Flag(0);
}
 
 
template <class Traits>
void ProblemStat<Traits>::
buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
{
  Dune::Timer t;
  Dune::Timer t2;
 
  // 0. initialize boundary condition and other constraints
  for (auto& bc : constraints_)
    bc.init();
 
  t2.reset();
 
  // 1. init matrix and rhs vector and initialize dirichlet boundary conditions
  if (asmMatrix)
    systemMatrix_->init();
  if (asmVector)
    rhs_->init(*globalBasis_, true);
 
  // statistic about system size
  if (gridView().comm().size() > 1)
    msg("{} local DOFs, {} global DOFs", rhs_->localSize(), rhs_->globalSize());
  else
    msg("{} local DOFs", rhs_->localSize());
 
  // 2. traverse grid and assemble operators on the elements
  auto localView = globalBasis_->localView();
  for (auto const& element : entitySet(*globalBasis_)) {
    localView.bind(element);
 
    if (asmMatrix)
      systemMatrix_->assemble(localView, localView);
    if (asmVector)
      rhs_->assemble(localView);
 
    localView.unbind();
  }
 
  // 3. finish matrix insertion and apply dirichlet boundary conditions
  if (asmMatrix)
    systemMatrix_->finish();
  if (asmVector)
    rhs_->finish();
 
  info(2,"  assemble operators needed {} seconds", t2.elapsed());
  t2.reset();
 
  solution_->resize(*globalBasis_);
 
  // 4. apply boundary condition and constraints to matrix, solution, and rhs
  for (auto& bc : constraints_)
    bc.apply(*systemMatrix_, *solution_, *rhs_);
 
  rhs_->finish();
  solution_->finish();
  info(2,"  assemble boundary conditions needed {} seconds", t2.elapsed());
 
  msg("fill-in of assembled matrix: {}", systemMatrix_->nnz());
  msg("assemble needed {} seconds", t.elapsed());
}
 
 
template <class Traits>
void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
  Dune::Timer t;
  for (auto writer : filewriter_)
    writer->write(adaptInfo, force);
  msg("writeFiles needed {} seconds", t.elapsed());
}
 
} // end namespace AMDiS