follow

lhy11009 · lhy11009 · commit 26f828cfba2e · 2026-03-05T17:57:42.000-08:00
diff --git a/.vscode/launch.json b/.vscode/launch.json
@@ -8,7 +8,7 @@
             // Resolved by CMake Tools:
             "program": "${workspaceFolder}/build/aspect-debug",
 	    // *** replace --test by the name of the .prm file you want to run: ***
-            "args": ["TwoDSubduction_skip_expensive_stokes.prm"],
+            "args": ["TwoDSubduction_skip_expensive_stokes_smaller.prm"],
             "stopAtEntry": true,
 	    // *** append to the path here if you want to run inside some directory like /cookbooks/ ***
             "cwd": "/home/lochy/Softwares/aspect/dtemp",
diff --git a/source/simulator/newton.cc b/source/simulator/newton.cc
@@ -185,7 +185,7 @@ namespace aspect
                              "true, the derivatives are slowly introduced by the following equation: $\\max(0.0, "
                              "(1.0-(residual/switch\\_initial\\_residual)))$, where switch\\_initial\\_residual is the "
                              "residual at the time when the Newton solver is switched on.");
-          
+
           prm.declare_entry ("Maximum linear Stokes solver tolerance", "1e-2",
                              Patterns::Double (0., 1.),
                              "The linear Stokes solver tolerance is dynamically chosen for the Newton solver, based "
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
@@ -283,14 +283,14 @@ namespace aspect
                        "This parameter is only relevant if the `Nonlinear solver scheme' does nonlinear "
                        "iterations, in other words, if it is set to something other than "
                        "`single Advection, single Stokes' or `single Advection, no Stokes'.");
-          
-                       // todo_solver
-          prm.declare_entry ("Continue nonlinear solver loop after linear solver failure", "false",
-                             Patterns::Bool (),
-                             "This method allows to continue nonlinear solver loop after linear solver fails "
-                             "If set to false, the linear solver failure is thrown as normally. When this is "
-                             "set to true, the failure is tolerated and the nonlinear solver continues with "
-                             "the next linear solve.");
+
+    // todo_solver
+    prm.declare_entry ("Continue nonlinear solver loop after linear solver failure", "false",
+                       Patterns::Bool (),
+                       "This method allows to continue nonlinear solver loop after linear solver fails "
+                       "If set to false, the linear solver failure is thrown as normally. When this is "
+                       "set to true, the failure is tolerated and the nonlinear solver continues with "
+                       "the next linear solve.");
 
     prm.declare_entry ("Pressure normalization", "surface",
                        Patterns::Selection ("surface|volume|no"),
diff --git a/source/simulator/solver.cc b/source/simulator/solver.cc
@@ -910,13 +910,23 @@ namespace aspect
                 if (parameters.n_expensive_stokes_solver_steps > 0)
                   solver_controls.push_back(solver_control_expensive);
 
-                // Exit with an exception that describes the underlying cause:
-                Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
-                                                                 "Simulator::solve_stokes",
-                                                                 solver_controls,
-                                                                 exc,
-                                                                 mpi_communicator,
-                                                                 parameters.output_directory+"solver_history.txt");
+                //todo_solver
+                // when linear failure is allowed, we need to have the following steps other than throwing failure
+                if (parameters.continue_nonlinear_solver_loop_after_linear_solver_failure)
+                  {
+                    pcout << "linear solver failed, but this is allowed and nonlinear solver loop continues (inside solve_stokes())" << std::endl;
+                  }
+                else
+                  {
+
+                    // Exit with an exception that describes the underlying cause:
+                    Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
+                                                                     "Simulator::solve_stokes",
+                                                                     solver_controls,
+                                                                     exc,
+                                                                     mpi_communicator,
+                                                                     parameters.output_directory+"solver_history.txt");
+                  }
               }
           }
 
diff --git a/source/simulator/solver/stokes_matrix_free_local_smoothing.cc b/source/simulator/solver/stokes_matrix_free_local_smoothing.cc
@@ -1522,12 +1522,21 @@ namespace aspect
             if (this->get_parameters().n_expensive_stokes_solver_steps > 0)
               solver_controls.push_back(solver_control_expensive);
 
-            Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
-                                                             "StokesMatrixFreeHandlerLocalSmoothingImplementation::solve",
-                                                             solver_controls,
-                                                             exc,
-                                                             this->get_mpi_communicator(),
-                                                             this->get_parameters().output_directory+"solver_history.txt");
+            //todo_solver
+            // when linear failure is allowed, we need to have the following steps other than throwing failure
+            if (this->get_parameters().continue_nonlinear_solver_loop_after_linear_solver_failure)
+              {
+                this->get_pcout() << "linear solver failed, but this is allowed and nonlinear solver loop continues (inside StokesMatrixFreeHandlerLocalSmoothingImplementation:solve())" << std::endl;
+              }
+            else
+              {
+                Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
+                                                                 "StokesMatrixFreeHandlerLocalSmoothingImplementation::solve",
+                                                                 solver_controls,
+                                                                 exc,
+                                                                 this->get_mpi_communicator(),
+                                                                 this->get_parameters().output_directory+"solver_history.txt");
+              }
           }
       }
 
diff --git a/source/simulator/solver_schemes.cc b/source/simulator/solver_schemes.cc
@@ -683,61 +683,69 @@ namespace aspect
         // why that happened:
 
         // todo_solver
-
-        // If the exception we got is not one of the two documented by
-        // throw_linear_solver_failure_exception(), then we have a genuine
-        // problem here, and will need to get outta here right away:
-        if ((dynamic_cast<const ExcMessage *>(&exc)==nullptr) &&
-            (dynamic_cast<const QuietException *>(&exc)==nullptr))
-          throw;
-
-        // Otherwise, the solver presumably failed because the Newton matrix
-        // is not symmetric or not definite. If we do not actually want to
-        // handle this case, then again leave:
-        if (newton_handler->parameters.use_Newton_failsafe == false)
-          throw;
-
-        // Otherwise, start the solve over again and try with a stabilized version:
-        pcout << "failed, trying again with stabilization" << std::endl;
-        newton_handler->parameters.preconditioner_stabilization = Newton::Parameters::Stabilization::SPD;
-        newton_handler->parameters.velocity_block_stabilization = Newton::Parameters::Stabilization::SPD;
-
-        // If the Stokes matrix depends on the solution, or we have active
-        // velocity boundary conditions, we need to re-assemble the system matrix
-        // (and preconditioner) every time. If we have active boundary conditions,
-        // they could a) depend on the solution, or b) be inhomogeneous. In both
-        // cases, just assembling the RHS will be incorrect.  If no active
-        // boundaries exist, we only have no-slip or free slip conditions, so we
-        // don't need to force assembly of the matrix.
-        if (stokes_matrix_depends_on_solution()
-            ||
-            (nonlinear_iteration == 0 && boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size() > 0))
-          rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true;
-        else if (parameters.enable_prescribed_dilation)
-          // The dilation requires the Stokes matrix (which is on the rhs
-          // in the Newton solver) to be updated.
-          rebuild_stokes_matrix = true;
-
-        assemble_stokes_system();
-
-        /**
-         * Eisenstat Walker method for determining the tolerance
-         */
-        if (nonlinear_iteration > 1)
+        // question, how the dcr.stokes_residuals got handled if solve_stokes throws some exception
+        if (parameters.continue_nonlinear_solver_loop_after_linear_solver_failure)
+          {
+            pcout << "linear solver failed, but this is allowed and nonlinear solver loop continues" << std::endl;
+          }
+        else
           {
-            update_residuals();
 
-            // Eisenstat Walker method for determining the linear solver tolerance
-            if (!use_picard)
+            // If the exception we got is not one of the two documented by
+            // throw_linear_solver_failure_exception(), then we have a genuine
+            // problem here, and will need to get outta here right away:
+            if ((dynamic_cast<const ExcMessage *>(&exc)==nullptr) &&
+                (dynamic_cast<const QuietException *>(&exc)==nullptr))
+              throw;
+
+            // Otherwise, the solver presumably failed because the Newton matrix
+            // is not symmetric or not definite. If we do not actually want to
+            // handle this case, then again leave:
+            if (newton_handler->parameters.use_Newton_failsafe == false)
+              throw;
+
+            // Otherwise, start the solve over again and try with a stabilized version:
+            pcout << "failed, trying again with stabilization" << std::endl;
+            newton_handler->parameters.preconditioner_stabilization = Newton::Parameters::Stabilization::SPD;
+            newton_handler->parameters.velocity_block_stabilization = Newton::Parameters::Stabilization::SPD;
+
+            // If the Stokes matrix depends on the solution, or we have active
+            // velocity boundary conditions, we need to re-assemble the system matrix
+            // (and preconditioner) every time. If we have active boundary conditions,
+            // they could a) depend on the solution, or b) be inhomogeneous. In both
+            // cases, just assembling the RHS will be incorrect.  If no active
+            // boundaries exist, we only have no-slip or free slip conditions, so we
+            // don't need to force assembly of the matrix.
+            if (stokes_matrix_depends_on_solution()
+                ||
+                (nonlinear_iteration == 0 && boundary_velocity_manager.get_prescribed_boundary_velocity_indicators().size() > 0))
+              rebuild_stokes_matrix = rebuild_stokes_preconditioner = assemble_newton_stokes_matrix = true;
+            else if (parameters.enable_prescribed_dilation)
+              // The dilation requires the Stokes matrix (which is on the rhs
+              // in the Newton solver) to be updated.
+              rebuild_stokes_matrix = true;
+
+            assemble_stokes_system();
+
+            /**
+             * Eisenstat Walker method for determining the tolerance
+             */
+            if (nonlinear_iteration > 1)
               {
-                update_Eisenstat_Walker_tolerance();
+                update_residuals();
+
+                // Eisenstat Walker method for determining the linear solver tolerance
+                if (!use_picard)
+                  {
+                    update_Eisenstat_Walker_tolerance();
+                  }
               }
-          }
 
-        build_preconditioner();
+            build_preconditioner();
 
-        // Give this another try:
-        dcr.stokes_residuals = solve_stokes(search_direction);
+            // Give this another try:
+            dcr.stokes_residuals = solve_stokes(search_direction);
+          }
       }
 
     // Apply the solution or the update of the solution
diff --git a/tests/tolerate_linear_solver_failure.prm b/tests/tolerate_linear_solver_failure.prm
diff --git a/tests/tolerate_linear_solver_failure/screen-output b/tests/tolerate_linear_solver_failure/screen-output
