Nonlinear convergence problem with static friction fault

PyLith
1

Hello, CIG team and users!

I’m running a 3d-mesh model with a static friction fault by using Pylith 2.2.2. I apply velocity boundary conditions in east-west and north-south directions and fix the bottom boundary. I want to model the situations with full locking fault or full unlocking fault. So I set the friction coefficient to 1000000000 (approximate infinity) or 0 for the above situations. I also set the solver settings according to Chapter 6.4.5.2 in Pylith manual. Unfortunately, it was stuck after linear friction solve converged. The model can run when we remove the static friction fault. So it must be the static friction problem. I also test the influence of the normal traction on the fault. The results show that the problem persists regardless of the normal traction. The pictures show the fault nodes in 3d-mesh model and the stuck situation, respectively. The BHsim.cfg and pylithapp.cfg show as follows. Please give me some advice to solve the problem. Thanks!

###################BHsim.cfg################
[pylithapp]
[pylithapp.problem.formulation]
time_step = pylith.problems.TimeStepAdapt ;
[pylithapp.problem.formulation.time_step]
total_time = 100.0year
max_dt = 1year
adapt_skip = 10 ; Default value
stability_factor = 2
[pylithapp.timedependent]
###################################
[pylithapp.timedependent]
bc=[modelBC,bott]
[pylithapp.timedependent.implicit]

Set the output to an array of 2 output managers.

We will output the solution over the domain and the ground surface.

output = [domain,subdomain]

Set subdomain component to OutputSolnSubset (subset of domain).

output.subdomain = pylith.meshio.OutputSolnSubset
solver = pylith.problems.SolverNonlinear
#-----------------------------------------------------------------------
[pylithapp.timedependent.bc.modelBC]
bc_dof = [0,1]
label = modelBC
db_initial.label = Dirichlet BC on modelBC
db_rate = spatialdata.spatialdb.SimpleDB
db_rate.label = Dirichlet rate BC on modelBC
db_rate.iohandler.filename = spatialdb/modelBCdV.spatialdb

[pylithapp.timedependent.bc.bott]
bc_dof = [2]
label = bott
db_initial.label = Dirichlet BC on bott

----------------------------------------------------------------------

[pylithapp.timedependent]
interfaces=[lmsf]

Set the type of fault interface condition.

[pylithapp.timedependent.interfaces]
lmsf = pylith.faults.FaultCohesiveDyn

[pylithapp.timedependent.interfaces.lmsf]
zero_tolerance = 1.0e-7
zero_tolerance_normal=1.0e-9
label =lmsf
id=21
edge=lmsf_edge
id=22
friction = pylith.friction.StaticFriction
friction.label = static friction
friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = static friction
friction.db_properties.values = [friction-coefficient, cohesion]
friction.db_properties.data = [1000000000, 0.0*Pa]
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2

output

----------------------------------------------------------------------

[pylithapp.problem.formulation.output.domain]
writer.filename = output/CDstress.h5

Ground surface

[pylithapp.problem.formulation.output.subdomain]
writer.filename = output/Sgroundsurf.h5

----------------------------------------------------------------------

PETSc

----------------------------------------------------------------------

[pylithapp.petsc]
friction_pc_type = asm
friction_sub_pc_factor_shift_type = nonzero
friction_ksp_max_it = 25
friction_ksp_gmres_restart = 30
friction_ksp_converged_reason = true
###################BHsim.cfg################
###################Pylithapp.cfg##############
[pylithapp]

----------------------------------------------------------------------

launcher

----------------------------------------------------------------------

[pylithapp.launcher]
command = mpirun -np 2

----------------------------------------------------------------------

journal

----------------------------------------------------------------------

[pylithapp.journal.info]
timedependent = 1
implicit = 1
petsc = 1
solverlinear = 1
meshiocubit = 1
implicitelasticity = 1
faultcohesivekin = 1
fiatsimplex = 1
pylithapp = 1
materials = 1

----------------------------------------------------------------------

mesh_generator

----------------------------------------------------------------------

[pylithapp.mesh_generator]
reader = pylith.meshio.MeshIOCubit
[pylithapp.mesh_generator.reader]
filename = mesh/mesh_tet.exo
coordsys.space_dim = 3

----------------------------------------------------------------------

problem

----------------------------------------------------------------------

[pylithapp.timedependent]
dimension = 3

----------------------------------------------------------------------

materials

----------------------------------------------------------------------

[pylithapp.timedependent]
materials=[bh1,bh2,bh3,qq1,qq2,qq3,qt1,qt2,qt3,sc1,sc2,sc3]

[pylithapp.timedependent.materials]

bh1 = pylith.materials.ElasticIsotropic3D
bh2 = pylith.materials.MaxwellIsotropic3D
bh3 = pylith.materials.MaxwellIsotropic3D

qq1 = pylith.materials.ElasticIsotropic3D
qq2 = pylith.materials.MaxwellIsotropic3D
qq3 = pylith.materials.MaxwellIsotropic3D

qt1 = pylith.materials.ElasticIsotropic3D
qt2 = pylith.materials.MaxwellIsotropic3D
qt3 = pylith.materials.MaxwellIsotropic3D

sc1 = pylith.materials.ElasticIsotropic3D
sc2 = pylith.materials.MaxwellIsotropic3D
sc3 = pylith.materials.MaxwellIsotropic3D

Continental crust -----------------------------

[pylithapp.timedependent.materials.bh1]
label = bh1
id = 4
db_properties.label = bh1 properties
db_properties.iohandler.filename = spatialdb/bhup.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

[pylithapp.timedependent.materials.bh2]
label = bh2
id = 7
db_properties.label = bh2 properties
db_properties.iohandler.filename = spatialdb/bhlow.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

[pylithapp.timedependent.materials.bh3]
label = bh3
id = 8
db_properties.label = bh3 properties
db_properties.iohandler.filename = spatialdb/bhmantle.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

[pylithapp.timedependent.materials.sc1]
label = sc1
id = 3
db_properties.label = sc1 properties
db_properties.iohandler.filename = spatialdb/scup.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

[pylithapp.timedependent.materials.sc2]
label = sc2
id = 9
db_properties.label = sc2 properties
db_properties.iohandler.filename = spatialdb/sclow.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

[pylithapp.timedependent.materials.sc3]
label = sc3
id = 10
db_properties.label = sc3 properties
db_properties.iohandler.filename = spatialdb/scmantle.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

[pylithapp.timedependent.materials.qq1]
label = qq1
id = 2
db_properties.label = qq1 properties
db_properties.iohandler.filename = spatialdb/qqup.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

[pylithapp.timedependent.materials.qq2]
label = qq2
id = 11
db_properties.label = qq2 properties
db_properties.iohandler.filename = spatialdb/qqlow.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

[pylithapp.timedependent.materials.qq3]
label = qq3
id = 12
db_properties.label = qq3 properties
db_properties.iohandler.filename = spatialdb/qqmantle.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

[pylithapp.timedependent.materials.qt1]
label = qt1
id = 1
db_properties.label = qt1 properties
db_properties.iohandler.filename = spatialdb/qtup.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3

[pylithapp.timedependent.materials.qt2]
label = qt2
id = 5
db_properties.label = qt2 properties
db_properties.iohandler.filename = spatialdb/qtlow.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

[pylithapp.timedependent.materials.qt3]
label = qt3
id = 6
db_properties.label = qt3 properties
db_properties.iohandler.filename = spatialdb/qtmantle.spatialdb
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_info_fields = [density,mu,lambda,maxwell_time]

----------------------------------------------------------------------

output

----------------------------------------------------------------------

[pylithapp.timedependent.formulation]
output = [domain,subdomain]
output.subdomain = pylith.meshio.OutputSolnSubset

Domain

[pylithapp.problem.formulation.output.domain]
output_freq = time_step
time_step = 10*year
writer = pylith.meshio.DataWriterHDF5
vertex_data_fields = [displacement,velocity]

Ground surface

[pylithapp.problem.formulation.output.subdomain]
label = groundface ; Name of CUBIT nodeset for ground surface.
writer = pylith.meshio.DataWriterHDF5
vertex_data_fields = [displacement,velocity]

----------------------------------------------------------------------

PETSc

----------------------------------------------------------------------

Set the solver options.

[pylithapp.petsc]

Preconditioner settings.

pc_type = asm
sub_pc_factor_shift_type = nonzero

Convergence parameters.

ksp_rtol = 1.0e-20
ksp_atol = 1.0e-10
ksp_max_it = 1000
ksp_gmres_restart = 50

Linear solver monitoring options.

ksp_monitor = true
#ksp_view = true
ksp_converged_reason = true
ksp_error_if_not_converged = true

Nonlinear solver monitoring options.

snes_rtol = 1.0e-20
snes_atol = 1.0e-2
snes_max_it = 100
snes_monitor = true
snes_linesearch_monitor = true
#snes_view = true
snes_converged_reason = true
snes_error_if_not_converged = true

PETSc summary – useful for performance information.

#log_view = true

Uncomment to launch gdb when starting PyLith.

start_in_debugger = true

###################Pylithapp.cfg##############

problem
problem886×1048 218 KB

Best,
Egg_chen

2

Is the figure of the mesh you show your entire domain? That is, do you have a U shaped domain or are the cells in the center filtered out? Can you provide a diagram/cartoon of your problem with all boundary conditions and faults labeled?

A SNES absolute tolerance of 1.0e-2 is much too large. It should be just 2-4 orders of magnitude larger than the KSP absolute tolerance. The zero tolerance should be between the KSP and SNES absolute tolerances.

If the fault normal traction is zero, then no matter what the static coefficient of friction is, then friction will be zero. You need to make sure the fault normal traction is compressive (<0) and large enough to generate sufficient friction to keep the fault locked.

3

Dear Brad,
Thanks for your reply. I provide a picture of my model. And I apply Dirichlet velocity boundary conditions on all vertical boundaries. I want to simulate the interseismic velocity with a locking fault model and the long-term fault slip rate with an unlocking fault model.

problem1
problem11534×893 204 KB

Thanks for your advice again. I impose a normal traction using the following input format and it works when the fault is locking. But it is stuck again when I set the friction coefficient to 0 or set the normal traction to 0 Pa in order to model an unlocking fault (i.e. the shear traction is almost 0). But when I change the friction coefficient and normal traction to 0.05 and -10 Mpa, the model works. Does it represent the unlocking fault? And I don’t know what is the lower bound of the shear traction that can make the fault slide free.

##############Initial fault tractions###############
traction_perturbation = pylith.faults.TractPerturbation
traction_perturbation.db_initial = spatialdata.spatialdb.UniformDB
traction_perturbation.db_initial.label = Initial fault tractions
traction_perturbation.db_initial.values = [traction-shear-leftlateral, traction-shear-updip, traction-normal]
traction_perturbation.db_initial.data = [0.0 * MPa,0.0 * MPa, -10000000*MPa]
up_dir = [0,1,0]
##############Initial fault tractions###############

No matter what the total simulation time is, the simulated velocity results on the ground surface are almost the same with the observed velocities under the situation with the locking fault or the unlocking fault (the shear traction is 0.5 MPa). I think the simulated veloctiy must be different between the above two situations. I can not figure out the problem. Please give me some advice. Thanks for your time.

Egg_chen

4

What do you mean by “it is stuck again” when the friction is 0? Do you mean the slip is zero, the simulation hangs, or something else? If the simulation does run to completion and the linear and nonlinear solvers converge, are the tractions on the shear tractions on the fault zero?

If your model is dominated by the deformation associated with the boundary conditions, then whether the fault is locked or free slipping will not have a significant impact. The slip on the fault might only be evident in the observations close to the fault. The lmsf fault is quite close to the edge of the domain, so the deformation field in that region are probably very sensitive to your boundary condition. Normally, we place the fault in the middle of the domain and several fault lengths away from boundaries to better match conditions on Earth.

5

Dear Brad,

I’ m sorry I didn’t make myself clear. “it is stuck” should mean that the simulation hangs. it seems to hang but not, because when I run the simulation with multiple core, no error is displayed on the terminal. but when I set the friction coefficient to 0 and use one core to run it, the error “RuntimeError: WARNING! Fault opening with nonzero traction., v_fault: 1716, opening: 1.52101e-10, normal traction: -3.33333e-05” occurs.

Even if I input the normal traction using the following format, the error also happens. I try to change the value of “zero_tolerance_normal” and specify the buried edges of the fault according to the advice from other topics, it doesn’t work either. How to input the normal traction on the lmsf fault correctly? Thank you.

In addition, the simulation can run to completion when the friction coefficient is not 0 and the normal traction is large enough. I check the slip and slip rate on the lmsf fault, both of them are 0.

##############Initial fault tractions###############
traction_perturbation = pylith.faults.TractPerturbation
traction_perturbation.db_initial = spatialdata.spatialdb.UniformDB
traction_perturbation.db_initial.label = Initial fault tractions
traction_perturbation.db_initial.values = [traction-shear-leftlateral, traction-shear-updip, traction-normal]
traction_perturbation.db_initial.data = [0.0 * MPa,0.0 * MPa, -10000*MPa]
up_dir = [0,1,0]
##############Initial fault tractions###############

Chen

6

I also try to set the friction to be small so that the fault friction can be overcome and the fault can slip at a steady rate like the step12 in section 7.9.7.4. The simulation can run at the first serval time step and then the nonlinear residual norm can’t converaged and become larger and larger until it fails.

7

This is most likely the underlying problem:

RuntimeError: WARNING! Fault opening with nonzero traction., v_fault: 1716, opening: 1.52101e-10, normal traction: -3.33333e-05

With nonplanar fault geometry, slip will tend to cause fault opening due to the geometric incompatibility. You need to increase the tolerance for detecting fault opening. Try setting zero_tolerance_normal to a value that is 3-4 orders of magnitude larger (for example 1.0e-6) than the detected opening (1.0e-10). Keep the friction coefficient set to zero, which is what you want.

Did you adjust the SNES absolute tolerance value like I suggested?

9

Yes, I did. I set the solver setting as follow:

zero_tolerance = 1.0e-9
zero_tolerance_normal = 1.0e-6
ksp_rtol = 1.0e-20
ksp_atol = 1.0e-12
snes_rtol = 1.0e-20
snes_atol = 1.0e-8

But, I got the same error. Should I use simpler model, such as reducing the curvature of the fault surface? Thank you.

10

zero_tolerance and zero_tolerance_normal are not PETSc settings. They are settings specific to FaultCohesiveDyn. See Section 6.4.5.2 Dynamic Rupture Parameters in the PyLith v2.2.2 manual.

11

Sorry, I didn’t make myself clear. I did set the zero_tolerance and zero_tolerance_normal under the

[pylithapp.timedependent.interfaces]
lmsf= pylith.faults.FaultCohesiveDyn
[pylithapp.timedependent.interfaces.lmsf]
12

I set the fault and petsc setting as follows:

[pylithapp.timedependent.interfaces.lmsf]
#open_free_surface=False
zero_tolerance = 1.0e-9
zero_tolerance_normal = 1.0e-6
label = lmsf
id=21
edge=lmsf_edge
id=22
friction = pylith.friction.StaticFriction
friction.label = Static friction
friction.db_properties = spatialdata.spatialdb.UniformDB
friction.db_properties.label = Static friction
friction.db_properties.values = [friction-coefficient, cohesion]
friction.db_properties.data = [0.00, 0.0*Pa]
# Initial fault tractions
traction_perturbation = pylith.faults.TractPerturbation
traction_perturbation.db_initial = spatialdata.spatialdb.UniformDB
traction_perturbation.db_initial.label = Initial fault tractions
traction_perturbation.db_initial.values = [traction-shear-leftlateral, traction-shear-updip, traction-normal]
traction_perturbation.db_initial.data = [0.0*MPa,0.0*MPa, -1000*MPa]
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2
up_dir = [0, 1, 0]

# ----------------------------------------------------------------------
# PETSc
# ----------------------------------------------------------------------
[pylithapp.petsc]
# Preconditioner settings.
#malloc_dump =
pc_type = asm
sub_pc_factor_shift_type = nonzero
# Convergence parameters.
ksp_rtol = 1.0e-20
ksp_atol = 1.0e-12
ksp_max_it = 1000
ksp_gmres_restart = 100
# Linear solver monitoring options.
ksp_monitor = true
#ksp_view = true
ksp_converged_reason = true
ksp_error_if_not_converged = true
# Nonlinear solver monitoring options.
snes_rtol = 1.0e-20
snes_atol = 1.0e-8
snes_max_it = 500
snes_monitor = true
snes_view = true
snes_linesearch_monitor = true
snes_converged_reason = true
snes_error_if_not_converged = true
# End of file

Is it reasonalbe? Please give me some advice. Thank you.

13

The ASM preconditioner will not give good performance. See Table 4.4 in the PyLith v2.2.2 manual for the recommended preconditioner for quasistatic fault problems. We provide them in src/pylith-2.2.2/share/settings/solver_fault_fieldsplit.cfg, so you can simply add this file to the command line when running PyLith.

Can you clarify what you mean by “I got the same error”? If you increased the zero_tolerancce_normal value, the error can’t be the same. Did the fault opening value increase? If so, then keep increasing the zero_tolerance_normal.

14

When I set the zero_tolerancce_normal value to 1.0e-6, the error about Fault opening with nonzero traction still occurred and the fault opening value did increase.

RuntimeError: WARNING! Fault opening with nonzero traction., v_fault: 1556, opening: 1.61471e-06, normal traction: -0.0333585

Then I replaced the ASM preconditioner with the text in the solver_fault_fieldsplit.cfg when running Pylith, but the error of Fault opening with nonzero traction still occurred. So I increased the zero_tolerancce_normal value until up to 1.0e-2, the above error disappeared. But the nonlinear residual norm can’t converaged and become larger and larger until Nonlinear solve did not converge due to DIVERGED_MAX_IT iterations. The problem with details shows as follows:

picture1
picture11082×614 21.4 KB

May I ask what factors can cause this kind of problem? Thank you.

15

The nonlinear solve may fail to converge when the linear approximation is far from the nonlinear system. In your case, this might be related to the amount of curvature on the fault and using a friction coefficient of 0. If the fault is very nonplanar, then there will be lots of little gaps forming, which are deviations from the linear approximation. Using a friction value of 0 means the entire fault is slipping, which also increases the nonlinearity.

Two ways to reduce the nonlinearity include:

  1. Smooth the fault surface so there is less curvature.
  2. Use smaller increments in loading so that there is less slip at each time step.

Addendum: The good news is that you now have better solver settings that are exposing the real underlying issue.

16

Hi Brad, thanks for your replying. How to smooth the fault surface? What I understand is to reduce the curvature of the fault. And about the second way, do you mean that reducing the velocity of the boundary condition? The boundary condition is refered from the realistic GPS velocity field, so it can’t be changed.

17

Most mesh generation tools can smooth surfaces.

If you take smaller time steps, then the slip at each time step will be smaller.

18

Thanks, Brad. I’ll try.

19

Hi, Brad. I did a simple test about replacing the complex curved fault with a planar fault. It can run to completion when I set the friction to 0. Then, I checked the slip on the fault, and found that there are no slip on the edge nodes of the fault. Is it reasonable? There are only dynamic fault in my model. How can I verify that the slip on the fault is correct? Thank you.

20

With buried fault edges, the slip will be zero at the edges. This is necessary for maintaining the correct topology at the fault edges and a stable numerical solution.

There is no way to verify the slip on the fault is correct for complex models with a single simulation. However, you can run multiple simulations using different discretization sizes (usually varying by a factor of 2) to verify that you achieve stable results for the output you are interested in.

21

Thank you for your reply.