You need to sign in or sign up before continuing.
Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
#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 = 0; 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_) {
fixedPointIteration++;
previousVelocityIterate = velocityIterate;
}
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
velocityUpdater->postProcess(velocityIterate);
return fixedPointIteration;
}
#include "fixedpointiterator_tmpl.cc"