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>
  template <class OperatorType>
void ProblemStat<Traits>::
addMatrixOperator(OperatorType const& op, int i, int j, double* factor, double* estFactor)
{
  matrixOperators[i][j].push_back({std::make_shared<OperatorType>(op), factor, estFactor});
  matrixChanging[i][j] = true;
}

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

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

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


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

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

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

template <class Traits>
void ProblemStat<Traits>::
addVectorOperator(std::map< int, std::pair<std::shared_ptr<ElementOperator>, bool> > ops)
{
  for (auto op : ops) {
    vectorOperators[op.first].push_back({op.second.first});
    vectorChanging[op.first] = vectorChanging[op.first] || 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>;
  dirichletBc[row][col].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;

  // update global feSpace, i.e. necessary after mesh change
  forEach(range_<0, nComponents>, [this](auto const _r) {
    this->getFeSpace(_r).update(this->leafGridView());
  });


  std::size_t nnz = 0;
  forEach(range_<0, nComponents>, [&,this](auto const _r)
  {
    static const int R = decltype(_r)::value;
    msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", R, "]");

    bool asmVector = asmVector_ && (!vectorAssembled[R] || vectorChanging[R]);

    if (asmVector) {
      rhs->compress(_r);
      rhs->getVector(_r) = 0.0;
    }

    forEach(range_<0, nComponents>, [&,this](auto const _c)
    {
      static const int C = decltype(_c)::value;

      using MatrixData = typename ProblemStat<Traits>::template MatrixData<R, C>;
      using VectorData = typename ProblemStat<Traits>::template VectorData<R>;

      // The DOFMatrix which should be assembled
      auto& dofmatrix    = (*systemMatrix)(_r, _c);
      auto& solution_vec = (*solution)[_c];
      auto& rhs_vec      = (*rhs)[_r];

      bool assembleMatrix = asmMatrix_ && (!matrixAssembled[R][C] || matrixChanging[R][C]);
      bool assembleVector = asmVector  && R == C;

      if (assembleMatrix) {
        // init boundary condition
        for (int c = 0; c < nComponents; ++c)
          for (auto bc : dirichletBc[R][c])
            bc->init(c == C, dofmatrix, solution_vec, rhs_vec);
      }

      auto mat = MatrixData{dofmatrix, matrixOperators[R][C], matrixBoundaryOperators[R][C], matrixIntersectionOperators[R][C], assembleMatrix};
      auto vec = VectorData{rhs_vec, vectorOperators[R], vectorBoundaryOperators[R], vectorIntersectionOperators[R], assembleVector};

      // assemble the DOFMatrix block and the corresponding rhs vector if r==c

      dofmatrix.init(assembleMatrix);
      this->assemble(mat, vec, this->leafGridView());
      dofmatrix.finish();

      if (assembleMatrix)
        matrixAssembled[R][C] = true;
      if (assembleVector)
        vectorAssembled[R] = true;

      if (assembleMatrix) {
        // finish boundary condition
        for (int c = 0; c < nComponents; ++c) {
          for (int r = 0; r < nComponents; ++r) {
            if (r != R && c != C)
              continue;
            for (auto bc : dirichletBc[r][c])
              bc->finish(r == R, c == C, dofmatrix, solution_vec, rhs_vec);
          }
        }

        nnz += dofmatrix.getNnz();
      }
    });
  });

  msg("fillin of assembled matrix: ", nnz);
  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");
}


template <class Traits>
  template <class LhsData, class RhsData, class GV>
void ProblemStat<Traits>::
assemble(LhsData lhs, RhsData rhs, GV const& gridView)
{
  auto const& rowFeSpace = lhs.matrix.getRowFeSpace();
  auto const& colFeSpace = lhs.matrix.getColFeSpace();

  if ((lhs.operators.empty() || !lhs.assemble) &&
      (rhs.operators.empty() || !rhs.assemble))
    return; // nothing to do

  for (auto scaledOp : lhs.operators)
    scaledOp.op->init(rowFeSpace, colFeSpace);
  for (auto scaledOp : rhs.operators)
    scaledOp.op->init(rowFeSpace, colFeSpace);

  auto rowLocalView = rowFeSpace.localView();
  auto rowIndexSet  = rowFeSpace.localIndexSet();

  auto colLocalView = colFeSpace.localView();
  auto colIndexSet  = colFeSpace.localIndexSet();

  for (auto const& element : elements(gridView)) {
    // TODO: use only one localView and localIndexSet if feSpaces are equal
    rowLocalView.bind(element);
    colLocalView.bind(element);
    rowIndexSet.bind(rowLocalView); // NOTE: expensive operation!
    colIndexSet.bind(colLocalView);

    if (lhs.assemble) {
      ElementMatrix elementMatrix;
      bool add = getElementMatrix(rowLocalView, colLocalView, elementMatrix,
                    lhs.operators, lhs.boundary_operators, lhs.intersection_operators);
      if (add)
        addElementMatrix(lhs.matrix, rowIndexSet, colIndexSet, elementMatrix);
    }

    if (rhs.assemble) {
      ElementVector elementVector;
      bool add = getElementVector(rowLocalView, elementVector,
                    rhs.operators, rhs.boundary_operators, rhs.intersection_operators);
      if (add)
        addElementVector(rhs.vector, rowIndexSet, elementVector);
    }

    rowIndexSet.unbind();
    rowLocalView.unbind();
    colLocalView.unbind();
    colIndexSet.unbind();
  }
}


template <class Traits>
  template <class RowView, class ColView>
bool ProblemStat<Traits>::
getElementMatrix(RowView const& rowLocalView,
                 ColView const& colLocalView,
                 ElementMatrix& elementMatrix,
                 std::list<Scaled<ElementOperator>>& operators,
                 std::list<Scaled<IntersectionOperator>>& boundary_operators,
                 std::list<Scaled<IntersectionOperator>>& intersection_operators)
{
  auto const nRows = rowLocalView.tree().finiteElement().size();
  auto const nCols = colLocalView.tree().finiteElement().size();

  auto const& element = rowLocalView.element();
  auto const& gridView = rowLocalView.globalBasis().gridView();

  // fills the entire matrix with zeroes
  elementMatrix.change_dim(nRows, nCols);
  set_to_zero(elementMatrix);

  bool add = false;

  auto assemble_operators = [&](auto& e, auto& operator_list) {
    for (auto scaled : operator_list) {
      bool add_op = scaled.op->getElementMatrix(e, rowLocalView, colLocalView, elementMatrix, scaled.factor);
      add = add || add_op;
    }
  };

  // assemble element operators
  assemble_operators(element, operators);

  // assemble intersection operators
  if (!intersection_operators.empty()
      || (!boundary_operators.empty() && element.hasBoundaryIntersections()))
  {
    for (auto const& intersection : intersections(gridView, element)) {
      if (intersection.boundary())
        assemble_operators(intersection, boundary_operators);
      else
        assemble_operators(intersection, intersection_operators); // TODO: assemble intersection operators also on the boundary?
    }
  }

  return add;
}


template <class Traits>
  template <class RowView>
bool ProblemStat<Traits>::
getElementVector(RowView const& rowLocalView,
                 ElementVector& elementVector,
                 std::list<Scaled<ElementOperator>>& operators,
                 std::list<Scaled<IntersectionOperator>>& boundary_operators,
                 std::list<Scaled<IntersectionOperator>>& intersection_operators)
{
  auto const nRows = rowLocalView.tree().finiteElement().size();

  auto const& element = rowLocalView.element();
  auto const& gridView = rowLocalView.globalBasis().gridView();

  // Set all entries to zero
  elementVector.change_dim(nRows);
  set_to_zero(elementVector);

  bool add = false;

  auto assemble_operators = [&](auto& e, auto& operator_list) {
    for (auto scaled : operator_list) {
      bool add_op = scaled.op->getElementVector(e, rowLocalView, elementVector, scaled.factor);
      add = add || add_op;
    }
  };

  // assemble element operators
  assemble_operators(element, operators);

  // assemble intersection operators
  if (!intersection_operators.empty()
      || (!boundary_operators.empty() && element.hasBoundaryIntersections()))
  {
    for (auto const& intersection : intersections(gridView, element)) {
      if (intersection.boundary())
        assemble_operators(intersection, boundary_operators);
      else
        assemble_operators(intersection, intersection_operators); // TODO: assemble intersection operators also on the boundary?
    }
  }

  return add;
}


template <class Traits>
  template <class Matrix, class RowIndexSet, class ColIndexSet>
void ProblemStat<Traits>::
addElementMatrix(Matrix& dofmatrix,
                 RowIndexSet const& rowIndexSet,
                 ColIndexSet const& colIndexSet,
                 ElementMatrix const& elementMatrix)
{
  for (std::size_t i = 0; i < num_rows(elementMatrix); ++i) {
    // The global index of the i−th vertex of the element
    auto const row = rowIndexSet.index(i);
    for (std::size_t j = 0; j < num_cols(elementMatrix); ++j) {
      // The global index of the j−th vertex of the element
      auto const col = colIndexSet.index(j);
      dofmatrix(row, col) += elementMatrix[i][j];
    }
  }
}


template <class Traits>
  template <class Vector, class IndexSet>
void ProblemStat<Traits>::
addElementVector(Vector& dofvector,
                 IndexSet const& indexSet,
                 ElementVector const& elementVector)
{
  for (std::size_t i = 0; i < size(elementVector); ++i) {
    // The global index of the i-th vertex
    auto const row = indexSet.index(i);
    dofvector[row] += elementVector[i];
  }
}

} // end namespace AMDiS