From 2780a619e17af5a13d9de2abb07b92c206891451 Mon Sep 17 00:00:00 2001 From: Brad Aagaard Date: Tue, 11 Feb 2025 14:46:06 -0700 Subject: [PATCH 1/3] FIX: Fix pythia import in pylith_eqinfo. --- applications/pylith_eqinfo | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/applications/pylith_eqinfo b/applications/pylith_eqinfo index 03737322e8..422f0fdc57 100755 --- a/applications/pylith_eqinfo +++ b/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 From 209f9d13bc58fd81e2b324dbca1731fc07501708 Mon Sep 17 00:00:00 2001 From: Brad Aagaard Date: Tue, 11 Feb 2025 14:45:43 -0700 Subject: [PATCH 2/3] DOCS: Improve documentation for pylith_eqinfo and illustrate use in examples/strikeslip-2d and examples/crustal-strikeslip-*. --- docs/references.bib | 12 ++++++ docs/user/run-pylith/utilities.md | 45 ++++++++++++++++++-- examples/crustal-strikeslip-2d/Makefile.am | 1 + examples/crustal-strikeslip-2d/eqinfoapp.cfg | 12 ++++++ examples/crustal-strikeslip-3d/Makefile.am | 1 + examples/crustal-strikeslip-3d/eqinfoapp.cfg | 12 ++++++ examples/strikeslip-2d/Makefile.am | 1 + examples/strikeslip-2d/eqinfoapp.cfg | 12 ++++++ 8 files changed, 92 insertions(+), 4 deletions(-) create mode 100644 examples/crustal-strikeslip-2d/eqinfoapp.cfg create mode 100644 examples/crustal-strikeslip-3d/eqinfoapp.cfg create mode 100644 examples/strikeslip-2d/eqinfoapp.cfg diff --git a/docs/references.bib b/docs/references.bib index 93a6dc5981..2f6f369e55 100644 --- a/docs/references.bib +++ b/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}, diff --git a/docs/user/run-pylith/utilities.md b/docs/user/run-pylith/utilities.md index 02c8296de8..a3b658b664 100644 --- a/docs/user/run-pylith/utilities.md +++ b/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 diff --git a/examples/crustal-strikeslip-2d/Makefile.am b/examples/crustal-strikeslip-2d/Makefile.am index 4472c12c21..880cc4c56b 100644 --- a/examples/crustal-strikeslip-2d/Makefile.am +++ b/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 \ diff --git a/examples/crustal-strikeslip-2d/eqinfoapp.cfg b/examples/crustal-strikeslip-2d/eqinfoapp.cfg new file mode 100644 index 0000000000..5c3811c0f6 --- /dev/null +++ b/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] diff --git a/examples/crustal-strikeslip-3d/Makefile.am b/examples/crustal-strikeslip-3d/Makefile.am index 52725d678f..ca1cab5d69 100644 --- a/examples/crustal-strikeslip-3d/Makefile.am +++ b/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 \ diff --git a/examples/crustal-strikeslip-3d/eqinfoapp.cfg b/examples/crustal-strikeslip-3d/eqinfoapp.cfg new file mode 100644 index 0000000000..5c3811c0f6 --- /dev/null +++ b/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] diff --git a/examples/strikeslip-2d/Makefile.am b/examples/strikeslip-2d/Makefile.am index b62d9e970f..141662a99b 100644 --- a/examples/strikeslip-2d/Makefile.am +++ b/examples/strikeslip-2d/Makefile.am @@ -25,6 +25,7 @@ dist_noinst_DATA = \ mesh_tri.exo \ mesh_quad.exo \ pylithapp.cfg \ + eqinfoapp.cfg \ step01a_slip.cfg \ step01b_slip.cfg \ step01c_slip.cfg \ diff --git a/examples/strikeslip-2d/eqinfoapp.cfg b/examples/strikeslip-2d/eqinfoapp.cfg new file mode 100644 index 0000000000..638e5b11fd --- /dev/null +++ b/examples/strikeslip-2d/eqinfoapp.cfg @@ -0,0 +1,12 @@ +[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] From 67d3d4ba69a3eec6eced0bd7135a72ef5595e593 Mon Sep 17 00:00:00 2001 From: Brad Aagaard Date: Tue, 11 Feb 2025 16:13:50 -0700 Subject: [PATCH 3/3] DOCS: Demonstrate use of pylith_eqinfo in examples/strikeslip-2d and examples/crustal-strikeslip-*. --- .../crustal-strikeslip-2d/step01-slip.md | 55 +++++++++++++++++++ .../crustal-strikeslip-3d/step01-slip.md | 52 ++++++++++++++++++ .../examples/strikeslip-2d/step04-varslip.md | 42 ++++++++++++++ 3 files changed, 149 insertions(+) diff --git a/docs/user/examples/crustal-strikeslip-2d/step01-slip.md b/docs/user/examples/crustal-strikeslip-2d/step01-slip.md index c2dbd6da97..c4318b3d3f 100644 --- a/docs/user/examples/crustal-strikeslip-2d/step01-slip.md +++ b/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. diff --git a/docs/user/examples/crustal-strikeslip-3d/step01-slip.md b/docs/user/examples/crustal-strikeslip-3d/step01-slip.md index 38e8f96373..b1593b1b29 100644 --- a/docs/user/examples/crustal-strikeslip-3d/step01-slip.md +++ b/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. diff --git a/docs/user/examples/strikeslip-2d/step04-varslip.md b/docs/user/examples/strikeslip-2d/step04-varslip.md index 1119f0cb34..c892b22317 100644 --- a/docs/user/examples/strikeslip-2d/step04-varslip.md +++ b/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.