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"
   ]
  },
  {
   "cell_type": "markdown",
   "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 msh2 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*](https://gmsh.info/doc/texinfo/gmsh.html#Command_002dline-options) PARAMETER VALUE 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"
   ]
  },
  {
   "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')"
   ]
  },
  {
   "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."
   ]
  },
  {
   "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.\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": [
    "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)\n",
    "\n",
    "def load_mesh(prefix):\n",
    "    # Read xdmf data\n",
    "    reader = vtk.vtkXdmfReader()\n",
    "    reader.SetFileName(f'output/{prefix}-surface.xdmf')\n",
    "    reader.Update()\n",
    "    reader.UpdateTimeStep(5.)\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",
    "\n",
    "    p.subplot(0, i)\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",
    "    \n",
    "# plot subsidence/uplift profile\n",
    "for i, pref in enumerate(prefixList):\n",
    "    load_mesh(pref).plot_over_line(\n",
    "        a,\n",
    "        b,\n",
    "        resolution=10000,\n",
    "        title=f\"subsidence/uplift profile for {pref}\",\n",
    "        ylabel=\"vectical displacement (m)\",\n",
    "        scalars='u3'\n",
    "    )"
   ]
  },
 
  {
   "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",
    "export my_docker=$(docker ps --format \"{{.Names}}\")",
    "\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"
    ]
  }
 
 ],
 "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
}