[Cleanup] Make adaptive timestepping comprehensible

src/time-stepping/adaptivetimestepper.cc

@@ -32,6 +32,7 @@ AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
       parset_(parset),
       globalFriction_(globalFriction),
       current_(current),
+      R1_(),
       externalForces_(externalForces),
       mustRefine_(mustRefine),
       errorNorm_(errorNorm) {}
@@ -42,73 +43,68 @@ bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::coarsen() {
-  bool didCoarsen = false;
-  while (relativeTime_ + relativeTau_ < 1.0) {
-    R2_ = {R1_.updaters.clone(), FixedPointIterationCounter()};
-    step(R2_, relativeTime_ + relativeTau_, relativeTau_);
-
-    UpdatersWithCount C{current_.clone(), FixedPointIterationCounter()};
-    step(C, relativeTime_, 2.0 * relativeTau_);
+IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
+  /*
+    |     C     | We check here if making the step R1 of size tau is a
+    |  R1 | R2  | good idea. To check if we can coarsen, we compare
+    |F1|F2|     | the result of (R1+R2) with C, i.e. two steps of size
+                  tau with one of size 2*tau. To check if we need to
+    refine, we compare the result of (F1+F2) with R1, i.e. two steps
+    of size tau/2 with one of size tau. The method makes multiple
+    coarsening/refining attempts, with coarsening coming first. */
+  if (R1_.updaters == Updaters())
+    R1_ = step(current_, relativeTime_, relativeTau_);
 
-    if (!mustRefine_(C.updaters, R2_.updaters)) {
-      R2_ = {};
-      R1_ = C;
-      relativeTau_ *= 2.0;
-      didCoarsen = true;
-    } else {
+  bool didCoarsen = false;
+  iterationRegister_.reset();
+  UpdatersWithCount R2;
+  UpdatersWithCount C;
+  while (relativeTime_ + relativeTau_ <= 1.0) {
+    R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
+    C = step(current_, relativeTime_, 2 * relativeTau_);
+    if (mustRefine_(R2.updaters, C.updaters))
       break;
-    }
-  }
-  return didCoarsen;
-}
 
-template <class Factory, class Updaters, class ErrorNorm>
-void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::refine() {
-  while (true) {
-    UpdatersWithCount F1{current_.clone(), FixedPointIterationCounter()};
-    step(F1, relativeTime_, relativeTau_ / 2.0);
-    UpdatersWithCount F2{F1.updaters.clone(), FixedPointIterationCounter()};
-    step(F2, relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
+    didCoarsen = true;
+    relativeTau_ *= 2;
+    R1_ = C;
+  }
+  UpdatersWithCount F1;
+  UpdatersWithCount F2;
+  if (!didCoarsen) {
+    while (true) {
+      F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
+      F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
+                relativeTau_ / 2.0);
+      if (!mustRefine_(F2.updaters, R1_.updaters))
+        break;
 
-    if (!mustRefine_(R1_.updaters, F2.updaters)) {
-      break;
-    } else {
-      R1_ = F1;
-      R2_ = F2;
       relativeTau_ /= 2.0;
+      R1_ = F1;
+      R2 = F2;
     }
   }
-}
 
-template <class Factory, class Updaters, class ErrorNorm>
-IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::advance() {
-  if (R2_.updaters == Updaters()) {
-    R1_ = {current_.clone(), FixedPointIterationCounter()};
-    step(R1_, relativeTime_, relativeTau_); // counting again upon restart
-  } else {
-    R1_ = R2_;
-  }
-
-  iterationRegister_.reset();
-  if (!coarsen())
-    refine();
-
-  current_ = R1_.updaters;
   iterationRegister_.registerFinalCount(R1_.count);
   relativeTime_ += relativeTau_;
+  current_ = R1_.updaters;
+  R1_ = R2;
 
   return iterationRegister_;
 }
 
 template <class Factory, class Updaters, class ErrorNorm>
-void AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
-    UpdatersWithCount &updatersAndCount, double rTime, double rTau) {
-  updatersAndCount.count =
+typename AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
+    Updaters const &oldUpdaters, double rTime, double rTau) {
+  UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
+  newUpdatersAndCount.count =
       MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
-                           updatersAndCount.updaters, errorNorm_,
-                           externalForces_).step(rTime, rTau);
-  iterationRegister_.registerCount(updatersAndCount.count);
+                           newUpdatersAndCount.updaters, errorNorm_,
+                           externalForces_)
+          .step(rTime, rTau);
+  iterationRegister_.registerCount(newUpdatersAndCount.count);
+  return newUpdatersAndCount;
 }
 
 #include "adaptivetimestepper_tmpl.cc"

src/time-stepping/adaptivetimestepper.hh

@@ -38,15 +38,14 @@ class AdaptiveTimeStepper {
                       std::function<bool(Updaters &, Updaters &)> mustRefine);
 
   bool reachedEnd();
-  bool coarsen();
-  void refine();
   IterationRegister advance();
 
   double relativeTime_;
   double relativeTau_;
 
 private:
-  void step(UpdatersWithCount &updatersWithCount, double rTime, double rTau);
+  UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
+                         double rTau);
 
   double finalTime_;
   Factory &factory_;
@@ -54,7 +53,6 @@ class AdaptiveTimeStepper {
   std::shared_ptr<Nonlinearity> globalFriction_;
   Updaters &current_;
   UpdatersWithCount R1_;
-  UpdatersWithCount R2_;
   std::function<void(double, Vector &)> externalForces_;
   std::function<bool(Updaters &, Updaters &)> mustRefine_;
   ErrorNorm const &errorNorm_;