[Algorit] Use adaptive time stepping

src/sand-wedge-data/parset.cfg

@@ -36,6 +36,7 @@ a               = 0.030 # [1]      ?
 b               = 0.015 # [1]      ?
 
 [timeSteps]
+refinementTolerance = 1e-5
 number = 100000
 scheme = newmark
 

src/sand-wedge.cc

@@ -488,17 +488,119 @@ int main(int argc, char *argv[]) {
         parset.get<double>("boundary.friction.V0"));
 
     auto const finalTime = parset.get<double>("problem.finalTime");
-    auto const relativeTau = 1.0 / parset.get<size_t>("timeSteps.number");
     double relativeTime = 0.0;
+
+    double relativeTau = 1e-6; // FIXME (not really important, though)
+
+    using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+    using MyTimeStepper = TimeSteppingScheme<Vector, Matrix, Function, dims>;
+
+    auto const refinementTolerance =
+        parset.get<double>("timeSteps.refinementTolerance");
+
     size_t timeStep = 1;
-    while (relativeTime < 1.0 - 1e-10) {
 
-      CoupledTimeStepper<NonlinearFactory, StateUpdater<ScalarVector, Vector>,
-                         TimeSteppingScheme<Vector, Matrix, Function, dims>>
+    auto stateUpdaterR1 = stateUpdater->clone();
+    auto velocityUpdaterR1 = velocityUpdater->clone();
+    {
+      CoupledTimeStepper<NonlinearFactory, MyStateUpdater, MyTimeStepper>
       coupledTimeStepper(finalTime, factory, parset, myGlobalFriction,
-                         stateUpdater, velocityUpdater, computeExternalForces);
+                         stateUpdaterR1, velocityUpdaterR1,
+                         computeExternalForces);
       coupledTimeStepper.step(relativeTime, relativeTau);
+    }
+
+    std::shared_ptr<MyStateUpdater> stateUpdaterR2 = nullptr;
+    std::shared_ptr<MyTimeStepper> velocityUpdaterR2 = nullptr;
+
+    while (relativeTime < 1.0 - 1e-10) {
+      bool didCoarsen = false;
+
+      while (true) {
+        stateUpdaterR2 = stateUpdaterR1->clone();
+        velocityUpdaterR2 = velocityUpdaterR1->clone();
+
+        ScalarVector alphaR2;
+        {
+          CoupledTimeStepper<NonlinearFactory, MyStateUpdater, MyTimeStepper>
+          coupledTimeStepper(finalTime, factory, parset, myGlobalFriction,
+                             stateUpdaterR2, velocityUpdaterR2,
+                             computeExternalForces);
+          coupledTimeStepper.step(relativeTime + relativeTau, relativeTau);
+        }
+        stateUpdaterR2->extractAlpha(alphaR2);
+
+        auto stateUpdaterC = stateUpdater->clone();
+        auto velocityUpdaterC = velocityUpdater->clone();
+
+        ScalarVector alphaC;
+        {
+          CoupledTimeStepper<NonlinearFactory, MyStateUpdater, MyTimeStepper>
+          coupledTimeStepper(finalTime, factory, parset, myGlobalFriction,
+                             stateUpdaterC, velocityUpdaterC,
+                             computeExternalForces);
+          coupledTimeStepper.step(relativeTime, 2.0 * relativeTau);
+        }
+        stateUpdaterC->extractAlpha(alphaC);
+
+        auto const coarseningError = stateEnergyNorm.diff(alphaC, alphaR2);
+
+        if (coarseningError < refinementTolerance) {
+          stateUpdaterR2 = nullptr;
+          velocityUpdaterR2 = nullptr;
+
+          stateUpdaterR1 = stateUpdaterC;
+          velocityUpdaterR1 = velocityUpdaterC;
+          relativeTau *= 2.0;
+          didCoarsen = true;
+        } else {
+          break;
+        }
+      }
+
+      if (!didCoarsen) {
+        ScalarVector alphaR1;
+        while (true) {
+          auto stateUpdaterF2 = stateUpdater->clone();
+          auto velocityUpdaterF2 = velocityUpdater->clone();
+
+          std::shared_ptr<MyStateUpdater> stateUpdaterF1;
+          std::shared_ptr<MyTimeStepper> velocityUpdaterF1;
+
+          ScalarVector alphaF2;
+          {
+            CoupledTimeStepper<NonlinearFactory, MyStateUpdater, MyTimeStepper>
+            coupledTimeStepper(finalTime, factory, parset, myGlobalFriction,
+                               stateUpdaterF2, velocityUpdaterF2,
+                               computeExternalForces);
+            coupledTimeStepper.step(relativeTime, relativeTau / 2.0);
+
+            stateUpdaterF1 = stateUpdaterF2->clone();
+            velocityUpdaterF1 = velocityUpdaterF2->clone();
+
+            coupledTimeStepper.step(relativeTime + relativeTau / 2.0,
+                                    relativeTau / 2.0);
+          }
+          stateUpdaterF2->extractAlpha(alphaF2);
+          stateUpdaterR1->extractAlpha(alphaR1);
+          auto const refinementError = stateEnergyNorm.diff(alphaR1, alphaF2);
+
+          if (refinementError < refinementTolerance) {
+            break;
+          } else {
+            stateUpdaterR1 = stateUpdaterF1;
+            velocityUpdaterR1 = velocityUpdaterF1;
+            stateUpdaterR2 = stateUpdaterF2;
+            velocityUpdaterR2 = velocityUpdaterF2;
+            relativeTau /= 2.0;
+          }
+        }
+      }
 
+      stateUpdater = stateUpdaterR1;
+      velocityUpdater = velocityUpdaterR1;
+      stateUpdaterR1 = stateUpdaterR2;
+      velocityUpdaterR1 = velocityUpdaterR2;
       relativeTime += relativeTau;
 
       Vector u, ur, vr;