The Problem of Convergence

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 :

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 :

– 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: SNESSolve has not converged
[0]PETSC ERROR: See 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/
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/”, 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/”, line 138, in main
File “/home/wuhang/Desktop/pylith/pylith/pylith-2.2.2-linux-x86_64/lib/python2.7/site-packages/pylith/problems/”, 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/”, 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/”, 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 (14.7 KB)
pylithapp.cfg (24.0 KB)

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.

ksp_atol = 1.0e-10
ksp_rtol = 1.0e-8

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

pc_type = ml

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 ?

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.