#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"