Seeking assistance with balancing the gravity field in Pylith 2.2.2

PyLith
1

Dear Pylith community,

I hope this email finds you well. I am currently facing an issue while using Pylith 2.2.2. In my model, which considers the effect of gravity, I am trying to apply boundary stress loading or an initial stress field that would balance out the initial strain caused by gravity, ensuring that the strain is solely controlled by tectonic loading.

However, the approaches I have attempted, such as applying static lithostatic pressure at both ends and an initial static lithostatic pressure field, have not successfully eliminated the gravity field. I am seeking assistance on how to incorporate the element of balancing the gravity field. Here are the details:

Boundary condition code:
[pylithapp.timedependent.bc.y_pos_1]
bc_dof = [0,1]
label = boundary_ypos
db_initial.label = Dirichlet BC on boundary_ypos

[pylithapp.timedependent.bc.y_neg_1]
bc_dof = [0,1]
label = boundary_yneg
db_initial.label = Dirichlet BC on boundary_yneg
db_rate = spatialdata.spatialdb.UniformDB
db_rate.label = Dirichlet rate BC on boundary_yneg
db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
db_rate.data = [5cm/year,0cm/year,0.0*year]

[pylithapp.timedependent.bc]
x_pos_2=pylith.bc.Neumann
x_neg_2=pylith.bc.Neumann

[pylithapp.timedependent.bc.x_pos_2]
label = boundary_xpos
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC on boundary_xpos
db_initial.iohandler.filename = spatialdb/xpos_srp.spatialdb
db_initial.query_type=linear
output.cell_info_fields = [initial_value]
output.writer.filename = output/step09-traction.vtk
output.cell_filter = pylith.meshio.CellFilterAvg
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2

[pylithapp.timedependent.bc.x_pos_1]
bc_dof = [1]
label = boundary_xpos
db_initial.label = Dirichlet BC on boundary_xpos

[pylithapp.timedependent.bc.x_neg_2]
label = boundary_xneg
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC on boundary_xneg
db_initial.iohandler.filename = spatialdb/xneg_srp.spatialdb
db_initial.query_type=linear
output.cell_info_fields = [initial_value]
output.writer.filename = output/step09-traction.vtk
output.cell_filter = pylith.meshio.CellFilterAvg
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2

[pylithapp.timedependent.bc.x_neg_1]
bc_dof = [1]
label = boundary_xneg
db_initial.label = Dirichlet BC on boundary_xneg

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

Initial Stress field code:
[pylithapp.timedependent.materials.upper_crust1]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear

[pylithapp.timedependent.materials.upper_crust2]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear

[pylithapp.timedependent.materials.fault]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear

[pylithapp.timedependent.materials.lower_crust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear

result picture:

重初边
重初边1330×787 59.8 KB

I would greatly appreciate it if you could provide any insights or suggestions on how to address this issue. Your expertise and guidance would be invaluable to me.

Thank you very much for your time and assistance.

Best regards

Stress field file:
xneg_srp.txt (795 Bytes)
initial_stress1.txt (959 Bytes)

2

It would be helpful if you could supply a diagram of the boundary value problem you are trying to solve, including boundary conditions. Include information about the spatial distribution for all initial conditions and boundary conditions.

3

Thank you for your reply. I have found the cause of this problem. There is a mistake in the pylith-2.2.2_manual.pdf file, and the boundary conditions I constructed with the gravity acceleration value (9.80655) are wrong, so the final displacement field results show anomalies. The correct value of gravity acceleration should be 9.80665.

4

Now, however, there is a new puzzle. The stress distribution on the fault plane that I solved with pylith2.2.2 is extremely discontinuous. The elastic-plastic material I use and the corresponding nonlinear solver, step09.cfg in example/3d, which is referenced in the configuration file, have a very long calculation time (up to one month). Could you help me check whether the parameter setting of the configuration file is improper, or is it caused by other reasons such as improper grid of the model? Thank you very much!
step09.cfg (13.6 KB)

fault_tracton
fault_tracton1842×773 83.5 KB

pylithapp.cfg (3.5 KB)

5

It would be helpful if you could supply a diagram of the boundary value problem you are trying to solve, including boundary conditions. Include information about the spatial distribution for all initial conditions and boundary conditions.

6

The boundary and initial stress fields are set as follows:

Boundary condition

bc = [x_pos_1,x_neg_1,x_pos_2,x_neg_2,y_pos_1,y_neg_1,z_neg_1,yz_1]

+y face

[pylithapp.timedependent.bc.y_pos_1]
bc_dof = [0,1]
label = boundary_ypos
db_initial.label = Dirichlet BC on boundary_ypos

-y face

[pylithapp.timedependent.bc.y_neg_1]
bc_dof = [0,1]
label = boundary_yneg
db_initial.label = Dirichlet BC on boundary_yneg
db_rate = spatialdata.spatialdb.UniformDB
db_rate.label = Dirichlet rate BC on boundary_yneg
db_rate.values = [displacement-rate-x,displacement-rate-y,rate-start-time]
db_rate.data = [5cm/year,0cm/year,0.0*year]
[pylithapp.timedependent.bc]
x_pos_2=pylith.bc.Neumann
x_neg_2=pylith.bc.Neumann

+x face

[pylithapp.timedependent.bc.x_pos_2]
label = boundary_xpos
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC on boundary_xpos
db_initial.iohandler.filename = spatialdb/xpos_srp.spatialdb
db_initial.query_type=linear
output.cell_info_fields = [initial_value]
output.writer.filename = output/step09-traction.vtk
output.cell_filter = pylith.meshio.CellFilterAvg
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2
[pylithapp.timedependent.bc.x_pos_1]
bc_dof = [1]
label = boundary_xpos_noyz
db_initial.label = Dirichlet BC on boundary_xpos_noyz

-x face

[pylithapp.timedependent.bc.x_neg_2]
label = boundary_xneg
db_initial = spatialdata.spatialdb.SimpleDB
db_initial.label = Neumann BC on boundary_xneg
db_initial.iohandler.filename = spatialdb/xneg_srp.spatialdb
db_initial.query_type=linear
output.cell_info_fields = [initial_value]
output.writer.filename = output/step09-traction.vtk
output.cell_filter = pylith.meshio.CellFilterAvg
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 2
quadrature.cell.quad_order = 2
[pylithapp.timedependent.bc.x_neg_1]
bc_dof = [1]
label = boundary_xneg_noyz
db_initial.label = Dirichlet BC on boundary_xneg_noyz

-z face

[pylithapp.timedependent.bc.z_neg_1]
bc_dof = [2]
label = boundary_zneg
db_initial.label = Dirichlet BC on boundary_zneg
[pylithapp.timedependent.bc.yz_1]
bc_dof = [1,2]
label = boundary_yz
db_initial.label = Dirichlet BC on boundary_yz
Initial stress field setting:

initial stresses

[pylithapp.timedependent.materials.upper_crust1]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear
[pylithapp.timedependent.materials.upper_crust2]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear
[pylithapp.timedependent.materials.fault]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear
[pylithapp.timedependent.materials.lower_crust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in model
db_initial_stress.iohandler.filename = spatialdb/initial_stress1.spatialdb
db_initial_stress.query_type = linear
xpos_srp.txt (800 Bytes)
xneg_sr.txt (791 Bytes)
initial_stress1.txt (989 Bytes)

7

These are all of the parameter settings you are using, but we want to make sure your parameter settings are matching the boundary value problem you are wanting to solve. Please provide a sketch or diagram of the boundary value problem you are trying to solve, including boundary conditions. Include information about the spatial distribution for all initial conditions and boundary conditions.

8

I’m sorry I didn’t understand you the first time. The attached picture is a sketch of my model and boundary conditions. The upper layer is the upper crust made of elastic-plastic material, and the lower layer is the lower crust and upper mantle made of viscoelastic material. The yellow line in the middle is two faults, which are set to dynamic static friction conditions. I want to know the stress evolution on the fault under this long-term shear loading state.

dcfbd5187da5596a08c1ed64753ce48
dcfbd5187da5596a08c1ed64753ce481504×846 153 KB

9

Here are some suggestions for improving the problem setup and solution for this type of simulation using PyLith v2.2.x:

  1. Start with a coarser mesh and refine only in areas where you need it. Using a uniform mesh with a fine discretization size will increase the computational cost.
  2. Use the custom field split preconditioner. The settings are in solver_fault_fieldsplit.cfg. This should help reduce the runtime.
  3. Your model has a high aspect ratio. That is, it is quite large in the horizontal direction compared with the vertical direction. In most cases, you will want a domain in which the lateral dimensions are about 2x the vertical dimension, so that all boundaries are roughly equal distant from the region of interest. For this simple geometry, I would use hexahedral cells. Hexahedral cells are slightly more efficient (fewer cells fill a given volume) and they are slightly more accurate.
  4. Start with a purely elastic model and make sure your solver settings are working correctly. This is especially important when using friction. Your SNES solver settings in step09.cfg are very likely contributing to the poor results. They are much too large. I would start with the solver tolerances you have in pylithapp.cfg. You also need to make sure the friction zero tolerance is between the KSP and SNES absolute solver tolerances. Refer to the PyLith v2.2 section 6.4.5.2. You may need to make small changes to these values so that they work for your case.
  5. Carefully examine the solution at each time step to ensure the physics are well resolved. For example, the solution should be smooth until slip starts. Once slip starts you want to make sure that the solution (slip, stresses, etc) are reasonable. If the solver tolerances are too large, then you will likely see spatial variations that do not make sense.