From e9e3fedddc94df1b81ca87ca396f3b318e3f458b Mon Sep 17 00:00:00 2001 From: aangniu Date: Thu, 16 Jan 2025 15:10:43 +0100 Subject: [PATCH] update notebook for cdb_tpv23 --- quakeworx/cdb_tpv23/cdb_tpv23.ipynb | 309 ++++++++++++++++++++++++---- quakeworx/tpv13/tpv13_qwx.ipynb | 6 +- 2 files changed, 267 insertions(+), 48 deletions(-) diff --git a/quakeworx/cdb_tpv23/cdb_tpv23.ipynb b/quakeworx/cdb_tpv23/cdb_tpv23.ipynb index 957e56a..5ebc461 100644 --- a/quakeworx/cdb_tpv23/cdb_tpv23.ipynb +++ b/quakeworx/cdb_tpv23/cdb_tpv23.ipynb @@ -2,11 +2,49 @@ "cells": [ { "cell_type": "markdown", - "id": "c8cc6a9b", + "id": "47d1a1c0-1cff-45d2-9a0d-d90c329c5d61", + "metadata": {}, + "source": [ + "Contact: Zihua Niu (zihua.niu@lmu.de)\n", + "\n", + "I this example, we will learn about dynamic rupture simulation with the continuum damage breakage model implemented in the off-fault material. We will explore together what is the influence of off-fault damage on the dynamic rupture process." + ] + }, + { + "cell_type": "markdown", + "id": "1e0eb4f7-da4b-42b2-a61b-4b809e15de15", + "metadata": {}, + "source": [ + "## Learning objectives\n", + "- Brief overview of the new parameters in the continuum damage breakage model\n", + "- Visualize the effect of off-fault damage on earthquake interaction.\n", + "- Understand what causes the difference in cases with and without off-fault damage." + ] + }, + { + "cell_type": "markdown", + "id": "e16090c3-13d1-46da-962d-7d48a71d129f", "metadata": {}, "source": [ "# Training example TPV 23-3D\n", "\n", + "We will use the TPV 23-3D configuration to achieve the above objectives.\n", + "\n", + "### Contents\n", + "\n", + "- [Model setup](#Model-setup)\n", + " - [Bulk Material Parameters of CDBM](#Bulk-Material-Parameters-of-CDBM)\n", + " - [Localized off-fault damage zone](#Localized-off-fault-damage-zone)\n", + "- [How the off-fault damage impacts on earthquake interaction](How-the-off-fault-damage-impacts-on-earthquake-interaction)" + ] + }, + { + "cell_type": "markdown", + "id": "c8cc6a9b", + "metadata": {}, + "source": [ + "## Model setup\n", + "\n", "TPV 23-3D is a benchmark exercise that is designed to test if computer codes that simulate dynamic earthquake rupture are working as intended [(Harris et al., SRL 2018)](https://pubs.geoscienceworld.org/ssa/srl/article/89/3/1146/530061/A-Suite-of-Exercises-for-Verifying-Dynamic). It was designed by the [SCEC/USGS Spontaneous Rupture Code Verification Project](https://strike.scec.org/cvws/) and features:\n", "\n", "* spontaneous rupture on two mutually parallel 2D vertical planar faults \n", @@ -24,7 +62,9 @@ "id": "e6367031", "metadata": {}, "source": [ - "### Bulk Material Parameters of CDBM [(Lyakhovsky et al., 2016)](https://academic.oup.com/gji/article-abstract/206/2/1126/2606007)\n", + "### Bulk Material Parameters of CDBM \n", + "\n", + "CDB model is proposed within the framework of the continuum mechanics. The mechanical response of rocks from their intact states to their failure is mathematically described with a scalar damage variable ($\\alpha$) that represents the density of distributed micro-cracks and a scalar breakage variable ($B$) describing the grain size distribution in the post-failure stage of rocks. In CDB, the latter stage is named as the granular phase. Both $\\alpha$ and $B$ are defined in the range of [0,1] [(Lyakhovsky et al., 2016)](https://academic.oup.com/gji/article-abstract/206/2/1126/2606007).\n", "\n", "| Parameters (Col 1) | Values (Col 1) | Units (Col 1) | Parameters (Col 2) | Values (Col 2) | Units (Col 2) |\n", "|--------------------|----------------|--------------------|--------------------|----------------|--------------------|\n", @@ -43,55 +83,185 @@ "source": [ "### Localized off-fault damage zone\n", "\n", - "We can see localized damage structures that extend from the fault plane into the bulk material at an angle of around 35.6 degrees. The theory of CDB describes the internal friction angle of the material in analogy to \\citet{byerlee1978friction}. The equivalent internal friction angle in CDB is determined from the nonlinear material properties $\\xi_0$ in the parameter table above. As listed in the above table, $\\xi_0$ is -0.75 \\citet{lyakhovsky1997distributed}. This corresponds to an internal friction angle of $43$ degrees. In the simulation of the following figure, the angle between the maximum compressive principal stress and the fault plane is 59.1 degrees. The two conjugate weak planes should take an angle of 45 - 43/2 = 23.5 degrees from the maximum compressive principal stress. This means the two conjugate weak planes should take angles of 59.1 - 23.5 = 35.6 degrees or 59.1 + 23.5 = 82.6 degrees. We find this theoretical angle of 35.6 degrees from the fault plane in the following figure (red dashed line) agrees with our numerical simulation results.\n", + "We can see localized damage structures that extend from the fault plane into the bulk material at an angle of around 35.6 degrees. The theory of CDB describes the internal friction angle of the material in analogy to [Byerlee (1978)](https://link.springer.com/book/10.1007/978-3-0348-7182-2). The equivalent internal friction angle in CDB is determined from the nonlinear material properties $\\xi_0$ in the parameter table above. As listed in the above table, $\\xi_0$ is -0.75 [(Lyakhovsky et al., 1997)](https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/97jb01896). This corresponds to an internal friction angle of $43$ degrees. From the modulus reduction in following figure, the angle between the maximum compressive principal stress and the fault plane is 59.1 degrees. The two conjugate weak planes should take an angle of 45 - 43/2 = 23.5 degrees from the maximum compressive principal stress. This means angles of 59.1 - 23.5 = 35.6 degrees or 59.1 + 23.5 = 82.6 degrees between the weak plane and the fault plane. We find this theoretical angle of 35.6 degrees from the fault plane in the following figure (red dashed line) agrees with our numerical simulation results.\n", "\n", "More detailed results are shown in our paper that is going to be submitted soon (Niu et al., 2025, in preparation).\n", "\n", "![](off_fault_damage.png)" ] }, + { + "cell_type": "markdown", + "id": "8125703a-3bde-4532-8372-7702cf242ee7", + "metadata": {}, + "source": [ + "Now we know the distribution of the off-fault damage structure around a single fault plane. Let's start to check how these structures influence earthquake interaction using our simulation on the Gateway." + ] + }, + { + "cell_type": "markdown", + "id": "5b16850b", + "metadata": {}, + "source": [ + "### How the off-fault damage impacts on earthquake interaction\n", + "\n", + "In this section, we will see how the off-fault damage impacts earthquake interaction, using the TPV23 configuration. \n", + "\n", + "In the linear (visco-)elastic case, the second fault will not be ruptured. It is ruptured because as the damage extends from the first fault, the stress field will be perturbed." + ] + }, + { + "cell_type": "markdown", + "id": "9c1ccdeb-8036-42b3-b919-6be58215d5c3", + "metadata": {}, + "source": [ + "We start from finding the location of our output files. " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "e48ebace-664e-4547-82f5-c75a3dc29340", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "/workspace/quakeworx/cdb_tpv23\n" + ] + } + ], + "source": [ + "!pwd # ! in the front keeps the command running in the terminal, instead of the notebook" + ] + }, + { + "cell_type": "markdown", + "id": "a2881c6a-e640-475d-882d-b1901a917fb9", + "metadata": {}, + "source": [ + "You will see you are at `/opt/notebooks/Training/quakeworx/cdb_tpv23`.\n", + "\n", + "`/opt/notebooks` is the system directory where all your jobs are stored and `Training/quakeworx/cdb_tpv23` is the structure of the Github repository you just cloned.\n", + "\n", + "You can check all the job names `{YOUR_JOB_NAME}` you have submitted:" + ] + }, { "cell_type": "code", "execution_count": null, - "id": "e989d406", + "id": "8bc36dab-a34a-4328-a55b-81455cacc7d9", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "!ls /opt/notebooks # list all job names" + ] + }, + { + "cell_type": "markdown", + "id": "a0816202-d21f-442d-9c04-443d8bade482", + "metadata": {}, + "source": [ + "You can then use this information to define the path to your solutions" + ] }, { "cell_type": "code", "execution_count": null, - "id": "78a1f71a", + "id": "999f4bab-2ace-4010-8ab3-89d22d41626e", "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# Replace {YOUR_JOB_NAME} with the job name you used to run your uniform shear modulus Tandem model\n", + "# result_dir = '/opt/notebooks/{YOUR_JOB_NAME}/'\n", + "result_dir = '/opt/notebooks/SeisSol_cdb_tpv23/'" + ] }, { "cell_type": "markdown", - "id": "5b16850b", + "id": "ef2f4b55-96ce-4562-9ae3-96b5b4738367", "metadata": {}, "source": [ - "### How the off-fault damage impacts on earthquake interaction?\n", - "\n", - "In this section, we will see how the off-fault damage impacts earthquake interaction, using the TPV23 configuration. \n", - "\n", - "In the linear (visco-)elastic case, the second fault will not be ruptured. It is ruptured because as the damage extends from the first fault, the stress field will be perturbed." + "First, load the necessary python packages that have been pre-installed." ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "id": "08b67ebc", - "metadata": {}, + "metadata": { + "scrolled": true + }, "outputs": [], "source": [ + "# For common data processing\n", "import numpy as np\n", - "from matplotlib import pyplot as plt\n", "\n", - "import seissolxdmf as seisx\n", - "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", + "# For visualization in 2D and 3D SeisSol mesh\n", + "from matplotlib import pyplot as plt\n", "from mpl_toolkits.mplot3d.art3d import Poly3DCollection\n", - "import matplotlib.colors as mcolors" + "import matplotlib.colors as mcolors\n", + "\n", + "# Dedicated library for post-processing HDF5 format SeisSol output\n", + "import seissolxdmf as seisx" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d96853a8-3cf0-4efa-a2b6-a9a495f335b2", + "metadata": {}, + "outputs": [], + "source": [ + "# Load data\n", + "xdmfFilename = \"output/cdb_1401CFL1p0_tpv23-fault.xdmf\"\n", + "sx = seisx.seissolxdmf(xdmfFilename)\n", + "ndt = sx.ReadNdt() - 1 #sx.ReadNdt() provides the total number of time steps that have been written out\n", + "xyz = sx.ReadGeometry()\n", + "connect = sx.ReadConnect()" + ] + }, + { + "cell_type": "markdown", + "id": "2ce81996-bebe-4965-a0a4-81f86d3c66f2", + "metadata": {}, + "source": [ + "Then you can define the time steps that you want to visualize. We first take a look at the fault output `6 seconds` after the nucleation of the rupture on the first fault plane. \n", + "\n", + "This is the time step when the **fault 1** has just ruptured. The **stress waves** reach the **fault 2**. But the damage zone just starts to accumulate from **fault 1** and does not reach **fault 2** yet." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2abf9fd-9a70-40c6-b140-afc62f9fe1d8", + "metadata": {}, + "outputs": [], + "source": [ + "ndt = 6\n", + "## Here you can load the physical field that you want to visualize. We list some of the variables that are often needed.\n", + "## More can be found here: https://seissol.readthedocs.io/en/latest/fault-output.html\n", + "## SRs and SRd: Slip rates in along-strike and along-dip directions\n", + "## T_s, T_d: Shear stress in strike and dip directions, P_n: Normal stress\n", + "## u_n: In the CDB model, u_n is used to store value of the damage variable $\\alpha$\n", + "sRate = sx.ReadData('SRs', ndt)\n", + "damage = sx.ReadData('u_n', ndt)\n", + "\n", + "# Define the colormap ranges\n", + "sRate_vmin, sRate_vmax = np.min(sRate), np.max(sRate)\n", + "sRate_norm = mcolors.Normalize(vmin=sRate_vmin, vmax=sRate_vmax)\n", + "## We manually set the damage range to keep it the same for all time steps, easier for comparison\n", + "damage_vmin, damage_vmax = 0, 0.25\n", + "damage_norm = mcolors.Normalize(vmin=damage_vmin, vmax=damage_vmax)" + ] + }, + { + "cell_type": "markdown", + "id": "b550daa8-880d-4611-a88e-3640fa541d33", + "metadata": {}, + "source": [ + "We here try to convert the data structure from SeisSol to those can be used for visualization with `matplotlib`." ] }, { @@ -112,23 +282,6 @@ } ], "source": [ - "# Load data\n", - "xdmfFilename = \"output_cdb/cdb_1401CFL1p0_tpv23-fault.xdmf\"\n", - "sx = seisx.seissolxdmf(xdmfFilename)\n", - "ndt = sx.ReadNdt() - 1\n", - "xyz = sx.ReadGeometry()\n", - "connect = sx.ReadConnect()\n", - "ndt = 6\n", - "sRate = sx.ReadData('SRs', ndt)\n", - "damage = sx.ReadData('u_n', ndt)\n", - "\n", - "# Define the colormap ranges\n", - "sRate_vmin, sRate_vmax = np.min(sRate), np.max(sRate)\n", - "sRate_norm = mcolors.Normalize(vmin=sRate_vmin, vmax=sRate_vmax)\n", - "\n", - "damage_vmin, damage_vmax = 0, 0.25\n", - "damage_norm = mcolors.Normalize(vmin=damage_vmin, vmax=damage_vmax)\n", - "\n", "# Extract vertices for each triangular face in the mesh\n", "triangles = [xyz[connect[i]] for i in range(connect.shape[0])]\n", "\n", @@ -180,6 +333,18 @@ "plt.show()\n" ] }, + { + "cell_type": "markdown", + "id": "f0ff1bec-c56d-4f8d-8a78-0e58106e9418", + "metadata": {}, + "source": [ + "The left pannel above shows the distribution of slip rate at 6 s. It shows the rupture front (high slip rate regions in darker red) has already reached the edge of the **fault 1** (the fault plane at the back). \n", + "\n", + "At the same time, the **off-fault damage** is still close to zero on **fault 2** (the fault plane at the front).\n", + "\n", + "It takes some longer time for the damage zone to extend to **fault 2**. At the end of the simulation (14 s), we can see the rupture also spontaneously grow on **fault 2**. We show the results in the follow cell." + ] + }, { "cell_type": "code", "execution_count": 28, @@ -260,6 +425,38 @@ "plt.show()" ] }, + { + "cell_type": "markdown", + "id": "4cadbd1a-3a5e-4871-bc22-c959dd8dba31", + "metadata": {}, + "source": [ + "From the left panel, we see the rupture front on **fault 2** now. \n", + "\n", + "This is a result of the damage zone reaching **fault 2**. We can tell the distribution of the **damage variable $\\alpha$** on the **fault 2** from the right panel.\n", + "\n", + "The solution at the nucleation region is not very smooth. This is a result of the low-resolution model that we run in the example to reduce the amount of time that we have to wait." + ] + }, + { + "cell_type": "markdown", + "id": "355b0d44-4010-4cc9-9b28-28a4c6b4abe1", + "metadata": {}, + "source": [ + "### Shear traction and strength evolution on the fault\n", + "\n", + "You may now wonder why **fault 2** can be ruptured with the accumulation of damage. In the following part, we will show how the shear traction and shear strength (frictional strength) on **fault 2** evolves in time.\n", + "\n", + "For this purpose, we need point-wise output on the fault fr a clearer quantitative visualization." + ] + }, + { + "cell_type": "markdown", + "id": "13e876d3-b4ba-4757-9040-a550bca8a2fa", + "metadata": {}, + "source": [ + "To read from the fault receiver files, we will need `pandas` to systematically read in different fields of data." + ] + }, { "cell_type": "code", "execution_count": 1, @@ -267,9 +464,11 @@ "metadata": {}, "outputs": [], "source": [ + "## Additional standard python library to access file names inside a folder\n", "import os\n", - "import pandas as pd\n", - "import re" + "import re\n", + "\n", + "import pandas as pd" ] }, { @@ -279,13 +478,14 @@ "metadata": {}, "outputs": [], "source": [ + "# Function to find the complete name of a receiver file for a given receiver {number}\n", "def find_receiver(directory, prefix, number):\n", " receiver_re = re.compile(f\"{prefix}-faultreceiver-{number:05d}-(\\\\d)+.dat\")\n", " for fn in os.listdir(directory):\n", " if receiver_re.match(fn):\n", " print(os.path.join(directory, fn))\n", " return os.path.join(directory, fn)\n", - "\n", + "# Function to read in time series with pandas once we know the specific receiver {file name}\n", "def read_receiver(filename):\n", " with open(filename) as receiver_file:\n", " receiver_file.readline()\n", @@ -300,6 +500,18 @@ " return receiver_df" ] }, + { + "cell_type": "markdown", + "id": "7e6f6b9d-5f0d-4624-b189-bda9dd6f8b4e", + "metadata": {}, + "source": [ + "Now let's plot the time series for shear traction and shear strength at a point on **fault 2** that locates inside the nucleation area.\n", + "\n", + "The point receiver number is **9**, located at $x$ = 5 km, $y$ = -3 km and $z$ = 5.5 km. \n", + "\n", + "The location of th fault can be found inside the input file `tpv23_3km_faultreceivers`." + ] + }, { "cell_type": "code", "execution_count": 4, @@ -381,12 +593,19 @@ ] }, { - "cell_type": "code", - "execution_count": null, - "id": "8895c950", + "cell_type": "markdown", + "id": "1d381632-6c65-4c3d-9fe5-4d23014da5f1", "metadata": {}, - "outputs": [], - "source": [] + "source": [ + "From the figure above, we can define **4 phases** in the evolution of **fault 2**\n", + "\n", + "- Phase 1: From 0 to 5 s, the stress wave from the rupturing of **fault 1** does not yet reach **fault 2**.\n", + "- Phase 2: From 5 to 8 s, the stress wave reaches **fault 2**, but the dynamic stress field is not contributing to a larger *shear traction* (the solid curve) that can exceed *shear strength* (the dashed curve).\n", + "- Phase 3: From 8 to 9.5 s, as the damage zone reaches **fault 2** and the *shear traction* starts to increase, getting closer to *shear strength*.\n", + "- Phase 4: After 9.5 s, **fault 2** rupture spontaneously and the *shear stress* drops again.\n", + "\n", + "More details will be soon available in our publication that is soon to be submitted. Please stay tuned :)" + ] }, { "cell_type": "code", @@ -413,7 +632,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.7" + "version": "3.12.7" } }, "nbformat": 4, diff --git a/quakeworx/tpv13/tpv13_qwx.ipynb b/quakeworx/tpv13/tpv13_qwx.ipynb index db37f52..9f40412 100644 --- a/quakeworx/tpv13/tpv13_qwx.ipynb +++ b/quakeworx/tpv13/tpv13_qwx.ipynb @@ -623,9 +623,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python (zihua_venv)", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "zihua_venv" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -637,7 +637,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.2" + "version": "3.9.7" } }, "nbformat": 4,