Skip to content

Commit a7c84dc

Browse files
lhy11009lhy11009
committed
Tolerate linear solver failure in non-linear schemes
1 parent 606854e commit a7c84dc

File tree

8 files changed

+661
-13
lines changed

8 files changed

+661
-13
lines changed

‎include/aspect/parameters.h‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -128,6 +128,32 @@ namespace aspect
128128
}
129129
};
130130

131+
/**
132+
* A struct that contains the enums to decide what to do when a linear solver fails.
133+
*/
134+
struct LinearSolverFailureStrategy
135+
{
136+
enum Kind
137+
{
138+
continue_with_nonlinear_solver,
139+
abort
140+
};
141+
/**
142+
* Parse the enum value from a string.
143+
*/
144+
static
145+
Kind
146+
parse(const std::string &input)
147+
{
148+
if (input == "continue with nonlinear solver")
149+
return continue_with_nonlinear_solver;
150+
else if (input == "abort")
151+
return abort;
152+
else
153+
AssertThrow(false, ExcNotImplemented());
154+
}
155+
};
156+
131157
/**
132158
* @brief The NullspaceRemoval struct
133159
*/
@@ -549,6 +575,7 @@ namespace aspect
549575
*/
550576
typename NonlinearSolver::Kind nonlinear_solver;
551577
typename NonlinearSolverFailureStrategy::Kind nonlinear_solver_failure_strategy;
578+
typename LinearSolverFailureStrategy::Kind linear_solver_failure_strategy;
552579

553580
typename AdvectionStabilizationMethod::Kind advection_stabilization_method;
554581
double nonlinear_tolerance;

‎source/simulator/parameters.cc‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -284,6 +284,13 @@ namespace aspect
284284
"iterations, in other words, if it is set to something other than "
285285
"`single Advection, single Stokes' or `single Advection, no Stokes'.");
286286

287+
prm.declare_entry ("Linear solver failure strategy", "abort",
288+
Patterns::Selection("continue with nonlinear solver|abort"),
289+
"Select the strategy on what to do if the linear solver scheme fails to "
290+
"converge. The options are:\n"
291+
"`continue with nonlinear solver`: ignore error and continue with the nonlinear solver\n"
292+
"`abort`: abort with an error message");
293+
287294
prm.declare_entry ("Pressure normalization", "surface",
288295
Patterns::Selection ("surface|volume|no"),
289296
"If and how to normalize the pressure after the solution step. "
@@ -1611,6 +1618,8 @@ namespace aspect
16111618
}
16121619
nonlinear_solver_failure_strategy = NonlinearSolverFailureStrategy::parse(
16131620
prm.get("Nonlinear solver failure strategy"));
1621+
linear_solver_failure_strategy = LinearSolverFailureStrategy::parse(
1622+
prm.get("Linear solver failure strategy"));
16141623

16151624
prm.enter_subsection ("Solver parameters");
16161625
{

‎source/simulator/solver.cc‎

Lines changed: 21 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -910,13 +910,27 @@ namespace aspect
910910
if (parameters.n_expensive_stokes_solver_steps > 0)
911911
solver_controls.push_back(solver_control_expensive);
912912

913-
// Exit with an exception that describes the underlying cause:
914-
Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
915-
"Simulator::solve_stokes",
916-
solver_controls,
917-
exc,
918-
mpi_communicator,
919-
parameters.output_directory+"solver_history.txt");
913+
// when linear failure is allowed, we need to have the following steps other than throwing failure
914+
switch (parameters.linear_solver_failure_strategy)
915+
{
916+
case Parameters<dim>::LinearSolverFailureStrategy::continue_with_nonlinear_solver:
917+
{
918+
pcout << " linear solver failed, continuing" << std::endl;
919+
break;
920+
}
921+
case Parameters<dim>::LinearSolverFailureStrategy::abort:
922+
{
923+
Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
924+
"Simulator::solve_stokes",
925+
solver_controls,
926+
exc,
927+
mpi_communicator,
928+
parameters.output_directory+"solver_history.txt");
929+
break;
930+
}
931+
default:
932+
AssertThrow(false, ExcNotImplemented());
933+
}
920934
}
921935
}
922936

‎source/simulator/solver/stokes_matrix_free_local_smoothing.cc‎

Lines changed: 21 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1522,12 +1522,27 @@ namespace aspect
15221522
if (this->get_parameters().n_expensive_stokes_solver_steps > 0)
15231523
solver_controls.push_back(solver_control_expensive);
15241524

1525-
Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
1526-
"StokesMatrixFreeHandlerLocalSmoothingImplementation::solve",
1527-
solver_controls,
1528-
exc,
1529-
this->get_mpi_communicator(),
1530-
this->get_parameters().output_directory+"solver_history.txt");
1525+
// when linear failure is allowed, we need to have the following steps other than throwing failure
1526+
switch (this->get_parameters().linear_solver_failure_strategy)
1527+
{
1528+
case Parameters<dim>::LinearSolverFailureStrategy::continue_with_nonlinear_solver:
1529+
{
1530+
this->get_pcout() << " linear solver failed (GMG), continuing" << std::endl;
1531+
break;
1532+
}
1533+
case Parameters<dim>::LinearSolverFailureStrategy::abort:
1534+
{
1535+
Utilities::throw_linear_solver_failure_exception("iterative Stokes solver",
1536+
"StokesMatrixFreeHandlerLocalSmoothingImplementation::solve",
1537+
solver_controls,
1538+
exc,
1539+
this->get_mpi_communicator(),
1540+
this->get_parameters().output_directory+"solver_history.txt");
1541+
break;
1542+
}
1543+
default:
1544+
AssertThrow(false, ExcNotImplemented());
1545+
}
15311546
}
15321547
}
15331548

‎tests/tolerate_linear_solver_failure.prm‎

Lines changed: 227 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,227 @@
1+
#### Modified from the continental extension cookbook to test the implementation of linear solver failure strategy
2+
#### This test looks into the linear solver failure strategy within a block GMG solver
3+
4+
#### Global parameters
5+
set Dimension = 2
6+
set Start time = 0
7+
set End time = 0
8+
set Use years in output instead of seconds = true
9+
set Nonlinear solver scheme = single Advection, iterated defect correction Stokes
10+
set Nonlinear solver tolerance = 1e-4
11+
set Max nonlinear iterations = 100
12+
set CFL number = 0.5
13+
set Maximum time step = 20e3
14+
set Pressure normalization = no
15+
set Linear solver failure strategy = continue with nonlinear solver
16+
17+
#### Parameters describing the model
18+
19+
# Governing equations
20+
subsection Formulation
21+
set Formulation = Boussinesq approximation
22+
end
23+
24+
# Model geometry (200x100 km, 20 km spacing)
25+
subsection Geometry model
26+
set Model name = box
27+
subsection Box
28+
set X repetitions = 4
29+
set Y repetitions = 2
30+
set X extent = 200e3
31+
set Y extent = 100e3
32+
end
33+
end
34+
35+
# Globally refinement only
36+
subsection Mesh refinement
37+
set Initial global refinement = 2
38+
set Time steps between mesh refinement = 0
39+
set Skip solvers on initial refinement = true
40+
end
41+
42+
# Use the same value of linear solver tolerance as the nonlinear
43+
# solver tolerance
44+
subsection Solver parameters
45+
subsection Stokes solver parameters
46+
set Stokes solver type = block GMG
47+
set Number of cheap Stokes solver steps = 10
48+
set Maximum number of expensive Stokes solver steps = 2
49+
set Linear solver tolerance = 1e-4
50+
end
51+
end
52+
53+
54+
subsection Mesh deformation
55+
set Mesh deformation boundary indicators = top: free surface, top: diffusion
56+
subsection Free surface
57+
set Surface velocity projection = normal
58+
end
59+
subsection Diffusion
60+
# Diffusivity term. Increasing this value will result
61+
# in a smoother free surface and lower topography
62+
# amplitudes.
63+
set Hillslope transport coefficient = 1.e-8
64+
end
65+
end
66+
67+
68+
subsection Boundary velocity model
69+
set Prescribed velocity boundary indicators = left x: function, right x:function, bottom y:function
70+
71+
subsection Function
72+
set Variable names = x,y
73+
set Function constants = v=0.0025, w=200.e3, d=100.e3
74+
set Function expression = if (x < w/2 , -v, v) ; v*2*d/w
75+
end
76+
end
77+
78+
79+
# Number and names of compositional fields
80+
# The five compositional fields represent:
81+
# 1. The plastic strain that accumualates over time, with the initial plastic strain removed
82+
# 2. The plastic strain that accumulated over time, including the initial plastic strain values
83+
# 3. The upper crust
84+
# 4. The lower crust
85+
# 5. The mantle lithosphere
86+
subsection Compositional fields
87+
set Number of fields = 5
88+
set Names of fields = noninitial_plastic_strain, plastic_strain, crust_upper, crust_lower, mantle_lithosphere
89+
end
90+
91+
92+
# Initial values of different compositional fields
93+
# The upper crust (20 km thick), lower crust (20 km thick)
94+
# and mantle (60 km thick) are continuous horizontal layers
95+
# of constant thickness. The non initial plastic strain is set
96+
# to 0 and the initial plastic strain is randomized between
97+
# 0.5 and 1.5.
98+
subsection Initial composition model
99+
set Model name = function
100+
subsection Function
101+
set Variable names = x,y
102+
set Function expression = 0; \
103+
if(x>50.e3 && x<150.e3 && y>50.e3, 0.5 + rand_seed(1), 0); \
104+
if(y>=80.e3, 1, 0); \
105+
if(y<80.e3 && y>=60.e3, 1, 0); \
106+
if(y<60.e3, 1, 0);
107+
end
108+
end
109+
110+
# Composition: fixed on bottom (inflow boundary), free on sides and top
111+
subsection Boundary composition model
112+
set Fixed composition boundary indicators = bottom
113+
set List of model names = initial composition
114+
end
115+
116+
# Temperature boundary conditions
117+
# Top and bottom (fixed) temperatures are consistent with the initial temperature field
118+
# Note that while temperatures are specified for the model sides, these values are
119+
# not used as the sides are not specified "Fixed temperature boundaries". Rather,
120+
# these boundaries are insulating (zero net heat flux).
121+
subsection Boundary temperature model
122+
set Fixed temperature boundary indicators = bottom, top
123+
set List of model names = box
124+
subsection Box
125+
set Bottom temperature = 1613
126+
set Top temperature = 273
127+
end
128+
end
129+
130+
# Initial temperature field
131+
subsection Initial temperature model
132+
set Model name = function
133+
subsection Function
134+
set Variable names = x,y
135+
set Function constants = h=100e3, ts1=273, ts2=633, ts3=893, \
136+
A1=1.e-6, A2=0.25e-6, A3=0.0, \
137+
k1=2.5, k2=2.5, k3=2.5, \
138+
qs1=0.055, qs2=0.035, qs3=0.030
139+
set Function expression = if( (h-y)<=20.e3, \
140+
ts1 + (qs1/k1)*(h-y) - (A1*(h-y)*(h-y))/(2.0*k1), \
141+
if( (h-y)>20.e3 && (h-y)<=40.e3, \
142+
ts2 + (qs2/k2)*(h-y-20.e3) - (A2*(h-y-20.e3)*(h-y-20.e3))/(2.0*k2), \
143+
ts3 + (qs3/k3)*(h-y-40.e3) - (A3*(h-y-40.e3)*(h-y-40.e3))/(2.0*k3) ) );
144+
end
145+
end
146+
147+
# Constant internal heat production values (W/m^3) for background material
148+
# and compositional fields.
149+
subsection Heating model
150+
set List of model names = compositional heating
151+
subsection Compositional heating
152+
set Use compositional field for heat production averaging = 1, 0, 0, 1, 1, 1
153+
set Compositional heating values = 0.0, 0.0, 0.0, 1.0e-6, 0.25e-6, 0.
154+
end
155+
end
156+
157+
# Material model
158+
# Rheology: Non-linear viscous flow and Drucker Prager Plasticity
159+
subsection Material model
160+
set Model name = visco plastic
161+
set Material averaging = harmonic average
162+
163+
subsection Visco Plastic
164+
165+
# Reference temperature and viscosity
166+
set Reference temperature = 273
167+
168+
# The minimum strain-rate helps limit large viscosities values that arise
169+
# as the strain-rate approaches zero.
170+
# The reference strain-rate is used on the first non-linear iteration
171+
# of the first time step when the velocity has not been determined yet.
172+
set Minimum strain rate = 1.e-20
173+
set Reference strain rate = 1.e-16
174+
175+
# Limit the viscosity with minimum and maximum values
176+
set Minimum viscosity = 1e18
177+
set Maximum viscosity = 1e26
178+
179+
# Thermal diffusivity is adjusted to match thermal conductivities
180+
# assumed in assigning the initial geotherm
181+
set Define thermal conductivities = true
182+
set Thermal conductivities = 2.5
183+
set Heat capacities = 750.
184+
185+
# Density values of 1 are assigned to "strain" fields, which are not taken into
186+
# account when computing material properties.
187+
set Densities = 3300, 2700, 2900, 3300
188+
set Thermal expansivities = 2e-5
189+
190+
# Harmonic viscosity averaging
191+
set Viscosity averaging scheme = harmonic
192+
193+
# Choose to have the viscosity (pre-yield) follow a dislocation
194+
# diffusion or composite flow law. Here, dislocation is selected
195+
# so no need to specify diffusion creep parameters below, which are
196+
# only used if "diffusion" or "composite" option is selected.
197+
set Viscous flow law = dislocation
198+
199+
# Dislocation creep parameters
200+
set Prefactors for dislocation creep = 6.52e-16, 8.57e-28, 7.13e-18, 6.52e-16
201+
set Stress exponents for dislocation creep = 3.5, 4.0, 3.0, 3.5
202+
set Activation energies for dislocation creep = 530.e3, 223.e3, 345.e3, 530.e3
203+
set Activation volumes for dislocation creep = 18.e-6, 0.0, 0.0, 18.e-6
204+
205+
# Plasticity parameters
206+
set Angles of internal friction = 30
207+
set Cohesions = 20.e6
208+
209+
# The parameters below weaken the friction and cohesion by a
210+
# a factor of 4 between plastic strain values of 0.5 and 1.5.
211+
set Strain weakening mechanism = plastic weakening with plastic strain only
212+
set Start plasticity strain weakening intervals = 0.5
213+
set End plasticity strain weakening intervals = 1.5
214+
set Cohesion strain weakening factors = 0.25
215+
set Friction strain weakening factors = 0.25
216+
217+
end
218+
end
219+
220+
# Gravity model
221+
subsection Gravity model
222+
set Model name = vertical
223+
224+
subsection Vertical
225+
set Magnitude = 9.81
226+
end
227+
end

0 commit comments

Comments
 (0)