Skip to content

Tolerate linear solver failure in non-linear schemes #6825

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
lhy11009 wants to merge 2 commits into
base: main
Choose a base branch
from
+661 −13
Open

Tolerate linear solver failure in non-linear schemes #6825

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
a7c84dc
Tolerate linear solver failure in non-linear schemes
lhy11009 Mar 6, 2026
cc69f40
follow
lhy11009 Mar 12, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 27 additions & 0 deletions include/aspect/parameters.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -128,6 +128,32 @@ namespace aspect
}
};

/**
* A struct that contains the enums to decide what to do when a linear solver fails.
*/
struct LinearSolverFailureStrategy
{
enum Kind
{
continue_with_nonlinear_solver,
abort
};
/**
Copy link
Member

@gassmoeller gassmoeller Mar 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

put an empty line here before the comment

* Parse the enum value from a string.
*/
static
Kind
parse(const std::string &input)
{
if (input == "continue with nonlinear solver")
return continue_with_nonlinear_solver;
else if (input == "abort")
return abort;
else
AssertThrow(false, ExcNotImplemented());
}
};

/**
* @brief The NullspaceRemoval struct
*/
Expand Down Expand Up @@ -549,6 +575,7 @@ namespace aspect
*/
typename NonlinearSolver::Kind nonlinear_solver;
typename NonlinearSolverFailureStrategy::Kind nonlinear_solver_failure_strategy;
typename LinearSolverFailureStrategy::Kind linear_solver_failure_strategy;

typename AdvectionStabilizationMethod::Kind advection_stabilization_method;
double nonlinear_tolerance;
Expand Down
9 changes: 9 additions & 0 deletions source/simulator/parameters.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -284,6 +284,13 @@ namespace aspect
"iterations, in other words, if it is set to something other than "
"`single Advection, single Stokes' or `single Advection, no Stokes'.");

prm.declare_entry ("Linear solver failure strategy", "abort",
Patterns::Selection("continue with nonlinear solver|abort"),
"Select the strategy on what to do if the linear solver scheme fails to "
"converge. The options are:\n"
"`continue with nonlinear solver`: ignore error and continue with the nonlinear solver\n"
"`abort`: abort with an error message");

prm.declare_entry ("Pressure normalization", "surface",
Patterns::Selection ("surface|volume|no"),
"If and how to normalize the pressure after the solution step. "
Expand Down Expand Up @@ -1611,6 +1618,8 @@ namespace aspect
}
nonlinear_solver_failure_strategy = NonlinearSolverFailureStrategy::parse(
prm.get("Nonlinear solver failure strategy"));
linear_solver_failure_strategy = LinearSolverFailureStrategy::parse(
prm.get("Linear solver failure strategy"));

prm.enter_subsection ("Solver parameters");
{
Expand Down
28 changes: 21 additions & 7 deletions source/simulator/solver.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -910,13 +910,27 @@ 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");
// Determine whether to warn or throw an exception due to linear solver failure
switch (parameters.linear_solver_failure_strategy)
{
case Parameters<dim>::LinearSolverFailureStrategy::continue_with_nonlinear_solver:
{
pcout << " linear solver failed, continuing" << std::endl;
break;
}
case Parameters<dim>::LinearSolverFailureStrategy::abort:
{
Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
"Simulator::solve_stokes",
solver_controls,
exc,
mpi_communicator,
parameters.output_directory+"solver_history.txt");
break;
Copy link
Contributor Author

@lhy11009 lhy11009 Mar 12, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here, I would get a warning if I don't have this break, but I think as long as the previous "throw" is executed, this break is never reached, is that correct?

Copy link
Member

@gassmoeller gassmoeller Mar 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes that is correct. It is good practice to end every case statement with break otherwise the other cases will be checked as well (unless of course that is the intended behavior). But throwing the exception will exit this function immediately, so we never reach the break statement.

}
default:
AssertThrow(false, ExcNotImplemented());
}
}
}

Expand Down
27 changes: 21 additions & 6 deletions source/simulator/solver/stokes_matrix_free_local_smoothing.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1522,12 +1522,27 @@ 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");
// Determine whether to warn or throw an exception due to linear solver failure
switch (this->get_parameters().linear_solver_failure_strategy)
{
case Parameters<dim>::LinearSolverFailureStrategy::continue_with_nonlinear_solver:
{
this->get_pcout() << " linear solver failed (GMG), continuing" << std::endl;
Copy link
Contributor Author

@lhy11009 lhy11009 Mar 12, 2026 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have put this both in the AMG and the BMG solver, and tweaked the outputs slightly to show these differences. I also have separate tests for these different linear solver schemes.

break;
}
case Parameters<dim>::LinearSolverFailureStrategy::abort:
{
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");
break;
}
default:
AssertThrow(false, ExcNotImplemented());
}
}
}

Expand Down
227 changes: 227 additions & 0 deletions tests/tolerate_linear_solver_failure.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,227 @@
#### Modified from the continental extension cookbook to test the implementation of linear solver failure strategy
#### This test looks into the linear solver failure strategy within a block GMG solver

#### Global parameters
set Dimension = 2
set Start time = 0
set End time = 0
set Use years in output instead of seconds = true
Copy link
Member

@gassmoeller gassmoeller Mar 17, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please use the new spelling to avoid the deprecation warnings in the output:

Suggested change
set Use years in output instead of seconds = true
set Use years instead of seconds = true

set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
set Nonlinear solver tolerance = 1e-4
set Max nonlinear iterations = 100
set CFL number = 0.5
set Maximum time step = 20e3
set Pressure normalization = no
set Linear solver failure strategy = continue with nonlinear solver

#### Parameters describing the model

# Governing equations
subsection Formulation
set Formulation = Boussinesq approximation
end

# Model geometry (200x100 km, 20 km spacing)
subsection Geometry model
set Model name = box
subsection Box
set X repetitions = 4
set Y repetitions = 2
set X extent = 200e3
set Y extent = 100e3
end
end

# Globally refinement only
subsection Mesh refinement
set Initial global refinement = 2
set Time steps between mesh refinement = 0
set Skip solvers on initial refinement = true
end

# Use the same value of linear solver tolerance as the nonlinear
# solver tolerance
subsection Solver parameters
subsection Stokes solver parameters
set Stokes solver type = block GMG
set Number of cheap Stokes solver steps = 10
set Maximum number of expensive Stokes solver steps = 2
set Linear solver tolerance = 1e-4
end
end


subsection Mesh deformation
set Mesh deformation boundary indicators = top: free surface, top: diffusion
subsection Free surface
set Surface velocity projection = normal
end
subsection Diffusion
# Diffusivity term. Increasing this value will result
# in a smoother free surface and lower topography
# amplitudes.
set Hillslope transport coefficient = 1.e-8
end
end


subsection Boundary velocity model
set Prescribed velocity boundary indicators = left x: function, right x:function, bottom y:function

subsection Function
set Variable names = x,y
set Function constants = v=0.0025, w=200.e3, d=100.e3
set Function expression = if (x < w/2 , -v, v) ; v*2*d/w
end
end


# Number and names of compositional fields
# The five compositional fields represent:
# 1. The plastic strain that accumualates over time, with the initial plastic strain removed
# 2. The plastic strain that accumulated over time, including the initial plastic strain values
# 3. The upper crust
# 4. The lower crust
# 5. The mantle lithosphere
subsection Compositional fields
set Number of fields = 5
set Names of fields = noninitial_plastic_strain, plastic_strain, crust_upper, crust_lower, mantle_lithosphere
end


# Initial values of different compositional fields
# The upper crust (20 km thick), lower crust (20 km thick)
# and mantle (60 km thick) are continuous horizontal layers
# of constant thickness. The non initial plastic strain is set
# to 0 and the initial plastic strain is randomized between
# 0.5 and 1.5.
subsection Initial composition model
set Model name = function
subsection Function
set Variable names = x,y
set Function expression = 0; \
if(x>50.e3 && x<150.e3 && y>50.e3, 0.5 + rand_seed(1), 0); \
if(y>=80.e3, 1, 0); \
if(y<80.e3 && y>=60.e3, 1, 0); \
if(y<60.e3, 1, 0);
end
end

# Composition: fixed on bottom (inflow boundary), free on sides and top
subsection Boundary composition model
set Fixed composition boundary indicators = bottom
set List of model names = initial composition
end

# Temperature boundary conditions
# Top and bottom (fixed) temperatures are consistent with the initial temperature field
# Note that while temperatures are specified for the model sides, these values are
# not used as the sides are not specified "Fixed temperature boundaries". Rather,
# these boundaries are insulating (zero net heat flux).
subsection Boundary temperature model
set Fixed temperature boundary indicators = bottom, top
set List of model names = box
subsection Box
set Bottom temperature = 1613
set Top temperature = 273
end
end

# Initial temperature field
subsection Initial temperature model
set Model name = function
subsection Function
set Variable names = x,y
set Function constants = h=100e3, ts1=273, ts2=633, ts3=893, \
A1=1.e-6, A2=0.25e-6, A3=0.0, \
k1=2.5, k2=2.5, k3=2.5, \
qs1=0.055, qs2=0.035, qs3=0.030
set Function expression = if( (h-y)<=20.e3, \
ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
if( (h-y)>20.e3 && (h-y)<=40.e3, \
ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3) ) );
end
end

# Constant internal heat production values (W/m^3) for background material
# and compositional fields.
subsection Heating model
set List of model names = compositional heating
subsection Compositional heating
set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1
set Compositional heating values = 0.0, 0.0, 0.0, 1.0e-6, 0.25e-6, 0.
end
end

# Material model
# Rheology: Non-linear viscous flow and Drucker Prager Plasticity
subsection Material model
set Model name = visco plastic
set Material averaging = harmonic average

subsection Visco Plastic

# Reference temperature and viscosity
set Reference temperature = 273

# The minimum strain-rate helps limit large viscosities values that arise
# as the strain-rate approaches zero.
# The reference strain-rate is used on the first non-linear iteration
# of the first time step when the velocity has not been determined yet.
set Minimum strain rate = 1.e-20
set Reference strain rate = 1.e-16

# Limit the viscosity with minimum and maximum values
set Minimum viscosity = 1e18
set Maximum viscosity = 1e26

# Thermal diffusivity is adjusted to match thermal conductivities
# assumed in assigning the initial geotherm
set Define thermal conductivities = true
set Thermal conductivities = 2.5
set Heat capacities = 750.

# Density values of 1 are assigned to "strain" fields, which are not taken into
# account when computing material properties.
set Densities = 3300, 2700, 2900, 3300
set Thermal expansivities = 2e-5

# Harmonic viscosity averaging
set Viscosity averaging scheme = harmonic

# Choose to have the viscosity (pre-yield) follow a dislocation
# diffusion or composite flow law. Here, dislocation is selected
# so no need to specify diffusion creep parameters below, which are
# only used if "diffusion" or "composite" option is selected.
set Viscous flow law = dislocation

# Dislocation creep parameters
set Prefactors for dislocation creep = 6.52e-16, 8.57e-28, 7.13e-18, 6.52e-16
set Stress exponents for dislocation creep = 3.5, 4.0, 3.0, 3.5
set Activation energies for dislocation creep = 530.e3, 223.e3, 345.e3, 530.e3
set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 18.e-6

# Plasticity parameters
set Angles of internal friction = 30
set Cohesions = 20.e6

# The parameters below weaken the friction and cohesion by a
# a factor of 4 between plastic strain values of 0.5 and 1.5.
set Strain weakening mechanism = plastic weakening with plastic strain only
set Start plasticity strain weakening intervals = 0.5
set End plasticity strain weakening intervals = 1.5
set Cohesion strain weakening factors = 0.25
set Friction strain weakening factors = 0.25

end
end

# Gravity model
subsection Gravity model
set Model name = vertical

subsection Vertical
set Magnitude = 9.81
end
end
Loading
Loading