Gravity and initial stress

[ShaozhuoLiu]

Dear Pylith developers and users,

Now I am trying to add the effects of gravity. I read the manual carefully and tested section 7.9.8 (Step15.cfg and Step16.cfg), section 7.17, and section 7.18.12. Those examples work well! So, I tried to test my case.

The domain is meshed by tetrahedral elements. The left-, right-, front-, and back surfaces are fixed in the normal directions of those faces. The bottom is fixed in the vertical direction. Body force caused by gravity is considered. The initial stress field was set according to the density structure.

[pylithapp]
[pylithapp.problem]
bc = [x_pos,x_neg,y_pos,y_neg,z_neg]
gravity_field = spatialdata.spatialdb.GravityField
implicit.solver = pylith.problems.SolverLinear

[pylithapp.problem.formulation.time_step]
total_time = 0.0year
dt = 0.05year

[pylithapp.problem.materials]
UpperCrust = pylith.materials.ElasticIsotropic3D
… …

#for simplicity, other inputs for the lower crust and upper mantle are not shown)

UpperCrust

[pylithapp.problem.materials.UpperCrust]
db_initial_stress = spatialdata.spatialdb.SimpleDB
db_initial_stress.label = Initial stress in UpperCrust
db_initial_stress.iohandler.filename = spatialdb/mat_initial_stress_grav.spatialdb
db_initial_stress.query_type = linear
db_properties = spatialdata.spatialdb.SimpleDB
db_properties.label = Properties for UpperCrust
db_properties.iohandler.filename = spatialdb/mat_elastic.spatialdb
… …

#-----------------------------------------------------------------------

+x face

[pylithapp.problem.bc.x_pos]
bc_dof = [0]
label = face_xpos
db_initial.label = Dirichlet BC on +x

-x face

[pylithapp.problem.bc.x_neg]
bc_dof = [0]
label = face_xneg
db_initial.label = Dirichlet BC on -x

+y face

[pylithapp.problem.bc.y_pos]
bc_dof = [1]
label = face_ypos
db_initial.label = Dirichlet BC on +y

-y face

[pylithapp.problem.bc.y_neg]
bc_dof = [1]
label = face_yneg
db_initial.label = Dirichlet BC on -y

-z face

[pylithapp.problem.bc.z_neg]
bc_dof = [2]
label = face_zneg
db_initial.label = Dirichlet BC on -z

[pylithapp.problem.formulation.output.domain]
writer.filename = output_prestress/test.h5
[pylithapp.problem.materials.UpperCrust.output]
writer.filename = output_prestress/UpperCrust.h5

The resultant simulation does not show subsidence but uplift! I checked the density distribution and found that my settings are right (The density increases with depth!). So does the initial stress setting:
#SPATIAL.ascii 1
SimpleDB {
num-values = 6 // number of stress components
value-names = stress-xx stress-yy stress-zz stress-xy stress-yz stress-xz
value-units = GPa GPa GPa GPa GPa GPa // units
num-locs = 20 // number of locations
data-dim = 1 // 1D variations
space-dim = 3
cs-data = cartesian {
to-meters = 1.0e+3 // specify coordinates in km
space-dim = 3
}
}
// Columns are
// (1) x coordinate (km)
// (2) y coordinate (km)
// (3) z coordinate (km)
// (4) stress-xx ¶
// (5) stress-yy ¶
// (6) stress-zz ¶
// (7) stress-xy ¶
// (8) stress-yz ¶
// (9) stress-xz ¶
0.0 0.0 0.0 -.000000E+00 -.000000E+00 -.000000E+00 0.0 0.0 0.0
0.0 0.0 -2.0 -0.052564 -0.052564 -0.052564 0.0 0.0 0.0
0.0 0.0 -10.0 -0.266741 -0.266741 -0.266741 0.0 0.0 0.0
0.0 0.0 -20.0 -0.539366 -0.539366 -0.539366 0.0 0.0 0.0
0.0 0.0 -35.0 -0.961542 -0.961542 -0.961542 0.0 0.0 0.0
0.0 0.0 -50.0 -1.392544 -1.392544 -1.392544 0.0 0.0 0.0
0.0 0.0 -68.0 -1.988651 -1.988651 -1.988651 0.0 0.0 0.0
0.0 0.0 -80.0 -2.386056 -2.386056 -2.386056 0.0 0.0 0.0
0.0 0.0 -115.0 -3.545153 -3.545153 -3.545153 0.0 0.0 0.0
0.0 0.0 -150.0 -4.704250 -4.704250 -4.704250 0.0 0.0 0.0
0.0 0.0 -155.0 -4.869835 -4.869835 -4.869835 0.0 0.0 0.0
0.0 0.0 -165.0 -5.201006 -5.201006 -5.201006 0.0 0.0 0.0
0.0 0.0 -175.0 -5.532176 -5.532176 -5.532176 0.0 0.0 0.0
0.0 0.0 -185.0 -5.863347 -5.863347 -5.863347 0.0 0.0 0.0
0.0 0.0 -200.0 -6.360103 -6.360103 -6.360103 0.0 0.0 0.0
0.0 0.0 -220.0 -7.022444 -7.022444 -7.022444 0.0 0.0 0.0
0.0 0.0 -265.0 -8.512711 -8.512711 -8.512711 0.0 0.0 0.0
0.0 0.0 -310.0 -10.002979 -10.002979 -10.002979 0.0 0.0 0.0
0.0 0.0 -355.0 -11.493247 -11.493247 -11.493247 0.0 0.0 0.0
0.0 0.0 -400.0 -12.983514 -12.983514 -12.983514 0.0 0.0 0.0

Using the model-produced stress field as a new initial stress field (and I chose the option “db_initial_stress.query_type = nearest”), I re-run the simulation. The resultant vertical displacement field looks like

It does not show that most regions have ZERO z-displacement. However, Figure 7.93 has that pattern.

I am puzzled by what I obtained. Does someone else meet similar issues?

Best,
SZ

PS:
in the pylithapp.cfg, some settings are as follows
… …

[pylithapp.problem.materials.UpperCrust]
label = UpperCrust
id = 1
quadrature.cell = pylith.feassemble.FIATSimplex
quadrature.cell.dimension = 3
output.cell_filter = pylith.meshio.CellFilterAvg
output.writer = pylith.meshio.DataWriterHDF5
db_properties = spatialdata.spatialdb.SimpleDB
db_properties.label = Properties for UpperCrust
db_properties.iohandler.filename = spatialdb/mat_elasticspatialdb
… …

[pylithapp.problem.materials.UpperCrust.output]
cell_filter = pylith.meshio.CellFilterAvg
output_freq = time_step
time_step = 1*year
writer = pylith.meshio.DataWriterHDF5
cell_data_fields = [stress,total_strain]

[baagaard]

Do you have uniform density over the entire domain?

Assuming you have uniform density over the entire domain…Gravitational body forces results in a linear variation in stress. Tetrahedral cells with linear basis functions only capture linear variations in displacement, which means the strain and stress are uniform within each tetrahedron. As your mesh size increases with depth, the discretization error will increase. This may be one of the reasons for the issues you see. Additionally, when you vary the mesh size laterally, the discretization error will also vary, which can introduce additional lateral variations in the displacement field that would not be present with a more uniform mesh. You may need to adjust your mesh if you want to resolve the gravitational body forces more accurately.

[ShaozhuoLiu]

Hi Brad,

Thank you for your answer!

In my case, the density has both lateral and vertical variations. Regarding the vertical variation, we can see how the density increases with depth based on the initial stress setting.

So, an intrinsic mismatch does exist between the spatial variations of the gravitation-induced body force and that of the strain/stress in each element. I don’t know such a fundamental issue.

When comparing my results with those in Figure 7.92-93, I suspected a possibility that the weird results I obtained arise from a bad mesh (Figure 7.92-93 has very nice meshes, with almost uniformly sized tetrahedrons). Now, with your professional explanations, I have more confidence in that and I will test a new case with a better mesh.

I will post new results later.

Best,
SZ

[ShaozhuoLiu]

Hi Brad,

I tested more using better meshes. However, I still don’t see subsidence in the shallow portion when gravity is turned on.

Here is the new mesh.

A fault is wedged into the model, but the rupture on that fault is not considered for my tests (my purpose is to check the effects of gravity and initial stress). The model is divided into several blocks. Each block has apsect ratios smaller than 5 (the previous mesh has aspect ratio up to 70!). The mesh quality is much better than the previous one.

Here is the Z-displacement after turning on the gravity.

On the topmost surface, the two major red regions indicate two major domains with different elastic properties.

Here is the Z-displacement using the initial stress generated by the previous step.

I tested my cfg file and later just used the modified step08a.cfg and step08b.cfg. All give me the same weird results.

For such simple elastic cases, why can’t I get nice results as shown in Figure 7.92 and 7.93?

Best,
SZ

[baagaard]

The displacements in your bottom figure show z displacements on the order of 1.0e-5. This is about what I would expect for stresses that are balancing the gravitational body forces.

The displacements in your middle figure appear to show uplift everywhere on the top surface, so your initial stresses appear to be consistent with a density that is larger than your average density at any given depth. This may not be a significant issue if your resulting stress field is smooth.

I recommend checking the stress field to verify that it is relatively smooth. This gives you the best indication that you have sufficient resolution to represent the stresses due to gravity.

[ShaozhuoLiu]

Hi Brad,

You have sharp eyes and a nice physical intuition!
The second paragraph you wrote gives the answer.

I checked carefully the internal inconsistency between the density and initial stress in my previous tests. In my last post, the density has two major jumps across the upper crust-lower crust and lower crust-upper mantle boundaries; the initial stress was only set for two surfaces: one on the topmost surface, the other at the bottom. The linear interpolation of the input initial stress for lower crust section gives sigma_zz ~20% larger than that consistent with the input density structure (N.B.: Here, the crust can be 65 km thick.).

After correcting such an internal incosistency, the results look much better!

My tests reveal that even ~3%-5% error in setting the initial stress gives me quite bad results (no subsidence at all on the shallow portion). So, once we consider using PREM-type model to constrain the density and initial stress, how we mesh the domain really matters!

Best regards,
SZ

[ShaozhuoLiu]

Hi Brad,

In example 7.17, we use generate_statedb.py to get the initial stress. Here, for the 2D quadrilateral element, we use a basis function of the form

Basis functions for quad4 cell evaluated at quadrature points. Use

to compute coordinate of quadrature points in each cell from

coordinates of vertices. Note the order must correspond to the order

of the data at the quadrature points in the output.

qpts = numpy.array([[ 0.62200847, 0.16666667, 0.0446582, 0.16666667],
[ 0.16666667, 0.62200847, 0.16666667, 0.0446582 ],
[ 0.16666667, 0.0446582, 0.16666667, 0.62200847],
[ 0.0446582, 0.16666667, 0.62200847, 0.16666667]], dtype=numpy.float64)

How can we get the corresponding basis function for 3D tetrahedron and hexahedron?

Best,
SZ

[baagaard]

The quadrature information is setup in the FIATLagrange and FIATSimplex objects. You can use the journals to turn on dumping the information to the screen in a PyLith simulation using

[pylithapp.journal.info]
fiatsimplex = 1
fiatlagrange = 1

Here is the answer for fiatlagrange and fiatsimplex in 3D:

 >> fiatlagrange(info)
 -- Cell geometry:
 -- <pylith.feassemble.CellGeometry.GeometryHex3D; proxy of <Swig Object of type 'pylith::feassemble::GeometryHex3D *' at 0x10581a360> >
 -- Vertices:
 -- [[-1. -1. -1.]
 [-1.  1. -1.]
 [ 1.  1. -1.]
 [ 1. -1. -1.]
 [-1. -1.  1.]
 [ 1. -1.  1.]
 [ 1.  1.  1.]
 [-1.  1.  1.]]
 -- Quad pts:
 -- [[-0.57735027 -0.57735027 -0.57735027]
 [ 0.57735027 -0.57735027 -0.57735027]
 [-0.57735027  0.57735027 -0.57735027]
 [ 0.57735027  0.57735027 -0.57735027]
 [-0.57735027 -0.57735027  0.57735027]
 [ 0.57735027 -0.57735027  0.57735027]
 [-0.57735027  0.57735027  0.57735027]
 [ 0.57735027  0.57735027  0.57735027]]

 >> fiatsimplex(info)
 -- Cell geometry:
 -- <pylith.feassemble.CellGeometry.GeometryTet3D; proxy of <Swig Object of type 'pylith::feassemble::GeometryTet3D *' at 0x108d3a360> >
 -- Vertices:
 -- [[ 1. -1. -1.]
 [-1. -1. -1.]
 [-1.  1. -1.]
 [-1. -1.  1.]]
 -- Quad pts:
 -- [[-0.5 -0.5 -0.5]]

This information can be extracted in a few lines of Python code and I will post a script that shows how to do this. I need to fix one issue before posting.

[baagaard]

Here is the Python script. You can save it to a file (cell_info.py), make the file executable (chmod +x cell_info.py) and then run it directly cell_info.py.

#!/usr/bin/env mpinemesis

from pylith.feassemble.FIATLagrange import FIATLagrange
from pylith.feassemble.FIATSimplex import FIATSimplex
import journal

# Turn on journals that will display the cell information
journal.info("fiatlagrange").activate()
journal.info("fiatsimplex").activate()

def initialize_cell(cell, dim):
    """Initialize the reference cell."""
    # Set the appropriate Pyre properties for the reference cell.
    # In a PyLith simulation this is normally done in a .cfg fle.
    cell.inventory.dimension = dim
    cell._configure()
    cell.initialize(spaceDim=dim)

initialize_cell(FIATLagrange(), 2) # Quad
initialize_cell(FIATLagrange(), 3) # Hex
initialize_cell(FIATSimplex(), 2) # Tri
initialize_cell(FIATSimplex(), 3) # Tet