# The Problem of Convergence

**URL:** <https://community.geodynamics.org/t/the-problem-of-convergence/3013>\
**Category:** PyLith\
**Created:** [June 18, 2023, 4:00pm UTC](https://community.geodynamics.org/t/the-problem-of-convergence/3013 "2023-06-18T16:00:19Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![HangWu](https://avatars.discourse-cdn.com/v4/letter/h/f17d59/32.png) [@HangWu](https://community.geodynamics.org/u/HangWu)\
**Post date:** [June 18, 2023, 4:00pm UTC](https://community.geodynamics.org/t/the-problem-of-convergence/3013/1 "2023-06-18T16:00:19Z")

</div>

Hello all,  
I successfully installed Pylith 2.22. But I had a little problem with the simulation .  
1、 I want to simulate the deformation of the interior of the lithosphere under gravity. I take the 2D model, and I take the quadrilateral mesh.  
2、 I set up two lithospheres, plateau 、transition zone and basin, and set their brittle layer as elastic material and ductile layer as viscoelastic material through rheological strength.  
3、 For the transition zone, the density and viscosity are given as a linear vary along the transverse direction and depth. For plateaus and basins, only vary in depth are considered .  
4、 No faults are defined within the lithosphere.  
5、 I set the sides and bottom of the model to free slip boundaries .  
6、 In order to test the correctness of the simulation, I have set the mesh to be rough. In previous simulations, it might have taken only a few minutes to compute, but now it takes so long to execute that it seems impossible to converge .

The petsc message is as follows ：  
[pylithapp.petsc]

# Preconditioner settings.

pc\_type = asm

# Convergence parameters.

ksp\_rtol = 1.0e-9  
ksp\_atol = 1.0e-9

# 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-9  
snes\_atol = 1.0e-9  
snes\_max\_it = 500  
snes\_monitor = true  
#snes\_view = true  
snes\_converged\_reason = true  
snes\_linesearch\_monitor = true  
snes\_error\_if\_not\_converged = true

The error message is as follows ：

> > /home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/lib/python2.7/site-packages/pylith/problems/Implicit.py:207:step  
> > – implicit(info)  
> > – Solving equations.  
> > 0 SNES Function norm 3.504722374800e-09  
> > 0 KSP Residual norm 9.999723312736e-10  
> > Linear solve converged due to CONVERGED\_ATOL iterations 0  
> > Line search: Initial direction and size is 0  
> > Nonlinear solve did not converge due to DIVERGED\_LINE\_SEARCH iterations 0  
> > [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------  
> > [0]PETSC ERROR:  
> > [0]PETSC ERROR: SNESSolve has not converged  
> > [0]PETSC ERROR: See [http://www.mcs.anl.gov/petsc/documentation/faq.html](http://www.mcs.anl.gov/petsc/documentation/faq.html) for trouble shooting.  
> > [0]PETSC ERROR: Petsc Release Version 3.10.2, Jul, 01, 2019  
> > [0]PETSC ERROR: /home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/bin/mpinemesis on a named ubuntu by wuhang Sun Jun 18 20:33:56 2023  
> > [0]PETSC ERROR: Configure options --prefix=/home/brad/pylith-binary/dist --with-c2html=0 --with-x=0 --with-clanguage=C --with-mpicompilers=1 --with-shared-libraries=1 --with-64-bit-points=1 --with-large-file-io=1 --download-chaco=1 --download-ml=1 --download-f2cblaslapack=1 --with-hwloc=0 --with-ssl=0 --with-x=0 --with-c2html=0 --with-lgrind=0 --with-hdf5=1 --with-hdf5-dir=/home/brad/pylith-binary/dist --with-zlib=1 --LIBS=-lz --with-debugging=0 --with-fc=0 CPPFLAGS=“-I/home/brad/pylith-binary/dist/include -I/home/brad/pylith-binary/dist/include " LDFLAGS=”-L/home/brad/pylith-binary/dist/lib -L/home/brad/pylith-binary/dist/lib64 -L/home/brad/pylith-binary/dist/lib -L/home/brad/pylith-binary/dist/lib64 " CFLAGS=“-g -O2” CXXFLAGS=“-g -O2 -DMPICH\_IGNORE\_CXX\_SEEK” FCFLAGS= PETSC\_DIR=/home/brad/pylith-binary/build/petsc-pylith PETSC\_ARCH=arch-pylith  
> > [0]PETSC ERROR: #1 SNESSolve() line 4408 in /home/brad/pylith-binary/build/petsc-pylith/src/snes/interface/snes.c  
> > [0]PETSC ERROR: #2 void pylith::problems::SolverNonlinear::solve(pylith::topology::Field\*, pylith::topology::Jacobian\*, const pylith::topology::Field&)() line 152 in ../../../pylith-2.2.2/libsrc/pylith/problems/SolverNonlinear.cc  
> > Fatal error. Calling MPI\_Abort() to abort PyLith application.  
> > Traceback (most recent call last):  
> > File “/home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/lib/python2.7/site-packages/pylith/apps/PetscApplication.py”, line 74, in onComputeNodes  
> > self.main(\*args, \*\*kwds)  
> > File “/home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/lib/python2.7/site-packages/pylith/apps/PyLithApp.py”, line 138, in main  
> > self.problem.run(self)  
> > File “/home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/lib/python2.7/site-packages/pylith/problems/TimeDependent.py”, line 203, in run  
> > self.formulation.step(t, dt)  
> > File “/home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/lib/python2.7/site-packages/pylith/problems/Implicit.py”, line 212, in step  
> > self.solver.solve(dispIncr, self.jacobian, residual)  
> > File “/home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/lib/python2.7/site-packages/pylith/problems/problems.py”, line 186, in solve  
> > def solve(self, \*args): return \_problems.SolverNonlinear\_solve(self, \*args)  
> > RuntimeError: Error detected while in PETSc function.  
> > application called MPI\_Abort(MPI\_COMM\_WORLD, -1) - process 0  
> > /home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/bin/nemesis: mpirun: exit 255  
> > /home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/bin/pylith: /home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86\_64/bin/nemesis: exit 1

Does anyone know what caused this and how I can fix it?  
Thank you

[gravity\_infstrain.cfg](https://community.geodynamics.org/uploads/short-url/9CCgrxK5bjySUCDMAlAOJOSVzUx.cfg) (14.7 KB)  
[pylithapp.cfg](https://community.geodynamics.org/uploads/short-url/tG80wBezjosxPU63Uvq34sz6X3g.cfg) (24.0 KB)

---

<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:** [June 19, 2023, 2:46pm UTC](https://community.geodynamics.org/t/the-problem-of-convergence/3013/2 "2023-06-19T14:46:51Z")

</div>

The linear and nonlinear solver tolerances are the same. In general, we recommend using nonlinear solver tolerances that are about 2 orders of magnitude larger than the linear solver tolerances. We also recommend using the ML (algebraic multigrid perconditioner) rather than the ASM preconditioner for elasticity problems without faults.

These settings will likely work better than what you have. If these do not work, then provide an update that includes a diagram of the boundary value problem you are trying to solve with all materials and boundary conditions labeled.

```auto
[pylithapp.petsc]
ksp_atol = 1.0e-10
ksp_rtol = 1.0e-8

snes_atol = 10e-8
snes_rtol = 1.0e-8

pc_type = ml

```

---

<div class="post-metadata">

**Author:** ![HangWu](https://avatars.discourse-cdn.com/v4/letter/h/f17d59/32.png) [@HangWu](https://community.geodynamics.org/u/HangWu)\
**Post date:** [June 20, 2023, 1:41am UTC](https://community.geodynamics.org/t/the-problem-of-convergence/3013/3 "2023-06-20T01:41:15Z")

</div>

Thank you Brad!  
After I changed the solver tolerance, the problem was solved .  
I have two more questions:  
1、 For maxwell materials, we need to give the viscosity, so can we output the viscosity to the info file like the density after the simulation, so that we can use Paraview to check whether the viscosity is correct.  
2、 For a mantle of 50km-100km, I give the viscosity change from 0km-100km, so does pylith read the viscosity start from 0km or from the 50km ？

---

<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:** [June 20, 2023, 2:25am UTC](https://community.geodynamics.org/t/the-problem-of-convergence/3013/4 "2023-06-20T02:25:02Z")

</div>

See `examples/2d/subduction` for an example of how to use a linear Maxwell model in PyLith v2.2.2. The input field is viscosity, but the Maxwell time is stored internally, so it is available as an output field. For example, for a material component `mantle`, you can use `mantle.cell_info_fields = [density, mu, lambda, maxwell_time]`. The Maxwell time = viscosity / shear\_modulus.

Yes, you can use a spatial database that specifies values over a larger coordinate range than that needed in a simulation.
