From 3550423d6d3be58c46af9df81a3b23c650392ec4 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Mon, 8 Sep 2014 01:17:14 +0200
Subject: [PATCH] [Algorit] Use state energy norm to terminate
---
src/adaptivetimestepper.hh | 34 +++++++++++++-----------
src/coupledtimestepper.cc | 29 ++++++++++++---------
src/coupledtimestepper.hh | 5 +++-
src/coupledtimestepper_tmpl.cc | 9 +++++--
src/explicitvectors.hh | 1 +
src/fixedpointiterator.cc | 47 ++++++++++++++++++----------------
src/fixedpointiterator.hh | 7 +++--
src/fixedpointiterator_tmpl.cc | 9 +++++--
src/sand-wedge-data/parset.cfg | 2 +-
src/sand-wedge.cc | 7 ++---
10 files changed, 90 insertions(+), 60 deletions(-)
diff --git a/src/adaptivetimestepper.hh b/src/adaptivetimestepper.hh
index 0a034172..48634e36 100644
--- a/src/adaptivetimestepper.hh
+++ b/src/adaptivetimestepper.hh
@@ -5,7 +5,8 @@ std::pair<T1, T2> clonePair(std::pair<T1, T2> in) {
return { in.first->clone(), in.second->clone() };
}
-template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
+template <class Factory, class UpdaterPair, class ErrorNorm>
+class AdaptiveTimeStepper {
using StateUpdater = typename UpdaterPair::first_type::element_type;
using VelocityUpdater = typename UpdaterPair::second_type::element_type;
using Vector = typename Factory::Vector;
@@ -13,13 +14,14 @@ template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
using Nonlinearity = typename ConvexProblem::NonlinearityType;
using MyCoupledTimeStepper =
- CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>;
+ CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>;
public:
AdaptiveTimeStepper(
Factory &factory, Dune::ParameterTree const &parset,
std::shared_ptr<Nonlinearity> globalFriction, UpdaterPair ¤t,
std::function<void(double, Vector &)> externalForces,
+ ErrorNorm const &errorNorm,
std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine)
: finalTime_(parset.get<double>("problem.finalTime")),
relativeTime_(0.0),
@@ -31,10 +33,11 @@ template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
R1_(clonePair(current_)),
externalForces_(externalForces),
mustRefine_(mustRefine),
- iterationWriter_("iterations", std::fstream::out) {
- MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
- globalFriction_, R1_.first,
- R1_.second, externalForces_);
+ iterationWriter_("iterations", std::fstream::out),
+ errorNorm_(errorNorm) {
+ MyCoupledTimeStepper coupledTimeStepper(
+ finalTime_, factory_, parset_, globalFriction_, R1_.first, R1_.second,
+ errorNorm, externalForces_);
stepAndReport("R1", coupledTimeStepper, relativeTime_, relativeTau_);
iterationWriter_ << std::endl;
}
@@ -47,18 +50,18 @@ template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
while (relativeTime_ + relativeTau_ < 1.0) {
R2_ = clonePair(R1_);
{
- MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
- globalFriction_, R2_.first,
- R2_.second, externalForces_);
+ MyCoupledTimeStepper coupledTimeStepper(
+ finalTime_, factory_, parset_, globalFriction_, R2_.first,
+ R2_.second, errorNorm_, externalForces_);
stepAndReport("R2", coupledTimeStepper, relativeTime_ + relativeTau_,
relativeTau_);
}
UpdaterPair C = clonePair(current_);
{
- MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
- globalFriction_, C.first,
- C.second, externalForces_);
+ MyCoupledTimeStepper coupledTimeStepper(
+ finalTime_, factory_, parset_, globalFriction_, C.first, C.second,
+ errorNorm_, externalForces_);
stepAndReport("C", coupledTimeStepper, relativeTime_,
2.0 * relativeTau_);
@@ -81,9 +84,9 @@ template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
UpdaterPair F2 = clonePair(current_);
UpdaterPair F1;
{
- MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
- globalFriction_, F2.first,
- F2.second, externalForces_);
+ MyCoupledTimeStepper coupledTimeStepper(
+ finalTime_, factory_, parset_, globalFriction_, F2.first, F2.second,
+ errorNorm_, externalForces_);
stepAndReport("F1", coupledTimeStepper, relativeTime_,
relativeTau_ / 2.0);
@@ -135,4 +138,5 @@ template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
std::function<void(double, Vector &)> externalForces_;
std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine_;
std::fstream iterationWriter_;
+ ErrorNorm const &errorNorm_;
};
diff --git a/src/coupledtimestepper.cc b/src/coupledtimestepper.cc
index f76b47e1..1dd286ff 100644
--- a/src/coupledtimestepper.cc
+++ b/src/coupledtimestepper.cc
@@ -5,23 +5,28 @@
#include "coupledtimestepper.hh"
#include "fixedpointiterator.hh"
-template <class Factory, class StateUpdater, class VelocityUpdater>
-CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>::CoupledTimeStepper(
- double finalTime, Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction,
- std::shared_ptr<StateUpdater> stateUpdater,
- std::shared_ptr<VelocityUpdater> velocityUpdater,
- std::function<void(double, Vector &)> externalForces)
+template <class Factory, class StateUpdater, class VelocityUpdater,
+ class ErrorNorm>
+CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
+ CoupledTimeStepper(double finalTime, Factory &factory,
+ Dune::ParameterTree const &parset,
+ std::shared_ptr<Nonlinearity> globalFriction,
+ std::shared_ptr<StateUpdater> stateUpdater,
+ std::shared_ptr<VelocityUpdater> velocityUpdater,
+ ErrorNorm const &errorNorm,
+ std::function<void(double, Vector &)> externalForces)
: finalTime_(finalTime),
factory_(factory),
parset_(parset),
globalFriction_(globalFriction),
stateUpdater_(stateUpdater),
velocityUpdater_(velocityUpdater),
- externalForces_(externalForces) {}
+ externalForces_(externalForces),
+ errorNorm_(errorNorm) {}
-template <class Factory, class StateUpdater, class VelocityUpdater>
-int CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>::step(
+template <class Factory, class StateUpdater, class VelocityUpdater,
+ class ErrorNorm>
+int CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::step(
double relativeTime, double relativeTau) {
stateUpdater_->nextTimeStep();
velocityUpdater_->nextTimeStep();
@@ -38,8 +43,8 @@ int CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>::step(
stateUpdater_->setup(tau);
velocityUpdater_->setup(ell, tau, newRelativeTime, velocityRHS,
velocityIterate, velocityMatrix);
- FixedPointIterator<Factory, StateUpdater, VelocityUpdater> fixedPointIterator(
- factory_, parset_, globalFriction_);
+ FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>
+ fixedPointIterator(factory_, parset_, globalFriction_, errorNorm_);
auto const iterations =
fixedPointIterator.run(stateUpdater_, velocityUpdater_, velocityMatrix,
velocityRHS, velocityIterate);
diff --git a/src/coupledtimestepper.hh b/src/coupledtimestepper.hh
index 52ece50c..f8122e62 100644
--- a/src/coupledtimestepper.hh
+++ b/src/coupledtimestepper.hh
@@ -6,7 +6,8 @@
#include <dune/common/parametertree.hh>
-template <class Factory, class StateUpdater, class VelocityUpdater>
+template <class Factory, class StateUpdater, class VelocityUpdater,
+ class ErrorNorm>
class CoupledTimeStepper {
using Vector = typename Factory::Vector;
using Matrix = typename Factory::Matrix;
@@ -19,6 +20,7 @@ class CoupledTimeStepper {
std::shared_ptr<Nonlinearity> globalFriction,
std::shared_ptr<StateUpdater> stateUpdater,
std::shared_ptr<VelocityUpdater> velocityUpdater,
+ ErrorNorm const &errorNorm,
std::function<void(double, Vector &)> externalForces);
int step(double relativeTime, double relativeTau);
@@ -31,5 +33,6 @@ class CoupledTimeStepper {
std::shared_ptr<StateUpdater> stateUpdater_;
std::shared_ptr<VelocityUpdater> velocityUpdater_;
std::function<void(double, Vector &)> externalForces_;
+ ErrorNorm const &errorNorm_;
};
#endif
diff --git a/src/coupledtimestepper_tmpl.cc b/src/coupledtimestepper_tmpl.cc
index a7831212..33f448bd 100644
--- a/src/coupledtimestepper_tmpl.cc
+++ b/src/coupledtimestepper_tmpl.cc
@@ -4,6 +4,7 @@
#include <dune/common/function.hh>
+#include <dune/solvers/norms/energynorm.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/globalfriction.hh>
@@ -18,9 +19,13 @@
using Function = Dune::VirtualFunction<double, double>;
using Factory = SolverFactory<
- MY_DIM, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+ MY_DIM,
+ MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyVelocityUpdater = TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>;
-template class CoupledTimeStepper<Factory, MyStateUpdater, MyVelocityUpdater>;
+using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+template class CoupledTimeStepper<Factory, MyStateUpdater, MyVelocityUpdater,
+ ErrorNorm>;
diff --git a/src/explicitvectors.hh b/src/explicitvectors.hh
index 923a8c35..5f2d887f 100644
--- a/src/explicitvectors.hh
+++ b/src/explicitvectors.hh
@@ -11,4 +11,5 @@ using LocalMatrix = Dune::FieldMatrix<double, MY_DIM, MY_DIM>;
using Vector = Dune::BlockVector<LocalVector>;
using Matrix = Dune::BCRSMatrix<LocalMatrix>;
using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
#endif
diff --git a/src/fixedpointiterator.cc b/src/fixedpointiterator.cc
index d868eee0..23569d9f 100644
--- a/src/fixedpointiterator.cc
+++ b/src/fixedpointiterator.cc
@@ -13,10 +13,12 @@
#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)
+template <class Factory, class StateUpdater, class VelocityUpdater,
+ class ErrorNorm>
+FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
+ FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
+ std::shared_ptr<Nonlinearity> globalFriction,
+ ErrorNorm const &errorNorm)
: factory_(factory),
parset_(parset),
globalFriction_(globalFriction),
@@ -25,10 +27,12 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::FixedPointIterator(
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")) {}
+ verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
+ errorNorm_(errorNorm) {}
-template <class Factory, class StateUpdater, class VelocityUpdater>
-int FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::run(
+template <class Factory, class StateUpdater, class VelocityUpdater,
+ class ErrorNorm>
+int FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
std::shared_ptr<StateUpdater> stateUpdater,
std::shared_ptr<VelocityUpdater> velocityUpdater,
Matrix const &velocityMatrix, Vector const &velocityRHS,
@@ -39,21 +43,12 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::run(
LoopSolver<Vector> velocityProblemSolver(
multigridStep, velocityMaxIterations_, velocityTolerance_, &energyNorm,
verbosity_, false); // absolute error
- Vector previousVelocityIterate = velocityIterate;
size_t fixedPointIteration;
+ ScalarVector alpha;
+ stateUpdater->extractAlpha(alpha);
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_,
@@ -63,13 +58,21 @@ int FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::run(
velocityProblemSolver.preprocess();
velocityProblemSolver.solve();
- if (energyNorm.diff(previousVelocityIterate, velocityIterate) <
- fixedPointTolerance_) {
+ 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 newAlpha;
+ stateUpdater->extractAlpha(newAlpha);
+
+ if (errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
fixedPointIteration++;
break;
}
-
- previousVelocityIterate = velocityIterate;
+ alpha = newAlpha;
}
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
diff --git a/src/fixedpointiterator.hh b/src/fixedpointiterator.hh
index 16b3eec5..34a7d42f 100644
--- a/src/fixedpointiterator.hh
+++ b/src/fixedpointiterator.hh
@@ -8,7 +8,8 @@
#include <dune/solvers/norms/norm.hh>
#include <dune/solvers/solvers/solver.hh>
-template <class Factory, class StateUpdater, class VelocityUpdater>
+template <class Factory, class StateUpdater, class VelocityUpdater,
+ class ErrorNorm>
class FixedPointIterator {
using ScalarVector = typename StateUpdater::ScalarVector;
using Vector = typename Factory::Vector;
@@ -19,7 +20,8 @@ class FixedPointIterator {
public:
FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
- std::shared_ptr<Nonlinearity> globalFriction);
+ std::shared_ptr<Nonlinearity> globalFriction,
+ ErrorNorm const &errorNorm_);
int run(std::shared_ptr<StateUpdater> stateUpdater,
std::shared_ptr<VelocityUpdater> velocityUpdater,
@@ -37,5 +39,6 @@ class FixedPointIterator {
size_t velocityMaxIterations_;
double velocityTolerance_;
Solver::VerbosityMode verbosity_;
+ ErrorNorm const &errorNorm_;
};
#endif
diff --git a/src/fixedpointiterator_tmpl.cc b/src/fixedpointiterator_tmpl.cc
index 41d4f910..50040293 100644
--- a/src/fixedpointiterator_tmpl.cc
+++ b/src/fixedpointiterator_tmpl.cc
@@ -4,6 +4,7 @@
#include <dune/common/function.hh>
+#include <dune/solvers/norms/energynorm.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
#include <dune/tectonic/globalfriction.hh>
@@ -18,9 +19,13 @@
using Function = Dune::VirtualFunction<double, double>;
using Factory = SolverFactory<
- MY_DIM, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+ MY_DIM,
+ MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
Grid>;
using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
using MyVelocityUpdater = TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>;
-template class FixedPointIterator<Factory, MyStateUpdater, MyVelocityUpdater>;
+using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
+
+template class FixedPointIterator<Factory, MyStateUpdater, MyVelocityUpdater,
+ ErrorNorm>;
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
index 45a42435..52a760fe 100644
--- a/src/sand-wedge-data/parset.cfg
+++ b/src/sand-wedge-data/parset.cfg
@@ -57,7 +57,7 @@ maximumIterations = 100000
verbosity = quiet
[v.fpi]
-tolerance = 1e-8
+tolerance = 1e-5
maximumIterations = 10000
lambda = 0.5
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index e7273b23..aa323951 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -436,9 +436,10 @@ int main(int argc, char *argv[]) {
size_t timeStep = 1;
- AdaptiveTimeStepper<NonlinearFactory, UpdaterPair> adaptiveTimeStepper(
- factory, parset, myGlobalFriction, current, computeExternalForces,
- mustRefine);
+ AdaptiveTimeStepper<NonlinearFactory, UpdaterPair,
+ EnergyNorm<ScalarMatrix, ScalarVector>>
+ adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
+ computeExternalForces, stateEnergyNorm, mustRefine);
while (!adaptiveTimeStepper.reachedEnd()) {
adaptiveTimeStepper.advance();
reportTimeStep(adaptiveTimeStepper.getRelativeTime(),
--
GitLab