From 1b2bad0816bb8290b73d7af8d098149be6173416 Mon Sep 17 00:00:00 2001 From: Thomas-Ulrich Date: Mon, 11 Nov 2024 21:53:57 +0100 Subject: [PATCH] use seissol v1.3/Rusanov fluxes/ vtkhdf/ extend tpv13 (#52) * use rusanov * fix Dockerfile_jupyterlab * small improvements on tpv13_training.geo * update to msh4 * switch to main * update both dockers * use latest releases of SeisSol and PUMGen * use vtk output for fault * add rupture speed exercise * combine subsidence plots and fix typo in geo * close plotter to avoid memory leak * precise exercise and input file * disable vtk fault output for now * add small exercice related with the initial stress * add sympy * fix dictionnary * more precise instructions * fix missing import module and format files naemes and variables with * add subtitle for ground deformation exercise * improve exercise with varying dip * add pandas (for moment rate release) * fix bug intoduced previously * change task title * fix clim for Vr --- Dockerfile | 19 +-- Dockerfile_jupyterlab | 14 +-- kaikoura/parametersLSW.par | 2 + kaikoura/parametersRS.par | 3 + sulawesi/parametersLSW.par | 2 + sulawesi/parametersRS.par | 2 + tpv13/parameters.par | 6 +- tpv13/tpv12_13_fault.yaml | 38 +++--- tpv13/tpv12_13_initial_stress.yaml | 5 +- tpv13/tpv13.ipynb | 187 +++++++++++++++++++++++------ tpv13/tpv13_training.geo | 9 +- 11 files changed, 201 insertions(+), 86 deletions(-) diff --git a/Dockerfile b/Dockerfile index e9c7a0f..3b50174 100644 --- a/Dockerfile +++ b/Dockerfile @@ -86,20 +86,15 @@ RUN git clone --recursive https://github.com/TUM-I5/ASAGI.git \ && CC=mpicc CXX=mpicxx cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DSHARED_LIB=off -DSTATIC_LIB=on -DNONUMA=on \ && make -j$(nproc) && make install -RUN git clone https://github.com/uphoffc/ImpalaJIT.git \ - && cd ImpalaJIT \ - && mkdir build && cd build \ - && cmake .. && make -j $(nproc) install - RUN git clone https://github.com/SeisSol/easi \ && cd easi \ && mkdir build && cd build \ - && CC=mpicc CXX=mpicxx cmake .. -DEASICUBE=OFF -DLUA=ON -DCMAKE_PREFIX_PATH=/home/tools -DCMAKE_INSTALL_PREFIX=/home/tools -DASAGI=ON -DIMPALAJIT=ON .. \ + && CC=mpicc CXX=mpicxx cmake .. -DEASICUBE=OFF -DLUA=ON -DCMAKE_PREFIX_PATH=/home/tools -DCMAKE_INSTALL_PREFIX=/home/tools -DASAGI=ON -DIMPALAJIT=OFF .. \ && make -j$(nproc) && make install RUN pip install numpy && docker-clean -RUN git clone --recursive --depth 1 --single-branch --branch v1.1.3 https://github.com/SeisSol/SeisSol.git \ +RUN git clone --recursive --depth 1 --single-branch --branch v1.3.0 https://github.com/SeisSol/SeisSol.git \ && cd SeisSol \ && mkdir build_hsw && cd build_hsw \ && export PATH=$PATH:/home/tools/bin \ @@ -116,13 +111,7 @@ RUN cd SeisSol/preprocessing/science/rconv \ && CC=mpicc CXX=mpicxx cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_PREFIX_PATH=/home/tools \ && make -j$(nproc) && cp rconv /home/tools/bin/ -RUN git clone --recursive --depth 1 --single-branch --branch v2.2.7 https://github.com/SCOREC/core.git \ - && cd core \ - && mkdir build && cd build \ - && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release -DSCOREC_CXX_FLAGS="-Wno-error=array-bounds" \ - && make -j$(nproc) && make install - -RUN git clone --recursive --branch v1.0.1 https://github.com/SeisSol/PUMGen.git \ +RUN git clone --recursive --branch v1.1.0 https://github.com/SeisSol/PUMGen.git \ && cd PUMGen \ && mkdir build && cd build \ && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release \ @@ -163,6 +152,8 @@ RUN conda install \ matplotlib \ gmsh \ python-gmsh \ + sympy \ + pandas \ && docker-clean ENV PATH=/home/tools/bin:$PATH ENV OMP_PLACES="cores" diff --git a/Dockerfile_jupyterlab b/Dockerfile_jupyterlab index 2d71c4f..fa05cfb 100644 --- a/Dockerfile_jupyterlab +++ b/Dockerfile_jupyterlab @@ -93,7 +93,7 @@ RUN git clone https://github.com/SeisSol/easi \ && CC=mpicc CXX=mpicxx cmake .. -DEASICUBE=OFF -DLUA=ON -DCMAKE_PREFIX_PATH=/home/tools -DCMAKE_INSTALL_PREFIX=/home/tools -DASAGI=ON -DIMPALAJIT=OFF .. \ && make -j$(nproc) && make install -RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git@davschneller/compile-fixes && docker-clean +RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git && docker-clean #RUN git clone --depth 1 --single-branch --branch main https://github.com/libxsmm/libxsmm.git \ # && cd libxsmm \ @@ -104,7 +104,7 @@ RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git@d # -j$(nproc) generator \ # && cp bin/libxsmm_gemm_generator /home/tools/bin -RUN git clone --recursive --depth 1 --single-branch --branch v1.1.3 https://github.com/SeisSol/SeisSol.git \ +RUN git clone --recursive --depth 1 --single-branch --branch v1.3.0 https://github.com/SeisSol/SeisSol.git \ && cd SeisSol \ && mkdir build_hsw && cd build_hsw \ && export PATH=$PATH:/home/tools/bin \ @@ -121,13 +121,7 @@ RUN cd SeisSol/preprocessing/science/rconv \ && CC=mpicc CXX=mpicxx cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_PREFIX_PATH=/home/tools \ && make -j$(nproc) && cp rconv /home/tools/bin/ -RUN git clone --recursive --depth 1 --single-branch --branch v2.2.7 https://github.com/SCOREC/core.git \ - && cd core \ - && mkdir build && cd build \ - && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release -DSCOREC_CXX_FLAGS="-Wno-error=array-bounds" \ - && make -j$(nproc) && make install - -RUN git clone --recursive https://github.com/SeisSol/PUMGen.git \ +RUN git clone --recursive --branch v1.1.0 https://github.com/SeisSol/PUMGen.git \ && cd PUMGen \ && mkdir build && cd build \ && cmake .. -DCMAKE_INSTALL_PREFIX=/home/tools -DCMAKE_C_COMPILER=mpicc -DCMAKE_CXX_COMPILER=mpicxx -DCMAKE_BUILD_TYPE=Release \ @@ -169,6 +163,8 @@ RUN conda install \ matplotlib \ gmsh \ python-gmsh \ + sympy \ + pandas \ && docker-clean ENV PATH=/home/tools/bin:$PATH ENV OMP_PLACES="cores" diff --git a/kaikoura/parametersLSW.par b/kaikoura/parametersLSW.par index 3949d68..4edf97b 100644 --- a/kaikoura/parametersLSW.par +++ b/kaikoura/parametersLSW.par @@ -7,6 +7,8 @@ Tv=0.05 !Attenuation parameters (ignored if not compiled with attenuation) FreqCentral=0.5 FreqRatio=100 +numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov' +numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov' / &IniCondition diff --git a/kaikoura/parametersRS.par b/kaikoura/parametersRS.par index 92fdbbb..9dd594e 100644 --- a/kaikoura/parametersRS.par +++ b/kaikoura/parametersRS.par @@ -7,6 +7,9 @@ Tv=0.05 !Attenuation parameters (ignored if not compiled with attenuation) FreqCentral=0.5 FreqRatio=100 + +numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov' +numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov' / &IniCondition diff --git a/sulawesi/parametersLSW.par b/sulawesi/parametersLSW.par index 8f0ca36..f2a5292 100644 --- a/sulawesi/parametersLSW.par +++ b/sulawesi/parametersLSW.par @@ -7,6 +7,8 @@ Tv=0.05 !Attenuation parameters (ignored if not compiled with attenuation) FreqCentral=0.3 FreqRatio=100 +numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov' +numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov' / &IniCondition diff --git a/sulawesi/parametersRS.par b/sulawesi/parametersRS.par index 6bed675..753e173 100644 --- a/sulawesi/parametersRS.par +++ b/sulawesi/parametersRS.par @@ -7,6 +7,8 @@ Tv=0.05 !Attenuation parameters (ignored if not compiled with attenuation) FreqCentral=0.3 FreqRatio=100 +numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov' +numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov' / &IniCondition diff --git a/tpv13/parameters.par b/tpv13/parameters.par index ee3be63..367c296 100644 --- a/tpv13/parameters.par +++ b/tpv13/parameters.par @@ -4,6 +4,8 @@ MaterialFileName = 'tpv12_13_material.yaml' !enable off-fault plasticity (ignored in Plasticity=0) Plasticity = 1 Tv = 0.03 +numflux = 'rusanov' ! The numerical flux. Values: 'godunov', 'rusanov' +numfluxnearfault = 'rusanov' ! The numerical flux for faces of cells adjacent to the fault (on cells which have at least one DR face), excluding the DR faces themselves. Values: 'godunov', 'rusanov' / &IniCondition @@ -34,6 +36,7 @@ printtimeinterval_sec = 1.0 ! Time interval at which output w OutputMask = 1 1 1 1 1 1 1 1 1 1 1 ! turn on and off fault outputs refinement_strategy = 2 refinement = 1 +!vtkorder = 3 / !parameterize ascii fault file outputs @@ -68,8 +71,9 @@ refinement = 1 ! Free surface output SurfaceOutput = 1 -SurfaceOutputRefinement = 1 +SurfaceOutputRefinement = 0 SurfaceOutputInterval = 2.0 +surfacevtkorder = 2 pickdt = 0.005 ! Pickpoint Sampling RFileName = 'tpv12_13_receivers.dat' ! Record Points in extra file diff --git a/tpv13/tpv12_13_fault.yaml b/tpv13/tpv12_13_fault.yaml index 1eb8f2a..d2141e3 100644 --- a/tpv13/tpv12_13_fault.yaml +++ b/tpv13/tpv12_13_fault.yaml @@ -1,26 +1,30 @@ !Switch [s_xx, s_yy, s_zz, s_xy, s_yz, s_xz]: !Include tpv12_13_initial_stress.yaml -[mu_s, mu_d, d_c, cohesion]: !IdentityMap - components: - # Inside nucleation patch - - !AxisAlignedCuboidalDomainFilter - limits: - x: [-1500, 1500] - y: [-.inf, .inf] - z: [-11691.34295108992, -9093.266739736605] - components: !ConstantMap - map: - mu_s: 0.54 - mu_d: 0.10 - d_c: 0.50 - cohesion: -200000 - # Outside nucleation patch - - !ConstantMap +[mu_d, d_c, cohesion]: !ConstantMap map: - mu_s: 0.70 mu_d: 0.10 d_c: 0.50 cohesion: -200000 +[mu_s]: !LuaMap + returns: [mu_s] + function: | + function f(x) + mu_s = 0.7 + nucleation_size = 3000.0 + along_width_nucleation_center = 12000.0 + sinDip = math.sin(60 * math.pi / 180) + x_in_nucleation = (math.abs(x["x"]) <= 0.5 * nucleation_size) + z_in_nucleation = ( + math.abs(x["z"] + along_width_nucleation_center * sinDip) <= + 0.5 * nucleation_size * sinDip + ) + if x_in_nucleation and z_in_nucleation then + mu_s = 0.54 + end + + return { mu_s = mu_s } + end + [Tnuc_n, Tnuc_s, Tnuc_d, forced_rupture_time]: !ConstantMap map: Tnuc_n: 0 diff --git a/tpv13/tpv12_13_initial_stress.yaml b/tpv13/tpv12_13_initial_stress.yaml index 77c517f..3490980 100644 --- a/tpv13/tpv12_13_initial_stress.yaml +++ b/tpv13/tpv12_13_initial_stress.yaml @@ -11,9 +11,10 @@ depth = math.abs(x["z"]) s_max_minus_Pf = 9.8 * (2700.0 - 1000.0); if (depth<=11951.15) then + s3_to_s1 = 0.3496 -- math.floor(a+0.5) is round - s_xx = -0.01 * math.floor(100.0 * (0.5 * (1.0 + 0.3496) * s_max_minus_Pf) + 0.5) * depth - s_yy = -0.01 * math.floor(100.0 * (0.3496 * s_max_minus_Pf) + 0.5) * depth + s_xx = -0.01 * math.floor(100.0 * (0.5 * (1.0 + s3_to_s1) * s_max_minus_Pf) + 0.5) * depth + s_yy = -0.01 * math.floor(100.0 * (s3_to_s1 * s_max_minus_Pf) + 0.5) * depth s_zz = -s_max_minus_Pf * depth else s_xx = -s_max_minus_Pf * depth diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb index 5fecdcc..3af0c25 100644 --- a/tpv13/tpv13.ipynb +++ b/tpv13/tpv13.ipynb @@ -46,16 +46,17 @@ "metadata": {}, "outputs": [], "source": [ - "!gmsh -3 tpv13_training.geo\n", + "!gmsh -3 tpv13_training.geo -setnumber dip 60\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", + "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)." @@ -67,16 +68,16 @@ "metadata": {}, "outputs": [], "source": [ - "!mpirun -n 1 pumgen -s msh2 tpv13_training.msh\n", + "!mpirun -n 1 pumgen -s msh4 tpv13_training.msh\n", "# on Frontera with apptainer, replace with:\n", - "# !mpirun -n 1 apptainer run {\"~/my-training.sif\"} pumgen -s msh2 tpv13_training.msh" + "# !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", + "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." ] }, @@ -108,8 +109,8 @@ "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", + "* 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." ] }, @@ -182,33 +183,78 @@ " 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", + " #fix color range for rupture speed\n", + " clim_arg = {\"clim\":(0, 6000)} if var=='Vr' else {}\n", + " plotter.add_mesh(mesh, cmap='Blues', scalars=var, **clim_arg)\n", " plotter.view_xz()\n", - " plotter.show(jupyter_backend='panel')" + " plotter.show(jupyter_backend='panel')\n", + " plotter.close()" ] }, { "cell_type": "markdown", + "id": "5291f06b-f8c9-473f-9361-eddfa993b0c1", "metadata": {}, "source": [ - "# Local-time stepping exercises" + "# Initial stress and fault strength\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": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "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", + "import sympy as sp\n", + "import numpy as np\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", + "# 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", - "* 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", + "# 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", - "*Hint:* You may reduce the *EndTime* parameter for this exercise." + "# 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 `s3_to_s1` 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`)." ] }, { @@ -219,9 +265,9 @@ "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", + "* 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" ] }, @@ -231,6 +277,9 @@ "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", @@ -238,21 +287,20 @@ "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", + "p = pv.Plotter(shape=(1, nSim), notebook=True, window_size=[1200, 500])\n", "\n", "def load_mesh(prefix):\n", - " # Read xdmf data\n", - " reader = vtk.vtkXdmfReader()\n", - " reader.SetFileName(f'output/{prefix}-surface.xdmf')\n", + " # Read free-surface output data\n", + " reader = vtk.vtkHDFReader()\n", + " reader.SetFileName(f'output/{prefix}-free-surface-4.vtkhdf')\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.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", @@ -261,17 +309,56 @@ " )\n", " p.camera.zoom(3)\n", "p.show(jupyter_backend='panel')\n", - " \n", - "# plot subsidence/uplift profile\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", - " 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", - " )" + " 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()" ] }, { @@ -322,6 +409,30 @@ "\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": { diff --git a/tpv13/tpv13_training.geo b/tpv13/tpv13_training.geo index b08e8cf..6865c03 100644 --- a/tpv13/tpv13_training.geo +++ b/tpv13/tpv13_training.geo @@ -10,7 +10,7 @@ SetFactory("OpenCASCADE"); // Length of the fault l_f = 30e3; -// Depth of the fault +// Width of the fault (along-dip) w_f = 15e3; dip_rad = dip*Pi/180.; @@ -18,7 +18,7 @@ dip_rad = dip*Pi/180.; // Hypocenter, z_n is here the depth along-dip x_n = 0e3; z_n = -12e3; -// size of the nucleaton patch +// Dimensions of the nucleation patch l_n = 3e3; w_n = 3e3; @@ -32,7 +32,7 @@ Z0 = -42e3; // Create the domain as a box domain = newv; Box(domain) = {X0, Y0, Z0, X1-X0, Y1-Y0, -Z0}; -// Create the fault as a vertically dipping rectangle, centered in x at the hypocenter +// Create the fault as a rectangle in the x-y plane, centered in x at the hypocenter fault = news; Rectangle(fault) = {-l_f/2, -w_f, 0, l_f, w_f}; // Create the nucleation patch as a smaller rectangle @@ -44,7 +44,7 @@ Rotate{ {1, 0, 0}, {0, 0, 0}, dip_rad } { Surface{fault, nucl}; } // Intersect the domain box with the fault rectangle at the free surface v() = BooleanFragments{ Volume{domain}; Delete; }{ Surface{fault, nucl}; Delete; }; -// Update all coordinates that define important surfaces within the mesh +// Extract surfaces within the mesh for defining boundary conditions eps = 1e-3; fault_final[] = Surface In BoundingBox{-l_f/2-eps, -w_f-eps, -w_f-eps, l_f/2+eps, w_f+eps, w_f+eps}; top[] = Surface In BoundingBox{X0-eps, Y0-eps, -eps, X1+eps, Y1+eps, eps}; @@ -66,4 +66,3 @@ Physical Surface(3) = {fault_final[]}; Physical Surface(5) = {other[]}; Physical Volume(1) = {domain}; -Mesh.MshFileVersion = 2.2;