From cc84e3025f4823432c0538ba33cf52d94dd2965b Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Mon, 4 Nov 2024 15:08:43 +0100
Subject: [PATCH 01/24] use rusanov

---
 Dockerfile                 | 9 ++-------
 kaikoura/parametersLSW.par | 2 ++
 kaikoura/parametersRS.par  | 3 +++
 sulawesi/parametersLSW.par | 2 ++
 sulawesi/parametersRS.par  | 2 ++
 tpv13/parameters.par       | 2 ++
 6 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/Dockerfile b/Dockerfile
index e9c7a0f..26c57a0 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.2.0 https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
     && export PATH=$PATH:/home/tools/bin \
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..be904c9 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

From e1cddd3a164652ace7be79e87e28a5ff6cebbf40 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Mon, 4 Nov 2024 15:21:52 +0100
Subject: [PATCH 02/24] fix Dockerfile_jupyterlab

---
 Dockerfile_jupyterlab | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Dockerfile_jupyterlab b/Dockerfile_jupyterlab
index 2d71c4f..64588c1 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.2.0 https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
     && export PATH=$PATH:/home/tools/bin \

From 4051c03771fe288bd21972d8f5a3585ffc47c0c6 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 10:24:47 +0100
Subject: [PATCH 03/24] small improvements on tpv13_training.geo

---
 tpv13/tpv13_training.geo | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/tpv13/tpv13_training.geo b/tpv13/tpv13_training.geo
index b08e8cf..48f859e 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 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;

From 010b1f013bc74da177d31391afe61e3e4a1dfa16 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 10:26:06 +0100
Subject: [PATCH 04/24] update to msh4

---
 tpv13/tpv13.ipynb | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 5fecdcc..80ccad3 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -67,9 +67,9 @@
    "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"
    ]
   },
   {

From 4e9b92e30fcce80b596a37e2ba1c5ed83f511b6f Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 11:06:47 +0100
Subject: [PATCH 05/24] switch to main

---
 Dockerfile        | 3 ++-
 tpv13/tpv13.ipynb | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/Dockerfile b/Dockerfile
index 26c57a0..ba9e451 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -94,7 +94,8 @@ RUN git clone https://github.com/SeisSol/easi \
 
 RUN pip install numpy && docker-clean
 
-RUN git clone --recursive --depth 1 --single-branch --branch v1.2.0 https://github.com/SeisSol/SeisSol.git \
+#RUN git clone --recursive --depth 1 --single-branch --branch v1.2.0 https://github.com/SeisSol/SeisSol.git \
+RUN git clone --recursive --depth 1 --single-branch --branch master https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
     && export PATH=$PATH:/home/tools/bin \
diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 80ccad3..3532576 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -269,7 +269,7 @@
     "        b,\n",
     "        resolution=10000,\n",
     "        title=f\"subsidence/uplift profile for {pref}\",\n",
-    "        ylabel=\"vectical displacement (m)\",\n",
+    "        ylabel=\"vertical displacement (m)\",\n",
     "        scalars='u3'\n",
     "    )"
    ]

From 2e007f3f468234933689de8f004da98d6e1be172 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 11:24:09 +0100
Subject: [PATCH 06/24] update both dockers

---
 Dockerfile            | 1 -
 Dockerfile_jupyterlab | 2 +-
 2 files changed, 1 insertion(+), 2 deletions(-)

diff --git a/Dockerfile b/Dockerfile
index ba9e451..9602861 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -94,7 +94,6 @@ RUN git clone https://github.com/SeisSol/easi \
 
 RUN pip install numpy && docker-clean
 
-#RUN git clone --recursive --depth 1 --single-branch --branch v1.2.0 https://github.com/SeisSol/SeisSol.git \
 RUN git clone --recursive --depth 1 --single-branch --branch master https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
diff --git a/Dockerfile_jupyterlab b/Dockerfile_jupyterlab
index 64588c1..e85aa0c 100644
--- a/Dockerfile_jupyterlab
+++ b/Dockerfile_jupyterlab
@@ -104,7 +104,7 @@ RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git &
 #       -j$(nproc) generator \
 #    && cp bin/libxsmm_gemm_generator /home/tools/bin
 
-RUN git clone --recursive --depth 1 --single-branch --branch v1.2.0 https://github.com/SeisSol/SeisSol.git \
+RUN git clone --recursive --depth 1 --single-branch --branch master https://github.com/SeisSol/SeisSol.git \
     && cd SeisSol \
     && mkdir build_hsw && cd build_hsw \
     && export PATH=$PATH:/home/tools/bin \

From 5abb3c77245b9cb6c746fc4540c066ed6846e5b4 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 14:39:01 +0100
Subject: [PATCH 07/24] use latest releases of SeisSol and PUMGen

---
 Dockerfile            | 10 ++--------
 Dockerfile_jupyterlab | 10 ++--------
 tpv13/parameters.par  |  3 ++-
 tpv13/tpv13.ipynb     |  7 +++----
 4 files changed, 9 insertions(+), 21 deletions(-)

diff --git a/Dockerfile b/Dockerfile
index 9602861..e59d226 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -94,7 +94,7 @@ RUN git clone https://github.com/SeisSol/easi \
 
 RUN pip install numpy && docker-clean
 
-RUN git clone --recursive --depth 1 --single-branch --branch master 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 \
@@ -111,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 \
diff --git a/Dockerfile_jupyterlab b/Dockerfile_jupyterlab
index e85aa0c..70e354a 100644
--- a/Dockerfile_jupyterlab
+++ b/Dockerfile_jupyterlab
@@ -104,7 +104,7 @@ RUN pip install numpy && pip install git+https://github.com/SeisSol/PSpaMM.git &
 #       -j$(nproc) generator \
 #    && cp bin/libxsmm_gemm_generator /home/tools/bin
 
-RUN git clone --recursive --depth 1 --single-branch --branch master 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 \
diff --git a/tpv13/parameters.par b/tpv13/parameters.par
index be904c9..3bc9dbd 100644
--- a/tpv13/parameters.par
+++ b/tpv13/parameters.par
@@ -70,8 +70,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/tpv13.ipynb b/tpv13/tpv13.ipynb
index 3532576..bdf0535 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -241,11 +241,10 @@
     "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",
+    "    # 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",

From 9e2929ca04fc978630c06264ebadb6c315406b7b Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 22:32:19 +0100
Subject: [PATCH 08/24] use vtk output for fault

---
 tpv13/parameters.par      |  3 +--
 tpv13/tpv12_13_fault.yaml | 37 ++++++++++++++++++++-----------------
 tpv13/tpv13.ipynb         | 32 ++++++++++++++++++++++++--------
 3 files changed, 45 insertions(+), 27 deletions(-)

diff --git a/tpv13/parameters.par b/tpv13/parameters.par
index 3bc9dbd..b1a164b 100644
--- a/tpv13/parameters.par
+++ b/tpv13/parameters.par
@@ -34,8 +34,7 @@ OutputPointType = 4
 &Elementwise
 printtimeinterval_sec = 1.0                    ! Time interval at which output will be written
 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
diff --git a/tpv13/tpv12_13_fault.yaml b/tpv13/tpv12_13_fault.yaml
index 1eb8f2a..07d073e 100644
--- a/tpv13/tpv12_13_fault.yaml
+++ b/tpv13/tpv12_13_fault.yaml
@@ -1,26 +1,29 @@
 !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 = 1500.0
+        along_width_nucleation_center = 12000.0
+        sinDip = math.sin(60 * math.pi / 180)
+        
+        x_in_nucleation = (math.abs(x["x"]) <= nucleation_size)
+        z_in_nucleation = (math.abs(x["z"] + along_width_nucleation_center * sinDip) <= 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/tpv13.ipynb b/tpv13/tpv13.ipynb
index bdf0535..05b73a3 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -170,16 +170,32 @@
    "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",
+    "def read_variable_names():\n",
+    "    # Initialize the reader and load the first file to extract variable names\n",
+    "    initial_reader = vtk.vtkHDFReader()\n",
+    "    initial_reader.SetFileName('output/tpv13-fault-elementwise-0.vtkhdf')\n",
+    "    initial_reader.Update()\n",
+    "    \n",
+    "    # Retrieve variables from PointData\n",
+    "    pd = initial_reader.GetOutput().GetPointData()\n",
+    "    return [pd.GetArrayName(i) for i in range(pd.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",
+    "variables = read_variable_names()\n",
+    "\n",
+    "@interact(t=(0, 8, 1), var=variables)\n",
+    "def plot(t=0, var='SRd'):\n",
+    "    # Construct the filename based on the selected time step\n",
+    "    filename = f'output/tpv13-fault-elementwise-{t}.vtkhdf'\n",
+    "    \n",
+    "    # Create and set up the reader for the selected file\n",
+    "    reader = vtk.vtkHDFReader()\n",
+    "    reader.SetFileName(filename)\n",
+    "    reader.Update()\n",
+    "    \n",
+    "    # Wrap the output in PyVista\n",
     "    mesh = pv.wrap(reader.GetOutput())\n",
+    "    \n",
+    "    # Set up the plotter and display the data\n",
     "    plotter = pv.Plotter(notebook=True)\n",
     "    plotter.set_background('white')\n",
     "    plotter.add_mesh(mesh, cmap='Blues', scalars=var)\n",

From 65cfc19049625f241f3af5c739a4e845b436af33 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Thu, 7 Nov 2024 23:28:01 +0100
Subject: [PATCH 09/24] add rupture speed exercise

---
 tpv13/tpv13.ipynb | 45 +++++++++++++++++++++++++++++----------------
 1 file changed, 29 insertions(+), 16 deletions(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 05b73a3..07cc0c1 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -207,24 +207,13 @@
    "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",
+    "# Rupture Speed Analysis\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",
+    "We analyze the mechanism causing supershear in this setup and disable it.\n",
     "\n",
-    "*Hint:* You may reduce the *EndTime* parameter for this exercise."
+    "* Examine the earthquake rupture speed and identify the conditions that trigger supershear rupture.\n",
+    "* Identify and modify the associated parameter to disable the current supershear triggering mechanism. \n",
+    "* You should now observe a delay in the supershear transition. This new supershear transition mechanism could be also suppressed by decreasing initial shear stress."
    ]
   },
   {
@@ -337,6 +326,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": {

From fbbf8b93040668c050a0d40dd603d59aaf52ce23 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sat, 9 Nov 2024 12:38:43 +0100
Subject: [PATCH 10/24] combine subsidence plots and fix typo in geo

---
 tpv13/tpv13.ipynb        | 26 ++++++++++++++++----------
 tpv13/tpv13_training.geo |  2 +-
 2 files changed, 17 insertions(+), 11 deletions(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 07cc0c1..55a1231 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -266,16 +266,22 @@
     "    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=\"vertical displacement (m)\",\n",
-    "        scalars='u3'\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.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",
+    "\n",
+    "plot_combined_profiles(prefixList, a, b)"
    ]
   },
   {
diff --git a/tpv13/tpv13_training.geo b/tpv13/tpv13_training.geo
index 48f859e..6865c03 100644
--- a/tpv13/tpv13_training.geo
+++ b/tpv13/tpv13_training.geo
@@ -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 in the x-y plane, 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

From 637d10ee8f85f61553e4d6d0da1dd5c48fad2b8d Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sat, 9 Nov 2024 14:28:41 +0100
Subject: [PATCH 11/24] close plotter to avoid memory leak

---
 tpv13/tpv13.ipynb | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 55a1231..e229b0f 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -200,7 +200,8 @@
     "    plotter.set_background('white')\n",
     "    plotter.add_mesh(mesh, cmap='Blues', scalars=var)\n",
     "    plotter.view_xz()\n",
-    "    plotter.show(jupyter_backend='panel')"
+    "    plotter.show(jupyter_backend='panel')\n",
+    "    plotter.close()"
    ]
   },
   {
@@ -265,7 +266,8 @@
     "    )\n",
     "    p.camera.zoom(3)\n",
     "p.show(jupyter_backend='panel')\n",
-    "    \n",
+    "p.close()\n",
+    "\n",
     "def plot_combined_profiles(prefixList, a, b, resolution=10000):\n",
     "    plt.figure(figsize=(6, 4))\n",
     "    for prefix in prefixList:\n",

From 05c2e87bc345c4f8d161ad5a9568391f39bc8e73 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sat, 9 Nov 2024 15:07:45 +0100
Subject: [PATCH 12/24] precise exercise and input file

---
 tpv13/tpv12_13_fault.yaml | 6 +++---
 tpv13/tpv13.ipynb         | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/tpv13/tpv12_13_fault.yaml b/tpv13/tpv12_13_fault.yaml
index 07d073e..b67aad2 100644
--- a/tpv13/tpv12_13_fault.yaml
+++ b/tpv13/tpv12_13_fault.yaml
@@ -10,12 +10,12 @@
   function: |
     function f(x)
         mu_s = 0.7
-        nucleation_size = 1500.0
+        nucleation_size = 3000.0
         along_width_nucleation_center = 12000.0
         sinDip = math.sin(60 * math.pi / 180)
         
-        x_in_nucleation = (math.abs(x["x"]) <= nucleation_size)
-        z_in_nucleation = (math.abs(x["z"] + along_width_nucleation_center * sinDip) <= nucleation_size * sinDip)
+        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
diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index e229b0f..98beecb 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -213,7 +213,7 @@
     "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",
-    "* Identify and modify the associated parameter to disable the current supershear triggering mechanism. \n",
+    "* Change the nucleation size from 3 to 1.2km in both tpv12_13_fault.yaml and tpv13_training.geo to disable the current supershear triggering mechanism. \n",
     "* You should now observe a delay in the supershear transition. This new supershear transition mechanism could be also suppressed by decreasing initial shear stress."
    ]
   },

From e7277370343cfa6278bfb2a3329a55a783d8ab79 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sat, 9 Nov 2024 16:45:44 +0100
Subject: [PATCH 13/24] disable vtk fault output for now

---
 tpv13/parameters.par               |  4 +++-
 tpv13/tpv12_13_initial_stress.yaml |  5 +++--
 tpv13/tpv13.ipynb                  | 36 +++++++++---------------------
 3 files changed, 16 insertions(+), 29 deletions(-)

diff --git a/tpv13/parameters.par b/tpv13/parameters.par
index b1a164b..367c296 100644
--- a/tpv13/parameters.par
+++ b/tpv13/parameters.par
@@ -34,7 +34,9 @@ OutputPointType = 4
 &Elementwise
 printtimeinterval_sec = 1.0                    ! Time interval at which output will be written
 OutputMask = 1 1 1 1 1 1 1 1 1 1 1             ! turn on and off fault outputs
-vtkorder = 3
+refinement_strategy = 2
+refinement = 1
+!vtkorder = 3
 /
 
 !parameterize ascii fault file outputs
diff --git a/tpv13/tpv12_13_initial_stress.yaml b/tpv13/tpv12_13_initial_stress.yaml
index 77c517f..758715c 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
+          sigma_3_to_sigma_1 = 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 + sigma_3_to_sigma_1) * s_max_minus_Pf) + 0.5) * depth
+          s_yy =  -0.01 * math.floor(100.0 *               (sigma_3_to_sigma_1 * 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 98beecb..0eb871f 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -170,32 +170,16 @@
    "source": [
     "from ipywidgets import interact\n",
     "\n",
-    "def read_variable_names():\n",
-    "    # Initialize the reader and load the first file to extract variable names\n",
-    "    initial_reader = vtk.vtkHDFReader()\n",
-    "    initial_reader.SetFileName('output/tpv13-fault-elementwise-0.vtkhdf')\n",
-    "    initial_reader.Update()\n",
-    "    \n",
-    "    # Retrieve variables from PointData\n",
-    "    pd = initial_reader.GetOutput().GetPointData()\n",
-    "    return [pd.GetArrayName(i) for i in range(pd.GetNumberOfArrays())]\n",
-    "\n",
-    "variables = read_variable_names()\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, 8, 1), var=variables)\n",
-    "def plot(t=0, var='SRd'):\n",
-    "    # Construct the filename based on the selected time step\n",
-    "    filename = f'output/tpv13-fault-elementwise-{t}.vtkhdf'\n",
-    "    \n",
-    "    # Create and set up the reader for the selected file\n",
-    "    reader = vtk.vtkHDFReader()\n",
-    "    reader.SetFileName(filename)\n",
-    "    reader.Update()\n",
-    "    \n",
-    "    # Wrap the output in PyVista\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",
-    "    \n",
-    "    # Set up the plotter and display the data\n",
     "    plotter = pv.Plotter(notebook=True)\n",
     "    plotter.set_background('white')\n",
     "    plotter.add_mesh(mesh, cmap='Blues', scalars=var)\n",
@@ -213,8 +197,8 @@
     "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 and tpv13_training.geo to disable the current supershear triggering mechanism. \n",
-    "* You should now observe a delay in the supershear transition. This new supershear transition mechanism could be also suppressed by decreasing initial shear stress."
+    "* Change the nucleation size from 3 to 1.2km in both tpv12_13_fault.yaml and tpv13_training.geo to disable the current supershear triggering mechanism. 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)."
    ]
   },
   {

From 8d7c653da9e03b574b1d9d41a5d40ee8fe56d8a6 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 11:43:25 +0100
Subject: [PATCH 14/24] add small exercice related with the initial stress

---
 tpv13/tpv13.ipynb | 54 ++++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 53 insertions(+), 1 deletion(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 0eb871f..1a77948 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -190,9 +190,61 @@
   },
   {
    "cell_type": "markdown",
+   "id": "5291f06b-f8c9-473f-9361-eddfa993b0c1",
    "metadata": {},
    "source": [
-    "# Rupture Speed Analysis\n",
+    "# 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",

From b366ca6f1539083691f8db23db97180b7d7c3c69 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 11:45:12 +0100
Subject: [PATCH 15/24] add sympy

---
 Dockerfile            | 1 +
 Dockerfile_jupyterlab | 1 +
 2 files changed, 2 insertions(+)

diff --git a/Dockerfile b/Dockerfile
index e59d226..ba1c8a2 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -152,6 +152,7 @@ RUN conda install \
     matplotlib \
     gmsh \
     python-gmsh \
+    sympy \
     && docker-clean
 ENV PATH=/home/tools/bin:$PATH
 ENV OMP_PLACES="cores"
diff --git a/Dockerfile_jupyterlab b/Dockerfile_jupyterlab
index 70e354a..c5918b8 100644
--- a/Dockerfile_jupyterlab
+++ b/Dockerfile_jupyterlab
@@ -163,6 +163,7 @@ RUN conda install \
     matplotlib \
     gmsh \
     python-gmsh \
+    sympy \
     && docker-clean
 ENV PATH=/home/tools/bin:$PATH
 ENV OMP_PLACES="cores"

From e73d4ed6f04c4b8e3a2f086323a6ae4d1355434a Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 14:24:33 +0100
Subject: [PATCH 16/24] fix dictionnary

---
 tpv13/tpv13.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 1a77948..38f6e06 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -229,7 +229,7 @@
     "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",
+    "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",

From 933bcaf9637c968c92de6313b5787bedc0b97122 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 14:28:29 +0100
Subject: [PATCH 17/24] more precise instructions

---
 tpv13/tpv13.ipynb | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 38f6e06..731fd71 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -249,7 +249,8 @@
     "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 and tpv13_training.geo to disable the current supershear triggering mechanism. You should now observe a delay in the supershear transition.\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)."
    ]
   },

From 358069b5fb3f586d9293d1d1bc6a9c8246ddbc14 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 22:14:33 +0100
Subject: [PATCH 18/24] fix missing import module and format files naemes and
 variables with

---
 tpv13/tpv13.ipynb | 23 +++++++++++++----------
 1 file changed, 13 insertions(+), 10 deletions(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 731fd71..858fbc7 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -53,9 +53,10 @@
   },
   {
    "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,7 +68,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "!mpirun -n 1 pumgen -s msh4 tpv13_training.msh\n",
+    "!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"
    ]
@@ -76,7 +77,7 @@
    "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."
    ]
   },
@@ -249,9 +250,9 @@
     "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",
+    "* 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)."
+    "* 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`)."
    ]
   },
   {
@@ -262,9 +263,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"
    ]
   },
@@ -274,6 +275,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "import matplotlib.pyplot as plt\n",
+    "\n",
     "prefixList = ['tpv13']\n",
     "nSim = len(prefixList)\n",
     "# coordinates of the line end points\n",

From a5ef4121700de487449568d9d9e11ad4d96a60c6 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 22:26:38 +0100
Subject: [PATCH 19/24] add subtitle for ground deformation exercise

---
 tpv13/tpv13.ipynb | 1 +
 1 file changed, 1 insertion(+)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 858fbc7..1db0686 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -298,6 +298,7 @@
     "    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",

From 03de6005897f4e4098e62314d3408cc78b1ec123 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 23:03:18 +0100
Subject: [PATCH 20/24] improve exercise with varying dip

---
 tpv13/tpv13.ipynb | 40 ++++++++++++++++++++++++++++++++++++----
 1 file changed, 36 insertions(+), 4 deletions(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 1db0686..67dde28 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -276,6 +276,7 @@
    "outputs": [],
    "source": [
     "import matplotlib.pyplot as plt\n",
+    "import pandas as pd\n",
     "\n",
     "prefixList = ['tpv13']\n",
     "nSim = len(prefixList)\n",
@@ -284,7 +285,7 @@
     "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 free-surface output data\n",
@@ -296,7 +297,6 @@
     "# 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",
@@ -309,6 +309,24 @@
     "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",
@@ -317,14 +335,28 @@
     "        distance = sampled_data['Distance']\n",
     "        displacement = sampled_data['u3']\n",
     "        plt.plot(distance/1000, displacement, label=prefix, linewidth=2)\n",
-    "    plt.title(\"Comparison of Subsidence/Uplift Profiles\")\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_combined_profiles(prefixList, a, b)"
+    "# 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()"
    ]
   },
   {

From d6cad509904144414fe65d5b35f081f76dd999cf Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Sun, 10 Nov 2024 23:04:30 +0100
Subject: [PATCH 21/24] add pandas (for moment rate release)

---
 Dockerfile            | 1 +
 Dockerfile_jupyterlab | 1 +
 2 files changed, 2 insertions(+)

diff --git a/Dockerfile b/Dockerfile
index ba1c8a2..3b50174 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -153,6 +153,7 @@ RUN conda install \
     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 c5918b8..fa05cfb 100644
--- a/Dockerfile_jupyterlab
+++ b/Dockerfile_jupyterlab
@@ -164,6 +164,7 @@ RUN conda install \
     gmsh \
     python-gmsh \
     sympy \
+    pandas \
     && docker-clean
 ENV PATH=/home/tools/bin:$PATH
 ENV OMP_PLACES="cores"

From cfc93eb1b27c5a06d1ed359d04b53aa697307a1d Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Mon, 11 Nov 2024 09:27:55 +0100
Subject: [PATCH 22/24] fix bug intoduced previously

---
 tpv13/tpv12_13_fault.yaml          | 7 ++++---
 tpv13/tpv12_13_initial_stress.yaml | 6 +++---
 tpv13/tpv13.ipynb                  | 6 +++---
 3 files changed, 10 insertions(+), 9 deletions(-)

diff --git a/tpv13/tpv12_13_fault.yaml b/tpv13/tpv12_13_fault.yaml
index b67aad2..d2141e3 100644
--- a/tpv13/tpv12_13_fault.yaml
+++ b/tpv13/tpv12_13_fault.yaml
@@ -13,10 +13,11 @@
         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)
-        
+        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
diff --git a/tpv13/tpv12_13_initial_stress.yaml b/tpv13/tpv12_13_initial_stress.yaml
index 758715c..3490980 100644
--- a/tpv13/tpv12_13_initial_stress.yaml
+++ b/tpv13/tpv12_13_initial_stress.yaml
@@ -11,10 +11,10 @@
         depth = math.abs(x["z"])
         s_max_minus_Pf = 9.8 * (2700.0 - 1000.0);
         if (depth<=11951.15) then
-          sigma_3_to_sigma_1 = 0.3496
+          s3_to_s1 = 0.3496
           -- math.floor(a+0.5) is round
-          s_xx =  -0.01 * math.floor(100.0 * (0.5 * (1.0 + sigma_3_to_sigma_1) * s_max_minus_Pf) + 0.5) * depth
-          s_yy =  -0.01 * math.floor(100.0 *               (sigma_3_to_sigma_1 * 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 67dde28..3e91c1d 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -46,7 +46,7 @@
    "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"
    ]
@@ -68,7 +68,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "!mpirun -n 1 pumgen -s msh4 tpv13_training.msh  -setnumber dip 60\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 msh4 tpv13_training.msh"
    ]
@@ -252,7 +252,7 @@
     "* 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`)."
+    "* 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`)."
    ]
   },
   {

From 7345b35785db63ab56ea6f116021d4b05f4bc262 Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Mon, 11 Nov 2024 10:09:33 +0100
Subject: [PATCH 23/24] change task title

---
 tpv13/tpv13.ipynb | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index 3e91c1d..e312d01 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -194,7 +194,7 @@
    "id": "5291f06b-f8c9-473f-9361-eddfa993b0c1",
    "metadata": {},
    "source": [
-    "# Initial stress\n",
+    "# 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",

From 459abea4a5f67a2eaeeb0a85e459dff94156d67f Mon Sep 17 00:00:00 2001
From: Thomas-Ulrich <ulrich@geophysik.uni-muenchen.de>
Date: Mon, 11 Nov 2024 10:19:27 +0100
Subject: [PATCH 24/24] fix clim for Vr

---
 tpv13/tpv13.ipynb | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/tpv13/tpv13.ipynb b/tpv13/tpv13.ipynb
index e312d01..3af0c25 100644
--- a/tpv13/tpv13.ipynb
+++ b/tpv13/tpv13.ipynb
@@ -183,7 +183,9 @@
     "    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')\n",
     "    plotter.close()"