ProblemStat.inc.hpp

#pragma once

#include <dune/amdis/AdaptInfo.hpp>
#include <dune/amdis/common/Loops.hpp>
#include <dune/amdis/common/Timer.hpp>

namespace AMDiS {

template <class Traits>
void ProblemStat<Traits>::initialize(
    Flag initFlag,
    Self* adoptProblem,
    Flag adoptFlag)
{
  // create meshes
  if (mesh) {
    warning("mesh already created");
  }
  else {
    if (initFlag.isSet(CREATE_MESH) ||
        (!adoptFlag.isSet(INIT_MESH) &&
        (initFlag.isSet(INIT_SYSTEM) || initFlag.isSet(INIT_FE_SPACE)))) {
      createMesh();
    }

    if (adoptProblem &&
        (adoptFlag.isSet(INIT_MESH) ||
        adoptFlag.isSet(INIT_SYSTEM) ||
        adoptFlag.isSet(INIT_FE_SPACE))) {
      mesh = adoptProblem->getMesh();
    }

    componentMeshes.resize(nComponents, mesh.get());
  }

  if (!mesh)
    warning("no mesh created");

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


  // create fespace
  if (feSpaces) {
    warning("feSpaces already created");
  }
  else {
    if (initFlag.isSet(INIT_FE_SPACE) ||
        (initFlag.isSet(INIT_SYSTEM) && !adoptFlag.isSet(INIT_FE_SPACE))) {
      createFeSpaces();
    }

    if (adoptProblem &&
        (adoptFlag.isSet(INIT_FE_SPACE) || adoptFlag.isSet(INIT_SYSTEM))) {
      feSpaces = adoptProblem->getFeSpaces();
    }
  }

  if (!feSpaces)
    warning("no feSpaces created\n");


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

  if (adoptProblem && adoptFlag.isSet(INIT_SYSTEM)) {
    solution = adoptProblem->getSolution();
    rhs = adoptProblem->getRhs();
    systemMatrix = adoptProblem->getSystemMatrix();
  }


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

  if (!linearSolver) {
    warning("no solver created\n");
  }


  // create file writer
  if (initFlag.isSet(INIT_FILEWRITER))
    createFileWriter();

  solution->compress();
}


// add matrix/vector operator terms

template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(ElementOperator const& op, int i, int j, double* factor, double* estFactor)
{
  matrixOperators[i][j].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
{
  matrixOperators[i][j].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
  matrixOperators[i][j].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(BoundaryType b, IntersectionOperator const& op, int i, int j, double* factor, double* estFactor)
{
  matrixOperators[i][j].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  matrixOperators[i][j].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(std::map< std::pair<int,int>, ElementOperator> ops)
{
  for (auto op : ops){
    int r = op.first.first;
    int c = op.first.second;
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second)});
    matrixOperators[r][c].changing = true;
  }
}

template <class Traits>
void ProblemStat<Traits>::
addMatrixOperator(std::map< std::pair<int,int>, std::pair<ElementOperator,bool> > ops)
{
  for (auto op : ops) {
    int r = op.first.first;
    int c = op.first.second;
    matrixOperators[r][c].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    matrixOperators[r][c].changing = matrixOperators[r][c].changing || op.second.second;
  }
}


template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(ElementOperator const& op, int i, double* factor, double* estFactor)
{
  rhsOperators[i].element.push_back({std::make_shared<ElementOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(IntersectionOperator const& op, int i, double* factor, double* estFactor)
{
  rhsOperators[i].intersection.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor});
  rhsOperators[i].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(BoundaryType b, IntersectionOperator const& op, int i, double* factor, double* estFactor)
{
  rhsOperators[i].boundary.push_back({std::make_shared<IntersectionOperator>(op), factor, estFactor, b});
  rhsOperators[i].changing = true;
}

template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(std::map< int, ElementOperator> ops)
{
  for (auto op : ops) {
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second)});
    rhsOperators[op.first].changing = true;
  }
}

template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(std::map< int, std::pair<ElementOperator, bool> > ops)
{
  for (auto op : ops) {
    rhsOperators[op.first].element.push_back({std::make_shared<ElementOperator>(op.second.first)});
    rhsOperators[op.first].changing = rhsOperators[op.first].changing || op.second.second;
  }
}


// Adds a Dirichlet boundary condition
template <class Traits>
  template <class Predicate, class Values>
void ProblemStat<Traits>::
addDirichletBC(Predicate const& predicate, int row, int 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");

  static_assert( Concepts::Callable<Values, WorldVector>,
    "Function passed to addDirichletBC for values does not model the Callable<WorldVector> concept");

  test_exit(row >= 0 && row < nComponents, "row number out of range: ", row);
  test_exit(col >= 0 && col < nComponents, "col number out of range: ", col);

//     boundaryConditionSet = true;

  using BcType = DirichletBC<WorldVector>;
  matrixOperators[row][col].dirichlet.emplace_back( new BcType{predicate, values} );
}


template <class Traits>
void ProblemStat<Traits>::
solve(AdaptInfo& adaptInfo, bool createMatrixData, bool storeMatrixData)
{
  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>::
buildAfterCoarsen(AdaptInfo& /*adaptInfo*/, Flag flag, bool asmMatrix, bool asmVector)
{
  Timer t;

  Assembler<FeSpaces> assembler(getFeSpaces(), asmMatrix, asmVector);

  auto gridView = leafGridView();
  assembler.assemble(gridView,
                     *systemMatrix, *solution, *rhs,
                     matrixOperators, rhsOperators);

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


template <class Traits>
void ProblemStat<Traits>::
writeFiles(AdaptInfo& adaptInfo, bool force)
{
  Timer t;
  filewriter->write(adaptInfo.getTime(), *solution);
  msg("writeFiles needed ", t.elapsed(), " seconds");
}

} // end namespace AMDiS