The problem of using pylith2.2.2 to simulate 2-dimensional plate subduction event

[Ultrafollow]

Dear Sir, I encountered some problems and doubts while simulating a 2-D plate subduction event.

Simulated event description

I first created a dive model using cubit, as shown in the image. My idea is to use the speed difference between the east-west boundary as the driving force behind the plate subduction. Set slab_top and slab_bot to FaultCohesiveDyn models, using StaticFriction. I named the parameter configuration script for the event simulation subduction2d.cfg.

I changed the driving force because of the mistake. I set slab_bot as the FaultCohesiveKin model, ConstRateSlipFn. A sliding rate of 4cm/year is specified. slab_top holds FaultCohesiveDyn unchanged.I named the parameter configuration script for the event simulation subduction2d_slip.cfg.

Error feedback

Error feedback from running subduction2d.cfg script:

Nonlinear solve did not converge due to DIVERGED_FUNCTION_COUNT iterations 3351
Fatal error. Calling MPI_Abort() to abort PyLith application.
[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 for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.10.2, Jul, 01, 2019 
[0]PETSC ERROR: /home/fallow/pylith-2.2.2/bin/mpinemesis on a  named fallow-virtual-machine by fallow Sun Feb 23 15:24:00 2025
[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
Traceback (most recent call last):
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/apps/PetscApplication.py", line 74, in onComputeNodes
    self.main(*args, **kwds)
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/apps/PyLithApp.py", line 138, in main
    self.problem.run(self)
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/problems/TimeDependent.py", line 203, in run
    self.formulation.step(t, dt)
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/problems/Implicit.py", line 212, in step
    self.solver.solve(dispIncr, self.jacobian, residual)
  File "/home/fallow/pylith-2.2.2/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/fallow/pylith-2.2.2/bin/nemesis: mpirun: exit 255
/home/fallow/pylith-2.2.2/bin/pylith: /home/fallow/pylith-2.2.2/bin/nemesis: exit 1

And error feedback from running subduction2d_slip.cfg script:

Nonlinear solve did not converge due to DIVERGED_FUNCTION_COUNT iterations 3367
[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 for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.10.2, Jul, 01, 2019 
[0]PETSC ERROR: /home/fallow/pylith-2.2.2/bin/mpinemesis on a  named fallow-virtual-machine by fallow Sun Feb 23 16:11:42 2025
[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/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/apps/PetscApplication.py", line 74, in onComputeNodes
    self.main(*args, **kwds)
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/apps/PyLithApp.py", line 138, in main
    self.problem.run(self)
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/problems/TimeDependent.py", line 203, in run
    self.formulation.step(t, dt)
  File "/home/fallow/pylith-2.2.2/lib/python2.7/site-packages/pylith/problems/Implicit.py", line 212, in step
    self.solver.solve(dispIncr, self.jacobian, residual)
  File "/home/fallow/pylith-2.2.2/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/fallow/pylith-2.2.2/bin/nemesis: mpirun: exit 255
/home/fallow/pylith-2.2.2/bin/pylith: /home/fallow/pylith-2.2.2/bin/nemesis: exit 1

I don’t know if there is something wrong with my parameter setting, can you give me some advice. I will attach my parameter configuration file and model file at the end of the article.
subduction.zip (71.7 KB)

[baagaard]

Your diagram and description do not mention what you are using for the normal traction on the fault. Do you have normal tractions consistent with gravitational body forces. You can impose an initial normal traction consistent with gravitational body forces without having to run a separate simulation. For example, you may want to start with overburden pressure equal to the lithostatic load assuming hydrostatic pore pressure and a uniform density (normal traction = - (density - density_water) * g * depth).

Here are the steps I suggest you take to help debug the error. Each step introduces slightly more complexity, working towards the model that you want.

1. Start with both faults having a prescribed slip of 0 (FaultCohesiveKin). This will help identify any issues with the boundary conditions.
2. Same simulation as #1 but change the slab_bot fault to having a constant slip rate.
3. Same as #2 but change the slab_top fault to have a static coefficient of friction of 1.0e+6. This should keep the fault locked. If there is no slip, then the nonlinear solve should converge in a single iteration. If there is slip, examine the normal and shear traction components.

[Ultrafollow]

I have seen step05 in pylith2.2.2subduction2d example, which simulates traction_change according to step01, and then calculates traction_normal according to gravity force. I have some questions about whether this traction_normal does not need to calculate the fault normal component of traction based on the fault Angle.
Other than that, I’m not very clear on the mechanics of elastic_prestep.
Also, why not enable gravity(ravityField) in the subduction2d example.
And why should the tractive force of the fault plane without traction_change be set separately artificially? Is this because the fault plane is separated from the material?
If I wanted to simulate a 2-dimensional event, I would describe it as follows: The model has two layers, the upper crust and the lower mantle. The crust on the right side of the model has a tear in it and is filled with mantle material. How to simulate the crust above the fill area just after the fill (the initial state) : The crust above the fill area has settled in 1000 years.
Do I have to turn on the gravity force option, or how do I set the initial state.
There is a simple diagram at the end of the article:

[baagaard]

For uniform or horizontally stratified materials, we usually assume that gravitational body forces create an overburden pressure that results in equal, compressive axial stresses (sigma_xx = sigma_yy.= sigma_zz). Consequently, the normal traction on a fault surface is independent of its orientation.

In the absence of fault friction, for a linear elastic material and infinitesimal strain, the deformation is independent of the absolute stress state. The deformation is driven by stress changes. That is, you can add a stress field that satisfies the elasticity equation, and no additional deformation will occur. This also means that when using fault friction, you can impose stress that satisfies equilibrium in the domain or just on the fault. In most cases, it is easiest to just impose it on the fault, especially when the stress is just the overburden pressure (as discussed above).

For a cross section of a subduction zone, we often have the oceanic crust subducting beneath the continental crust over the mantle (Subduction Zone (3D) — PyLith 4.2.0 documentation). The stress state in such a cross-section depends on the history of the deformation. You can do a literature search for geodynamic models that simulate the long-term evolution of subduction zones to get estimates of the stress state. Alternatively, you can either assume an initial stress state or estimate it using a separate static model. In PyLith v2.2.2, the examples/2d/gravity/gravity_initstress.cfg example demonstrates this approach.

[Ultrafollow]

Thank you for your reply.

As you said, I want to get an initial gravitational field for my model, but since the model coordinates are difficult to determine, I want to simulate the initial stress conditions.

For example pylith-2.2.2\examples\3d\subduction\step08a.cfg. It sets the simulation time to 0year and then outputs the adjusted stress as the initial stress condition for 08b.

I have a problem with this, because the 08a assumes that the full model uses the density of the mantle, so the change in stress caused by this density difference will cause the model to deform. Then, why can this deformed stress condition be applied to 08b as the initial condition?

Also in the joint simulation of steps 3 and 8 of pylith-2.2.2\examples\2d\gravity, elastic_prestep is turned off, while 3d\subduction\step08a.cfg is not. But they can also generate initial stress fields. I’m confused about the difference between them. I tried to close the elastic_prestep in subduction\step08a.cfg, and the result was that there was a large displacement in the domain. If the elastic_prestep was opened, the displacement could be basically ignored.

[Ultrafollow]

I’m sorry to bother you again. I now want to simulate a gravitational field based on your suggestion, and I have noticed that there is an incompressible material in the newer version of pylith.

So I modify the configuration file in the example pylith-4.1.3/src/pylith-4.1.3/examples/reverse-2d/step03_gravity_incompressible.cfg to get the gravity field. However, in the simulation results, there are still some displacements that cannot be ignored.

And the deformation is exaggerated by a factor of 1000.

[pylithapp.problem.materials]
upper_crust = pylith.materials.IncompressibleElasticity
lower_crust = pylith.materials.IncompressibleElasticity
slab = pylith.materials.IncompressibleElasticity
A_mantle = pylith.materials.IncompressibleElasticity
slab_split = pylith.materials.IncompressibleElasticity
eclogite_top = pylith.materials.IncompressibleElasticity
eclogite_bot = pylith.materials.IncompressibleElasticity

[pylithapp.problem.materials.upper_crust]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_upper_crust.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s): 2.600000e+03  3.600000e+03  1.000000e+16

[pylithapp.problem.materials.lower_crust]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_lower_crust.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s): 2.650000e+03  3.820000e+03  1.000000e+16

[pylithapp.problem.materials.slab]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_slab.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s): 3.250000e+03  4.630000e+03  1.000000e+16

[pylithapp.problem.materials.Asthenosphere_mantle]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_Amantle.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s):3.200000e+03  4.400000e+03  1.000000e+16

[pylithapp.problem.materials.slab_tear]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_slab.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s):3.250000e+03  4.400000e+03  1.000000e+16

[pylithapp.problem.materials.eclogite_top]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_eclogite.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s):3.400000e+03  4.500000e+03  1.000000e+16

[pylithapp.problem.materials.eclogite_bot]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_eclogite.spatialdb
# density  vs  vp(kg/m**3  m/s  m/s):3.400000e+03  4.500000e+03  1.000000e+16

[baagaard]

We strongly recommend that you use PyLith v4.2.0 rather than v4.1.3 as v4.2.0 contains some important bug fixes.

If the model contains lateral variations in the density, then there will be some deformation even with incompressible elasticity. I suggest looking at the stress and strain fields.

If you are unsure if you are getting reasonable displacements, then I suggest starting with uniform material properties and transitioning to your heterogeneous material properties step-by-step. You could adjust your the material properties one material at a time and adjust the material properties in steps from the uniform values to the desired values.

[Ultrafollow]

I’m sorry to bother you again.

For this simulation, I took your advice and used version 4.2 of pylith to simulate the initial force field. Then transfer the output file to pylith2.2.2 and extract the data into the spatialdb folder with the python script. This allows the fault to rupture spontaneously.
Besides, I have a question. As shown in the picture below:

This figure is from the paper Diverse Styles of Lithospheric Dripping: Synthesizing Gravitational Instability Models, Continental Tectonics, and Geologic Observations - McMillan - 2023 - Geochemistry, Geophysics, Geosystems - Wiley Online Library.

The figure shows the surface uplift and subsidence caused by delamination and dripping, in which the factors affecting the surface response are – η 'or wavelength.

Note: η '= ηc/ηm denotes the viscosity of the crust ηc scaled to that of the mantle lithosphere ηm.
The principle driving the motion is Rayleigh Taylor instability.

Can pylith simulate this kind of surface uplift? I make η '=1, but it seems like the ground is still settling.

[baagaard]

The types of models you show in the figure are long-term deformation models that involve substantial changes in geometry. Consequently, PyLith is not suitable for those types of models. Other codes, such as ASPECT, are designed for those types of models.

[Ultrafollow]

Sorry to disturb you again, I have a question about cubit modeling.

# -------------------------
# Domain and Fault Parameters
# -------------------------

km = 1000.0

MODEL_CENTER_UTM = (362208.36, 3330507.04, 0)

DOMAIN_X = 450*km 
DOMAIN_Y = 230*km 
DOMAIN_Z = 40*km  
FAULT_DEPTH1 = 15.0*km
FAULT_DEPTH2 = 25.0*km

FAULT_LENGTH = 50*km
FAULT_WIDTH = 50*km
DX_FAULT = 0.5*km
DX_BIAS = 1.07
CELL = "tet"

# -------------------------------------------------------------------------------------------------
# Start cubit and setup domain
# -------------------------------------------------------------------------------------------------
setup_cubit()
cubit.reset()

domain = cubit.brick(DOMAIN_X, DOMAIN_Y, DOMAIN_Z)
cubit.move(domain, (MODEL_CENTER_UTM[0], MODEL_CENTER_UTM[1], -DOMAIN_Z/2))


cubit.create_vertex(337223.59, 3329634.55, -13000)
cubit.create_vertex(387193.13, 3331379.53, -13000)
cubit.create_vertex(387891.826291, 3311372.276970, 0)
cubit.create_vertex(337922.286291, 3309627.296970, 0)
cubit.cmd(f"create curve spline vertex 9 10")
faulttrace = get_last_entity("curve")
cubit.cmd(f"create curve spline vertex 9 12")
cubit.cmd(f"create curve spline vertex 11 12")
cubit.cmd(f"create curve spline vertex 11 10")
cubit.cmd(f"create surface curve 13 14 15 16")

cubit.cmd(f"webcut volume 1 with sheet extended from surface 7")
cubit.cmd(f"imprint volume 1 3 with curve 13 to 26")
cubit.cmd(f"delete surface 7")
cubit.cmd("delete vertex all")
cubit.cmd("imprint all")
cubit.cmd(f"merge all")

Just like the above-mentioned imprint step, the imprint of surface 7 has hidden errors.
And if I set the surface 7 as a 90° vertical fault, then the imprint can be successful.

cubit.create_vertex(337223.59, 3329634.55, -13000)
cubit.create_vertex(387193.13, 3331379.53, -13000)
cubit.create_vertex337223.59, 3329634.55, 0)
cubit.create_vertex(387193.13, 3331379.53, 0)

[MingqianShi]

I used my own script to create this model and encountered the same error.
This is the first time I’ve run into an issue with the imprint command.

[baagaard]

I think the commands below do what you want. The comments note the changes from your lines.

setup_cubit()
cubit.reset()

domain = cubit.brick(DOMAIN_X, DOMAIN_Y, DOMAIN_Z)
cubit.move(domain, (MODEL_CENTER_UTM[0], MODEL_CENTER_UTM[1], -DOMAIN_Z/2))


cubit.create_vertex(337223.59, 3329634.55, -13000)
cubit.create_vertex(387193.13, 3331379.53, -13000)
cubit.create_vertex(387891.826291, 3311372.276970, 0)
cubit.create_vertex(337922.286291, 3309627.296970, 0)
cubit.cmd(f"create curve spline vertex 9 10")
cubit.cmd(f"create curve spline vertex 9 12")
cubit.cmd(f"create curve spline vertex 11 12")
cubit.cmd(f"create curve spline vertex 11 10")
cubit.cmd(f"create surface curve 13 14 15 16")

cubit.cmd(f"webcut volume 1 with sheet extended from surface 7")

# Imprinting the volumes doesn't seem to work, but imprinting the curves from surface 7 onto surface 13 (surface in voiume) does
cubit.cmd(f"imprint surface 13 with curve 13 to 16")

# Delete the body containing surface 7
cubit.cmd(f"delete body 2")

cubit.cmd("delete vertex all")
cubit.cmd("imprint all")
cubit.cmd(f"merge all")

[Ultrafollow]

Dear sir. As shown in the picture (fig1), I want to apply an additional force to the area marked with the black diagonal line. The source of this force is the gravitational gap of the phase transition in this region.

I have attempted to independently simulate the gravitational field of an incompressible material using pylith4.2. Keep the material properties of the slab at the values before the phase transformation.:

[pylithapp.problem]
gravity_field = spatialdata.spatialdb.GravityField
gravity_field.gravity_dir = [0.0, -1.0, 0.0]

solution = pylith.problems.SolnDispPres

[pylithapp.problem.materials]
in_up_crust = pylith.materials.IncompressibleElasticity
zang_up_crust = pylith.materials.IncompressibleElasticity
......(more materials)
[pylithapp.problem.materials.in_up_crust]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_in_up_crust.spatialdb

[pylithapp.problem.materials.slab]
db_auxiliary_field.iohandler.filename = ./spatialdb/mat_in_low_crust.spatialdb
......(more materials)

Then apply the simulation result of 4.2 to pylith2.2.2. Finally, set the property of the slab material in the pylith2.2.2 configuration file to the value of eclogite.

[pylithapp.timedependent.materials.in_up_crust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in in_up_crust
db_initial_stress.iohandler.filename = spatialdb/initial_field-in_up_crust.spatialdb
db_initial_stress.query_type = nearest

[pylithapp.timedependent.materials.in_low_crust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in in_low_crust
db_initial_stress.iohandler.filename = spatialdb/initial_field-in_low_crust.spatialdb
db_initial_stress.query_type = nearest

db_initial_state = spatialdata.spatialdb.SimpleDB
db_initial_state.label = Initial state in in_low_crust
db_initial_state.iohandler.filename = spatialdb/initial_field-in_low_crust.spatialdb
db_initial_state.query_type = nearest
......(more materials)

However, such a method has several drawbacks: Firstly, the model is not laterally uniform, so the simulation results of 4.2 will have a slight displacement; Secondly, if there are open boundaries in the 2.2.2 simulation, the materials will be squeezed out (fig2). In this case, I can only expand the model to reduce the impact caused by the boundaries.

At your convenience, would you please inform me of other alternative solutions. Thank you very much.

[baagaard]

To estimate the gravitational stresses using the incompressible elasticity formulation in PyLith v4.2, you should ensure that the shear modulus is the same for the both the incompressible and compressible cases.

I am not sure what you are trying to simulate after applying the gravitational stresses. It is true that you need to confine the material for consistency with the gravitational loading.

I recommend using a very simple example to determine the workflow for doing the type of simulations you ate trying to do. You may want to start with a simple small box and calculate the analytical solution for each step and then implementing the problems in PyLith. This should make it easier to debug the workflow as you can check the displacements, stresses, and strains between the numerical model and the analytical solution at each step.

[Ultrafollow]

Sorry, I might not have made it clear.

What I want to obtain is the gravitational difference of a subducting slab and see the influence of this gravitational difference on the fault area of the upper crust.

The driving force of this simulation is the fault slip rate specified by the lower boundary of the slab, and the remaining faults are spontaneous rapture.

For this, I have taken the following steps:

1. First, simulate the gravity of state a of the subducting slab material using version 4.2.
2. Put the simulation result of step 1 into 2.2.2.
3. Modify the material parameters of the subducting plate to state b.

In this way, state b does not match the gravitational field of the imported state a. Because state b has a greater density and the S-wave velocity is faster, the subducting slab will be pulled down by the excess gravity.