[WIP]Nonlinear solver: deal with failure

This introduces a new setting to decide what should happen if the
nonlinear solver fails to converge.

The old behavior was to "just continue", which is an unsafe choice.
Therefor, I am changing the default. This is not backwards compatible,
of course.

include/aspect/parameters.h

@@ -97,6 +97,38 @@ namespace aspect
              false;
     }
 
+    /**
+     * A struct that contains the enums to decide what to do when a nonlinear solver fails.
+     */
+    struct NonlinearSolverFailureStrategy
+    {
+      enum Kind
+      {
+        continue_with_next_timestep,
+        cut_timestep_size,
+        abort_program
+      };
+
+      /**
+       * parse the enum value from a string
+       */
+      static
+      Kind
+      parse(const std::string &input)
+      {
+        if (input == "continue with next timestep")
+          return continue_with_next_timestep;
+        else if (input == "cut timestep size")
+          return cut_timestep_size;
+        else if (input == "abort program")
+          return abort_program;
+        else
+          AssertThrow(false, ExcNotImplemented());
+
+        return Kind();
+      }
+    };
+
     /**
      * @brief The NullspaceRemoval struct
      */
@@ -450,6 +482,7 @@ namespace aspect
      * @{
      */
     typename NonlinearSolver::Kind nonlinear_solver;
+    typename NonlinearSolverFailureStrategy::Kind nonlinear_solver_failure_strategy;
 
     typename AdvectionStabilizationMethod::Kind advection_stabilization_method;
     double                         nonlinear_tolerance;

include/aspect/simulator.h

@@ -168,6 +168,14 @@ namespace aspect
     Tensor<1,dim> tensor_rotation;
   };
 
+  /**
+   * Exception to be thrown when the nonlinear solver needs to many iterations to converge.
+   */
+  DeclExceptionMsg(ExcNonlinearSolverNoConvergence,
+                   "Nonlinear solver failed to converge in the prescribed number of steps. "
+                   "Consider changing `Max nonlinear iterations` or `Nonlinear solver failure "
+                   "strategy`.");
+
   /**
    * This is the main class of ASPECT. It implements the overall simulation
    * algorithm using the numerical methods discussed in the papers and manuals

source/simulator/core.cc

@@ -2018,15 +2018,53 @@ namespace aspect
     simulator_is_past_initialization = true;
     do
       {
-        // Only solve if we are not in pre-refinement, or we do not want to skip
-        // solving in pre-refinement.
+        // During pre-refinement, do not solve if we are asked to skip it:
         if (! (parameters.skip_solvers_on_initial_refinement
                && pre_refinement_step < parameters.initial_adaptive_refinement))
           {
             start_timestep ();
 
             // then do the core work: assemble systems and solve
-            solve_timestep ();
+            try
+              {
+                solve_timestep ();
+              }
+            catch (...)
+              {
+                pcout << "ERROR: the solver in the current timestep failed to converge." << std::endl;
+
+                if (parameters.nonlinear_solver_failure_strategy
+                    == Parameters<dim>::NonlinearSolverFailureStrategy::continue_with_next_timestep)
+                  {
+                    // do nothing and continue
+                  }
+                else if (parameters.nonlinear_solver_failure_strategy
+                         == Parameters<dim>::NonlinearSolverFailureStrategy::cut_timestep_size)
+                  {
+                    if (timestep_number == 0)
+                      {
+                        pcout << "Note: we can not cut the timestep in step 0, so we are aborting."
+                              << std::endl;
+                        throw;
+                      }
+
+                    // reduce timestep size and update time to "go back":
+                    const double new_time_step = 0.5*time_step;
+                    pcout << "\tReducing timestep to " << new_time_step << '\n' << std::endl;
+
+                    time -= time_step - new_time_step;
+                    time_step = new_time_step;
+                    continue; // skip the rest of the main loop that would advance in time
+                  }
+                else if (parameters.nonlinear_solver_failure_strategy
+                         == Parameters<dim>::NonlinearSolverFailureStrategy::abort_program)
+                  {
+                    // rethrow the current exception
+                    throw;
+                  }
+                else
+                  AssertThrow(false, ExcNotImplemented());
+              }
           }
 
         // See if we have to start over with a new adaptive refinement cycle

source/simulator/parameters.cc

@@ -265,6 +265,14 @@ namespace aspect
                        "The `Newton Stokes' scheme is deprecated and only allowed for reasons of "
                        "backwards compatibility. It is the same as `iterated Advection and Newton Stokes'.");
 
+    prm.declare_entry ("Nonlinear solver failure strategy", "abort program",
+                       Patterns::Selection("continue with next timestep|cut timestep size|abort program"),
+                       "Select the strategy on what to do if the nonlinear solver scheme fails to "
+                       "converge. The options are:\n"
+                       "`continue with next timestep`: ignore error and continue to the next timestep\n"
+                       "`cut timestep size`: reduce the current timestep size and redo the timestep\n"
+                       "`abort program`: abort the program with an error message.");
+
     prm.declare_entry ("Nonlinear solver tolerance", "1e-5",
                        Patterns::Double(0., 1.),
                        "A relative tolerance up to which the nonlinear solver will iterate. "
@@ -1417,6 +1425,8 @@ namespace aspect
       else
         AssertThrow (false, ExcNotImplemented());
     }
+    nonlinear_solver_failure_strategy = NonlinearSolverFailureStrategy::parse(
+                                          prm.get("Nonlinear solver failure strategy"));
 
     prm.enter_subsection ("Solver parameters");
     {

source/simulator/solver_schemes.cc

@@ -814,6 +814,7 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
     signals.post_nonlinear_solver(nonlinear_solver_control);
   }
 
@@ -872,6 +873,8 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
+
     // Reset the linear tolerance to what it was at the beginning of the time step.
     parameters.linear_stokes_solver_tolerance = begin_linear_tolerance;
 
@@ -942,6 +945,8 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
+
     // Reset the linear tolerance to what it was at the beginning of the time step.
     parameters.linear_stokes_solver_tolerance = begin_linear_tolerance;
 
@@ -1049,6 +1054,8 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
+
     // Reset the linear tolerance to what it was at the beginning of the time step.
     parameters.linear_stokes_solver_tolerance = begin_linear_tolerance;
 
@@ -1169,6 +1176,7 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
     signals.post_nonlinear_solver(nonlinear_solver_control);
   }
 
@@ -1216,6 +1224,7 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
     signals.post_nonlinear_solver(nonlinear_solver_control);
   }
 
@@ -1315,6 +1324,8 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
+
     // Reset the Newton stabilization at the end of the timestep.
     newton_handler->parameters.preconditioner_stabilization = starting_preconditioner_stabilization;
     newton_handler->parameters.velocity_block_stabilization = starting_velocity_block_stabilization;
@@ -1419,6 +1430,8 @@ namespace aspect
       }
     while (nonlinear_solver_control.check(nonlinear_iteration, relative_residual) == SolverControl::iterate);
 
+    AssertThrow(nonlinear_solver_control.last_check() != SolverControl::failure, ExcNonlinearSolverNoConvergence());
+
     // Reset the Newton stabilization at the end of the timestep.
     newton_handler->parameters.preconditioner_stabilization = starting_preconditioner_stabilization;
     newton_handler->parameters.velocity_block_stabilization = starting_velocity_block_stabilization;

tests/nonlinear_failure_strategy_cut.cc

@@ -0,0 +1 @@
+#include "../benchmarks/nonlinear_channel_flow/simple_nonlinear.cc"

tests/nonlinear_failure_strategy_cut.prm

@@ -0,0 +1,11 @@
+# Check that nonlinear solver strategy "cut timestep size" works
+#
+# Like iterated_advection_and_stokes_residual.prm
+
+include $ASPECT_SOURCE_DIR/tests/iterated_advection_and_stokes_residual.prm
+
+set Nonlinear solver failure strategy = cut timestep size
+set Max nonlinear iterations = 7
+
+set Dimension = 2
+set End time = 20000