Merge pull request #811 from baagaard-usgs/fix-eqinfo-docs

Improve documentation for pylith_eqinfo and illustrate its use in examples.

Author: baagaard-usgs

applications/pylith_eqinfo

@@ -28,7 +28,7 @@ NOTE: Works with HDF5 files, not VTK files.
 if __name__ == "__main__":
 
     from pylith.apps.EqInfoApp import EqInfoApp
-    from pyre.applications import start
+    from pythia.pyre.applications import start
     start(applicationClass=EqInfoApp)
 
 # End of file

docs/references.bib

@@ -82,6 +82,18 @@ @article{Brune:1970
   pages =	 {4997--5009},
 }
 
+@article{Wu:Chen:2003,
+  author = {Wu, Z.~L. and Chen, Y.~T.},
+  title = {Definition of seismic moment at a discontinuity interface},
+  journal = {Bulleting of the Seismological Society of America},
+  volume = {93},
+  number = {4},
+  year = 2003,
+  month = aug,
+  pages = {1832--1834},
+  doi = {10.1785/0120020234}
+}
+
 @article{Courant:etal:1967,
   author = {Courant, R. and Friedrichs, K. and Lewy, H.},
   title = {On the Partial Difference Equations of Mathematical Physics},

docs/user/examples/crustal-strikeslip-2d/step01-slip.md

@@ -142,6 +142,61 @@ At the end of the output written to the terminal, we see that the solver advance
 The linear solve converged after 19 iterations and the norm of the residual met the absolute convergence tolerance (`ksp_atol`) .
 The nonlinear solve converged in 1 iteration, which we expect because this is a linear problem, and the residual met the absolute convergence tolerance (`snes_atol`).
 
+### Earthquake rupture parameters
+
+We use the `pylith_eqinfo` utility to compute rupture information, such as earthquake magnitude, seismic moment, seismic potency, and average slip.
+For 2D simulations, the average slip provides the most useful information.
+[`pylith_eqinfo`](../../run-pylith/utilities.md) is a Pyre application and you specify parameters using `cfg` files and the command line.
+It writes results to a Python script for use in post-processing of simulation results.
+
+The file `eqinfoapp.cfg` holds the parameters for `pylith_eqinfo` in this example.
+By default, `pylith_eqinfo` extracts information for the final time step in the output file.
+In order to compute the seismic moment and moment magnitude, we need the shear modulus, which in this case is uniform over the domain.
+Consquently, we specify the density and shear wave speed using a `UniformDB` in `eqinfoapp.cfg`.
+
+```{code-block} console
+---
+caption: Run `pylith_eqinfo` for the Step 1 simulation.
+---
+$ pylith_eqinfo
+```
+
+```{code-block} python
+---
+caption: Contents of `output/step01_slip-eqinfo.py` generated by `pylith_eqinfo`. The `all` object includes earthquake rupture information combined from the three the individual faults. The average slip matches the uniform slip prescribed on each fault.
+---
+class RuptureStats(object):
+    pass
+all = RuptureStats()
+all.timestamp = [  3.155760e+07]
+all.ruparea = [  5.387992e+04]
+all.potency = [  1.733080e+05]
+all.moment = [  3.899430e+15]
+all.avgslip = [  3.216560e+00]
+all.mommag = [  4.360667e+00]
+main_fault = RuptureStats()
+main_fault.timestamp = [  3.155760e+07]
+main_fault.ruparea = [  3.577375e+04]
+main_fault.potency = [  1.430950e+05]
+main_fault.moment = [  3.219637e+15]
+main_fault.avgslip = [  4.000000e+00]
+main_fault.mommag = [  4.305205e+00]
+east_branch = RuptureStats()
+east_branch.timestamp = [  3.155760e+07]
+east_branch.ruparea = [  5.999357e+03]
+east_branch.potency = [  5.999357e+03]
+east_branch.moment = [  1.349855e+14]
+east_branch.avgslip = [  1.000000e+00]
+east_branch.mommag = [  3.386858e+00]
+west_branch = RuptureStats()
+west_branch.timestamp = [  3.155760e+07]
+west_branch.ruparea = [  1.210682e+04]
+west_branch.potency = [  2.421364e+04]
+west_branch.moment = [  5.448068e+14]
+west_branch.avgslip = [  2.000000e+00]
+west_branch.mommag = [  3.790828e+00]
+```
+
 ## Visualizing the results
 
 The `output` directory contains the simulation output.

docs/user/examples/crustal-strikeslip-3d/step01-slip.md

@@ -111,6 +111,58 @@ At the end of the output written to the terminal, we see that the solver advance
 The linear solve converged after 88 iterations and the norm of the residual met the absolute convergence tolerance (`ksp_atol`) .
 The nonlinear solve converged in 1 iteration, which we expect because this is a linear problem, and the residual met the absolute convergence tolerance (`snes_atol`).
 
+### Earthquake rupture parameters
+
+We use the `pylith_eqinfo` utility to compute rupture information, such as earthquake magnitude, seismic moment, seismic potency, and average slip.
+The file `eqinfoapp.cfg` holds the parameters for `pylith_eqinfo` in this example.
+By default, `pylith_eqinfo` extracts information for the final time step in the output file.
+In order to compute the seismic moment and moment magnitude, we need the shear modulus, which in this case is uniform over the domain.
+Consquently, we specify the density and shear wave speed using a `UniformDB` in `eqinfoapp.cfg`.
+The average slip, rupture area, seismic moment, and seismic potency are all given in SI units.
+
+```{code-block} console
+---
+caption: Run `pylith_eqinfo` for the Step 1 simulation.
+---
+$ pylith_eqinfo
+```
+
+```{code-block} python
+---
+caption: Contents of `output/step01_slip-eqinfo.py` generated by `pylith_eqinfo`. The `all` object includes earthquake rupture information combined from the three the individual faults. The average slip matches the uniform slip prescribed on each fault.
+---
+class RuptureStats(object):
+    pass
+all = RuptureStats()
+all.timestamp = [  3.155760e+07]
+all.ruparea = [  8.140700e+08]
+all.potency = [  2.623065e+09]
+all.moment = [  5.901897e+19]
+all.avgslip = [  3.222162e+00]
+all.mommag = [  7.147328e+00]
+main_fault = RuptureStats()
+main_fault.timestamp = [  3.155760e+07]
+main_fault.ruparea = [  5.424291e+08]
+main_fault.potency = [  2.169716e+09]
+main_fault.moment = [  4.881862e+19]
+main_fault.avgslip = [  4.000000e+00]
+main_fault.mommag = [  7.092390e+00]
+east_branch = RuptureStats()
+east_branch.timestamp = [  3.155760e+07]
+east_branch.ruparea = [  8.993288e+07]
+east_branch.potency = [  8.993288e+07]
+east_branch.moment = [  2.023490e+18]
+east_branch.avgslip = [  1.000000e+00]
+east_branch.mommag = [  6.170734e+00]
+west_branch = RuptureStats()
+west_branch.timestamp = [  3.155760e+07]
+west_branch.ruparea = [  1.817081e+08]
+west_branch.potency = [  3.634162e+08]
+west_branch.moment = [  8.176864e+18]
+west_branch.avgslip = [  2.000000e+00]
+west_branch.mommag = [  6.575058e+00]
+```
+
 ## Visualizing the results
 
 The `output` directory contains the simulation output.

docs/user/examples/strikeslip-2d/step04-varslip.md

@@ -151,6 +151,48 @@ At the end of the output written to the terminal, we see that the solver advance
 The linear solve converged after 27 iterations and the norm of the residual met the absolute convergence tolerance (`ksp_atol`).
 The nonlinear solve converged in 1 iteration, which we expect because this is a linear problem, and the residual met the absolute convergence tolerance (`snes_atol`).
 
+## Earthquake rupture parameters
+
+We use the `pylith_eqinfo` utility to compute rupture information, such as earthquake magnitude, seismic moment, seismic potency, and average slip.
+For 2D simulations, the average slip provides the most useful information.
+[`pylith_eqinfo`](../../run-pylith/utilities.md) is a Pyre application and you specify parameters using `cfg` files and the command line.
+It writes results to a Python script for use in post-processing of simulation results.
+
+The file `eqinfoapp.cfg` holds the parameters for `pylith_eqinfo` in this example.
+By default, `pylith_eqinfo` extracts information for the final time step in the output file.
+In order to compute the seismic moment and moment magnitude, we need the shear modulus for the fault.
+We account for the contrast in rigidity across the fault by using the effective shear modulus from {cite:t}`Wu:Chen:2003` and set the shear wave speed to 3.46 km/s.
+
+
+```{code-block} console
+---
+caption: Run `pylith_eqinfo` for the Step 4 simulation.
+---
+$ pylith_eqinfo
+```
+
+```{code-block} python
+---
+caption: Contents of `output/step04_varslip-eqinfo.py` generated by `pylith_eqinfo`. The `all` object includes earthquake rupture information combined from all of the individual faults.
+---
+class RuptureStats(object):
+    pass
+all = RuptureStats()
+all.timestamp = [  3.155760e+07]
+all.ruparea = [  6.300870e+04]
+all.potency = [  1.306180e+04]
+all.moment = [  3.909267e+14]
+all.avgslip = [  2.073016e-01]
+all.mommag = [  3.694730e+00]
+fault = RuptureStats()
+fault.timestamp = [  3.155760e+07]
+fault.ruparea = [  6.300870e+04]
+fault.potency = [  1.306180e+04]
+fault.moment = [  3.909267e+14]
+fault.avgslip = [  2.073016e-01]
+fault.mommag = [  3.694730e+00]
+```
+
 ## Visualizing the results
 
 In {numref}`fig:example:strikeslip:2d:step04:solution` we use the `pylith_viz` utility to visualize the y displacement field.

docs/user/run-pylith/utilities.md

@@ -337,19 +337,56 @@ pylith_dumpparameters [--quiet] [--format=ascii|json] [--filnemame=FILENAME]
 ## pylith_eqinfo
 
 This utility computes the moment magnitude, seismic moment, seismic potency, and average slip at user-specified time snapshots from PyLith fault HDF5 output.
-The utility works with output from simulations with either prescribed slip and/or spontaneous rupture.
-Currently, we compute the shear modulus from a user-specified spatial database at the centroid of the fault cells. In the future we plan to account for lateral variations in shear modulus across the fault when calculating the seismic moment.
-The Python script is a Pyre application, so its parameters can be specified using `cfg` and command line arguments just like PyLith.
+The utility works with output from simulations with either prescribed slip or spontaneous rupture.
+The moment magnitudes, seismic moment, and seismic potentency for 2D simulation have limited value, because they use on a 1D fault.
+Currently, we compute the shear modulus from a user-specified spatial database at the centroid of the fault cells.
+In the future we plan to automatically account for lateral variations in shear modulus across the fault when calculating the seismic moment.
+The average slip, rupture area, seismic moment, and seismic potency are all given in SI units.
+
+```{tip}
+From {cite:t}`Wu:Chen:2003` the seismic moment, $M_0$, when considering a laterial variation in the shear modulus across the fault is
+
+\begin{equation}
+M_0 = \mu_\mathit{eff} A D
+\end{equation}
+
+where $A$ is the rupture area, $D$ is the average slip, and $\mu_\mathit{eff}$ is the effective shear modulus given by
+
+\begin{equation}
+\mu_\mathit{eff} = \frac{1}{2} \left( \mu^+ + \mu^- \right) \left(1 - \left(\frac{\mu^+ - \mu^-}{\mu^+ + \mu^-} \right)^2 \right)
+\end{equation}
+
+where $\mu^+$ and $\mu^-$ are the shear modulus on each side of the fault.
+```
+
+The Python script is a Pyre application, so its parameters can be specified using `cfg` files (`eqinfoapp.cfg` will be read if it exists) or command line arguments just like PyLith.
 
 The Pyre properties and facilities include:
 
 :output_filename: Filename for output of slip information.
+:coordsys: Coordinate system associated with mesh in simulation.
 :faults: Array of fault names.
 :filename_pattern: Filename pattern in C/Python format for creating filename for each fault. Default is `output/fault_\%s.h5`.
 :snapshots: Array of timestamps for slip snapshosts ([-1] means use last time step in file, which is the default).
 :snapshot_units: Units for timestamps in array of snapshots.
 :db_properties: Spatial database for elastic properties.
-:coordsys: Coordinate system associated with mesh in simulation.
+```{code-block} cfg
+---
+caption: Example `cfg` file for running `pylith_eqinfo` in `examples/strikeslip-2d` to compute the earthquake magnitude, seismic moment, and average slip for step04_varslip.
+---
+[eqinfoapp]
+output_filename = output/step04_varslip-eqinfo.py
+
+coordsys.space_dim = 2
+
+faults = [fault]
+filename_pattern = output/step04_varslip-%s.h5
+
+db_properties = spatialdata.spatialdb.UniformDB
+db_properties.description = Fault properties
+db_properties.values = [density, Vs]
+db_properties.data = [2500.0*kg/m**3, 3.46*km/s]
+```
 
 (sec-user-run-pylith-pylith-genxdmf)=
 ## pylith_genxdmf

examples/crustal-strikeslip-2d/Makefile.am

@@ -17,6 +17,7 @@ dist_noinst_DATA = \
 	mesh_tri.msh \
 	mesh_tri.exo \
 	pylithapp.cfg \
+    eqinfoapp.cfg \
 	step01_slip.cfg \
 	step01_slip_cubit.cfg \
 	step02_varslip.cfg \

examples/crustal-strikeslip-2d/eqinfoapp.cfg

@@ -0,0 +1,12 @@
+[eqinfoapp]
+output_filename = output/step01_slip-eqinfo.py
+
+coordsys.space_dim = 3
+
+faults = [main_fault, east_branch, west_branch]
+filename_pattern = output/step01_slip-%s.h5
+
+db_properties = spatialdata.spatialdb.UniformDB
+db_properties.description = Fault properties
+db_properties.values = [density, Vs]
+db_properties.data = [2500.0*kg/m**3, 3.00*km/s]

examples/crustal-strikeslip-3d/Makefile.am

@@ -16,6 +16,7 @@ dist_noinst_DATA = \
 	README.md \
 	mesh_tet.msh \
 	pylithapp.cfg \
+    eqinfoapp.cfg \
 	step01_slip.cfg \
 	step01_slip_cubit.cfg \
 	step02_varslip.cfg \

examples/crustal-strikeslip-3d/eqinfoapp.cfg

@@ -0,0 +1,12 @@
+[eqinfoapp]
+output_filename = output/step01_slip-eqinfo.py
+
+coordsys.space_dim = 3
+
+faults = [main_fault, east_branch, west_branch]
+filename_pattern = output/step01_slip-%s.h5
+
+db_properties = spatialdata.spatialdb.UniformDB
+db_properties.description = Fault properties
+db_properties.values = [density, Vs]
+db_properties.data = [2500.0*kg/m**3, 3.00*km/s]