Step 2: No faults with flexure#
Features
Quasi-static problem
triangular cells
LU preconditioner
pylith.materials.Poroelasticity
pylith.meshio.MeshIOPetsc
pylith.problems.TimeDependent
pylith.problems.SolnDispPresTracStrainVelPdotTdot
pylith.problems.InitialConditionDomain
pylith.bc.DirichletTimeDependent
pylith.bc.NeumannTimeDependent
pylith.meshio.DataWriterHDF5
spatialdata.spatialdb.SimpleGridDB
Simulation parameters#
This example uses poroelasticity to model the infiltration of seawater through a slab of oceanic lithosphere. The permeability field is depth dependent, decreasing with depth but does not vary laterally. The lithosphere is now subject to deformation, over 300 kyr the slab bends to simulate extensional stresses in the outer-rise of a subduction zone. A fluid pressure is applied to the top boundary that is equivalent to the pressure exerted on the seafloor by the water column. This simulates what the hydration state of the oceanic lithosphere as it is about to enter a convergent margin.
Fig. 171 shows the boundary conditions on the domain.
The parameters specific to this example are in step02_no_faults_flexure.cfg.
Fig. 171 Boundary and initial conditions for Step 2.
We fix the left boundary, but we now apply a spatially varying velocity condition on the top boundary using a SimpleDB file, while leaving the right and bottom boundaries unconstrained.
We impose a fluid pressure on the +y boundary equal to the weight of the water column to generate fluid flow.#
[pylithapp.problem.bc.boundary_top]
use_initial = False
use_rate = True
db_auxiliary_field = spatialdata.spatialdb.SimpleDB
db_auxiliary_field.description = Dirichlet BC +y boundary
db_auxiliary_field.iohandler.filename = top_velocity_boundary.spatialdb
SimpleGridDB file that does not contain enhanced permeability due to outer rise faults.#[pylithapp.problem]
[pylithapp.problem.materials.slab]
db_auxiliary_field.filename = no_faultzone_permeability.spatialdb
Running the simulation#
$ pylith step02_no_faults_flexure.cfg
software/pylith-debug/lib/python3.11/site-packages/pylith/apps/PyLithApp.py:77:main
-- pylithapp(info)
-- Running on 1 process(es).
software/pylith-debug/lib/python3.11/site-packages/pylith/meshio/MeshIOObj.py:38:read
-- meshiopetsc(info)
-- Reading finite-element mesh
>> /pylith-main/libsrc/pylith/meshio/MeshIO.cc:85:void pylith::meshio::MeshIO::read(pylith::topology::Mesh *, const bool)
-- meshiopetsc(info)
-- Component 'reader': Domain bounding box:
(0, 150000)
(-30000, 0)
# -- many lines omitted --
>> /software/pylith-debug/lib/python3.11/site-packages/pylith/problems/TimeDependent.py:132:run
-- timedependent(info)
-- Solving problem.
0 TS dt 6000. time -6000.
0 SNES Function norm 3.746797722103e+02
Linear solve converged due to CONVERGED_RTOL iterations 25
1 SNES Function norm 2.600147960846e-10
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
1 TS dt 6000. time 0.
0 SNES Function norm 1.697747712084e+02
Linear solve converged due to CONVERGED_RTOL iterations 23
1 SNES Function norm 5.641202104179e-10
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
# -- many lines omitted --
50 TS dt 6000. time 294000.
0 SNES Function norm 1.697753481522e+02
Linear solve converged due to CONVERGED_RTOL iterations 14
1 SNES Function norm 8.568690135687e-10
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
51 TS dt 6000. time 300000.
>> /software/pylith-debug/lib/python3.11/site-packages/pylith/problems/Problem.py:199:finalize
-- timedependent(info)
-- Finalizing problem.
Visualizing the results#
In Fig. 172 we use the pylith_viz utility to visualize the porosity field.
pylith_viz --filenames=output/step02_no_faults_flexure-slab.h5 warp_grid --field=porosity --exaggeration=1 --hide-edges
Fig. 172 Porosity field at the end of the simulation for Step 2.#