fixedpointiterator.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/common/exceptions.hh>

#include <dune/solvers/common/arithmetic.hh>
#include <dune/solvers/solvers/loopsolver.hh>

#include "enums.hh"
#include "enumparser.hh"

#include "fixedpointiterator.hh"

template <class Factory, class StateUpdater, class VelocityUpdater>
FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::FixedPointIterator(
    Factory &factory, Dune::ParameterTree const &parset,
    std::shared_ptr<Nonlinearity> globalFriction)
    : factory_(factory),
      parset_(parset),
      globalFriction_(globalFriction),
      fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
      fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
      lambda_(parset.get<double>("v.fpi.lambda")),
      velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
      velocityTolerance_(parset.get<double>("v.solver.tolerance")),
      verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")) {}

template <class Factory, class StateUpdater, class VelocityUpdater>
int FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::run(
    std::shared_ptr<StateUpdater> stateUpdater,
    std::shared_ptr<VelocityUpdater> velocityUpdater,
    Matrix const &velocityMatrix, Norm<Vector> const &velocityMatrixNorm,
    Vector const &velocityRHS, Vector &velocityIterate) {
  auto multigridStep = factory_.getSolver();

  LoopSolver<Vector> velocityProblemSolver(
      multigridStep, velocityMaxIterations_, velocityTolerance_,
      &velocityMatrixNorm, verbosity_, false); // absolute error

  Vector previousVelocityIterate = velocityIterate;

  size_t fixedPointIteration;
  for (fixedPointIteration = 1; fixedPointIteration <= fixedPointMaxIterations_;
       ++fixedPointIteration) {
    Vector v_m;
    velocityUpdater->extractOldVelocity(v_m);
    v_m *= 1.0 - lambda_;
    Arithmetic::addProduct(v_m, lambda_, velocityIterate);

    // solve a state problem
    stateUpdater->solve(v_m);
    ScalarVector alpha;
    stateUpdater->extractAlpha(alpha);

    // solve a velocity problem
    globalFriction_->updateAlpha(alpha);
    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                velocityRHS, velocityIterate);
    BlockProblem velocityProblem(parset_, convexProblem);
    multigridStep->setProblem(velocityIterate, velocityProblem);
    velocityProblemSolver.preprocess();
    velocityProblemSolver.solve();

    if (velocityMatrixNorm.diff(previousVelocityIterate, velocityIterate) <
        fixedPointTolerance_)
      break;

    previousVelocityIterate = velocityIterate;
  }
  if (fixedPointIteration == fixedPointMaxIterations_)
    DUNE_THROW(Dune::Exception, "FPI failed to converge");

  velocityUpdater->postProcess(velocityIterate);

  return fixedPointIteration;
}

#include "fixedpointiterator_tmpl.cc"