Step 2: Magma inflation with evolution of porosity#
New in v4.0.0
Features
Quasi-static problem
LU preconditioner
pylith.materials.Poroelasticity
pylith.meshio.MeshIOCubit
pylith.problems.TimeDependent
pylith.problems.SolnDispPresTracStrain
pylith.problems.InitialConditionDomain
pylith.bc.DirichletTimeDependent
pylith.bc.NeumannTimeDependent
pylith.meshio.DataWriterHDF5
spatialdata.spatialdb.SimpleGridDB
spatialdata.spatialdb.UniformDB
Poroelasticity with porosity state variable
Isotropic linear poroelasticity with reference state
Simulation parameters#
We extend the simulation in Step 1 by including evolution of the porosity, which depends on the time derivative of the pressure and trace strain. We also compute the deformation relative to a uniform reference compressive pressure of 5 MPa to illustrate how to use a reference state with poroelasticity. We use the same initial conditions and boundary conditions as in Step 1.
Because the evolution of porosity depends on the time derivative of the solution subfields, we need to include the time derivatives in the solution field. As a result, we have 6 subfields in our solution field.
# Poroelasticity with porosity state variable requires solution with time derivatives
solution = pylith.problems.SolnDispPresTracStrainVelPdotTdot
# Set basis order for all solution subfields
[pylithapp.problem.solution.subfields]
displacement.basis_order = 2
pressure.basis_order = 1
trace_strain.basis_order = 1
velocity.basis_order = 2
pressure_t.basis_order = 1
trace_strain_t.basis_order = 1
[pylithapp.problem.materials.crust]
use_state_variables = True
db_auxiliary_field.values = [
solid_density, fluid_density, fluid_viscosity, porosity, shear_modulus, drained_bulk_modulus, biot_coefficient, fluid_bulk_modulus, solid_bulk_modulus, isotropic_permeability,
reference_stress_xx, reference_stress_yy, reference_stress_zz, reference_stress_xy,
reference_strain_xx, reference_strain_yy, reference_strain_zz, reference_strain_xy
]
db_auxiliary_field.data = [
2500*kg/m**3, 1000*kg/m**3, 0.001*Pa*s, 0.01, 6.0*GPa, 10.0*GPa, 1.0, 2.0*GPa, 20.0*GPa, 1e-15*m**2,
-5.0*MPa, -5.0*MPa, -5.0*MPa, 0.0*MPa,
0.0, 0.0, 0.0, 0.0
]
auxiliary_subfields.porosity.basis_order = 1
[pylithapp.problem.materials.crust.bulk_rheology]
use_reference_state = True
[pylithapp.problem.materials.intrusion]
use_state_variables = True
db_auxiliary_field.values = [
solid_density, fluid_density, fluid_viscosity, porosity, shear_modulus, drained_bulk_modulus, biot_coefficient, fluid_bulk_modulus, solid_bulk_modulus, isotropic_permeability,
reference_stress_xx, reference_stress_yy, reference_stress_zz, reference_stress_xy,
reference_strain_xx, reference_strain_yy, reference_strain_zz, reference_strain_xy
]
db_auxiliary_field.data = [
2500*kg/m**3, 1000*kg/m**3, 0.001*Pa*s, 0.1, 6.0*GPa, 10.0*GPa, 0.8, 2.0*GPa, 20.0*GPa, 1e-13*m**2,
-5.0*MPa, -5.0*MPa, -5.0*MPa, 0.0*Pa,
0.0, 0.0, 0.0, 0.0
]
auxiliary_subfields.porosity.basis_order = 1
[pylithapp.problem.materials.intrusion.bulk_rheology]
use_reference_state = True
The changes in the physics as well as the default solver settings that impact the initial guess leads to a large initial residual at the second time step. Consequently, we increase the divergence tolerance for the linear solver.
[pylithapp.petsc]
# Increase divergence tolerance. Initial guess at second time step is not accurate.
ksp_divtol = 1.0e+5
Running the simulation#
$ pylith step02_inflation.cfg
# The output should look something like the following.
>> /home/pylith-user/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/apps/PyLithApp.py:84:main
-- pylithapp(info)
-- Running on 1 process(es).
>> /home/pylith-user/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/meshio/MeshIOObj.py:43:read
-- meshiocubit(info)
-- Reading finite-element mesh
>> /home/pylith-user/src/cig/pylith/libsrc/pylith/meshio/MeshIOCubit.cc:156:void pylith::meshio::MeshIOCubit::_readVertices(ExodusII &, scalar_array *, int *, int *) const
-- meshiocubit(info)
-- Component 'reader': Reading 747 vertices.
# -- many lines omitted --
-- Setting PETSc options:
fieldsplit_displacement_pc_type = lu
fieldsplit_pressure_pc_type = ilu
fieldsplit_pressure_t_pc_type = ilu
fieldsplit_trace_strain_pc_type = ilu
fieldsplit_trace_strain_t_pc_type = ilu
fieldsplit_velocity_pc_type = ilu
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_0_fields = 2
pc_fieldsplit_1_fields = 1
pc_fieldsplit_2_fields = 0
pc_fieldsplit_3_fields = 3
pc_fieldsplit_4_fields = 4
pc_fieldsplit_5_fields = 5
pc_fieldsplit_type = multiplicative
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
>> /home/pylith-user/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/problems/TimeDependent.py:132:run
-- timedependent(info)
-- Solving problem.
0 TS dt 1. time -1.
0 SNES Function norm 7.521665654021e-01
Linear solve converged due to CONVERGED_RTOL iterations 42
1 SNES Function norm 1.098964504117e-10
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
# -- many lines omitted --
50 TS dt 1. time 49.
0 SNES Function norm 1.583815770370e-01
Linear solve converged due to CONVERGED_ATOL iterations 10
1 SNES Function norm 7.897032897534e-12
Nonlinear solve converged due to CONVERGED_FNORM_ABS iterations 1
51 TS dt 1. time 50.
>> /home/pylith-user/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/problems/Problem.py:199:finalize
-- timedependent(info)
-- Finalizing problem.
The linear solver exhibits similar performance with less than 50 iterations at most time steps. Furthermore, the problem is still linear, so the nonlinear solver converges in one iteration.
Visualizing the results#
In Fig. 114 we use ParaView to visualize the evolution of the y displacement component using the viz/plot_dispwarp.py Python script.
First, we start ParaView from the examples/magma-2d directory.
$ PATH_TO_PARAVIEW/paraview
# For macOS, it will be something like
$ /Applications/ParaView-5.10.1.app/Contents/MacOS/paraview
Next, we override the default name of the simulation file with the name of the current simulation.
>>> SIM = "step02_inflation_statevars"
Next we run the viz/plot_dispwarp.py Python script as described in ParaView Python Scripts.
Fig. 114 Solution for Step 2 at t=100 yr. The colors of the shaded surface indicate the fluid pressure, and the deformation is exaggerated by a factor of 1000. The reference state gives rise to greater vertical deformation compared to Step 1. The choice of material properties does not lead to significant changes in porosity in either material during the simulation (not shown in the figure).#