Step 8: Slip on Two Faults and Power-law Viscoelastic Materials#
Features
Triangular cells
pylith.meshio.MeshIOPetsc
pylith.problems.TimeDependent
pylith.bc.DirichletTimeDependent
spatialdata.spatialdb.SimpleDB
spatialdata.spatialdb.ZeroDB
pylith.meshio.OutputSolnBoundary
pylith.meshio.DataWriterHDF5
Quasi-static simulation
pylith.materials.Elasticity
pylith.materials.IsotropicPowerLaw
pylith.faults.FaultCohesiveKin
pylith.faults.KinSrcStep
spatialdata.spatialdb.UniformDB
spatialdata.spatialdb.CompositeDB
Simulation parameters#
In this example we replace the linear, isotropic Maxwell viscoelastic bulk rheology in the slab in Step 7 with an isotropic powerlaw viscoelastic rheology.
The other parameters remain the same as those in Step 7.
The parameters specific to this example are in step08_twofaults_powerlaw.cfg.
[pylithapp.problem.materials]
slab.bulk_rheology = pylith.materials.IsotropicPowerLaw
[pylithapp.problem.materials.slab]
db_auxiliary_field = spatialdata.spatialdb.CompositeDB
db_auxiliary_field.description = Power law material properties
bulk_rheology.auxiliary_subfields.power_law_reference_strain_rate.basis_order = 0
bulk_rheology.auxiliary_subfields.power_law_reference_stress.basis_order = 0
bulk_rheology.auxiliary_subfields.power_law_exponent.basis_order = 0
[pylithapp.problem.materials.slab.db_auxiliary_field]
# Elastic properties
values_A = [density, vs, vp]
db_A = spatialdata.spatialdb.SimpleDB
db_A.description = Elastic properties for slab
db_A.iohandler.filename = mat_elastic.spatialdb
# Power law properties
values_B = [
power_law_reference_stress, power_law_reference_strain_rate, power_law_exponent,
viscous_strain_xx, viscous_strain_yy, viscous_strain_zz, viscous_strain_xy,
reference_stress_xx, reference_stress_yy, reference_stress_zz, reference_stress_xy,
reference_strain_xx, reference_strain_yy, reference_strain_zz, reference_strain_xy,
deviatoric_stress_xx, deviatoric_stress_yy, deviatoric_stress_zz, deviatoric_stress_xy
]
db_B = spatialdata.spatialdb.SimpleDB
db_B.description = Material properties specific to power law bulk rheology for the slab
db_B.iohandler.filename = mat_powerlaw.spatialdb
db_B.query_type = linear
Power-law spatial database#
New in v4.0.0.
We use the utility script pylith_powerlaw_gendb (see pylith_powerlaw_gendb) to generate the spatial database mat_powerlaw.spatialdb with the power-law bulk rheology parameters.
We provide the parameters for pylith_powerlaw_gendb in powerlaw_gendb.cfg, which follows the same formatting conventions as the PyLith parameter files.
pylith_powerlaw_gendb powerlaw_gendb.cfg
Running the simulation#
$ pylith step08_twofaults_powerlaw.cfg
# The output should look something like the following.
>> software/pylith-debug/lib/python3.12/site-packages/pylith/apps/PyLithApp.py:79:main
-- info (application-flow)
-- Running on 1 process(es).
# -- many lines omitted --
>> src/cig/pylith/libsrc/pylith/problems/TimeDependent.cc:473:void pylith::problems::TimeDependent::solve()
-- info (application-flow)
-- Component 'timedependent.problem': Solving equations.
0 TS dt 0.2 time -0.2
0 SNES Function norm 6.969252188583e+00
Linear solve converged due to CONVERGED_ATOL iterations 23
1 SNES Function norm 3.410200001998e-01
Linear solve converged due to CONVERGED_ATOL iterations 15
2 SNES Function norm 1.086370548021e-02
Linear solve converged due to CONVERGED_ATOL iterations 11
3 SNES Function norm 3.546734950165e-04
Linear solve converged due to CONVERGED_ATOL iterations 7
4 SNES Function norm 1.156742268560e-05
Linear solve converged due to CONVERGED_ATOL iterations 3
5 SNES Function norm 3.793910664255e-07
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 5
# -- many lines omitted --
25 TS dt 0.2 time 4.8
0 SNES Function norm 8.706014661840e-02
Linear solve converged due to CONVERGED_ATOL iterations 15
1 SNES Function norm 1.154191339686e-03
Linear solve converged due to CONVERGED_ATOL iterations 9
2 SNES Function norm 2.319736370610e-05
Linear solve converged due to CONVERGED_ATOL iterations 4
3 SNES Function norm 5.416653473355e-07
Linear solve converged due to CONVERGED_ATOL iterations 1
4 SNES Function norm 3.211036236790e-08
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 4
26 TS dt 0.2 time 5.
>> software/pylith-debug/lib/python3.12/site-packages/pylith/problems/Problem.py:222:finalize
-- info (application-flow)
-- Finalizing problem.
As in Step 7, the simulation advances 26 time steps. With a nonlinear bulk rheology, the nonlinear solver now requires a few iterations to converge at each time step.
Visualizing the results#
In Fig. 111 we use the pylith_viz utility to visualize the x displacement field.
You can move the slider or use the p and n keys to change the increment or decrement time.
pylith_viz --filename=output/step08_twofaults_powerlaw-domain.h5 warp_grid --component=x
Fig. 111 Solution for Step 8 at t=100 years. The colors of the shaded surface indicate the x displacement, and the deformation is exaggerated by a factor of 1000. The undeformed configuration is shown by the gray wireframe. Our parameters for the power-law bulk rheology result in much less viscoelastic relaxation in this case compared to Step 7.#
Step 8b: Adaptive Time Stepping#
In Step 8b we demonstrate the use adaptive time stepping. We start with an initial time step of 0.2 years and let the adaptive time stepping algorithm increase the time step as the rate of deformation decreases.
$ pylith step08_twofaults_maxwell.cfg step08b_twofaults_maxwell.cfg
# The output should look something like the following.
>> software/pylith-debug/lib/python3.12/site-packages/pylith/apps/PyLithApp.py:79:main
-- info (application-flow)
-- Running on 1 process(es).
# -- many lines omitted --
>> src/cig/pylith/libsrc/pylith/problems/TimeDependent.cc:473:void pylith::problems::TimeDependent::solve()
-- info (application-flow)
-- Component 'timedependent.problem': Solving equations.
0 TS dt 0.01 time -0.01
0 SNES Function norm 6.969252188583e+00
Linear solve converged due to CONVERGED_ATOL iterations 23
1 SNES Function norm 2.536646382544e-02
Linear solve converged due to CONVERGED_ATOL iterations 11
2 SNES Function norm 8.123236908769e-05
Linear solve converged due to CONVERGED_ATOL iterations 5
3 SNES Function norm 2.954288172419e-07
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 3
TSAdapt basic beuler 0: step 0 accepted t=-0.01 + 1.000e-02 dt=1.000e-02
# -- many lines omitted --
19 TS dt 0.776784 time 4.22322
0 SNES Function norm 3.212286568880e-01
Linear solve converged due to CONVERGED_ATOL iterations 17
1 SNES Function norm 8.434433734286e-03
Linear solve converged due to CONVERGED_ATOL iterations 12
2 SNES Function norm 2.503978690373e-04
Linear solve converged due to CONVERGED_ATOL iterations 7
3 SNES Function norm 9.325515431121e-06
Linear solve converged due to CONVERGED_ATOL iterations 3
4 SNES Function norm 3.575321328214e-07
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 4
TSAdapt basic beuler 0: step 19 accepted t=4.22322 + 7.768e-01 dt=2.634e+00 wlte=0.00348 wltea= -1 wlter= -1
20 TS dt 2.63361 time 5.
>> software/pylith-debug/lib/python3.12/site-packages/pylith/problems/Problem.py:222:finalize
-- info (application-flow)
-- Finalizing problem.
The deformation rate is smoother with the power-law rheology yielding more efficient results with the adaptive time stepper which uses 6 fewer time steps compared with the uniform time step of 4 years.
./plot_compare.py --step=8
Fig. 112 Displacement and shear stress time histories for Step 8 and 8b at a location in the viscoelastic slab below the bottom ends of the main fault and splay fault. The adaptive time stepping does a better job of capturing the variable rate of deformation associated with the power-law viscoelastic rheology; however, the simulation steps over the time of the earthquake on the splay fault before reducing the time step.#