Skip to content

Commit 57359a4

Browse files
baagaard-usgsbaagaard-usgs
authored
Merge pull request #666 from baagaard-usgs/example-poroelasticity-statevars
Add Step 2 to examples/magma-2d demonstrating use of porosity state variable
2 parents 48f264a + 42ce815 commit 57359a4

31 files changed

+1249
-582
lines changed

‎CHANGES.md‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,21 @@
11
See <https://github.com/geodynamics/pylith/commits/main> for the complete log of changes made to PyLith.
22

3-
## Version 3.1.0
3+
## Version 4.0.0
44

55
* Removed support for importing meshes from LaGriT.
66
* Add output of change in fault tractions for prescribed slip.
7-
* By default use PETSc proper orthogonal decomposition (POD) methodology for initial guess of solutions to improve convergence.
7+
* State variables are now included in the default `data_fields` for simulation output.
8+
* The default solver settings use the PETSc proper orthogonal decomposition (POD) methodology for initial guess of solutions to improve convergence.
9+
* Changed name of fault Lagrange multiplier field for solution component in Python from `lagrange_fault` to `lagrange_multiplier_fault` to match name of solution field in C++.
810
* Add demonstration of `pylith_powerlaw_gendb` in Step 8 of `examples/reverse-2d`.
11+
* Add demonstration of using poroelasticity with porosity as a state variable to `examples/magma-2d`.
912
* Switched from CppUnit to Catch2 as the C++ testing framework.
1013
* Update to PETSc 3.19.5
1114
* Improve integration with VSCode for testing and debugging (see Developer Guide)
1215
* Bug fixes
1316
* Fix errors in KinSrcTimeHistory.py
1417
* Fix creation of PETSc label for edges when importing Gmsh files. This fixes creation of faults with buried edges for 3D meshes imported from Gmsh.
18+
* Add containers for solution fields for poroelasticity with faults.
1519

1620
## Version 3.0.3
1721

‎docs/Makefile.am‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -267,6 +267,8 @@ dist_noinst_DATA = \
267267
user/examples/magma-2d/figs/geometry.pdf \
268268
user/examples/magma-2d/figs/step01-diagram.pdf \
269269
user/examples/magma-2d/figs/step01-diagram.svg \
270+
user/examples/magma-2d/figs/step01-solution.jpg \
271+
user/examples/magma-2d/figs/step01-solution.jpg \
270272
user/examples/reverse-2d/index.md \
271273
user/examples/reverse-2d/common-information.md \
272274
user/examples/reverse-2d/exercises.md \

‎docs/user/examples/magma-2d/figs/step02-solution.jpg‎

54.6 KB
Loading

‎docs/user/examples/magma-2d/index.md‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,7 @@ The files and directories for this set of examples includes:
1616
This example demonstrates the use of poroelasticity to model magma flow in a conduit and reservoir in the crust.
1717

1818
:Step 1: Magma influx with displacement and pressure boundary conditions.
19+
:Step 2: Same as Step 1 with evolution of porosity state variable.
1920

2021
:::{figure-md} fig:example:magma:2d:overview
2122
<img src="figs/geometry.*" alt="Diagram of geometry for magma reservoir." scale="100%"/>
@@ -31,5 +32,6 @@ We refer to the domain boundaries using the names shown in the diagram.
3132
meshing-cubit.md
3233
common-information.md
3334
step01-inflation.md
35+
step02-inflation-statevars.md
3436
exercises.md
3537
:::

‎docs/user/examples/magma-2d/step01-inflation.md‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -3,10 +3,10 @@
33
```{include} step01_inflation-synopsis.md
44
```
55

6-
## Simulation parameter
6+
## Simulation parameters
77

88
This example uses poroelasticity to model flow of magma up through a conduit and into a magma reservoir.
9-
The magma reseroir and conduit have a higher permeability than the surrounding crust.
9+
The magma reservoir and conduit have a higher permeability than the surrounding crust.
1010
We generate flow by imposing a pressure on the external boundary of the conduit that is higher than the uniform initial pressure in the domain.
1111
{numref}`fig:example:magma:2d:step01:diagram` shows the boundary conditions on the domain.
1212
The parameters specific to this example are in `step01_inflation.cfg`.
@@ -176,7 +176,7 @@ ts_error_if_step_fails = true
176176
ts_monitor = true
177177
ts_type = beuler
178178
179-
# -- many lines ommitted --
179+
# -- many lines omitted --
180180
181181
50 TS dt 1. time 49.
182182
0 SNES Function norm 3.049429649018e-03

‎docs/user/examples/magma-2d/step01_inflation-synopsis.md‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11

22
**Features**
33

4-
* Quasistatic problem
4+
* Quasi-static problem
55
* LU preconditioner
66
* pylith.materials.Poroelasticity
77
* pylith.meshio.MeshIOCubit

‎docs/user/examples/magma-2d/step02-inflation-statevars.md‎

Lines changed: 259 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,259 @@
1+
# Step 2: Magma inflation with evolution of porosity
2+
3+
```{include} step02_inflation_statevars-synopsis.md
4+
```
5+
6+
## Simulation parameters
7+
8+
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.
9+
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.
10+
We use the same initial conditions and boundary conditions as in Step 1.
11+
12+
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.
13+
As a result, we have 6 subfields in our solution field.
14+
15+
```{code-block} cfg
16+
---
17+
caption: Solution subfields for Step 2.
18+
---
19+
# Poroelasticity with porosity state variable requires solution with time derivatives
20+
solution = pylith.problems.SolnDispPresTracStrainVelPdotTdot
21+
22+
# Set basis order for all solution subfields
23+
[pylithapp.problem.solution.subfields]
24+
displacement.basis_order = 2
25+
pressure.basis_order = 1
26+
trace_strain.basis_order = 1
27+
velocity.basis_order = 2
28+
pressure_t.basis_order = 1
29+
trace_strain_t.basis_order = 1
30+
```
31+
32+
```{code-block} cfg
33+
---
34+
caption: Material parameters for poroelasticity with state variables and reference state for Step 2.
35+
---
36+
[pylithapp.problem.materials.crust]
37+
use_state_variables = True
38+
39+
db_auxiliary_field.values = [
40+
solid_density, fluid_density, fluid_viscosity, porosity, shear_modulus, drained_bulk_modulus, biot_coefficient, fluid_bulk_modulus, solid_bulk_modulus, isotropic_permeability,
41+
reference_stress_xx, reference_stress_yy, reference_stress_zz, reference_stress_xy,
42+
reference_strain_xx, reference_strain_yy, reference_strain_zz, reference_strain_xy
43+
]
44+
db_auxiliary_field.data = [
45+
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,
46+
-5.0*MPa, -5.0*MPa, -5.0*MPa, 0.0*MPa,
47+
0.0, 0.0, 0.0, 0.0
48+
]
49+
50+
auxiliary_subfields.porosity.basis_order = 1
51+
52+
[pylithapp.problem.materials.crust.bulk_rheology]
53+
use_reference_state = True
54+
55+
56+
[pylithapp.problem.materials.intrusion]
57+
use_state_variables = True
58+
59+
db_auxiliary_field.values = [
60+
solid_density, fluid_density, fluid_viscosity, porosity, shear_modulus, drained_bulk_modulus, biot_coefficient, fluid_bulk_modulus, solid_bulk_modulus, isotropic_permeability,
61+
reference_stress_xx, reference_stress_yy, reference_stress_zz, reference_stress_xy,
62+
reference_strain_xx, reference_strain_yy, reference_strain_zz, reference_strain_xy
63+
]
64+
db_auxiliary_field.data = [
65+
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,
66+
-5.0*MPa, -5.0*MPa, -5.0*MPa, 0.0*Pa,
67+
0.0, 0.0, 0.0, 0.0
68+
]
69+
70+
auxiliary_subfields.porosity.basis_order = 1
71+
72+
[pylithapp.problem.materials.intrusion.bulk_rheology]
73+
use_reference_state = True
74+
```
75+
76+
## Running the simulation
77+
78+
```{code-block} console
79+
---
80+
caption: Run Step 2 simulation
81+
---
82+
$ pylith step02_inflation.cfg
83+
84+
# The output should look something like the following.
85+
86+
>> /Users/baagaard/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/apps/PyLithApp.py:84:main
87+
-- pylithapp(info)
88+
-- Running on 1 process(es).
89+
>> /Users/baagaard/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/meshio/MeshIOObj.py:43:read
90+
-- meshiocubit(info)
91+
-- Reading finite-element mesh
92+
>> /Users/baagaard/src/cig/pylith/libsrc/pylith/meshio/MeshIOCubit.cc:156:void pylith::meshio::MeshIOCubit::_readVertices(ExodusII &, scalar_array *, int *, int *) const
93+
-- meshiocubit(info)
94+
-- Component 'reader': Reading 747 vertices.
95+
96+
# -- many lines omitted --
97+
98+
>> /Users/baagaard/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/problems/TimeDependent.py:137:run
99+
-- timedependent(info)
100+
-- Solving problem.
101+
0 TS dt 1. time -1.
102+
0 SNES Function norm 7.521665654021e-01
103+
Linear solve converged due to CONVERGED_ATOL iterations 1
104+
1 SNES Function norm 5.516056137722e-02
105+
Linear solve converged due to CONVERGED_ATOL iterations 1
106+
2 SNES Function norm 1.166275103468e-02
107+
Linear solve converged due to CONVERGED_ATOL iterations 1
108+
3 SNES Function norm 3.829709526109e-03
109+
Linear solve converged due to CONVERGED_ATOL iterations 1
110+
4 SNES Function norm 1.461994548486e-03
111+
Linear solve converged due to CONVERGED_ATOL iterations 1
112+
5 SNES Function norm 6.066575419934e-04
113+
Linear solve converged due to CONVERGED_ATOL iterations 1
114+
6 SNES Function norm 2.655403615568e-04
115+
Linear solve converged due to CONVERGED_ATOL iterations 1
116+
7 SNES Function norm 1.202717080189e-04
117+
Linear solve converged due to CONVERGED_ATOL iterations 1
118+
8 SNES Function norm 5.567347859811e-05
119+
Linear solve converged due to CONVERGED_ATOL iterations 1
120+
9 SNES Function norm 2.613778728633e-05
121+
Linear solve converged due to CONVERGED_ATOL iterations 1
122+
10 SNES Function norm 1.238875375025e-05
123+
Linear solve converged due to CONVERGED_ATOL iterations 1
124+
11 SNES Function norm 5.911640181274e-06
125+
Linear solve converged due to CONVERGED_ATOL iterations 1
126+
12 SNES Function norm 2.834946322421e-06
127+
Linear solve converged due to CONVERGED_ATOL iterations 1
128+
13 SNES Function norm 1.364700102970e-06
129+
Linear solve converged due to CONVERGED_ATOL iterations 1
130+
14 SNES Function norm 6.589360471852e-07
131+
Linear solve converged due to CONVERGED_ATOL iterations 1
132+
15 SNES Function norm 3.189484959867e-07
133+
Linear solve converged due to CONVERGED_ATOL iterations 1
134+
16 SNES Function norm 1.547003895237e-07
135+
Linear solve converged due to CONVERGED_ATOL iterations 1
136+
17 SNES Function norm 7.516606485367e-08
137+
Linear solve converged due to CONVERGED_ATOL iterations 0
138+
18 SNES Function norm 3.657704494472e-08
139+
Linear solve converged due to CONVERGED_ATOL iterations 0
140+
19 SNES Function norm 1.782258437327e-08
141+
Linear solve converged due to CONVERGED_ATOL iterations 1
142+
20 SNES Function norm 8.694405655955e-09
143+
Linear solve converged due to CONVERGED_ATOL iterations 0
144+
21 SNES Function norm 4.245852416264e-09
145+
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 21
146+
1 TS dt 1. time 0.
147+
0 SNES Function norm 1.826079838384e+01
148+
Linear solve converged due to CONVERGED_ATOL iterations 31
149+
1 SNES Function norm 4.299962863814e-02
150+
Linear solve converged due to CONVERGED_ATOL iterations 1
151+
2 SNES Function norm 6.431697063491e-03
152+
Linear solve converged due to CONVERGED_ATOL iterations 1
153+
3 SNES Function norm 1.639359669638e-03
154+
Linear solve converged due to CONVERGED_ATOL iterations 1
155+
4 SNES Function norm 5.026774903429e-04
156+
Linear solve converged due to CONVERGED_ATOL iterations 1
157+
5 SNES Function norm 1.680648206351e-04
158+
Linear solve converged due to CONVERGED_ATOL iterations 1
159+
6 SNES Function norm 5.967136123324e-05
160+
Linear solve converged due to CONVERGED_ATOL iterations 1
161+
7 SNES Function norm 2.230420331383e-05
162+
Linear solve converged due to CONVERGED_ATOL iterations 1
163+
8 SNES Function norm 8.721830601750e-06
164+
Linear solve converged due to CONVERGED_ATOL iterations 1
165+
9 SNES Function norm 3.542878294948e-06
166+
Linear solve converged due to CONVERGED_ATOL iterations 1
167+
10 SNES Function norm 1.484027205445e-06
168+
Linear solve converged due to CONVERGED_ATOL iterations 1
169+
11 SNES Function norm 6.368036188603e-07
170+
Linear solve converged due to CONVERGED_ATOL iterations 1
171+
12 SNES Function norm 2.784493002713e-07
172+
Linear solve converged due to CONVERGED_ATOL iterations 1
173+
13 SNES Function norm 1.235724894148e-07
174+
Linear solve converged due to CONVERGED_ATOL iterations 1
175+
14 SNES Function norm 5.549400015107e-08
176+
Linear solve converged due to CONVERGED_ATOL iterations 0
177+
15 SNES Function norm 2.516326072475e-08
178+
Linear solve converged due to CONVERGED_ATOL iterations 0
179+
16 SNES Function norm 1.150182762543e-08
180+
Linear solve converged due to CONVERGED_ATOL iterations 1
181+
17 SNES Function norm 5.293129922422e-09
182+
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 17
183+
184+
# -- many lines omitted --
185+
186+
50 TS dt 1. time 49.
187+
0 SNES Function norm 1.583816680399e-01
188+
Linear solve converged due to CONVERGED_ATOL iterations 1
189+
1 SNES Function norm 3.297785796277e-05
190+
Linear solve converged due to CONVERGED_ATOL iterations 1
191+
2 SNES Function norm 1.625267288996e-05
192+
Linear solve converged due to CONVERGED_ATOL iterations 1
193+
3 SNES Function norm 8.011050419800e-06
194+
Linear solve converged due to CONVERGED_ATOL iterations 1
195+
4 SNES Function norm 3.949253118601e-06
196+
Linear solve converged due to CONVERGED_ATOL iterations 1
197+
5 SNES Function norm 1.947152025331e-06
198+
Linear solve converged due to CONVERGED_ATOL iterations 1
199+
6 SNES Function norm 9.601580478216e-07
200+
Linear solve converged due to CONVERGED_ATOL iterations 1
201+
7 SNES Function norm 4.735242133215e-07
202+
Linear solve converged due to CONVERGED_ATOL iterations 1
203+
8 SNES Function norm 2.335591654437e-07
204+
Linear solve converged due to CONVERGED_ATOL iterations 1
205+
9 SNES Function norm 1.152140892580e-07
206+
Linear solve converged due to CONVERGED_ATOL iterations 0
207+
10 SNES Function norm 5.684167635325e-08
208+
Linear solve converged due to CONVERGED_ATOL iterations 0
209+
11 SNES Function norm 2.804655474534e-08
210+
Linear solve converged due to CONVERGED_ATOL iterations 0
211+
12 SNES Function norm 1.384019346908e-08
212+
Linear solve converged due to CONVERGED_ATOL iterations 0
213+
13 SNES Function norm 6.830511841192e-09
214+
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 13
215+
51 TS dt 1. time 50.
216+
>> /Users/baagaard/software/unix/py310-venv/pylith-debug/lib/python3.10/site-packages/pylith/problems/Problem.py:204:finalize
217+
-- timedependent(info)
218+
-- Finalizing problem.
219+
220+
221+
```
222+
223+
In contrast with Step 1, the evolution of porosity in Step 2 results in nonlinear governing equations.
224+
This requires multiple iterations of the linear solve for the nonlinear solver to converge.
225+
226+
## Visualizing the results
227+
228+
In {numref}`fig:example:magma:2d:step02:solution` we use ParaView to visualize the evolution of the y displacement component using the `viz/plot_dispwarp.py` Python script.
229+
First, we start ParaView from the `examples/magma-2d` directory.
230+
231+
```{code-block} console
232+
---
233+
caption: Open ParaView using the command line.
234+
---
235+
$ PATH_TO_PARAVIEW/paraview
236+
237+
# For macOS, it will be something like
238+
$ /Applications/ParaView-5.10.1.app/Contents/MacOS/paraview
239+
```
240+
241+
Next, we override the default name of the simulation file with the name of the current simulation.
242+
243+
```{code-block} python
244+
---
245+
caption: Set the simulation in the ParaView Python Shell.
246+
---
247+
>>> SIM = "step02_inflation_statevars"
248+
```
249+
250+
Next we run the `viz/plot_dispwarp.py` Python script as described in {ref}`sec-paraview-python-scripts`.
251+
252+
:::{figure-md} fig:example:magma:2d:step02:solution
253+
<img src="figs/step02-solution.*" alt="Solution for Step 2 at t=100 yr. The colors indicate the fluid pressure, and the deformation is exaggerated by a factor of 1000." width="75%"/>
254+
255+
Solution for Step 2 at t=100 yr.
256+
The colors of the shaded surface indicate the fluid pressure, and the deformation is exaggerated by a factor of 1000.
257+
The reference state gives rise to greater vertical deformation compared to Step 1.
258+
The choice of material properties does not lead to significant changes in porosity in either material during the simulation (not shown in the figure).
259+
:::

‎docs/user/examples/magma-2d/step02_inflation_statevars-synopsis.md‎

Lines changed: 17 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,17 @@
1+
2+
**Features**
3+
4+
* Quasi-static problem
5+
* LU preconditioner
6+
* pylith.materials.Poroelasticity
7+
* pylith.meshio.MeshIOCubit
8+
* pylith.problems.TimeDependent
9+
* pylith.problems.SolnDispPresTracStrain
10+
* pylith.problems.InitialConditionDomain
11+
* pylith.bc.DirichletTimeDependent
12+
* pylith.bc.NeumannTimeDependent
13+
* pylith.meshio.DataWriterHDF5
14+
* spatialdata.spatialdb.SimpleGridDB
15+
* spatialdata.spatialdb.UniformDB
16+
* Poroelasticity with porosity state variable
17+
* Isotropic linear poroelasticity with reference state

‎docs/user/examples/reverse-2d/step08_twofaults_powerlaw-synopsis.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,3 +15,4 @@
1515
* pylith.faults.FaultCohesiveKin
1616
* pylith.faults.KinSrcStep
1717
* spatialdata.spatialdb.UniformDB
18+
* spatialdata.spatialdb.CompositeDB

‎docs/user/physics/materials/poroelasticity.md‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -42,12 +42,16 @@ L: isotropic, linear poroelasticity
4242

4343
```{table} Derived subfields that are available for output for poroelasticity bulk rheologies.
4444
:name: tab:poroelasticity:derived:subfields
45-
| Subfield | L | LM | GM | PL | DP | Components |
46-
|:--------------|:---:|:---:|:---:|:---:|:---:|:-----------------------|
47-
| `cauchy_stress` | ✓ | ✓ | ✓ | ✓ | ✓ | xx, yy, zz, xy, yz, xz |
48-
| `cauchy_strain` | ✓ | ✓ | ✓ | ✓ | ✓ | xx, yy, zz, xy, yz, xz |
45+
| Subfield | L | Components |
46+
|:--------------|:---:|:-----------------------|
47+
| `cauchy_stress` | X | xx, yy, zz, xy, yz, xz |
48+
| `cauchy_strain` | X | xx, yy, zz, xy, yz, xz |
49+
| `bulk_density` | X | |
50+
| `water_content` | X | |
4951
```
5052

53+
When porosity is enabled as a state variable, it will be included in the output along with the derived subfields.
54+
5155
:::{seealso}
5256
See [`Poroelasticity` Component](../../components/materials/Poroelasticity.md) for the Pyre properties and facilities and configuration examples.
5357
:::

0 commit comments

Comments
 (0)