Step 8: Use of Gravitational Body Forces#

This example demonstrates the use of gravitational body forces as well as the use of initial stresses to balance the body forces. This involves enabling gravity within our domain with Dirichlet roller boundary conditions on the lateral and bottom boundaries; we do not include faults in this example. We also demonstrate what happens when the initial stresses are not in balance with the gravitational stresses, and show how viscoelastic problems with gravitational stresses will in general not reach a steady-state solution. This example consists of three sub-problems:

  • Step 8a: Gravitational body forces with 3-D density variations in elastic materials and initial stresses for a uniform density.

  • Step 8b: Gravitational body forces with 3-D density variations in incompressible elastic materials.

  • Step 8c: Gravitational body forces with 3-D density variations in elastic and viscoelastic materials and initial stresses from Step 8a plus finite strain formulation (does not reach a steady-state solution). This example is not yet implemented for PyLith v3+.

Warning

Step 8c is still being updated for use with PyLith v3+.

Danger

Step 8b takes an extremely long time to run and does not presently run in parallel. We are working to select better preconditioners to improve the solvers and reduce the runtime.

All these simulations use the default ZeroDB displacement boundary conditions on the lateral and bottom boundaries, and they all use spatialdata.spatialdb.GravityField to apply gravitational body forces. Step 8a uses reference stresses (provided in the associated spatial databases for each material) to balance the deformation induced by turning on body forces. The reference stresses assume a constant density, while there is actually a density contrast between materials, meaning that the stresses don’t balance and there is some resulting deformation. Step 8b uses incompressible elasticity (using pylith.problems.SolnDispPres) to prevent the vertical deformation that would otherwise occur. In this case it is necessary to provide an additional Dirichlet zero pressure boundary condition on the ground surface.

Fig. 156 shows the boundary conditions on the domain.

Fig. 156 Boundary conditions for gravitational body forces examples. Body forces are applied along with roller boundary conditions on the lateral sides and bottom of the domain. For Step 8b, we also apply zero pressure BC along the upper surface of the mesh.#

Features

  • Tetrahedral cells

  • pylith.meshio.MeshIOCubit

  • pylith.problems.TimeDependent

  • pylith.meshio.OutputSolnBoundary

  • pylith.meshio.DataWriterHDF5

  • pylith.bc.DirichletTimeDependent

  • pylith.bc.ZeroDB

  • spatialdata.geocoords.CSGeo

  • Static simulation

  • spatialdata.spatialdb.GravityField

  • Reference state

Simulation parameters#

The parameters specific to this example are in step08a_gravity_refstate.cfg and step08b_gravity_incompressible.cfg and include:

  • pylithapp.metadata Metadata for this simulation. Even when the author and version are the same for all simulations in a directory, we prefer to keep that metadata in each simulation file as a reminder to keep it up-to-date for each simulation.

  • pylithapp Parameters defining where to write the output.

  • pylithapp.problem Parameters to set quadrature and basis order, as well as provide a gravity_field. For Step 8b, we also set the solution to pylith.problems.SolnDispPres.

  • pylithapp.problem.bc For Step 8b, we add a zero pressure BC on the ground surface.

  • pylithapp.problem.materials For Step 8a, we set use_reference_state and provide spatial databases that include reference stresses.

For all of the problems involving gravitational body forces we use spatialdata.spatialdb.GravityField to define the gravity_field. For Step 8a, we use higher order quadrature_order and basis_order to accurately capture the displacement field. This is not needed for Step 8b, which has much smaller displacements. Note that Step 8b also requires mat_elastic_incompressible.cfg, which provides all of the necessary properties using a UniformDB.

Listing 218 Run Step 8a and 8b simulations#
$ pylith step08a_gravity_refstate.cfg
$ pylith step08b_gravity_incompressible.cfg mat_elastic_incompressible.cfg

# The output should look something like the following. This is for Step 8b.
>> /work/charlesw/virtualenv/python3.12/lib/python3.12/site-packages/pylith/apps/PyLithApp.py:77:main
 -- pylithapp(info)
 -- Running on 1 process(es).
 >> /work/charlesw/virtualenv/python3.12/lib/python3.12/site-packages/pylith/meshio/MeshIOObj.py:38:read
 -- meshiocubit(info)
 -- Reading finite-element mesh
 >> /home/charlesw/cig/pylith3/source/pylith-fork/libsrc/pylith/meshio/MeshIOCubit.cc:148:void pylith::meshio::MeshIOCubit::_readVertices(pylith::meshio::ExodusII&, pylith::scalar_array*, int*, int*) const
 -- meshiocubit(info)
 -- Component 'reader': Reading 24824 vertices.
 >> /home/charlesw/cig/pylith3/source/pylith-fork/libsrc/pylith/meshio/MeshIOCubit.cc:208:void pylith::meshio::MeshIOCubit::_readCells(pylith::meshio::ExodusII&, pylith::int_array*, pylith::int_array*, int*, int*) const
 -- meshiocubit(info)
 -- Component 'reader': Reading 134381 cells in 4 blocks.

# -- many lines omitted --

 >> /home/charlesw/cig/pylith3/source/pylith-fork/libsrc/pylith/utils/PetscOptions.cc:239:static void pylith::utils::_PetscOptions::write(pythia::journal::info_t&, const char*, const pylith::utils::PetscOptions&)
 -- petscoptions(info)
 -- Setting PETSc options:
fieldsplit_displacement_pc_type = lu
fieldsplit_pressure_pc_type = lu
ksp_atol = 1.0e-12
ksp_converged_reason = true
ksp_error_if_not_converged = true
ksp_guess_pod_size = 8
ksp_guess_type = pod
ksp_rtol = 1.0e-12
pc_fieldsplit_schur_factorization_type = full
pc_fieldsplit_schur_precondition = full
pc_fieldsplit_type = schur
pc_type = fieldsplit
snes_atol = 1.0e-9
snes_converged_reason = true
snes_error_if_not_converged = true
snes_monitor = true
snes_rtol = 1.0e-12
ts_error_if_step_fails = true
ts_monitor = true
ts_type = beuler

 >> /work/charlesw/virtualenv/python3.12/lib/python3.12/site-packages/pylith/problems/TimeDependent.py:132:run
 -- timedependent(info)
 -- Solving problem.
0 TS dt 3.15576e+07 time 0.
    0 SNES Function norm 1.847563517818e+03
      Linear solve converged due to CONVERGED_ATOL iterations 2
    1 SNES Function norm 1.534546536971e-09
    Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 1
1 TS dt 3.15576e+07 time 3.15576e+07
 >> /work/charlesw/virtualenv/python3.12/lib/python3.12/site-packages/pylith/problems/Problem.py:199:finalize
 -- timedependent(info)
 -- Finalizing problem.

The beginning of the output is nearly the same as in several previous examples. Both of these problems are static, however, so there is only a single time step.

Visualizing the results#

In Fig. 157 we use the pylith_viz utility to visualize the z displacement field.

Listing 219 Visualize PyLith output using pylith_viz.#
pylith_viz --filename=output/step08a_gravity_refstate-domain.h5 warp_grid --component=z --exaggeration=50
Solution for Step 8a. The colors indicate the z displacement, and the deformation is exaggerated by a factor of 50.

Fig. 157 Solution for Step 8a. The colors of the shaded surface indicate the z displacement, and the deformation is exaggerated by a factor of 50.#

Note that the vertical displacements for Step 8b are very close to zero, because the material is nearly incompressible.