ProblemStat.inc.hpp

#pragma once

#include <map>
#include <string>
#include <utility>

#include <dune/common/timer.hh>
#include <dune/typetree/childextraction.hh>

#include <amdis/AdaptInfo.hpp>
#include <amdis/Assembler.hpp>
#include <amdis/FileWriter.hpp>
#include <amdis/LocalAssembler.hpp>
#include <amdis/GridFunctionOperator.hpp>
#include <amdis/common/Loops.hpp>
// #include <amdis/Estimator.hpp>

namespace AMDiS {

template <class Traits>
void ProblemStat<Traits>::initialize(
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
{
  // create grids
  if (grid_) {
    warning("grid already created");
  }
  else {
    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))) {
      grid_ = adoptProblem->grid_;
    }
  }

  if (!grid_)
    warning("no grid created");

  int globalRefinements = 0;
  Parameters::get(gridName_ + "->global refinements", globalRefinements);
  if (globalRefinements > 0) {
    grid_->globalRefine(globalRefinements);
  }


  // create fespace
  if (globalBasis_) {
    warning("globalBasis already created");
  }
  else {
    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))) {
      globalBasis_ = adoptProblem->globalBasis_;
      initGlobalBasis(*globalBasis_);
    }
  }

  if (!globalBasis_)
    warning("no globalBasis created\n");


  // create system
  if (initFlag.isSet(INIT_SYSTEM))
    createMatricesAndVectors();

  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    solution_ = adoptProblem->solution_;
    estimates = adoptProblem->estimates;
    rhs_ = adoptProblem->rhs_;
    systemMatrix_ = adoptProblem->systemMatrix_;
  }


  // create solver
  if (linearSolver_) {
    warning("solver already created\n");
  }
  else {
    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 estimator
  estimator.resize(nComponents, NULL);

  if (initFlag.isSet(INIT_ESTIMATOR))
    createEstimator();

  if (adoptProblem && adoptFlag.isSet(INIT_ESTIMATOR))
    estimator = adoptProblem->estimator;
*/

  // 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_->compress();
}


template <class Traits>
void ProblemStat<Traits>::createMarker()
{
  auto localView = globalBasis->localView();
  forEachNode(localView.tree(), [&,this](auto const& node, auto treePath)
  {
    std::string componentName = name + "->marker[" + to_string(treePath) + "]";

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

    int i = std::stoi(to_string(treePath)); // TODO: To be removed
    // replace i with treePath once supported
    auto newMarker = Marker<Traits>::createMarker(
      componentName, std::to_string(i), estimates[std::to_string(i)], grid);
    if (newMarker)
    {
      marker.push_back(std::move(newMarker));
    }
    // If there is more than one marker, and all components are defined
    // on the same grid, the maximum marking has to be enabled.
    // TODO: What about two markers each for two grids?
    if (marker.size() > 1 && nGrids == 1)
 	    marker.back()->setMaximumMarking(true);
  });
}


template <class Traits>
void ProblemStat<Traits>::createFileWriter()
{
  AMDiS::forEachNode_(localView_->tree(), [&,this](auto const& node, auto treePath)
  {
    std::string componentName = name_ + "->output[" + to_string(treePath) + "]";

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

    auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
    filewriter_.push_back(std::move(writer));
  });
}


// add matrix/vector operator terms

template <class Traits>
  template <class Operator, class RowTreePath, class ColTreePath>
void ProblemStat<Traits>::addMatrixOperator(
    Operator const& preOp,
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
{
  static_assert( Concepts::PreTreePath<RowTreePath>,
      "row must be a valid treepath, or an integer/index-constant");
  static_assert( Concepts::PreTreePath<ColTreePath>,
      "col must be a valid treepath, or an integer/index-constant");

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

  auto op = makeGridOperator(preOp, globalBasis_->gridView());
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);

  matrixOperators_[i][j].element.push_back({localAssembler, factor, estFactor});
  matrixOperators_[i][j].changing = true;
}


template <class Traits>
  template <class Operator, class RowTreePath, class ColTreePath>
void ProblemStat<Traits>::addMatrixOperator(
    BoundaryType b,
    Operator const& preOp,
    RowTreePath row, ColTreePath col,
    double* factor, double* estFactor)
{
  static_assert( Concepts::PreTreePath<RowTreePath>,
      "row must be a valid treepath, or an integer/index-constant");
  static_assert( Concepts::PreTreePath<ColTreePath>,
      "col must be a valid treepath, or an integer/index-constant");

  auto i = child(localView_->tree(), makeTreePath(row));
  auto j = child(localView_->tree(), makeTreePath(col));
  using Intersection = typename GridView::Intersection;

  auto op = makeGridOperator(preOp, globalBasis_->gridView());
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);

  matrixOperators_[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
  matrixOperators_[i][j].changing = true;
}


template <class Traits>
  template <class Operator, class TreePath>
void ProblemStat<Traits>::addVectorOperator(
    Operator const& preOp,
    TreePath path,
    double* factor, double* estFactor)
{
  static_assert( Concepts::PreTreePath<TreePath>,
      "path must be a valid treepath, or an integer/index-constant");

  auto i = child(localView_->tree(), makeTreePath(path));

  auto op = makeGridOperator(preOp, globalBasis_->gridView());
  auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);

  rhsOperators_[i].element.push_back({localAssembler, factor, estFactor});
  rhsOperators_[i].changing = true;
}


template <class Traits>
  template <class Operator, class TreePath>
void ProblemStat<Traits>::addVectorOperator(
    BoundaryType b,
    Operator const& preOp,
    TreePath path,
    double* factor, double* estFactor)
{
  static_assert( Concepts::PreTreePath<TreePath>,
      "path must be a valid treepath, or an integer/index-constant");

  auto i = child(localView_->tree(), makeTreePath(path));
  using Intersection = typename GridView::Intersection;

  auto op = makeGridOperator(preOp, globalBasis_->gridView());
  auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);

  rhsOperators_[i].boundary.push_back({localAssembler, factor, estFactor, b});
  rhsOperators_[i].changing = true;
}


// 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_assert( Concepts::Functor<Predicate, bool(WorldVector)>,
    "Function passed to addDirichletBC for `predicate` does not model the Functor<bool(WorldVector)> concept");

  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_shared<BcType>(predicate, valueGridFct);
  constraints_[i][j].push_back(bc);
  // TODO: make DirichletBC an abstract class and add specialization with gridfunction type
}


template <class Traits>
void ProblemStat<Traits>::
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
  Dune::Timer t;

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

  solution_->compress();

  linearSolver_->solve(systemMatrix_->getMatrix(), solution_->getVector(), rhs_->getVector(),
                      solverInfo);

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

    if (solverInfo.getAbsResidual() >= 0.0) {
      if (solverInfo.getRelResidual() >= 0.0)
        msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual(),
                          ", ||b-Ax||/||b|| = ", solverInfo.getRelResidual());
      else
        msg("Residual norm: ||b-Ax|| = ", solverInfo.getAbsResidual());
    }
  }

  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() << " ";
    error_exit("Tolerance ", tol_str.str(), " could not be reached!");
  }
}


template <class Traits>
void ProblemStat<Traits>::
estimate(AdaptInfo& adaptInfo)
{
  Dune::Timer t;
/*
  for (std::size_t i = 0; i < nComponents; i++) {
    Estimator *scalEstimator = estimator[i];

    if (scalEstimator) {
      scalEstimator->estimate(adaptInfo->getTimestep());

      adaptInfo->setEstSum(scalEstimator->getErrorSum(), i);
      adaptInfo->setEstMax(scalEstimator->getErrorMax(), i);

      if (adaptInfo->getRosenbrockMode() == false) {
        adaptInfo->setTimeEstSum(scalEstimator->getTimeEst(), i);
        adaptInfo->setTimeEstMax(scalEstimator->getTimeEstMax(), i);
      }
	  }
  }

//#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
//    MPI::COMM_WORLD.Barrier();
//#endif
*/
  msg("estimation of the error needed ", t.elapsed(), "seconds\n");
}


template <class Traits>
Flag ProblemStat<Traits>::markElements(AdaptInfo& adaptInfo)
{
  Dune::Timer t;

  Flag markFlag = 0;
  for (auto currentMarker : marker)
    markFlag |= currentMarker->markGrid(adaptInfo);

  msg("markElements needed ", t.elapsed(), " seconds");

  return markFlag;
}


template <class Traits>
Flag ProblemStat<Traits>::refineMesh(AdaptInfo& adaptInfo)
{
  // TODO: actual element data, other components

  grid->preAdapt();

  std::map<typename Grid::LocalIdSet::IdType, double> vertexData;
//  std::map<typename Grid::LocalIdSet::IdType, double> elementData;
  const auto& gridView = globalBasis->gridView();
  const auto& indexSet = gridView.indexSet();
  const auto& idSet = grid->localIdSet();

  // Save data in container during adaptation
  for (const auto& v : vertices(gridView)) {
    vertexData[idSet.id(v)] = (*solution)[indexSet.index(v)];
  }/*
  for (const auto& e : elements(gridView)) {
    elementData[idSet.id(e)] = something[indexSet.index(e)];
  }*/

  Flag flag = grid->adapt();
/*
  // Unpack data
  solution->compress();
  for (const auto& v : vertices(gridView)) {
    (*solution)[indexSet.index(v)] = vertexData[idSet.id(v)];
  }/*
  for (const auto& e : elements(gridView)) {
    something[indexSet.index(e)] = elementData[idSet.id(e)];
  }*/

  grid->postAdapt();

  return flag;
}


template <class Traits>
void ProblemStat<Traits>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool asmVector)
{
  Dune::Timer t;

  Assembler<Traits> assembler(*globalBasis_, matrixOperators_, rhsOperators_, constraints_);
  auto gv = leafGridView();
  globalBasis_->update(gv);
  assembler.assemble(*systemMatrix_, *solution_, *rhs_, asmMatrix, asmVector);

  msg("buildAfterCoarsen needed ", t.elapsed(), " seconds");
}


template <class Traits>
void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
  Dune::Timer t;
  for (auto writer : filewriter_)
    writer->writeFiles(adaptInfo, force);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}


} // end namespace AMDiS