tpv13/tpv13.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training example TPV 13-3D\n",
    "\n",
    "TPV 13-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 a 2D planar 60-degree dipping normal fault  \n",
    "* homogeneous half-space: $V_p=5716\\,m/s, V_s=3300\\,m/s, \\rho=2700\\,kg/m^3$\n",
    "* non-associative Drucker-Prager plasticity with yielding in shear \n",
    "* linear-slip weakening (LSW) friction \n",
    "* initial stress conditions are depth-dependent and specified everywhere throughout the entire model volume\n",
    "* prescribed nucleation zone, square, 3 x 3 km in size, locally lower static coefficient of friction\n",
    "* strongly supershear rupture conditions\n",
    "\n",
    "![](tpv13.jpg)\n",
    "\n",
    "[Detailed benchmark description (SCEC)](https://strike.scec.org/cvws/tpv12_13docs.html).\n",
    "\n",
    "Note, that for this training we slightly reduced two material constants, the bulk friction and the bulk cohesion, to increase the effects of off-fault yielding."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mesh generation\n",
    "\n",
    "We use the open-source software Gmsh ([http://gmsh.info](http://gmsh.info)) for mesh generation. The following steps are required:\n",
    "\n",
    "* create a half-space box-shaped domain (box)\n",
    "* create a 60° dipping rectangular fault plane\n",
    "* ensure that the mesh explicitly accounts for the nucleation patch at the hypocenter (why? see [Pelties et al., 2014](https://gmd.copernicus.org/articles/7/847/2014/) and [Galis et al., 2015](https://academic.oup.com/gji/article/200/2/890/609776))\n",
    "* set all required boundary conditions\n",
    "\n",
    "Gmsh comes with its own scripting language, and provides a [manual](http://gmsh.info/doc/texinfo/gmsh.html) and [tutorial examples](http://gmsh.info/doc/texinfo/gmsh.html#Tutorial). The first step is to build a geometrical (CAD) model, which is then used by the automatic mesh generator.\n",
    "We use the gmsh script tpv13_training.geo to generate the 3D unstructured tetrahedral mesh in the next step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gmsh -3 tpv13_training.geo\n",
    "# on Frontera with apptainer, replace with:\n",
    "# !mpirun apptainer run {\"~/my-training.sif\"} gmsh -3 tpv13_training.geo"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d67515aa",
   "metadata": {},
   "source": [
    "You should now have a file called `tpv13_training.msh`.\n",
    "In the next step, we translate this file into the efficient HDF5 format, which is going to be read in by SeisSol. HDF5 stands for the Hierarchical Data Format version 5, and is an open-source file format that supports large, complex, heterogeneous data. HDF5 uses a \"file directory\" like structure that allows you to organize data within the file in many different structured ways, as you might do with files on your computer, and is built for fast I/O processing and storage.\n",
    "\n",
    "We use our tool pumgen for mesh generation for SeisSol, which is [open-source software](  https://github.com/SeisSol/PUMGen), see also [SeisSol's documentation]( https://seissol.readthedocs.io/en/latest/meshing-with-pumgen.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mpirun -n 1 pumgen -s msh4 tpv13_training.msh  -setnumber dip 60\n",
    "# on Frontera with apptainer, replace with:\n",
    "# !mpirun -n 1 apptainer run {\"~/my-training.sif\"} pumgen -s msh4 tpv13_training.msh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The files `tpv13_training.puml.h5` and `tpv13_training.xdmf` were created.\n",
    "While the .h5 file is read by SeisSol, the .xdmf file can be used to visualize the mesh, as in the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vtk\n",
    "import pyvista as pv\n",
    "\n",
    "reader = vtk.vtkXdmfReader()\n",
    "reader.SetFileName('tpv13_training.xdmf')\n",
    "reader.Update()\n",
    "mesh = pv.wrap(reader.GetOutput())\n",
    "pv.plot(mesh, cmap='BuGn', background='white', show_edges=True, jupyter_backend='panel')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the mesh is refined close to the fault."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises\n",
    "* Inspect the `tpv13_training.geo` file and try to understand the meshing script. Consult the [Gmsh documentation](https://gmsh.info/doc/texinfo/gmsh.html) about the BooleanFragments operation (which may not be intuitive but is a powerful operation gmsh provides).\n",
    "* Create a finer resolved mesh by changing the `h_fault` parameter. Do not change the .geo file but use the [`-setnumber PARAMETER VALUE`](https://gmsh.info/doc/texinfo/gmsh.html#Command_002dline-options) command line argument of above *gmsh* command. Try, for example, to mesh the nucleation patch with 300 m and the remaining fault with 600 m element edge length.\n",
    "* Change the dip of the fault by changing the dip parameter."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run SeisSol for a 3D dynamic rupture earthquake scenario with off-fault plastic deformation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the following, we start SeisSol using the provided parameter file.\n",
    "You should adjust the OMP_NUM_THREADS variable to match the number of cores on your PC or laptop.\n",
    "Note, that on larger clusters SeisSol uses an optimized hybrid MPI + OpenMP parallelization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!OMP_NUM_THREADS=4 mpirun -n 1 SeisSol_Release_dhsw_4_elastic parameters.par\n",
    "# on Frontera with apptainer, replace with:\n",
    "# !SEISSOL_COMMTHREAD=0 OMP_NUM_THREADS=28 mpirun -n 2 apptainer run {\"~/my-training.sif\"} SeisSol_Release_dhsw_4_elastic parameters.par"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While SeisSol is running, you might want to checkout the [documentation of the input files](https://seissol.readthedocs.io/en/latest/parameter-file.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now visualize the fault output generated by SeisSol.\n",
    "Check out the [documentation](https://seissol.readthedocs.io/en/latest/fault-output.html#outputmask) for an explanation of which output field each abbreviation means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "reader = vtk.vtkXdmfReader()\n",
    "reader.SetFileName('output/tpv13-fault.xdmf')\n",
    "reader.Update()\n",
    "cd = reader.GetOutput().GetCellData()\n",
    "variables = [cd.GetArrayName(i) for i in range(cd.GetNumberOfArrays())]\n",
    "\n",
    "@interact(t=(0.0, 8.0, 1.0), var=variables)\n",
    "def plot(t=0.0, var='SRd'):\n",
    "    reader.UpdateTimeStep(t)\n",
    "    mesh = pv.wrap(reader.GetOutput())\n",
    "    plotter = pv.Plotter(notebook=True)\n",
    "    plotter.set_background('white')\n",
    "    plotter.add_mesh(mesh, cmap='Blues', scalars=var)\n",
    "    plotter.view_xz()\n",
    "    plotter.show(jupyter_backend='panel')\n",
    "    plotter.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5291f06b-f8c9-473f-9361-eddfa993b0c1",
   "metadata": {},
   "source": [
    "# Initial stress\n",
    "\n",
    "The loading in tpv13 is defined as a stress tensor, which is converted into fault tractions by SeisSol.\n",
    "In the next cell, we compute analytically the normal and shear tractions, and compute the relative prestress ratio, a key parameter for rupture dynamics, characterizing the initial stress.\n",
    "By modifying params with the nucleation static friction (mu_s=0.54), verify that this area is critically stressed (R>1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp\n",
    "import numpy as np\n",
    "\n",
    "# Define symbolic variables\n",
    "s3_to_s1, sigma_max, dip, mu_s, mu_d = sp.symbols(\"s3_to_s1 sigma_max dip mu_s mu_d\")\n",
    "\n",
    "# stress tensor\n",
    "stress = sp.diag((1 + s3_to_s1) * sigma_max / 2, s3_to_s1 * sigma_max, sigma_max)\n",
    "# normal and along dip unit vectors\n",
    "u_n = sp.Matrix([0, -sp.sin(dip), sp.cos(dip)])\n",
    "udip = sp.Matrix([0, sp.cos(dip), sp.sin(dip)])\n",
    "\n",
    "# tractions\n",
    "traction = stress * u_n\n",
    "\n",
    "sigma_n = traction.dot(u_n)\n",
    "print(\"sigma_n:\", sigma_n)\n",
    "sigma_d = traction.dot(udip)\n",
    "print(\"sigma_d:\", sigma_d)\n",
    "\n",
    "R = (sigma_d - mu_d * sigma_n) / ((mu_s - mu_d) * sigma_n)\n",
    "print(\"relative prestress ratio R:\", R)\n",
    "\n",
    "params = {\"dip\": 60 * np.pi / 180, \"s3_to_s1\": 0.3495, \"mu_s\": 0.7, \"mu_d\": 0.1}\n",
    "print(\n",
    "    f\"relative prestress ration R evaluated with {params}:\",\n",
    "    R.subs(params),\n",
    ")\n",
    "\n",
    "print(\"sigma_n:\", sigma_n.subs(params))\n",
    "print(\"sigma_d:\", sigma_d.subs(params))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de67024b-6854-4fa3-9fe4-13562d828801",
   "metadata": {},
   "source": [
    "# Rupture speed analysis\n",
    "\n",
    "We analyze the mechanism causing supershear in this setup and disable it.\n",
    "\n",
    "* Examine the earthquake rupture speed and identify the conditions that trigger supershear rupture.\n",
    "* Change the nucleation size from 3 to 1.2km in both `tpv12_13_fault.yaml` (update parameters nucleation_size) and `tpv13_training.geo` (update parameters `l_n` and `w_n`) to disable the current supershear triggering mechanism.\n",
    "* You should now observe a delay in the supershear transition.\n",
    "* Verify that you can suppress the supershear transition mechanism by decreasing the initial shear stress, that is changing `sigma_3_to_sigma_1` from 0.3495 to 0.4 in `tpv12_13_initial_stress.yaml` and reenabling nucleation by decreasing `mu_s` in the nucleation area from 0.54 to 0.43 and increasing nucleation size from 1.2 to 1.8 km (in `tpv12_13_fault.yaml` and `tpv13_training.geo`)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fault-dip dependence of rupture dynamics and ground deformation exercise\n",
    "Now we compare the subsidence and uplift of dynamic rupture earthquake scenarios varying the dip of the fault\n",
    "\n",
    "To do so, perform the following steps:  \n",
    "* Generate a new mesh with an alternative dip angle of your choice. In addition to the mesh, don't forget to update the dip value in `tpv12_13_fault.yaml`, and adapt `mu_s` for the nucleation to remain overstressed (use the python code in the section \"Initial stress\").\n",
    "* Run a new simulation using the same command as above but with its output saved in a different folder, and change the prefix that all files will have (`OutputFile`) in the `&Output` namelist \n",
    "* Add the new prefix to the `prefixList` of the code block below\n",
    "* Analyse differences using the visualisation widget above and the cross-sections of surface displacements below\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "\n",
    "prefixList = ['tpv13']\n",
    "nSim = len(prefixList)\n",
    "# coordinates of the line end points\n",
    "a = [0, -10e3, 0]\n",
    "b = [0, 10e3, 0]\n",
    "# create line polydata for vizualizing\n",
    "line = pv.Line(a, b)\n",
    "p = pv.Plotter(shape=(1, nSim), notebook=True, window_size=[1200, 500])\n",
    "\n",
    "def load_mesh(prefix):\n",
    "    # Read free-surface output data\n",
    "    reader = vtk.vtkHDFReader()\n",
    "    reader.SetFileName(f'output/{prefix}-free-surface-4.vtkhdf')\n",
    "    reader.Update()\n",
    "    return pv.wrap(reader.GetOutput())\n",
    "\n",
    "# plot the line of the elevation profile over the surface displacement \n",
    "for i, pref in enumerate(prefixList):\n",
    "    mesh = load_mesh(pref)\n",
    "    p.subplot(0, i)\n",
    "    p.add_text(pref, font_size=10, position=\"upper_left\")\n",
    "    p.set_background('lightgrey')\n",
    "    p.add_mesh(mesh, scalars='u3')\n",
    "    p.add_mesh(line, color=\"white\", line_width=10)\n",
    "    p.add_point_labels(\n",
    "        [a, b], [\"A\", \"B\"], font_size=20, point_color=\"red\", text_color=\"red\"\n",
    "    )\n",
    "    p.camera.zoom(3)\n",
    "p.show(jupyter_backend='panel')\n",
    "p.close()\n",
    "\n",
    "# plot fault output\n",
    "p = pv.Plotter(shape=(1, nSim), notebook=True, window_size=[1200, 350])\n",
    "for i, pref in enumerate(prefixList):\n",
    "    reader = vtk.vtkXdmfReader()\n",
    "    reader.SetFileName(f'output/{pref}-fault.xdmf')\n",
    "    reader.Update()\n",
    "    reader.UpdateTimeStep(9)\n",
    "    mesh = pv.wrap(reader.GetOutput())\n",
    "    p.subplot(0, i)\n",
    "    p.add_text(pref, font_size=10, position=\"upper_left\")\n",
    "    p.set_background('white')\n",
    "    p.add_mesh(mesh, cmap='plasma', scalars=\"ASl\")\n",
    "    p.view_xz()\n",
    "    p.camera.zoom(1.6)  \n",
    "p.show(jupyter_backend='panel')\n",
    "p.close()\n",
    "\n",
    "# plot subsidence/uplift profiles\n",
    "def plot_combined_profiles(prefixList, a, b, resolution=10000):\n",
    "    plt.figure(figsize=(6, 4))\n",
    "    for prefix in prefixList:\n",
    "        mesh = load_mesh(prefix)\n",
    "        sampled_data = mesh.sample_over_line(a, b, resolution=resolution)\n",
    "        distance = sampled_data['Distance']\n",
    "        displacement = sampled_data['u3']\n",
    "        plt.plot(distance/1000, displacement, label=prefix, linewidth=2)\n",
    "    plt.axhline(y=0,color='k', linestyle='--', linewidth=0.5)\n",
    "    plt.title(\"Comparison of subsidence/uplift profiles\")\n",
    "    plt.xlabel(\"Distance (km)\")\n",
    "    plt.ylabel(\"Vertical Displacement (m)\")\n",
    "    plt.legend()\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "plot_combined_profiles(prefixList, a, b)\n",
    "\n",
    "# plot moment rate release\n",
    "fig, ax = plt.subplots(figsize=(6, 4))\n",
    "for i, pref in enumerate(prefixList):\n",
    "    df = pd.read_csv(f'output/{pref}-energy.csv')\n",
    "    df = df.pivot_table(index=\"time\", columns=\"variable\", values=\"measurement\")\n",
    "    df[\"seismic_moment_rate\"] = np.gradient(df[\"seismic_moment\"], df.index[1])\n",
    "    ax.plot(df.index, df[\"seismic_moment_rate\"], label=pref)\n",
    "\n",
    "ax.set_xlabel('Time (s)')\n",
    "ax.set_ylabel('Seismic moment rate (Nm/s)')\n",
    "ax.set_title('Seismic moment rate comparison')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization with Paraview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also visualize the SeisSol output in Paraview on your local machine. For that, you have to copy the files from your Docker container. \n",
    "\n",
    "The following commands allow you to copy the simulation output from docker to the native filesystem of your machine. Type the commands in the terminal of your local machine: \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Execute on your local machine, not within this notebook!\n",
    "\n",
    "export my_docker=$(docker ps --format \"{{.Names}}\")\n",
    "\n",
    "docker cp ${my_docker}:/home/training/tpv13/output ."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To visualize the copied SeisSol output, open your Paraview and follow these steps: \n",
    "\n",
    "* Click on 1) in the upper left corner\n",
    "* Choose the xdmf file you want to open and click *ok*\n",
    "* Choose *XDMF Reader* and click the *ok* button\n",
    "* The output file will appear in the left bar below *Pipeline Browser*. \n",
    "Further down, in *Properties*, the *Apply-button* 2) should have turned green. Click on that button to read in the data.\n",
    "* With 3) you can choose the property you want to display. If you choose *ASL*, the fault slip will be shown.\n",
    "* With 4) you can switch to the last timestep. With the arrows you can switch back and forth, or you can choose a specific timestep by typing a *time*. \n",
    "* 5) adjusts the colorscale automatically to the data range. \n",
    "* With 6) you can change the properties of your color scale (i.e. color, range, log scale). \n",
    "\n",
    "![](paraview.jpg)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local-time stepping exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One of the key features of SeisSol is an efficient local-time stepping scheme, which clusters elements with similar time steps. This method pays off for applications featuring big differences in element sizes, caused by geometric complexity ([Ulrich et al., 2021](https://eartharxiv.org/repository/view/39/)) and/or multi-physics such as fully coupled seismic-acoustic-tsunami simulations ([Krenz et al., 2021](https://arxiv.org/abs/2107.06640)) or poroelastic wave propagation ([Wolf et al., 2021](https://arxiv.org/abs/2108.10565)).  \n",
    "\n",
    "A time-step cluster $c \\in \\mathbb{N}_0$ has the time-step\n",
    "$$\\Delta t_c = r^c \\Delta t_{\\text{min}},$$\n",
    "where the rate $r \\in \\mathbb{N}$. In the parameter file the rate $r$ can be set with the *ClusteredLTS* parameter.\n",
    "\n",
    "* Run SeisSol once with rate-2 local time-stepping (LTS) and once with global time-stepping (GTS).\n",
    "* In the log, look for \"Elapsed time (via clock_gettime)\" for both LTS and GTS.\n",
    "* Compute the speed-up and compare it to the theoretical speed-up due to LTS.\n",
    "\n",
    "*Hint:* You may reduce the *EndTime* parameter for this exercise."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}