Be smarter about damping

src/one-body-sample.cc

@@ -134,6 +134,13 @@ template <class FunctionMap> void initPython(FunctionMap &functions) {
       .toC<typename FunctionMap::Base>(functions);
 }
 
+template <class VectorType>
+double diff_two_norm(VectorType const &v1, VectorType const &v2) {
+  VectorType tmp = v1;
+  tmp -= v2;
+  return tmp.two_norm();
+}
+
 int main(int argc, char *argv[]) {
   try {
     Dune::Timer timer;
@@ -343,11 +350,12 @@ int main(int argc, char *argv[]) {
     };
 
     VectorType ud = ud_initial;
+    SingletonVectorType alpha = alpha_initial;
     auto const state_fpi_max =
         parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
     for (size_t run = 1; run <= timesteps; ++run) {
       VectorType u;
-      SingletonVectorType alpha;
+      double lastCorrection;
 
       stateUpdater->nextTimeStep();
       timeSteppingScheme->nextTimeStep();
@@ -393,15 +401,23 @@ int main(int argc, char *argv[]) {
       double const fixedPointTolerance =
           parset.get<double>("solver.tnnmg.fixed_point_tolerance");
       double const damping = parset.get<double>("solver.damping");
+      double const minimalCorrectionReduction =
+          parset.get<double>("solver.minimal_correction_reduction");
+
       for (size_t state_fpi = 1; state_fpi <= state_fpi_max; ++state_fpi) {
         stateUpdater->solve(ud);
-        if (state_fpi == 1)
-          stateUpdater->extractState(alpha);
-        else {
+        {
           SingletonVectorType computed_state;
           stateUpdater->extractState(computed_state);
-          alpha *= damping;
-          alpha.axpy(1.0 - damping, computed_state);
+          double const correction = diff_two_norm(computed_state, alpha);
+          if (state_fpi <= 2 // Let the first two steps pass through unchanged
+              || correction < minimalCorrectionReduction * lastCorrection)
+            alpha = computed_state;
+          else {
+            alpha *= damping;
+            alpha.axpy(1.0 - damping, computed_state);
+          }
+          lastCorrection = correction;
         }
 
         solveDisplacementProblem(problem_iterate, alpha);

src/one-body-sample.parset

@@ -23,7 +23,8 @@ width = 5
 
 [solver]
 tolerance = 1e-10
-damping = 0.0
+damping = 0.5
+minimal_correction_reduction = 0.5
 
 [solver.tnnmg]
 maxiterations = 1000000