Seeking assistance with balancing the gravity field in Pylith 2.2.2

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:

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)

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.