Add a two-step implicit time stepping scheme

src/one-body-sample.cc

@@ -101,6 +101,119 @@ void setup_boundary(GridView const &gridView,
   }
 }
 
+// Implicit Euler: Solve the problem
+//
+// a(Delta u_new, v - Delta u_new) + j_tau(v) - j_tau(Delta u_new)
+// >= l(w - Delta u_new) - a(u_new, v - Delta u_new)
+//
+// Setup: Substract a(u_new, .) from rhs
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void implicitEulerSetup(VectorType const &ell, MatrixType const &A,
+                        VectorType const &u_old, VectorType const *u_old_old,
+                        VectorType &problem_rhs, VectorType &problem_iterate,
+                        MatrixType &problem_A,
+                        Dune::BitSetVector<dim> const &dirichletNodes,
+                        FunctionType const &dirichletFunction, double time) {
+  problem_A = A;
+  problem_rhs = ell;
+  problem_A.mmv(u_old, problem_rhs);
+
+  problem_iterate = u_old;
+  if (u_old_old)
+    problem_iterate -= *u_old_old;
+
+  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+    if (dirichletNodes[i].count() == dim) {
+      double val;
+      dirichletFunction.evaluate(time, val);
+      problem_iterate[i] = 0;                    // Everything prescribed
+      problem_iterate[i][0] = val - u_old[i][0]; // Time-dependent X direction
+    } else if (dirichletNodes[i][1])
+      problem_iterate[i][1] = 0; // Y direction described
+}
+
+// Extraction: Add previous solution
+template <class VectorType>
+void implicitEulerExtract(VectorType const &u_old, VectorType const *u_old_old,
+                          VectorType const &problem_iterate,
+                          VectorType &solution) {
+  solution = u_old;
+  solution += problem_iterate;
+}
+
+template <class VectorType>
+void implicitEulerExtractVelocity(VectorType const &u_old,
+                                  VectorType const *u_old_old,
+                                  VectorType const &problem_iterate,
+                                  VectorType &diff) {
+  diff = problem_iterate;
+}
+
+// two-Stage implicit algorithm: Solve the problem
+//
+// a(Delta u_new, v - Delta u_new) + j_tau(v) - j_tau(Delta u_new)
+// >= l(w - Delta u_new) - a(u_new, v - Delta u_new)
+//
+// Setup: Substract a(u_new, .) from rhs
+template <class VectorType, class MatrixType, class FunctionType, int dim>
+void twoStageImplicitSetup(VectorType const &ell, MatrixType const &A,
+                           VectorType const &u_old, VectorType const *u_old_old,
+                           VectorType &problem_rhs, VectorType &problem_iterate,
+                           MatrixType &problem_A,
+                           Dune::BitSetVector<dim> const &dirichletNodes,
+                           FunctionType const &dirichletFunction, double time) {
+  problem_A = A;
+  problem_A /= 1.5;
+  problem_rhs = ell;
+  problem_A.usmv(-2, u_old, problem_rhs);
+  problem_A.usmv(.5, *u_old_old, problem_rhs);
+
+  // The finite difference makes a good start
+  problem_iterate = u_old;
+  problem_iterate -= *u_old_old;
+
+  for (size_t i = 0; i < dirichletNodes.size(); ++i)
+    if (dirichletNodes[i].count() == dim) {
+      double val;
+      dirichletFunction.evaluate(time, val);
+      problem_iterate[i] = 0;
+      problem_iterate[i].axpy(-2, u_old[i]);
+      problem_iterate[i].axpy(.5, (*u_old_old)[i]);
+      problem_iterate[i][0] =
+          1.5 * val - 2 * u_old[i][0] + .5 * (*u_old_old)[i][0];
+    } else if (dirichletNodes[i][1])
+      // Y direction described
+      problem_iterate[i][1] = -2 * u_old[i][1] + .5 * (*u_old_old)[i][1];
+}
+
+template <class VectorType>
+void twoStageImplicitExtract(VectorType const &u_old,
+                             VectorType const *u_old_old,
+                             VectorType const &problem_iterate,
+                             VectorType &solution) {
+  solution = problem_iterate;
+  solution.axpy(2, u_old);
+  solution.axpy(-.5, *u_old_old);
+  solution *= 2.0 / 3.0;
+
+  // Check if we split correctly
+  {
+    VectorType test = problem_iterate;
+    test.axpy(-1.5, solution);
+    test.axpy(+2, u_old);
+    test.axpy(-.5, *u_old_old);
+    assert(test.two_norm() < 1e-10);
+  }
+}
+
+template <class VectorType>
+void twoStageImplicitExtractVelocity(VectorType const &u_old,
+                                     VectorType const *u_old_old,
+                                     VectorType const &problem_iterate,
+                                     VectorType &diff) {
+  diff = problem_iterate;
+}
+
 int main(int argc, char *argv[]) {
   try {
     typedef SharedPointerMap<std::string, Dune::VirtualFunction<double, double>>
@@ -191,15 +304,18 @@ int main(int argc, char *argv[]) {
     // {{{ Initialise vectors
     VectorType u(finestSize);
     u = 0.0; // Has to be zero!
+    VectorType u_previous;
     VectorType u_diff(finestSize);
     u_diff = 0.0; // Has to be zero!
 
+    // temporary storage for u
+    VectorType solution(finestSize);
+
     SingletonVectorType alpha_old(finestSize);
     alpha_old = parset.get<double>("boundary.friction.state.initial");
     SingletonVectorType alpha(alpha_old);
 
     SingletonVectorType vonMisesStress;
-    VectorType rhs;
     // }}}
 
     typedef MyConvexProblem<MatrixType, VectorType> MyConvexProblemType;
@@ -244,36 +360,61 @@ int main(int argc, char *argv[]) {
     for (size_t run = 1; run <= timesteps; ++run) {
       double const time = tau * run;
       {
+        VectorType ell(finestSize);
         assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
-            leafView, p1Basis, neumannNodes, rhs, neumannFunction, time);
-        stiffnessMatrix.mmv(u, rhs);
-        // Apply Dirichlet condition
-        for (size_t i = 0; i < finestSize; ++i)
-          if (ignoreNodes[i].count() == dim) {
-            dirichletFunction.evaluate(time, u_diff[i][0]);
-            u_diff[i][0] /= timesteps;
-          }
+            leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
+
+        MatrixType problem_A;
+        VectorType problem_rhs(finestSize);
+        VectorType problem_iterate(finestSize);
+
+        auto setupFunc =
+            (run == 1 || !parset.get<bool>("twoStageImplicit"))
+                ? &implicitEulerSetup<VectorType, MatrixType,
+                                      decltype(dirichletFunction), dim>
+                : &twoStageImplicitSetup<VectorType, MatrixType,
+                                         decltype(dirichletFunction), dim>;
 
+        VectorType *u_old_old_ptr = (run == 1) ? nullptr : &u_previous;
+        setupFunc(ell, stiffnessMatrix, u, u_old_old_ptr, problem_rhs,
+                  problem_iterate, problem_A, ignoreNodes, dirichletFunction,
+                  time);
+
+        VectorType solution_saved = u;
         auto const state_fpi_max =
             parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
         for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
           auto myGlobalNonlinearity =
               assemble_nonlinearity<MatrixType, VectorType>(
                   parset.sub("boundary.friction"), *nodalIntegrals, alpha, tau);
-          MyConvexProblemType const myConvexProblem(stiffnessMatrix,
-                                                    *myGlobalNonlinearity, rhs);
 
+          MyConvexProblemType const myConvexProblem(
+              problem_A, *myGlobalNonlinearity, problem_rhs);
           MyBlockProblemType myBlockProblem(parset, myConvexProblem);
           auto multigridStep = mySolver.getSolver();
-          multigridStep->setProblem(u_diff, myBlockProblem);
+          multigridStep->setProblem(problem_iterate, myBlockProblem);
 
-          VectorType const u_diff_saved = u_diff;
           LoopSolver<VectorType> overallSolver(
               multigridStep, parset.get<size_t>("solver.tnnmg.maxiterations"),
               solver_tolerance, &energyNorm, verbosity,
               false); // absolute error
           overallSolver.solve();
 
+          auto extractFunc = (run == 1 || !parset.get<bool>("twoStageImplicit"))
+                                 ? implicitEulerExtract<VectorType>
+                                 : twoStageImplicitExtract<VectorType>;
+
+          // Extract solution from solver
+          extractFunc(u, u_old_old_ptr, problem_iterate, solution);
+
+          auto extractDiffFunc =
+              (run == 1 || !parset.get<bool>("twoStageImplicit"))
+                  ? implicitEulerExtractVelocity<VectorType>
+                  : twoStageImplicitExtractVelocity<VectorType>;
+
+          // Extract difference from solver
+          extractDiffFunc(u, u_old_old_ptr, problem_iterate, u_diff);
+
           // Update the state
           for (size_t i = 0; i < frictionalNodes.size(); ++i) {
             if (frictionalNodes[i][0]) {
@@ -299,9 +440,11 @@ int main(int argc, char *argv[]) {
             std::cerr << '.';
             std::cerr.flush();
           }
-          if (energyNorm.diff(u_diff_saved, u_diff) <
+          if (energyNorm.diff(solution_saved, solution) <
               parset.get<double>("solver.tnnmg.fixed_point_tolerance"))
             break;
+          else
+            solution_saved = solution;
 
           if (state_fpi == state_fpi_max)
             std::cerr << "[ref = " << refinements
@@ -363,7 +506,8 @@ int main(int argc, char *argv[]) {
             std::cout << u_diff[i][0] * timesteps << " ";
         std::cout << std::endl;
       }
-      u += u_diff;
+      u_previous = u;
+      u = solution;
       alpha_old = alpha;
 
       // Compute von Mises stress and write everything to a file

src/one-body-sample.parset

@@ -13,6 +13,8 @@ printVelocitySteppingComparison = false
 
 enable_timer = false
 
+twoStageImplicit = false
+
 [grid]
 refinements = 4
 

src/one-body-sample.py

@@ -13,17 +13,7 @@ class neumannCondition:
 
 class dirichletCondition:
     def __call__(self, x):
-        return .005
-        # return 0
-        fst = 3e-1
-        snd = 5e-1
-        trd = 3e-1
-        if x < 1.0/5:
-            return fst
-        elif x < 3.0/5:
-            return snd
-        else:
-            return trd
+        return x * .005
 
 Functions = {
     'neumannCondition' : neumannCondition(),