[Algorit] Use state energy norm to terminate

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 &current,
       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_;
 };

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);

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

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>;

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

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");

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

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>;

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
 

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(),