# Nonlinear convergence problem with static friction fault

**URL:** <https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672>\
**Category:** PyLith\
**Created:** [November 13, 2022, 11:23am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672 "2022-11-13T11:23:56Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 13, 2022, 11:23am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/1 "2022-11-13T11:23:56Z")

</div>

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.0_year  
max\_dt = 1_year  
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](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/f/fb0c265ca93996e77d1742c2363867f882b4ea17.jpeg)

Best,  
Egg\_chen

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 14, 2022, 3:13am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/2 "2022-11-14T03:13:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 14, 2022, 7:30am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/3 "2022-11-14T07:30:54Z")

</div>

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](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/0/02b5e1d092c9a34fa1fa68f9415ad53a1fe2c379.jpeg)  
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

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 14, 2022, 9:50pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/4 "2022-11-14T21:50:38Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 15, 2022, 7:11am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/5 "2022-11-15T07:11:01Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 15, 2022, 1:39pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/6 "2022-11-15T13:39:32Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 15, 2022, 3:44pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/7 "2022-11-15T15:44:48Z")

</div>

This is most likely the underlying problem:

```auto
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?

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 16, 2022, 2:26am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/9 "2022-11-16T02:26:07Z")

</div>

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

```auto
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.

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 16, 2022, 2:54am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/10 "2022-11-16T02:54:42Z")

</div>

`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.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 16, 2022, 2:59am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/11 "2022-11-16T02:59:12Z")

</div>

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

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

```

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 16, 2022, 3:08am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/12 "2022-11-16T03:08:28Z")

</div>

I set the fault and petsc setting as follows:

```auto
[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]

```

```auto
# ----------------------------------------------------------------------
# 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.

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 16, 2022, 3:17am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/13 "2022-11-16T03:17:03Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 16, 2022, 3:53am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/14 "2022-11-16T03:53:21Z")

</div>

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.

```auto
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](https://canada1.discourse-cdn.com/flex036/uploads/geodynamics/original/2X/0/0ee9cfe86f3502f88008c15cf51ba3bb7def319b.png)

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

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 16, 2022, 3:16pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/15 "2022-11-16T15:16:15Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 17, 2022, 1:23am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/16 "2022-11-17T01:23:53Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 17, 2022, 2:37am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/17 "2022-11-17T02:37:27Z")

</div>

Most mesh generation tools can smooth surfaces.

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

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 17, 2022, 2:50am UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/18 "2022-11-17T02:50:10Z")

</div>

Thanks, Brad. I’ll try.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 17, 2022, 3:09pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/19 "2022-11-17T15:09:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![baagaard](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/baagaard/32/615_2.png) [@baagaard](https://community.geodynamics.org/u/baagaard)\
**Post date:** [November 17, 2022, 3:25pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/20 "2022-11-17T15:25:48Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![egg\_chen](https://avatars.discourse-cdn.com/v4/letter/e/d07c76/32.png) [@egg\_chen](https://community.geodynamics.org/u/egg_chen)\
**Post date:** [November 18, 2022, 12:48pm UTC](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672/21 "2022-11-18T12:48:03Z")

</div>

Thank you for your reply.

[Next page](https://community.geodynamics.org/t/nonlinear-convergence-problem-with-static-friction-fault/2672.md?page=2)