diff --git a/ProcessLib/ComponentTransport/Tests.cmake b/ProcessLib/ComponentTransport/Tests.cmake
index d966b69f208f82258133d841c2ac2ed276f84715..3d1f7e05d9b23892accfabe70c2c288d92fa9bc0 100644
--- a/ProcessLib/ComponentTransport/Tests.cmake
+++ b/ProcessLib/ComponentTransport/Tests.cmake
@@ -952,7 +952,7 @@ if(NOT OGS_USE_PETSC)
NotebookTest(NOTEBOOKFILE Parabolic/LiquidFlow/AxiSymTheis/axisym_theis.py RUNTIME 10)
NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/ReactiveTransport/DecayChain/DecayChain.py RUNTIME 160)
NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/ReactiveTransport/RadionuclidesMigration/RadionuclidesMigration.py RUNTIME 55)
- NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.md RUNTIME 5)
+ NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.py RUNTIME 5)
NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/MultiLayerDiffusion/MultiLayerDiffusion.py RUNTIME 25)
NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/DiffusionSorptionDecay/DiffusionSorptionDecay.py RUNTIME 16)
NotebookTest(NOTEBOOKFILE Parabolic/ComponentTransport/elder_jupyter/elder_jupyter.py RUNTIME 25)
diff --git a/ProcessLib/HeatConduction/Tests.cmake b/ProcessLib/HeatConduction/Tests.cmake
index a2fa990c28c9542c22c1ee9e54532deaea4d377f..df2aa4f98c5f6ad8e3d0c57049c9efe1bf62302e 100644
--- a/ProcessLib/HeatConduction/Tests.cmake
+++ b/ProcessLib/HeatConduction/Tests.cmake
@@ -635,5 +635,5 @@ AddTest(
if(NOT OGS_USE_PETSC)
# Both tests above are executed in this notebook (without diff check). Maybe remove regular tests later?
- NotebookTest(NOTEBOOKFILE Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.md RUNTIME 10)
+ NotebookTest(NOTEBOOKFILE Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.py RUNTIME 10)
endif()
diff --git a/ProcessLib/HeatTransportBHE/Tests.cmake b/ProcessLib/HeatTransportBHE/Tests.cmake
index ddaa11ce26d82c80d1c9cfb4cddff222fd306dd8..028d97295a5d76783ba2507fbd74daed07696b38 100644
--- a/ProcessLib/HeatTransportBHE/Tests.cmake
+++ b/ProcessLib/HeatTransportBHE/Tests.cmake
@@ -430,5 +430,5 @@ AddTest(
)
if(NOT OGS_USE_PETSC)
- NotebookTest(NOTEBOOKFILE Parabolic/T/BHE_1P/pipe_flow_ebhe.md RUNTIME 200)
+ NotebookTest(NOTEBOOKFILE Parabolic/T/BHE_1P/pipe_flow_ebhe.py RUNTIME 200)
endif()
diff --git a/ProcessLib/HydroMechanics/Tests.cmake b/ProcessLib/HydroMechanics/Tests.cmake
index 255bb862e05d5c581611555b405620abe00c02d0..5cf817c8b268ef0aafd961d2e2fc9c0c6f6e6f7a 100644
--- a/ProcessLib/HydroMechanics/Tests.cmake
+++ b/ProcessLib/HydroMechanics/Tests.cmake
@@ -706,7 +706,7 @@ AddTest(
reference_results_MandelCryerStaggered_ts_1_t_0.010000.vtu results_MandelCryerStaggered_ts_1_t_0.010000.vtu sigma sigma 1e-1 0
)
if(NOT OGS_USE_MPI)
- NotebookTest(NOTEBOOKFILE HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.md
+ NotebookTest(NOTEBOOKFILE HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.py
RUNTIME 800)
endif()
diff --git a/ProcessLib/LargeDeformation/Tests.cmake b/ProcessLib/LargeDeformation/Tests.cmake
index d3e61f913b887eb651c53f90a1088300e04b23a5..6a214907d1c1c66522486311d8ef1b1994daa35e 100644
--- a/ProcessLib/LargeDeformation/Tests.cmake
+++ b/ProcessLib/LargeDeformation/Tests.cmake
@@ -1,8 +1,8 @@
if (NOT OGS_USE_MPI)
OgsTest(PROJECTFILE LargeDeformation/RigidBody/square_1e0.prj RUNTIME 1)
if(OGS_BUILD_PROCESS_SmallDeformation)
- NotebookTest(NOTEBOOKFILE LargeDeformation/RigidBody/RigidBody.md RUNTIME 15)
+ NotebookTest(NOTEBOOKFILE LargeDeformation/RigidBody/RigidBody.py RUNTIME 15)
endif()
OgsTest(PROJECTFILE LargeDeformation/Torsion/bar1to6_torsion.prj RUNTIME 5)
- NotebookTest(NOTEBOOKFILE LargeDeformation/Torsion/Torsion_robustness.md RUNTIME 15)
+ NotebookTest(NOTEBOOKFILE LargeDeformation/Torsion/Torsion_robustness.py RUNTIME 15)
endif()
diff --git a/ProcessLib/SmallDeformation/Tests.cmake b/ProcessLib/SmallDeformation/Tests.cmake
index ef59fa740a51741f90557a3207348b62504baa64..60d64477467273c3b5d244ad6c617ecff70f8eba 100644
--- a/ProcessLib/SmallDeformation/Tests.cmake
+++ b/ProcessLib/SmallDeformation/Tests.cmake
@@ -354,14 +354,14 @@ if(NOT OGS_USE_PETSC)
NotebookTest(NOTEBOOKFILE Mechanics/CooksMembrane/CooksMembraneBbar.py RUNTIME 1)
NotebookTest(NOTEBOOKFILE Mechanics/Linear/SimpleMechanics.py RUNTIME 5)
NotebookTest(NOTEBOOKFILE Mechanics/EvaluatingBbarWithSimpleExamples/evaluating_bbbar_with_simple_examples.py RUNTIME 5)
- NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md RUNTIME 15)
+ NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.py RUNTIME 15)
if(NOT WIN32)
NotebookTest(NOTEBOOKFILE Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole_convergence_analysis.py RUNTIME 40)
endif()
if (OGS_USE_MFRONT)
NotebookTest(NOTEBOOKFILE Mechanics/PLLC/PLLC.py RUNTIME 7)
if(TFEL_WITH_PYTHON)
- NotebookTest(NOTEBOOKFILE Mechanics/HoekBrown/HoekBrownYieldCriterion.md RUNTIME 20)
+ NotebookTest(NOTEBOOKFILE Mechanics/HoekBrown/HoekBrownYieldCriterion.py RUNTIME 20)
endif()
endif()
endif()
diff --git a/Tests/Data/HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.md b/Tests/Data/HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.md
deleted file mode 100644
index 8acdbbdf4f85e8cd2cc454024097089a36c48c0b..0000000000000000000000000000000000000000
--- a/Tests/Data/HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.md
+++ /dev/null
@@ -1,130 +0,0 @@
-+++
-date = "2022-01-25"
-title = "Mandel-Cryer Effect"
-weight = 151
-project = ["HydroMechanics/StaggeredScheme/MandelCryer/MandelCryerStaggered.prj"]
-author = "Dominik Kern, Frieder Loer"
-image = "./figures/MandelCryer_mesh.png"
-web_subsection = "hydro-mechanics"
-+++
-
-
-## Mandel-Cryer Effect
-
-This is a classical example to demonstrate the effect of hydromechanical coupling in a poroelastic medium.
-For more details we refer to a textbook (Verruijt, 2009), in which also the analytical solution is derived.
-As domain we consider a sphere, by symmetry we need to simulate only an octant.
-
-
-
-The boundary conditions of hydraulics are atmospheric pressure on the sphere surface and impermeable elsewhere.
-The boundary conditions of mechanics are an overburden (Neumann) applied as step load on the sphere surface at initial time $t=0$.
-The remaining sides are fixed in normal direction (Dirichlet).
-
-
-
-The material is isotropic, linear elastic.
-Solid and fluid are assumed to be incompressible.
-In its initial state the sphere is not deformed and there is ambient pressure everywhere.
-A sudden load increase on the surface is instantly transferred on the pore pressure, whereas the solid needs time to deform, until it carries the load.
-Since the pore fluid is squeezed out of the outer layers first, they act like a tightening belt and consequently the pressure in the center rises, it may even exceed the value of the applied load.
-Finally the pore pressure approaches to ambient pressure.
-
-All parameters are concluded in the following tables.
-
-### Material Properties
-
-| Property | Value | Unit |
-| ------------------------- | -------------- | ------------ |
-| Fluid density | $10^3$ | kg/m$^3$ |
-| Viscosity | $10^{-3}$ | Pa$\cdot$s |
-| Porosity | $0.2$ | - |
-| Permeability | $10\cdot 10^{-12}$ | m$^2$ |
-| Young’s modulus (bulk) | $10\cdot 10^6$ | Pa |
-| Poisson ratio (bulk) | $0.1$ | - |
-| Biot coefficient | $1.0$ | - |
-| Solid density | $2.5\cdot 10^3$| kg/m$^3$ |
-| Overburden | $1000$ | Pa |
-| Atmospheric pressure | $0$ | Pa |
-
-### Dimensions and Discretization
-
-| Property | Value | Unit |
-| -------------------------- | -------- | ----------------------- |
-| Radius | $0.4$ | m |
-| Finite Elements | $8741$ | Taylor-Hood tetrahedra |
-| Time step | $10^{-2}$| s |
-| Coupling scheme parameter | $0.7774$ | - |
-
-## Numerical Simulation
-
-```python
-# Create output path if it doesn't exist yet
-import os
-from pathlib import Path
-
-out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
-if not out_dir.exists():
- out_dir.mkdir(parents=True)
-
-import ogstools as ot
-```
-
-```python
-# Initiate an OGS-object
-# Pass it the project file and set an output file
-model = ot.Project(input_file="MandelCryerStaggered.prj", output_file=f"{out_dir}/MandelCryerStaggered_modified.prj")
-
-# Increase end time
-t_end = 1.5
-model.replace_text(t_end, xpath="./time_loop/processes/process/time_stepping/t_end")
-model.write_input()
-
-# Run OGS
-model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
-```
-
-## Results
-
-```python
-import vtuIO
-import numpy as np
-
-# Load the result
-pvdfile = vtuIO.PVDIO(f"{out_dir}/results_MandelCryerStaggered.pvd", dim=3)
-
-# Get point field names
-fields = pvdfile.get_point_field_names()
-print(fields)
-
-# Get all written timesteps
-time = pvdfile.timesteps
-```
-
-```python
-import matplotlib.pyplot as plt
-
-
-# Define observation points
-observation_points = {'center':(0, 0, 0)}
-pressure = pvdfile.read_time_series('pressure', observation_points)
-
-# Plot soil temperature
-fig, ax1 = plt.subplots(figsize=(10, 8))
-ax1.plot(time, pressure['center'], color="tab:orange")
-ax1.set_ylabel('Pressure (Pa)', fontsize=20)
-ax1.set_xlabel('Time (s)', fontsize=20)
-ax1.set_xlim(0, t_end)
-ax1.set_ylim(0, 1500)
-ax1.grid()
-fig.supxlabel('Pressure at center of sphere')
-plt.show()
-```
-
-As predicted, the pressure in the center exceeds the applied load and then levels out.
-
-For more details about the staggered scheme we refer to the [user guide - conventions]({{< ref "conventions" >}}#staggered-scheme).
-
-## References
-
-[1] Verruijt, A. (2009): _An introduction to soil dynamics_. Springer Science and Business Media, DOI: <https://doi.org/10.1007/978-90-481-3441-0>, <https://link.springer.com/book/10.1007/978-90-481-3441-0>
diff --git a/Tests/Data/HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.py b/Tests/Data/HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.py
new file mode 100644
index 0000000000000000000000000000000000000000..f9c8319a54ae3ab486a493fea4bc8bb2c4c77314
--- /dev/null
+++ b/Tests/Data/HydroMechanics/StaggeredScheme/MandelCryer/mandelcryer.py
@@ -0,0 +1,128 @@
+# %% [markdown]
+# +++
+# date = "2022-01-25"
+# title = "Mandel-Cryer Effect"
+# weight = 151
+# project = ["HydroMechanics/StaggeredScheme/MandelCryer/MandelCryerStaggered.prj"]
+# author = "Dominik Kern, Frieder Loer"
+# image = "./figures/MandelCryer_mesh.png"
+# web_subsection = "hydro-mechanics"
+# +++
+
+# %% [markdown]
+# ## Mandel-Cryer Effect
+#
+# This is a classical example to demonstrate the effect of hydromechanical coupling in a poroelastic medium.
+# For more details we refer to a textbook (Verruijt, 2009), in which also the analytical solution is derived.
+# As domain we consider a sphere, by symmetry we need to simulate only an octant.
+#
+# 
+#
+# The boundary conditions of hydraulics are atmospheric pressure on the sphere surface and impermeable elsewhere.
+# The boundary conditions of mechanics are an overburden (Neumann) applied as step load on the sphere surface at initial time $t=0$.
+# The remaining sides are fixed in normal direction (Dirichlet).
+#
+# 
+#
+# The material is isotropic, linear elastic.
+# Solid and fluid are assumed to be incompressible.
+# In its initial state the sphere is not deformed and there is ambient pressure everywhere.
+# A sudden load increase on the surface is instantly transferred on the pore pressure, whereas the solid needs time to deform, until it carries the load.
+# Since the pore fluid is squeezed out of the outer layers first, they act like a tightening belt and consequently the pressure in the center rises, it may even exceed the value of the applied load.
+# Finally the pore pressure approaches to ambient pressure.
+#
+# All parameters are concluded in the following tables.
+#
+# ### Material Properties
+#
+# | Property | Value | Unit |
+# | ------------------------- | -------------- | ------------ |
+# | Fluid density | $10^3$ | kg/m$^3$ |
+# | Viscosity | $10^{-3}$ | Pa$\cdot$s |
+# | Porosity | $0.2$ | - |
+# | Permeability | $10\cdot 10^{-12}$ | m$^2$ |
+# | Young`s modulus (bulk) | $10\cdot 10^6$ | Pa |
+# | Poisson ratio (bulk) | $0.1$ | - |
+# | Biot coefficient | $1.0$ | - |
+# | Solid density | $2.5\cdot 10^3$| kg/m$^3$ |
+# | Overburden | $1000$ | Pa |
+# | Atmospheric pressure | $0$ | Pa |
+#
+# ### Dimensions and Discretization
+#
+# | Property | Value | Unit |
+# | -------------------------- | -------- | ----------------------- |
+# | Radius | $0.4$ | m |
+# | Finite Elements | $8741$ | Taylor-Hood tetrahedra |
+# | Time step | $10^{-2}$| s |
+# | Coupling scheme parameter | $0.7774$ | - |
+#
+# ## Numerical Simulation
+
+# %%
+# Create output path if it doesn't exist yet
+import os
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import ogstools as ot
+import vtuIO
+
+out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
+if not out_dir.exists():
+ out_dir.mkdir(parents=True)
+
+# %%
+# Initiate an OGS-object
+# Pass it the project file and set an output file
+model = ot.Project(
+ input_file="MandelCryerStaggered.prj",
+ output_file=f"{out_dir}/MandelCryerStaggered_modified.prj",
+)
+
+# Increase end time
+t_end = 1.5
+model.replace_text(t_end, xpath="./time_loop/processes/process/time_stepping/t_end")
+model.write_input()
+
+# Run OGS
+model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
+
+# %% [markdown]
+# ## Results
+
+# %%
+# Load the result
+pvdfile = vtuIO.PVDIO(f"{out_dir}/results_MandelCryerStaggered.pvd", dim=3)
+
+# Get point field names
+fields = pvdfile.get_point_field_names()
+print(fields)
+
+# Get all written timesteps
+time = pvdfile.timesteps
+
+# %%
+# Define observation points
+observation_points = {"center": (0, 0, 0)}
+pressure = pvdfile.read_time_series("pressure", observation_points)
+
+# Plot soil temperature
+fig, ax1 = plt.subplots(figsize=(10, 8))
+ax1.plot(time, pressure["center"], color="tab:orange")
+ax1.set_ylabel("Pressure (Pa)", fontsize=20)
+ax1.set_xlabel("Time (s)", fontsize=20)
+ax1.set_xlim(0, t_end)
+ax1.set_ylim(0, 1500)
+ax1.grid()
+fig.supxlabel("Pressure at center of sphere")
+plt.show()
+
+# %% [markdown]
+# As predicted, the pressure in the center exceeds the applied load and then levels out.
+#
+# For more details about the staggered scheme we refer to the [user guide - conventions]({{< ref "conventions" >}}#staggered-scheme).
+#
+# ## References
+#
+# [1] Verruijt, A. (2009): _An introduction to soil dynamics_. Springer Science and Business Media, DOI: <https://doi.org/10.1007/978-90-481-3441-0>, <https://link.springer.com/book/10.1007/978-90-481-3441-0>
diff --git a/Tests/Data/LargeDeformation/RigidBody/RigidBody.md b/Tests/Data/LargeDeformation/RigidBody/RigidBody.md
deleted file mode 100644
index 5ee1e132436e90a7b6113cc9b1b76d357b921eb7..0000000000000000000000000000000000000000
--- a/Tests/Data/LargeDeformation/RigidBody/RigidBody.md
+++ /dev/null
@@ -1,232 +0,0 @@
-+++
-title = "LargeDeformations: Rigid Body Rotation"
-date = "2023-03-01"
-author = "Thomas Nagel"
-web_subsection = "large-deformations"
-+++
-
-
-|<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
-|---|---|--:|
-
-In this benchmark we test a basic kinematic feature of the finite strain / large
-deformations implementation in OpenGeoSys.
-An element is subjected to a rigid body rotation. The expected result is a
-stress- and strain-free motion.
-This is confirmed for the finite-strain implementation while the small-strain
-implementation shows phantom strains and stresses.
-
-## Basic definitions
-
-The small deformation code uses the linearized strain tensor:
-
-\begin{align}
- \boldsymbol{\epsilon} = \frac{1}{2} \left( \text{grad}\, \mathbf{u} + \text{grad}\,^\text{T} \mathbf{u} \right)
-\end{align}
-
-while the large deformation code is set up in a Total Lagrangian formulation and
-rests on Green-Lagrange strains:
-
-\begin{align}
- \mathbf{E} = \frac{1}{2} \left( \text{Grad}\, \mathbf{U} + \text{Grad}^\text{T} \mathbf{U} + \text{Grad}^\text{T} \mathbf{U}\,\text{Grad}\, \mathbf{U}\right)
-\end{align}
-
-A rigid body rotation in 2D can be described by
-
-\begin{align}
- \mathbf{u} = [X_1 (\cos \vartheta - 1) - X_2 \sin \vartheta] \mathbf{E}_1 + [X_1 \sin \vartheta + X_2 (\cos \vartheta - 1)]\mathbf{E}_2
-\end{align}
-
-While this yields $\mathbf{E} = \mathbf{0}$ in finite deformation kinematics, we
-obtain a linearized strain tensor coordinate matrix with
-
-\begin{align}
- [\boldsymbol{\epsilon}]_{ij} = \left(
- \begin{array}{cc}
- \cos \vartheta - 1 & 0
- \\
- 0 & \cos \vartheta - 1
- \end{array}
- \right)
-\end{align}
-
-For both cases the OGS's linear elastic model is used to compute stresses. In
-the TL formulation this amounts to a Saint-Venant-Kirchhoff model.
-It thus suffices to illustrate the behaviour of strain values.
-
-```python jupyter={"source_hidden": true}
-# HIDDEN
-import numpy as np
-import matplotlib.pyplot as plt
-from matplotlib import cm
-import ogstools as ot
-import os
-import pyvista as pv
-
-# Some plot settings
-plt.style.use("seaborn-v0_8-deep")
-plt.rcParams["lines.linewidth"] = 2.0
-plt.rcParams["lines.color"] = "black"
-plt.rcParams["legend.frameon"] = True
-plt.rcParams["font.family"] = "serif"
-plt.rcParams["legend.fontsize"] = 14
-plt.rcParams["font.size"] = 14
-plt.rcParams["axes.spines.right"] = False
-plt.rcParams["axes.spines.top"] = False
-plt.rcParams["axes.spines.left"] = True
-plt.rcParams["axes.spines.bottom"] = True
-plt.rcParams["axes.axisbelow"] = True
-plt.rcParams["figure.figsize"] = (8, 6)
-```
-
-```python jupyter={"source_hidden": true}
-from pathlib import Path
-
-out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
-if not out_dir.exists():
- out_dir.mkdir(parents=True)
-```
-
-```python jupyter={"source_hidden": true}
-model_s = ot.Project(input_file="square_1e0.prj", output_file=f"{out_dir}/square_1e0_small.prj")
-model = ot.Project(input_file="square_1e0.prj", output_file="square_1e0.prj")
-```
-
-```python jupyter={"source_hidden": true}
-model_s.replace_text("SMALL_DEFORMATION", xpath="./processes/process/type")
-model_s.replace_text("StandardElasticityBrick", xpath="./processes/process/constitutive_relation/behaviour")
-model_s.replace_text("square_1e0_small", xpath="./time_loop/output/prefix")
-model_s.remove_element(xpath="./processes/process/secondary_variables/secondary_variable[@internal_name='deformation_gradient']")
-model_s.remove_element(xpath="./processes/process/secondary_variables/secondary_variable[@internal_name='volume_ratio']")
-model_s.remove_element(xpath=".//vtkdiff[field='deformation_gradient']")
-model_s.remove_element(xpath=".//vtkdiff[field='volume_ratio']")
-model_s.write_input()
-```
-
-```python jupyter={"source_hidden": true}
-model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")
-model_s.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
-```
-
-```python jupyter={"source_hidden": true}
-pv.set_plot_theme("document")
-pv.set_jupyter_backend("static") # comment out for interactive graphics
-```
-
-```python jupyter={"source_hidden": true}
-reader = pv.get_reader(f"{out_dir}/square_1e0.pvd")
-reader_s = pv.get_reader(f"{out_dir}/square_1e0_small.pvd")
-```
-
-```python jupyter={"source_hidden": true}
-reader.set_active_time_value(0.0)
-reader_s.set_active_time_value(0.0)
-mesh = reader.read()[0] # nulltes Gitter lesen
-mesh_s = reader_s.read()[0] # nulltes Gitter lesen
-```
-
-```python
-points = mesh.point_data["epsilon"].shape[0]
-```
-
-```python jupyter={"source_hidden": true}
-xs = mesh.points[:, 0]
-ys = mesh.points[:, 1]
-```
-
-```python jupyter={"source_hidden": true}
-def ploteps(time, angle, quantity="epsilon"):
- reader.set_active_time_value(time)
- reader_s.set_active_time_value(time)
- mesh = reader.read()[0] # nulltes Gitter lesen
- mesh_s = reader_s.read()[0] # nulltes Gitter lesen
- eps_vec = mesh.point_data[quantity][:, 0]
- eps_vec_s = mesh_s.point_data[quantity][:, 0]
-
- print(
- "Expected: %.2f (small strain) and 0 (large strain)"
- % (np.cos(np.deg2rad(angle)) - 1)
- )
-
- sargs = dict(
- title="small deformation, " + str(angle) + "°",
- title_font_size=15,
- label_font_size=15,
- n_labels=2,
- position_x=0.2,
- position_y=0.85,
- fmt="%.1f",
- width=0.6,
- )
-
- p = pv.Plotter(shape=(1, 2), border=False)
- p.subplot(0, 0)
- p.add_mesh(
- mesh,
- scalars=eps_vec_s,
- show_edges=False,
- show_scalar_bar=True,
- colormap="RdBu_r",
- scalar_bar_args=sargs,
- )
- # p.show_bounds(ticks="outside", xlabel="", ylabel="")
- # p.add_axes()
- p.view_xy()
- p.camera.zoom(1.2)
-
- sargs1 = dict(
- title="large deformation, " + str(angle) + "°",
- title_font_size=15,
- label_font_size=15,
- n_labels=2,
- position_x=0.2,
- position_y=0.85,
- fmt="%.1f",
- width=0.6,
- )
-
- p.subplot(0, 1)
- p.add_mesh(
- mesh,
- scalars=eps_vec,
- show_edges=False,
- show_scalar_bar=True,
- colormap="RdBu_r",
- scalar_bar_args=sargs1,
- )
- # p.show_bounds(ticks="outside", xlabel="", ylabel="")
- # p.add_axes()
- p.view_xy()
- p.camera.zoom(1.2)
-
- p.window_size = [800, 400]
- p.show();
-```
-
-We plot the normal strain in the $x$-direction for both kinematic formulations
-on the *undeformed* configuration as it undergoes a 360° rotation and find our
-expectations confirmed.
-
-```python
-ploteps(0, 0)
-```
-
-```python
-ploteps(0.5, 45)
-```
-
-```python
-ploteps(1, 90)
-```
-
-```python
-ploteps(2, 180)
-```
-
-```python
-ploteps(3, 270)
-```
-
-```python
-ploteps(4, 360)
-```
diff --git a/Tests/Data/LargeDeformation/RigidBody/RigidBody.py b/Tests/Data/LargeDeformation/RigidBody/RigidBody.py
new file mode 100644
index 0000000000000000000000000000000000000000..4bb1be69375ef5df38770c29114be64d53bc8626
--- /dev/null
+++ b/Tests/Data/LargeDeformation/RigidBody/RigidBody.py
@@ -0,0 +1,226 @@
+# %% [markdown]
+# +++
+# title = "LargeDeformations: Rigid Body Rotation"
+# date = "2023-03-01"
+# author = "Thomas Nagel"
+# web_subsection = "large-deformations"
+# +++
+
+# %% [markdown]
+# |<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
+# |---|---|--:|
+#
+# In this benchmark we test a basic kinematic feature of the finite strain / large
+# deformations implementation in OpenGeoSys.
+# An element is subjected to a rigid body rotation. The expected result is a
+# stress- and strain-free motion.
+# This is confirmed for the finite-strain implementation while the small-strain
+# implementation shows phantom strains and stresses.
+#
+# ## Basic definitions
+#
+# The small deformation code uses the linearized strain tensor:
+#
+# \begin{align}
+# \boldsymbol{\epsilon} = \frac{1}{2} \left( \text{grad}\, \mathbf{u} + \text{grad}\,^\text{T} \mathbf{u} \right)
+# \end{align}
+#
+# while the large deformation code is set up in a Total Lagrangian formulation and
+# rests on Green-Lagrange strains:
+#
+# \begin{align}
+# \mathbf{E} = \frac{1}{2} \left( \text{Grad}\, \mathbf{U} + \text{Grad}^\text{T} \mathbf{U} + \text{Grad}^\text{T} \mathbf{U}\,\text{Grad}\, \mathbf{U}\right)
+# \end{align}
+#
+# A rigid body rotation in 2D can be described by
+#
+# \begin{align}
+# \mathbf{u} = [X_1 (\cos \vartheta - 1) - X_2 \sin \vartheta] \mathbf{E}_1 + [X_1 \sin \vartheta + X_2 (\cos \vartheta - 1)]\mathbf{E}_2
+# \end{align}
+#
+# While this yields $\mathbf{E} = \mathbf{0}$ in finite deformation kinematics, we
+# obtain a linearized strain tensor coordinate matrix with
+#
+# \begin{align}
+# [\boldsymbol{\epsilon}]_{ij} = \left(
+# \begin{array}{cc}
+# \cos \vartheta - 1 & 0
+# \\
+# 0 & \cos \vartheta - 1
+# \end{array}
+# \right)
+# \end{align}
+#
+# For both cases the OGS's linear elastic model is used to compute stresses. In
+# the TL formulation this amounts to a Saint-Venant-Kirchhoff model.
+# It thus suffices to illustrate the behaviour of strain values.
+
+# %% jupyter={"source_hidden": true}
+import os
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import ogstools as ot
+import pyvista as pv
+
+# Some plot settings
+plt.style.use("seaborn-v0_8-deep")
+plt.rcParams["lines.linewidth"] = 2.0
+plt.rcParams["lines.color"] = "black"
+plt.rcParams["legend.frameon"] = True
+plt.rcParams["font.family"] = "serif"
+plt.rcParams["legend.fontsize"] = 14
+plt.rcParams["font.size"] = 14
+plt.rcParams["axes.spines.right"] = False
+plt.rcParams["axes.spines.top"] = False
+plt.rcParams["axes.spines.left"] = True
+plt.rcParams["axes.spines.bottom"] = True
+plt.rcParams["axes.axisbelow"] = True
+plt.rcParams["figure.figsize"] = (8, 6)
+
+# %% jupyter={"source_hidden": true}
+out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
+if not out_dir.exists():
+ out_dir.mkdir(parents=True)
+
+# %% jupyter={"source_hidden": true}
+model_s = ot.Project(
+ input_file="square_1e0.prj", output_file=f"{out_dir}/square_1e0_small.prj"
+)
+model = ot.Project(input_file="square_1e0.prj", output_file="square_1e0.prj")
+
+# %% jupyter={"source_hidden": true}
+model_s.replace_text("SMALL_DEFORMATION", xpath="./processes/process/type")
+model_s.replace_text(
+ "StandardElasticityBrick",
+ xpath="./processes/process/constitutive_relation/behaviour",
+)
+model_s.replace_text("square_1e0_small", xpath="./time_loop/output/prefix")
+model_s.remove_element(
+ xpath="./processes/process/secondary_variables/secondary_variable[@internal_name='deformation_gradient']"
+)
+model_s.remove_element(
+ xpath="./processes/process/secondary_variables/secondary_variable[@internal_name='volume_ratio']"
+)
+model_s.remove_element(xpath=".//vtkdiff[field='deformation_gradient']")
+model_s.remove_element(xpath=".//vtkdiff[field='volume_ratio']")
+model_s.write_input()
+
+# %% jupyter={"source_hidden": true}
+model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")
+model_s.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
+
+# %% jupyter={"source_hidden": true}
+pv.set_plot_theme("document")
+pv.set_jupyter_backend("static") # comment out for interactive graphics
+
+# %% jupyter={"source_hidden": true}
+reader = pv.get_reader(f"{out_dir}/square_1e0.pvd")
+reader_s = pv.get_reader(f"{out_dir}/square_1e0_small.pvd")
+
+# %% jupyter={"source_hidden": true}
+reader.set_active_time_value(0.0)
+reader_s.set_active_time_value(0.0)
+mesh = reader.read()[0] # nulltes Gitter lesen
+mesh_s = reader_s.read()[0] # nulltes Gitter lesen
+
+# %%
+points = mesh.point_data["epsilon"].shape[0]
+
+# %% jupyter={"source_hidden": true}
+xs = mesh.points[:, 0]
+ys = mesh.points[:, 1]
+
+
+# %% jupyter={"source_hidden": true}
+def ploteps(time, angle, quantity="epsilon"):
+ reader.set_active_time_value(time)
+ reader_s.set_active_time_value(time)
+ mesh = reader.read()[0] # nulltes Gitter lesen
+ mesh_s = reader_s.read()[0] # nulltes Gitter lesen
+ eps_vec = mesh.point_data[quantity][:, 0]
+ eps_vec_s = mesh_s.point_data[quantity][:, 0]
+
+ print(
+ "Expected: %.2f (small strain) and 0 (large strain)"
+ % (np.cos(np.deg2rad(angle)) - 1)
+ )
+
+ sargs = {
+ "title": "small deformation, " + str(angle) + "°",
+ "title_font_size": 15,
+ "label_font_size": 15,
+ "n_labels": 2,
+ "position_x": 0.2,
+ "position_y": 0.85,
+ "fmt": "%.1f",
+ "width": 0.6,
+ }
+
+ p = pv.Plotter(shape=(1, 2), border=False)
+ p.subplot(0, 0)
+ p.add_mesh(
+ mesh,
+ scalars=eps_vec_s,
+ show_edges=False,
+ show_scalar_bar=True,
+ colormap="RdBu_r",
+ scalar_bar_args=sargs,
+ )
+ # p.show_bounds(ticks="outside", xlabel="", ylabel="")
+ # p.add_axes()
+ p.view_xy()
+ p.camera.zoom(1.2)
+
+ sargs1 = {
+ "title": "large deformation, " + str(angle) + "°",
+ "title_font_size": 15,
+ "label_font_size": 15,
+ "n_labels": 2,
+ "position_x": 0.2,
+ "position_y": 0.85,
+ "fmt": "%.1f",
+ "width": 0.6,
+ }
+
+ p.subplot(0, 1)
+ p.add_mesh(
+ mesh,
+ scalars=eps_vec,
+ show_edges=False,
+ show_scalar_bar=True,
+ colormap="RdBu_r",
+ scalar_bar_args=sargs1,
+ )
+ # p.show_bounds(ticks="outside", xlabel="", ylabel="")
+ # p.add_axes()
+ p.view_xy()
+ p.camera.zoom(1.2)
+
+ p.window_size = [800, 400]
+ p.show()
+
+
+# %% [markdown]
+# We plot the normal strain in the $x$-direction for both kinematic formulations
+# on the *undeformed* configuration as it undergoes a 360° rotation and find our
+# expectations confirmed.
+
+# %%
+ploteps(0, 0)
+
+# %%
+ploteps(0.5, 45)
+
+# %%
+ploteps(1, 90)
+
+# %%
+ploteps(2, 180)
+
+# %%
+ploteps(3, 270)
+
+# %%
+ploteps(4, 360)
diff --git a/Tests/Data/LargeDeformation/Torsion/Torsion_robustness.md b/Tests/Data/LargeDeformation/Torsion/Torsion_robustness.py
similarity index 54%
rename from Tests/Data/LargeDeformation/Torsion/Torsion_robustness.md
rename to Tests/Data/LargeDeformation/Torsion/Torsion_robustness.py
index 7c26309ab727086b714fff191697294f6921b9d5..d6a3ac91924aad03881f68a208a442fb61a3d7f9 100644
--- a/Tests/Data/LargeDeformation/Torsion/Torsion_robustness.md
+++ b/Tests/Data/LargeDeformation/Torsion/Torsion_robustness.py
@@ -1,31 +1,33 @@
-+++
-title = "LargeDeformations: Torsion of a hyperelastic bar"
-date = "2023-03-02"
-author = "Thomas Nagel"
-web_subsection = "large-deformations"
-+++
-
-
-|<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
-|---|---|--:|
-
-## Convergence under large deformations (preliminary)
-
-In this benchmark we impose large tensile and torsional deformations on a
-hyperelastic prismatic bar.
-The material model is a Saint-Venant-Kirchhoff material law for initial testing.
+# %% [markdown]
+# +++
+# title = "LargeDeformations: Torsion of a hyperelastic bar"
+# date = "2023-03-02"
+# author = "Thomas Nagel"
+# web_subsection = "large-deformations"
+# +++
+
+# %% [markdown]
+# |<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
+# |---|---|--:|
+#
+# ## Convergence under large deformations (preliminary)
+#
+# In this benchmark we impose large tensile and torsional deformations on a
+# hyperelastic prismatic bar.
+# The material model is a Saint-Venant-Kirchhoff material law for initial testing.
+
+# %% jupyter={"source_hidden": true}
+import os
+from pathlib import Path
-```python jupyter={"source_hidden": true}
-# HIDDEN
-import numpy as np
import matplotlib.pyplot as plt
-from matplotlib import cm
+import numpy as np
import ogstools as ot
-from ogstools.logparser import parse_file, fill_ogs_context
-import os
-import pyvista as pv
import pandas as pd
+import pyvista as pv
+from ogstools.logparser import fill_ogs_context, parse_file
+# %% jupyter={"source_hidden": true}
pv.set_plot_theme("document")
pv.set_jupyter_backend("static")
@@ -42,16 +44,18 @@ plt.rcParams["axes.spines.left"] = True
plt.rcParams["axes.spines.bottom"] = True
plt.rcParams["axes.axisbelow"] = True
plt.rcParams["figure.figsize"] = (8, 6)
-```
-## Boundary conditions
+# %% [markdown]
+# ## Boundary conditions
+#
+# The bar of dimensions $1 \times 1 \times 6$ m³ is stretched by $\lambda = 1.5$
+# in the axial direction and twisted by 180°.
+# The following graph illustrates the torsion bc.
-The bar of dimensions $1 \times 1 \times 6$ m³ is stretched by $\lambda = 1.5$
-in the axial direction and twisted by 180°.
-The following graph illustrates the torsion bc.
-```python jupyter={"source_hidden": true}
-R = lambda x, y: np.sqrt(x**2 + y**2)
+# %% jupyter={"source_hidden": true}
+def R(x, y):
+ return np.sqrt(x**2 + y**2)
def phi(x, y):
@@ -63,9 +67,9 @@ def u(x, y):
ux = R(x, y) * np.cos(phi(x, y) + np.pi / 20) - x
uy = R(x, y) * np.sin(phi(x, y) + np.pi / 20) - y
return [ux, uy]
-```
-```python jupyter={"source_hidden": true}
+
+# %% jupyter={"source_hidden": true}
x = np.linspace(-0.5, 0.5, 10)
y = x.copy()
X, Y = np.meshgrid(x, y)
@@ -76,39 +80,29 @@ ax.set_aspect("equal")
ax.set_xlabel("$x$ / m")
ax.set_ylabel("$y$ / m")
ax.set_title("Displacement-controlled torsion on face of prismatic bar")
-fig.tight_layout();
-```
-
-```python jupyter={"source_hidden": true}
-import os
-from pathlib import Path
+fig.tight_layout()
+# %% jupyter={"source_hidden": true}
out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
out_dir.mkdir(parents=True)
-```
-```python jupyter={"source_hidden": true}
-model = ot.Project(
- input_file="bar1to6_torsion.prj",
- output_file="bar1to6_torsion.prj"
- )
-```
+# %% jupyter={"source_hidden": true}
+model = ot.Project(input_file="bar1to6_torsion.prj", output_file="bar1to6_torsion.prj")
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")
-```
-Let's plot the result.
+# %% [markdown]
+# Let's plot the result.
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
reader = pv.get_reader(f"{out_dir}/bar1to6_torsion.pvd")
print(reader.time_values)
reader.set_active_time_value(0.05)
mesh = reader.read()[0]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
plotter = pv.Plotter(shape=(1, 2), window_size=[1000, 500])
warped = mesh.warp_by_vector("displacement")
@@ -128,25 +122,24 @@ plotter.add_mesh(
plotter.view_isometric()
plotter.add_text("deformed", font_size=10)
plotter.show()
-```
-
-## Convergence
-What we're more interested in, is the convergence behaviour over time.
-The large deformations were applied in only 5 load steps to challenge the
-algorithm.
+# %% [markdown]
+# ## Convergence
+#
+# What we're more interested in, is the convergence behaviour over time.
+# The large deformations were applied in only 5 load steps to challenge the
+# algorithm.
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
log_file_raw = parse_file(f"{out_dir}/out.txt")
log_file_df = pd.DataFrame(log_file_raw)
-df = fill_ogs_context(log_file_df)
-```
+ogs_log_context = fill_ogs_context(log_file_df)
-```python jupyter={"source_hidden": true}
-last_ts = df["time_step"].values[-1]
+# %% jupyter={"source_hidden": true}
+last_ts = ogs_log_context["time_step"].to_numpy()[-1]
fig, ax = plt.subplots()
for i in range(1, last_ts + 1):
- timestep = df[df["time_step"] == i]
+ timestep = ogs_log_context[ogs_log_context["time_step"] == i]
ax.plot(
timestep["iteration_number"],
timestep["dx"],
@@ -161,10 +154,10 @@ ax.text(7.3, 2e-12, "1")
ax.text(6.7, 1e-9, "2")
ax.set_yscale("log")
ax.set_xlabel("iteration number")
-ax.set_ylabel("$|| \Delta \\vec{u}$ || / m")
+ax.set_ylabel("$|| \\Delta \\vec{u}$ || / m")
ax.legend()
-fig.tight_layout();
-```
+fig.tight_layout()
-We observe quadratic convergence in the proximity of the solution supporting the
-implementation of the geometric stiffness matrix.
+# %% [markdown]
+# We observe quadratic convergence in the proximity of the solution supporting the
+# implementation of the geometric stiffness matrix.
diff --git a/Tests/Data/Mechanics/HoekBrown/HoekBrownYieldCriterion.md b/Tests/Data/Mechanics/HoekBrown/HoekBrownYieldCriterion.py
similarity index 58%
rename from Tests/Data/Mechanics/HoekBrown/HoekBrownYieldCriterion.md
rename to Tests/Data/Mechanics/HoekBrown/HoekBrownYieldCriterion.py
index 0b50c692aefeb751dc6e8ef5f11220871823ced8..aa62f58c0429bd3525d7668f08e7a3aeeda672e9 100644
--- a/Tests/Data/Mechanics/HoekBrown/HoekBrownYieldCriterion.md
+++ b/Tests/Data/Mechanics/HoekBrown/HoekBrownYieldCriterion.py
@@ -1,119 +1,119 @@
-+++
-title = "Implementation of Hoek-Brown Failure Criterion in MFront"
-date = "2023-11-15"
-author = "Mehran Ghasabeh, Dmitri Naumov, Thomas Nagel"
-web_subsection = "small-deformations"
-+++
-
-
-|<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
-|---|---|--:|
-
-<!-- #region -->
-The generalized Hoek-Brown model, proposed in the work of Hoek and Brown in 1980 [1, 2], is derived based on the results of the brittle failure of intact rock [3] and the model studies of jointed rock mass behaviour [4].
-One of the remarkable priorities of the Hoek-Brown failure criterion over other criteria in the literature is that the strength and the deformation characteristics of heavily jointed rock masses can be evaluated regarding the geological strength index (GSI) as the information representing the field observation effects on the mechanical properties of a rock mass, such as the structure (or blockiness) and the condition of the joints [5].
-The GSI was introduced based on Bieniawski’s rock mass rating (RMR) system [6] and the Q-system [7].
-The GSI ranges from about $10$ for extremely poor rock masses to $100$ for intact rock.
-
-The generalized Hoek-Brown criterion for the estimation of rock mass strength is expressed as
-
-$$
- F_{\text{HB}} = \sigma_{1} - \sigma_{3} \ - \sigma_{\text{ci}} \ \left(s - m_b \cfrac{\sigma_{1}}{\sigma_{\text{ci}}}\right)^{a} = 0
- \tag{1}
-$$
-
-where $\sigma_{1} \geq \sigma_{2} \geq \sigma_{3}$ are the effective principal stresses, and the material constants $m_b$, $s$, and $a$ are determined as a function of the GSI:
-
-$$
- m_b = m_i \exp \left( \cfrac{\text{GSI}-100}{28-14D} \right), \quad s=\exp\left(\cfrac{\text{GSI}-100}{9-3D}\right) \quad \text{and} \quad a=\cfrac{1}{2}+\cfrac{1}{6}\left( \exp \left(-\cfrac{\text{GSI}}{15}\right) - \exp\left(-\cfrac{20}{3} \right)\right).
- \tag{2}
-$$
-
-The parameter $D$ is a factor that represents the degree of disturbance in the rock masses.
-The suggested value of the disturbance factor is $D=0$ for undisturbed in situ rock masses and $D=1$ for disturbed rock mass properties.
-
-In order to comprehensively analyze the behaviour of a slope, foundation, or tunnel, it is necessary to estimate not only the strength of intact rock and rock masses but also the deformation modulus of the rock mass in which these structures are excavated.
-The following equation is proposed for evaluating rock mass modulus based on the database of rock mass deformation modulus measurements [8]
-
-$$
- E_{\text{rm}} = E_i \left(0.02 + \cfrac{1- D / 2}{1 + \exp\left(\left(60+15D-\text{GSI}\right)/11\right)} \right) \quad (\text{MPa}),
- \tag{3}
-$$
-
-where $E_i$ denotes the intact rock deformation modulus (MPa).
-However, when the laboratory-measured values for $E_i$ are not available, the following alternative expression is used to calculate the rock mass modulus $E_{\text{rm}}$ [8]
-
-$$
- E_{\text{rm}} = \left(\cfrac{1- D / 2}{1 + \exp\left(\left(75+25D-\text{GSI}\right)/11\right)} \right) \quad (\text{MPa}).
- \tag{4}
-$$
-
-The generalized Hoek-Brown criterion is expressed in terms of stress invariants ($I_1$ and $J_2$) is given by
-
-$$
- F_{\text{HB}}= - m_b \ p \ \sigma_{\text{ci}} ^ {1/a -1} + (2 \ \sqrt {J_2} \ \cos \theta)^{1/a} + m_b \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1} \left( \cos \theta - \cfrac{\sin \theta}{\sqrt{3}} \right) - s \ \sigma_{\text{ci}}^{1/a} = 0 \quad \text{with} \quad p = -\dfrac{1}{3} \ I_1.
- \tag{5}
-$$
-
-However, this model presents computational challenges as a result of gradient discontinuities that show up at both the edges and apex of the hexagonal yield surface pyramid.
-These singularities frequently result in inefficient or failed stress integration schemes.
-In the current work, a straightforward hyperbolic yield surface is introduced, eliminating the singular apex in the meridian plane.
-Furthermore, a modified yield surface is developed based on a hyperbolic approximation with the octahedral rounding technique proposed in the work of Sloan and Booker [9], which is further adapted to the Mohr-Coulomb yield surface in the work of Abbo et al. [10], and Nagel et al. [11].
-To this, following the $C_2$ continuous smoothing approach explained in [12], the modified expression of the Hoek-Brown yield function is proposed as follows
-
-$$
- F_{\text{HB}}^{C_2}
- = - m_b \ p \ \sigma_{\text{ci}} ^ {1/a -1}
- + m_b \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1}
- \left( A + B \sin{3 \theta} + C \sin^2{3 \theta}\right)- s \ \sigma_{\text{ci}}^{1/a} = 0.
- \tag{6}
-$$
-
-This modified yield function is applied when the Lode angle is greater than or equal to the transition angle ($\theta \ge \theta_{T}$).
-The modified yield function is continuous and differentiable with respect to the stresses if the following conditions hold.
-These conditions are used to derive the coefficients $A$, $B$ and $C$ in (6):
-
-i) At a transition point, the $J_2$ invariant for the modified surface is identical to the $J_2$ invariant for the Mohr-Coulomb surface.
-
-ii) At a transition point, the derivative of the $J_2$ invariant with respect to the Lode angle $\theta$ for the modified surface is identical to one for the generalized Hoek-Brown surface.
-
-iii) The modified surface must be convex and definite, so at $\theta = \pm \pi / 6$, the derivative of the $J_2$ invariant with respect to the Lode angle is $0$.
-
-Moreover, the singularities and the corresponding difficulties of the derivation at the apex of the yield surface are resolved by adopting a quasi-hyperbolic approximation [11].
-To this end, the smoothing approach is applied by permuting $J_2$ with a small term $\epsilon$ according to $\hat{J}_2 = \sqrt{J_2 + \epsilon^2}$ in the yield surface function.
-The permutation parameter $\epsilon$ is determined by the following rule at $\theta=0$
-
-$$
- \epsilon = \text{min} \left( \delta, \mu \sqrt{J_2} \ | \ (2 \ \sqrt {J_2} \ )^{1/a} + m_b \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1} - s \ \sigma_{\text{ci}}^{1/a} = 0 \right)
- \quad \text{with} \quad \delta = 10^{-6} \quad \text{and} \quad \mu = 10^{-1}.
- \tag{7}
-$$
-
-For the nonassociated material behaviour, a plastic potential with an analogous structure to the yield function but different dilation-controlling parameters is defined by [13]
-
-$$
- G_{\text{HB}}=
- - m_g \ p \ \sigma_{\text{ci}} ^ {1/a_g -1} + (2 \ \sqrt {J_2} \ \cos \theta)^{1/a_g} + m_g \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1} \left( \cos \theta - \cfrac{\sin \theta}{\sqrt{3}} \right) - s_g \ \sigma_{\text{ci}}^{1/a_g} = 0,
- \tag{8}
-$$
-
-if $\theta \lt \pi/6$, and
-
-$$
- G_{\text{HB}}=
- - m_g \ p \ \sigma_{\text{ci}} ^ {1/a_g -1} + m_g \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a_g - 1} \left( A + B \sin{3 \theta} + C \sin^2{3 \theta}\right)- s_g \ \sigma_{\text{ci}}^{1/a_g} = 0,
- \tag{9}
-$$
-
-if $\theta \ge \pi/6$.
-
-According to the potential function, the rock mass dilation and its rate are mainly controlled by the curvature parameter $a_g$ and the dilation parameter $m_g$.
-
-The proposed formulation is then validated by capturing the $\pi$-plane for a set of material properties for an intact or massive rock mass with few widely spaced discontinuities with excellent quality control blasting or excavation by tunnel boring machine, i.e., $D=0$.
-Then the model is applied to an average-quality rock mass, as explained in detail in [14].
-<!-- #endregion -->
-
-```python
+# %% [markdown]
+# +++
+# title = "Implementation of Hoek-Brown Failure Criterion in MFront"
+# date = "2023-11-15"
+# author = "Mehran Ghasabeh, Dmitri Naumov, Thomas Nagel"
+# web_subsection = "small-deformations"
+# +++
+
+# %% [markdown]
+# |<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
+# |---|---|--:|
+
+# %% [markdown]
+# The generalized Hoek-Brown model, proposed in the work of Hoek and Brown in 1980 [1, 2], is derived based on the results of the brittle failure of intact rock [3] and the model studies of jointed rock mass behaviour [4].
+# One of the remarkable priorities of the Hoek-Brown failure criterion over other criteria in the literature is that the strength and the deformation characteristics of heavily jointed rock masses can be evaluated regarding the geological strength index (GSI) as the information representing the field observation effects on the mechanical properties of a rock mass, such as the structure (or blockiness) and the condition of the joints [5].
+# The GSI was introduced based on Bieniawski`s rock mass rating (RMR) system [6] and the Q-system [7].
+# The GSI ranges from about $10$ for extremely poor rock masses to $100$ for intact rock.
+#
+# The generalized Hoek-Brown criterion for the estimation of rock mass strength is expressed as
+#
+# $$
+# F_{\text{HB}} = \sigma_{1} - \sigma_{3} \ - \sigma_{\text{ci}} \ \left(s - m_b \cfrac{\sigma_{1}}{\sigma_{\text{ci}}}\right)^{a} = 0
+# \tag{1}
+# $$
+#
+# where $\sigma_{1} \geq \sigma_{2} \geq \sigma_{3}$ are the effective principal stresses, and the material constants $m_b$, $s$, and $a$ are determined as a function of the GSI:
+#
+# $$
+# m_b = m_i \exp \left( \cfrac{\text{GSI}-100}{28-14D} \right), \quad s=\exp\left(\cfrac{\text{GSI}-100}{9-3D}\right) \quad \text{and} \quad a=\cfrac{1}{2}+\cfrac{1}{6}\left( \exp \left(-\cfrac{\text{GSI}}{15}\right) - \exp\left(-\cfrac{20}{3} \right)\right).
+# \tag{2}
+# $$
+#
+# The parameter $D$ is a factor that represents the degree of disturbance in the rock masses.
+# The suggested value of the disturbance factor is $D=0$ for undisturbed in situ rock masses and $D=1$ for disturbed rock mass properties.
+#
+# In order to comprehensively analyze the behaviour of a slope, foundation, or tunnel, it is necessary to estimate not only the strength of intact rock and rock masses but also the deformation modulus of the rock mass in which these structures are excavated.
+# The following equation is proposed for evaluating rock mass modulus based on the database of rock mass deformation modulus measurements [8]
+#
+# $$
+# E_{\text{rm}} = E_i \left(0.02 + \cfrac{1- D / 2}{1 + \exp\left(\left(60+15D-\text{GSI}\right)/11\right)} \right) \quad (\text{MPa}),
+# \tag{3}
+# $$
+#
+# where $E_i$ denotes the intact rock deformation modulus (MPa).
+# However, when the laboratory-measured values for $E_i$ are not available, the following alternative expression is used to calculate the rock mass modulus $E_{\text{rm}}$ [8]
+#
+# $$
+# E_{\text{rm}} = \left(\cfrac{1- D / 2}{1 + \exp\left(\left(75+25D-\text{GSI}\right)/11\right)} \right) \quad (\text{MPa}).
+# \tag{4}
+# $$
+#
+# The generalized Hoek-Brown criterion is expressed in terms of stress invariants ($I_1$ and $J_2$) is given by
+#
+# $$
+# F_{\text{HB}}= - m_b \ p \ \sigma_{\text{ci}} ^ {1/a -1} + (2 \ \sqrt {J_2} \ \cos \theta)^{1/a} + m_b \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1} \left( \cos \theta - \cfrac{\sin \theta}{\sqrt{3}} \right) - s \ \sigma_{\text{ci}}^{1/a} = 0 \quad \text{with} \quad p = -\dfrac{1}{3} \ I_1.
+# \tag{5}
+# $$
+#
+# However, this model presents computational challenges as a result of gradient discontinuities that show up at both the edges and apex of the hexagonal yield surface pyramid.
+# These singularities frequently result in inefficient or failed stress integration schemes.
+# In the current work, a straightforward hyperbolic yield surface is introduced, eliminating the singular apex in the meridian plane.
+# Furthermore, a modified yield surface is developed based on a hyperbolic approximation with the octahedral rounding technique proposed in the work of Sloan and Booker [9], which is further adapted to the Mohr-Coulomb yield surface in the work of Abbo et al. [10], and Nagel et al. [11].
+# To this, following the $C_2$ continuous smoothing approach explained in [12], the modified expression of the Hoek-Brown yield function is proposed as follows
+#
+# $$
+# F_{\text{HB}}^{C_2}
+# = - m_b \ p \ \sigma_{\text{ci}} ^ {1/a -1}
+# + m_b \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1}
+# \left( A + B \sin{3 \theta} + C \sin^2{3 \theta}\right)- s \ \sigma_{\text{ci}}^{1/a} = 0.
+# \tag{6}
+# $$
+#
+# This modified yield function is applied when the Lode angle is greater than or equal to the transition angle ($\theta \ge \theta_{T}$).
+# The modified yield function is continuous and differentiable with respect to the stresses if the following conditions hold.
+# These conditions are used to derive the coefficients $A$, $B$ and $C$ in (6):
+#
+# i) At a transition point, the $J_2$ invariant for the modified surface is identical to the $J_2$ invariant for the Mohr-Coulomb surface.
+#
+# ii) At a transition point, the derivative of the $J_2$ invariant with respect to the Lode angle $\theta$ for the modified surface is identical to one for the generalized Hoek-Brown surface.
+#
+# iii) The modified surface must be convex and definite, so at $\theta = \pm \pi / 6$, the derivative of the $J_2$ invariant with respect to the Lode angle is $0$.
+#
+# Moreover, the singularities and the corresponding difficulties of the derivation at the apex of the yield surface are resolved by adopting a quasi-hyperbolic approximation [11].
+# To this end, the smoothing approach is applied by permuting $J_2$ with a small term $\epsilon$ according to $\hat{J}_2 = \sqrt{J_2 + \epsilon^2}$ in the yield surface function.
+# The permutation parameter $\epsilon$ is determined by the following rule at $\theta=0$
+#
+# $$
+# \epsilon = \text{min} \left( \delta, \mu \sqrt{J_2} \ | \ (2 \ \sqrt {J_2} \ )^{1/a} + m_b \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1} - s \ \sigma_{\text{ci}}^{1/a} = 0 \right)
+# \quad \text{with} \quad \delta = 10^{-6} \quad \text{and} \quad \mu = 10^{-1}.
+# \tag{7}
+# $$
+#
+# For the nonassociated material behaviour, a plastic potential with an analogous structure to the yield function but different dilation-controlling parameters is defined by [13]
+#
+# $$
+# G_{\text{HB}}=
+# - m_g \ p \ \sigma_{\text{ci}} ^ {1/a_g -1} + (2 \ \sqrt {J_2} \ \cos \theta)^{1/a_g} + m_g \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a - 1} \left( \cos \theta - \cfrac{\sin \theta}{\sqrt{3}} \right) - s_g \ \sigma_{\text{ci}}^{1/a_g} = 0,
+# \tag{8}
+# $$
+#
+# if $\theta \lt \pi/6$, and
+#
+# $$
+# G_{\text{HB}}=
+# - m_g \ p \ \sigma_{\text{ci}} ^ {1/a_g -1} + m_g \ \sqrt{J_2} \ \sigma_{\text{ci}}^{1/a_g - 1} \left( A + B \sin{3 \theta} + C \sin^2{3 \theta}\right)- s_g \ \sigma_{\text{ci}}^{1/a_g} = 0,
+# \tag{9}
+# $$
+#
+# if $\theta \ge \pi/6$.
+#
+# According to the potential function, the rock mass dilation and its rate are mainly controlled by the curvature parameter $a_g$ and the dilation parameter $m_g$.
+#
+# The proposed formulation is then validated by capturing the $\pi$-plane for a set of material properties for an intact or massive rock mass with few widely spaced discontinuities with excellent quality control blasting or excavation by tunnel boring machine, i.e., $D=0$.
+# Then the model is applied to an average-quality rock mass, as explained in detail in [14].
+
+# %%
import math as mt
import os
import shutil
@@ -123,16 +123,15 @@ from pathlib import Path
import matplotlib.pyplot as plt
import mtest
import numpy as np
+import ogstools as ot
import pandas as pd
import vtuIO
from IPython.display import Markdown, display
from matplotlib.ticker import MultipleLocator
-import ogstools as ot
from scipy.optimize import fsolve
from tfel.material import projectOnPiPlane
-```
-```python
+# %%
# Specify the folder names
calFolder_name = Path("calculatedData")
refFolder_name = Path("refData")
@@ -149,9 +148,8 @@ except Exception as e:
# Create the folder if it doesn't exist
if not calFolder_name.exists():
calFolder_name.mkdir(parents=True)
-```
-```python
+# %%
file_path = [
"calculated_piplane_data_Ex1.csv",
"calculated_piplane_data_Ex2.csv",
@@ -165,14 +163,12 @@ for file in map(Path, file_path):
if file.exists():
# Delete the file
file.unlink()
-```
-```python
-mfront_behaviour_lib = os.environ['VIRTUAL_ENV'] + '/../lib/libOgsMFrontBehaviour.so'
+# %%
+mfront_behaviour_lib = os.environ["VIRTUAL_ENV"] + "/../lib/libOgsMFrontBehaviour.so"
print(mfront_behaviour_lib)
-```
-```python
+# %%
plt.rcParams["lines.linewidth"] = 2.0
plt.rcParams["lines.color"] = "black"
plt.rcParams["legend.frameon"] = True
@@ -186,13 +182,14 @@ plt.rcParams["axes.spines.left"] = True
plt.rcParams["axes.spines.bottom"] = True
plt.rcParams["axes.axisbelow"] = True
mtest.setVerboseMode(mtest.VerboseLevel.VERBOSE_QUIET)
-```
-## Error Analysis Function
-```python
-def ErrorAnalysis(refFile, calFile, input1, input2, labelx, labely,
- tol):
+# %% [markdown]
+# ## Error Analysis Function
+
+
+# %%
+def ErrorAnalysis(refFile, calFile, input1, input2, labelx, labely, tol):
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(16, 5.5))
file_path1 = refFolder_name / refFile
file_path2 = calFolder_name / calFile
@@ -219,15 +216,15 @@ def ErrorAnalysis(refFile, calFile, input1, input2, labelx, labely,
mean_error = np.mean([mae0, mae1])
- if (abs(mae0) > 0.0):
+ if abs(mae0) > 0.0:
ax[0].set_ylim([errors0.min(), errors0.max()])
ax[0].scatter(range(len(errors0)), errors0, alpha=0.7, label=labelx)
- ax[0].axhline(y=0, color='red', linestyle='--', label='Zero Error Line')
+ ax[0].axhline(y=0, color="red", linestyle="--", label="Zero Error Line")
- if (abs(mae1) > 0.0):
+ if abs(mae1) > 0.0:
ax[1].set_ylim([errors1.min(), errors1.max()])
ax[1].scatter(range(len(errors1)), errors1, alpha=0.7, label=labely)
- ax[1].axhline(y=0, color='red', linestyle='--', label='Zero Error Line')
+ ax[1].axhline(y=0, color="red", linestyle="--", label="Zero Error Line")
# Customize the plot
ax[0].set_xlabel("points")
@@ -244,16 +241,18 @@ def ErrorAnalysis(refFile, calFile, input1, input2, labelx, labely,
# Check if the error exceeds the threshold
if mean_error > tolerance:
- msg = f"Error exceeds the threshold. Stopping the code. Mean Error: {mean_error}"
- raise ValueError(
- msg
+ msg = (
+ f"Error exceeds the threshold. Stopping the code. Mean Error: {mean_error}"
)
+ raise ValueError(msg)
print("Continue with the rest of the code.", "Mean Error:", mean_error)
-```
-## Verifications
-```python
+# %% [markdown]
+# ## Verifications
+
+
+# %%
def run_pi_plane(theta, material_set, data, calFile):
m_b = material_set["m_i"] * np.exp(
(material_set["gsi"] - 100.0) / (28.0 - 14.0 * material_set["Df"])
@@ -358,14 +357,13 @@ def run_pi_plane(theta, material_set, data, calFile):
pd.DataFrame(data).to_csv(calFolder_name / calFile, index=False)
return results_set
-```
-```python
+
+# %%
dataEx1 = {"sigma0": [], "sigma1": []}
calFile = "calculated_piplane_data_Ex1.csv"
-```
-```python
+# %%
piResults_set1 = {}
material_set = {}
fig, ax = plt.subplots(figsize=(8, 5))
@@ -399,9 +397,8 @@ ax.tick_params(width=2, length=5)
ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5))
ax.tick_params(which="minor", width=2, length=3, color="k")
-```
-```python
+# %%
ErrorAnalysis(
"reference_piplane_data_Ex1.csv",
"calculated_piplane_data_Ex1.csv",
@@ -411,14 +408,12 @@ ErrorAnalysis(
r"Errors of $\sigma_{yy}$",
tol=1e-6,
)
-```
-```python
+# %%
dataEx2 = {"sigma0": [], "sigma1": []}
calFile = "calculated_piplane_data_Ex2.csv"
-```
-```python
+# %%
piResults_set2 = {}
material_set = {}
fig, ax = plt.subplots(figsize=(8, 5))
@@ -451,9 +446,8 @@ ax.tick_params(width=2, length=5)
ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5))
ax.tick_params(which="minor", width=2, length=3, color="k")
-```
-```python
+# %%
ErrorAnalysis(
"reference_piplane_data_Ex2.csv",
"calculated_piplane_data_Ex2.csv",
@@ -463,50 +457,52 @@ ErrorAnalysis(
r"Errors of $\sigma_{yy}$",
tol=1e-6,
)
-```
-
-The yield function is approached along different stress paths within the $\pi$-plane.
-This shows that yield is correctly detected, and the stress state is correctly pulled back onto the yield surface.
-
-The tensile strength is calculated by setting $\sigma_1 = \sigma_2 = \sigma_3 = \sigma_t$ in (1) and calculated by
-$$
- \sigma_t = s \cfrac{\sigma_{ci}}{m_b}.
- \tag{10}
-$$
-## Calculation of Tensile Strength
-
-```python
+# %% [markdown]
+# The yield function is approached along different stress paths within the $\pi$-plane.
+# This shows that yield is correctly detected, and the stress state is correctly pulled back onto the yield surface.
+#
+# The tensile strength is calculated by setting $\sigma_1 = \sigma_2 = \sigma_3 = \sigma_t$ in (1) and calculated by
+# $$
+# \sigma_t = s \cfrac{\sigma_{ci}}{m_b}.
+# \tag{10}
+# $$
+#
+# ## Calculation of Tensile Strength
+
+# %%
sigma_t = round(
piResults_set2["sF"] / piResults_set2["m_b"] * (material_set["sigma_ci"]), 3
)
-stress_t = "$\\sigma_t = %2.3f$ MPa" % sigma_t
+stress_t = f"$\\sigma_t = {sigma_t:2.3f}$ MPa"
display(Markdown(stress_t))
-```
-
-The compressive strength is obtained by setting $\sigma_1 = 0$ in (1), then we can calculate it by
-$$
- \sigma_c = -\sigma_{ci} \ s^{a}.
- \tag{11}
-$$
-## Calculation of Compressive Strength
+# %% [markdown]
+# The compressive strength is obtained by setting $\sigma_1 = 0$ in (1), then we can calculate it by
+# $$
+# \sigma_c = -\sigma_{ci} \ s^{a}.
+# \tag{11}
+# $$
+#
+# ## Calculation of Compressive Strength
-```python
+# %%
sigma_c = round(
material_set["sigma_ci"] * piResults_set2["sF"] ** piResults_set2["alphF"], 3
)
-stress_c = "$\\sigma_c = %2.3f$ MPa" % sigma_c
+stress_c = f"$\\sigma_c = {sigma_c:2.3f}$ MPa"
display(Markdown(stress_c))
-```
-The validity of the Hoek-Brown model implemented in MFront is confirmed by the conduct of compressive and hydrostatic loading tests, which aim to ascertain the analytically derived compressive strength ($\sigma_c$) and tensile strength ($\sigma_t$).
-## Numerical Evaluation of Compressive and Hydrostatic Loading Tests
+# %% [markdown]
+# The validity of the Hoek-Brown model implemented in MFront is confirmed by the conduct of compressive and hydrostatic loading tests, which aim to ascertain the analytically derived compressive strength ($\sigma_c$) and tensile strength ($\sigma_t$).
+#
+# ## Numerical Evaluation of Compressive and Hydrostatic Loading Tests
+
-```python
+# %%
def runTest(material_set, data, calFile):
m_b = material_set["m_i"] * np.exp(
(material_set["gsi"] - 100.0) / (28.0 - 14.0 * material_set["Df"])
@@ -609,16 +605,16 @@ def runTest(material_set, data, calFile):
if material_set["Test"] == "Hydrostatic Loading Test":
pd.DataFrame(data).to_csv(calFolder_name / calFile, index=False)
return results_set
-```
-## Compressive Loading Test
-```python
+# %% [markdown]
+# ## Compressive Loading Test
+
+# %%
dataEx3 = {"epsilon2": [], "sigma2": []}
calFile = "calculated_CLT_data_Ex3.csv"
-```
-```python
+# %%
results_set1 = {}
material_set = {}
fig, ax = plt.subplots(figsize=(8, 6))
@@ -640,7 +636,10 @@ material_set = {
flgt = 1
results_set1 = runTest(material_set, dataEx3, calFile)
ax.plot(
- results_set1["eps_zz"] * 100.0, results_set1["sig_zz"], "-", label="Numerical Result"
+ results_set1["eps_zz"] * 100.0,
+ results_set1["sig_zz"],
+ "-",
+ label="Numerical Result",
)
ax.set_xlabel("$\\epsilon_{zz}$ / % ")
ax.set_ylabel("$\\sigma_{zz}$ / MPa")
@@ -651,12 +650,11 @@ ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5))
ax.tick_params(which="minor", width=2, length=3, color="k")
ax.axhline(
- y=-sigma_c, color="r", linestyle="--", label="$\\sigma_c = %2.3f$ MPa" % sigma_c
+ y=-sigma_c, color="r", linestyle="--", label=f"$\\sigma_c = {sigma_c:2.3f}$ MPa"
)
ax.legend()
-```
-```python
+# %%
ErrorAnalysis(
"reference_CLT_data_Ex3.csv",
"calculated_CLT_data_Ex3.csv",
@@ -666,14 +664,16 @@ ErrorAnalysis(
r"Errors of $\epsilon_{zz}$ / %",
tol=1e-6,
)
-```
-## Hydrostatic Loading Test
-This part validates the hydrostatic loading test by using analytically derived solutions.
-To this end, it is necessary to determine the bulk modulus and then the elastic limit, which indicates the point at which the material starts to undergo plastic deformation.
+# %% [markdown]
+# ## Hydrostatic Loading Test
+#
+# This part validates the hydrostatic loading test by using analytically derived solutions.
+# To this end, it is necessary to determine the bulk modulus and then the elastic limit, which indicates the point at which the material starts to undergo plastic deformation.
+
-```python
+# %%
# Calculate the analytical solution to the hydrostatic loading test
def analytical_solution(epst, sigma_t, bulkModulus):
def stressStrain(strain, bulkModulus):
@@ -693,34 +693,33 @@ def analytical_solution(epst, sigma_t, bulkModulus):
)
if strain.any() > epst:
ax.axhline(y=sigma_t, color="r", linestyle="--")
-```
-```python
+
+# %%
# Calculate the bulk modulus
def BulkModulus(results_set, material_set):
return results_set["Ec"] / (3.0 * (1.0 - 2.0 * material_set["nu"]))
-```
-```python
+
+# %%
# Calculate the pressure
def pressure(results_set):
return (
-(results_set["sig_xx"] + results_set["sig_yy"] + results_set["sig_zz"]) / 3.0
)
-```
-```python
+
+# %%
# Calculate the elastic limit
def elasticLimit(sigma_t, bulkModulus):
return sigma_t / bulkModulus * 100.0
-```
-```python
+
+# %%
dataEx4 = {"epsvol": [], "p": []}
calFile = "calculated_HLT_data_Ex4.csv"
-```
-```python
+# %%
results_set2 = {}
material_set = {}
fig, ax = plt.subplots(figsize=(8, 5))
@@ -753,14 +752,14 @@ epst = elasticLimit(sigma_t, bulkModulus)
analytical_solution(epst, sigma_t, bulkModulus)
ax.axhline(
- y=sigma_t, color="r", linestyle="--", label="$\\sigma_t = %2.3f$ MPa" % sigma_t
+ y=sigma_t, color="r", linestyle="--", label=f"$\\sigma_t = {sigma_t:2.3f}$ MPa"
)
ax.axvline(
x=epst,
color="g",
linestyle="--",
- label="$\\epsilon_t = %2.5f$ %%" % epst,
+ label=f"$\\epsilon_t = {epst:2.5f}$ %",
)
ax.plot(epst, sigma_t, "ko", label="Elastic Limit")
ax.set_xlabel("$\\epsilon_{vol}$ / %")
@@ -772,9 +771,8 @@ ax.xaxis.set_minor_locator(MultipleLocator(5))
ax.yaxis.set_minor_locator(MultipleLocator(5))
ax.tick_params(which="minor", width=2, length=3, color="k")
ax.legend()
-```
-```python
+# %%
ErrorAnalysis(
"reference_HLT_data_Ex4.csv",
"calculated_HLT_data_Ex4.csv",
@@ -784,97 +782,96 @@ ErrorAnalysis(
"Error of $\\epsilon_{\\text{vol}}$",
tol=1e-6,
)
-```
-
-## Elastoplastic Response of a Circular Hole to Cyclic Loading
-
-In this example, the Hoek-Brown model implemented based on the $C_2$ continuous approach is applied to a circular hole under cyclic loading and compare the results in terms of the radial and tangential stresses with the analytical solutions.
-The analytical solution for this problem is provided by the work of Carranza-Torres and Fairhurst [15].
-We consider a cylindrical hole with a radius of $r$, which is subjected to an in-situ uniform stress field of $\sigma_0$, and an internal pressure of $p_i$.
-The scaled far-field stresses, $S_0$, and scaled internal pressure, $P_0$, are determined by the following two equations, respectively
-
-$$
- S_0 = \frac{\sigma_0}{m_b \sigma_0} + \frac{s}{m_b ^ 2}, \quad \text{and} \quad P_i = \frac{P_0}{m_b \sigma_0} + \frac{s}{m_b ^ 2}.
- \tag{12}
-$$
-
-The elastic limit of rock mass is determined by the scaled far-field stress ($S_0$) and the corresponding scaled (dimensionless) critical internal pressure ($P^{\text{cr}}_i$) is calculated as
-
-$$
- P_{i}^{\text{cr}} = \frac{1}{16} \left(1 - \sqrt{1 + 16 S_0} \right)^2.
- \tag{13}
-$$
-Then, the critical internal pressure, $p^{\text{cr}}_i$ , is given by
-
-$$
- p_{i}^{\text{cr}} = \left( P_{i}^{\text{cr}}
- - \frac{s}{m_{b} ^ {2}} \right) m_b \sigma_{\text{ci}}.
- \tag{14}
-$$
-
-The plastic region grows equally around the hole. The failure zone reaches to
-
-$$
- b_{\text{pl}} = b \exp \left( 2 \sqrt{P_{i}^{\text{cr}}} - \sqrt{P_i} \right).
- \tag{15}
-$$
-
-In the plastic region where $r < b_{\text{pl}}$, the solutions for the radial and the tangential stresses are calculated as follows, respectively
-
-$$
- \sigma_r(r) = \left( S_r (r) - \frac{s}{m_b ^2} \right) m_b \sigma_{\text{ci}},
- \quad \text{and} \quad
- \sigma_{\theta}(r) = \left( S_{\theta} (r) - \frac{s}{m_b ^2} \right) m_b \sigma_{\text{ci}},
- \tag{16}
-$$
-
-where $S_r$ and $S_{\theta}$ are correspondingly given by
-
-$$
- S_r(r) = \left( \sqrt{P_{i}^{\text{cr}}} + \frac{1}{2} \ln \left( \frac{r}{b_{\text{pl}}} \right) \right)^2,
- \quad \text{and} \quad
- S_{\theta}(\theta) = \sqrt{S_{r}(r)} + S_r(r).
- \tag{17}
-$$
-
-The solutions for the radial and the tangential stresses in the elastic region are calculated by
-
-$$
- \sigma_r(r) = \sigma_0 - (\sigma_0 - p_{i}^{\text{cr}}) \left( \frac{b_{\text{pl}}}{r} \right)^2,
- \quad \text{and} \quad
- \sigma_{\theta}(r) = \sigma_0 + (\sigma_0 - p_{i}^{\text{cr}}) \left( \frac{b_{\text{pl}}}{r} \right)^2.
- \tag{18}
-$$
-
-```python
+# %% [markdown]
+# ## Elastoplastic Response of a Circular Hole to Cyclic Loading
+#
+# In this example, the Hoek-Brown model implemented based on the $C_2$ continuous approach is applied to a circular hole under cyclic loading and compare the results in terms of the radial and tangential stresses with the analytical solutions.
+# The analytical solution for this problem is provided by the work of Carranza-Torres and Fairhurst [15].
+# We consider a cylindrical hole with a radius of $r$, which is subjected to an in-situ uniform stress field of $\sigma_0$, and an internal pressure of $p_i$.
+# The scaled far-field stresses, $S_0$, and scaled internal pressure, $P_0$, are determined by the following two equations, respectively
+#
+# $$
+# S_0 = \frac{\sigma_0}{m_b \sigma_0} + \frac{s}{m_b ^ 2}, \quad \text{and} \quad P_i = \frac{P_0}{m_b \sigma_0} + \frac{s}{m_b ^ 2}.
+# \tag{12}
+# $$
+#
+# The elastic limit of rock mass is determined by the scaled far-field stress ($S_0$) and the corresponding scaled (dimensionless) critical internal pressure ($P^{\text{cr}}_i$) is calculated as
+#
+# $$
+# P_{i}^{\text{cr}} = \frac{1}{16} \left(1 - \sqrt{1 + 16 S_0} \right)^2.
+# \tag{13}
+# $$
+#
+# Then, the critical internal pressure, $p^{\text{cr}}_i$ , is given by
+#
+# $$
+# p_{i}^{\text{cr}} = \left( P_{i}^{\text{cr}}
+# - \frac{s}{m_{b} ^ {2}} \right) m_b \sigma_{\text{ci}}.
+# \tag{14}
+# $$
+#
+# The plastic region grows equally around the hole. The failure zone reaches to
+#
+# $$
+# b_{\text{pl}} = b \exp \left( 2 \sqrt{P_{i}^{\text{cr}}} - \sqrt{P_i} \right).
+# \tag{15}
+# $$
+#
+# In the plastic region where $r < b_{\text{pl}}$, the solutions for the radial and the tangential stresses are calculated as follows, respectively
+#
+# $$
+# \sigma_r(r) = \left( S_r (r) - \frac{s}{m_b ^2} \right) m_b \sigma_{\text{ci}},
+# \quad \text{and} \quad
+# \sigma_{\theta}(r) = \left( S_{\theta} (r) - \frac{s}{m_b ^2} \right) m_b \sigma_{\text{ci}},
+# \tag{16}
+# $$
+#
+# where $S_r$ and $S_{\theta}$ are correspondingly given by
+#
+# $$
+# S_r(r) = \left( \sqrt{P_{i}^{\text{cr}}} + \frac{1}{2} \ln \left( \frac{r}{b_{\text{pl}}} \right) \right)^2,
+# \quad \text{and} \quad
+# S_{\theta}(\theta) = \sqrt{S_{r}(r)} + S_r(r).
+# \tag{17}
+# $$
+#
+# The solutions for the radial and the tangential stresses in the elastic region are calculated by
+#
+# $$
+# \sigma_r(r) = \sigma_0 - (\sigma_0 - p_{i}^{\text{cr}}) \left( \frac{b_{\text{pl}}}{r} \right)^2,
+# \quad \text{and} \quad
+# \sigma_{\theta}(r) = \sigma_0 + (\sigma_0 - p_{i}^{\text{cr}}) \left( \frac{b_{\text{pl}}}{r} \right)^2.
+# \tag{18}
+# $$
+
+# %%
out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
out_dir.mkdir(parents=True)
-```
-```python
+# %%
model_hb = ot.Project(
- input_file="load_test_hb_nonassociated.prj", output_file="load_test_hb_nonassociated.prj"
- )
-```
+ input_file="load_test_hb_nonassociated.prj",
+ output_file="load_test_hb_nonassociated.prj",
+)
-```python
-model_hb.run_model(logfile=str(out_dir / "out.txt"), args=f"-o {out_dir}", write_prj_to_pvd=False)
-```
+# %%
+model_hb.run_model(
+ logfile=str(out_dir / "out.txt"), args=f"-o {out_dir}", write_prj_to_pvd=False
+)
-```python
+# %%
pvd_hb = vtuIO.PVDIO(f"{out_dir}/load_test_hb.pvd", dim=2)
-```
-```python
+# %%
b = 1000.0
rmax = 6000.0
npts = 120
raxis = [(i, 500, 0) for i in np.linspace(start=b, stop=rmax, num=npts)]
-```
-```python
+
+# %%
def analyticalSolution(material_set, data, analytFile):
S_0 = material_set["sigma_0"] / (
material_set["m_b"] * material_set["sigma_ci"]
@@ -934,9 +931,9 @@ def analyticalSolution(material_set, data, analytFile):
pd.DataFrame(data).to_csv(calFolder_name / analytFile, index=False)
return results_set
-```
-```python
+
+# %%
# Material properties
material_set_hb = {}
m_b = 0.55
@@ -958,16 +955,14 @@ material_set_hb = {
"rmax": rmax,
"npas": npas,
}
-```
-```python
+# %%
analytResults_set = {}
dataEx5 = {"rad": [], "sigma_rr": [], "sigma_tt": []}
analytFile = "analytical_data_Ex5.csv"
analytResults_set = analyticalSolution(material_set_hb, dataEx5, analytFile)
-```
-```python
+# %%
fig, ax = plt.subplots(figsize=(12, 6))
for i in [0, 10, 20, 30, 40]:
ax.plot(
@@ -977,23 +972,27 @@ for i in [0, 10, 20, 30, 40]:
ls="-",
)
-ax.axvline(x=analytResults_set["pl_range"], color="black", linestyle='--', label="$b_{\\text{pl}} = %2.3f$ m" % analytResults_set["pl_range"])
+ax.axvline(
+ x=analytResults_set["pl_range"],
+ color="black",
+ linestyle="--",
+ label="$b_{{\\text{{pl}}}} = {:2.3f}$ m".format(analytResults_set["pl_range"]),
+)
ax.set_xlabel("$r$ / m")
ax.set_ylabel("$\\epsilon_\\mathrm{p,eff}$")
ax.legend()
ax.grid("both")
fig.tight_layout()
-```
-```python
+# %%
fig, ax = plt.subplots(figsize=(14, 7))
radData = np.array(raxis).T[0] / 1000
calData_stress_rr = pvd_hb.read_set_data(i, "sigma", pointsetarray=raxis).T[0]
calData_stress_tt = pvd_hb.read_set_data(i, "sigma", pointsetarray=raxis).T[2]
pd.DataFrame(
- {"rad": radData, "sigma_rr": calData_stress_rr, "sigma_tt": calData_stress_tt}
- ).to_csv(calFolder_name/ "calculated_data_Ex5.csv", index=False)
+ {"rad": radData, "sigma_rr": calData_stress_rr, "sigma_tt": calData_stress_tt}
+).to_csv(calFolder_name / "calculated_data_Ex5.csv", index=False)
ax.plot(
np.array(raxis).T[0] / 1000,
@@ -1028,9 +1027,8 @@ ax.set_ylabel("$\\sigma$ / MPa")
ax.legend()
fig.tight_layout()
ax.grid("both")
-```
-```python
+# %%
ErrorAnalysis(
"reference_data_Ex5.csv",
"calculated_data_Ex5.csv",
@@ -1040,9 +1038,8 @@ ErrorAnalysis(
"Errors of $\\sigma_{\\theta \\theta}$",
tol=1e-6,
)
-```
-```python
+# %%
ErrorAnalysis(
"reference_data_Ex5.csv",
"analytical_data_Ex5.csv",
@@ -1052,38 +1049,38 @@ ErrorAnalysis(
"Errors of $\\sigma_{\\theta \\theta}$",
tol=5e-1,
)
-```
-
-## References
-
-[1] E. Hoek and E. T. Brown, “Empirical strength criterion for rock masses,†_Journal of the geotechnical engineering division_, vol. 106, no. 9, pp. 1013–1035, 1980.
-
-[2] E. T. Brown and E. Hoek, _Underground excavations in rock_. CRC Press, 1980.
-
-[3] E. Hoek, “Rock fracture under static stress conditions,†1965.
-
-[4] E. T. Brown, “Strength of models of rock with intermittent joints,†_Journal of the Soil Mechanics and Foundations Division_, vol. 96, no. 6, pp. 1935–1949, 1970.
-
-[5] E. Hoek and E. Brown, “The Hoek–brown failure criterion and GSI – 2018 edition,†_Journal of Rock Mechanics and Geotechnical Engineering_, vol. 11, no. 3, pp. 445–463, 2019.
-
-[6] Z. Bieniawski, “Classification of rock masses for engineering: the RMR system and future trends,†_Rock Testing and Site Characterization_, pp. 553–573, Elsevier, 1993.
-
-[7] N. Barton, “Some new Q-value correlations to assist in site characterisation and tunnel design,†_International Journal of Rock Mechanics and Mining Sciences_, vol. 39, no. 2, pp. 185–216, 2002.
-
-[8] E. Hoek, and Mark S Diederichs. "Empirical estimation of rock mass modulus." _International Journal of Rock Mechanics and Mining Sciences_ vol. 43, no. 2, pp. 203-215, 2006
-
-[9] S. Sloan and J. Booker, “Removal of singularities in Tresca and Mohr-Coulomb yield functions,†_Communications in Applied Numerical Methods_, vol. 2, no. 2, pp. 173–179, 1986.
-
-[10] A. Abbo, A. Lyamin, S. Sloan, and J. Hambleton, “A C2 continuous approximation to the Mohr–Coulomb yield surface,†_International Journal of Solids and Structures_, vol. 48, no. 21, pp. 3001–3010, 2011.
-
-[11] T. Nagel, W. Minkley, N. Böttcher, and D. Naumov, “Implicit numerical integration and consistent linearization of inelastic constitutive models of rock salt, _Computers Structures_, vol. 182, pp. 87-103, 2017.
-
-[12] R. Merifield, A. Lyamin, and S. Sloan, “Limit analysis solutions for the bearing capacity of rock masses using the generalised Hoek-Brown criterion,†_International Journal of Rock Mechanics and
-Mining Sciences_, vol. 43, no. 6, pp. 920–937, 2006
-
-[13] J. Clausen and L. Damkilde, “An exact implementation of the Hoek–Brown criterion for elasto-plastic finite element calculations,†_International Journal of Rock Mechanics and Mining Sciences_, vol. 45, no. 6, pp. 831–847, 2008.
-
-[14] E. Hoek, "Practical Rock Engineering. 2007." _Online. ed. Rocscience_ (<https://www.rocscience.com/assets/resources/learning/hoek/Practical-Rock-Engineering-Chapter-11-Rock-Mass-Properties.pdf>), 2007.
-[15] C. Carranza-Torres, C. Fairhurst, The elasto-plastic response of underground excavations in rock
-masses that satisfy the Hoek-Brown failure criterion, _International Journal of Rock Mechanics and Mining Sciences_, vol. 36, pp. 77-809, 1999
+# %% [markdown]
+# ## References
+#
+# [1] E. Hoek and E. T. Brown, “Empirical strength criterion for rock masses,†_Journal of the geotechnical engineering division_, vol. 106, no. 9, pp. 1013–1035, 1980.
+#
+# [2] E. T. Brown and E. Hoek, _Underground excavations in rock_. CRC Press, 1980.
+#
+# [3] E. Hoek, “Rock fracture under static stress conditions,†1965.
+#
+# [4] E. T. Brown, “Strength of models of rock with intermittent joints,†_Journal of the Soil Mechanics and Foundations Division_, vol. 96, no. 6, pp. 1935–1949, 1970.
+#
+# [5] E. Hoek and E. Brown, “The Hoek–brown failure criterion and GSI – 2018 edition,†_Journal of Rock Mechanics and Geotechnical Engineering_, vol. 11, no. 3, pp. 445–463, 2019.
+#
+# [6] Z. Bieniawski, “Classification of rock masses for engineering: the RMR system and future trends,†_Rock Testing and Site Characterization_, pp. 553–573, Elsevier, 1993.
+#
+# [7] N. Barton, “Some new Q-value correlations to assist in site characterisation and tunnel design,†_International Journal of Rock Mechanics and Mining Sciences_, vol. 39, no. 2, pp. 185–216, 2002.
+#
+# [8] E. Hoek, and Mark S Diederichs. "Empirical estimation of rock mass modulus." _International Journal of Rock Mechanics and Mining Sciences_ vol. 43, no. 2, pp. 203-215, 2006
+#
+# [9] S. Sloan and J. Booker, “Removal of singularities in Tresca and Mohr-Coulomb yield functions,†_Communications in Applied Numerical Methods_, vol. 2, no. 2, pp. 173–179, 1986.
+#
+# [10] A. Abbo, A. Lyamin, S. Sloan, and J. Hambleton, “A C2 continuous approximation to the Mohr–Coulomb yield surface,†_International Journal of Solids and Structures_, vol. 48, no. 21, pp. 3001–3010, 2011.
+#
+# [11] T. Nagel, W. Minkley, N. Böttcher, and D. Naumov, “Implicit numerical integration and consistent linearization of inelastic constitutive models of rock salt, _Computers Structures_, vol. 182, pp. 87-103, 2017.
+#
+# [12] R. Merifield, A. Lyamin, and S. Sloan, “Limit analysis solutions for the bearing capacity of rock masses using the generalised Hoek-Brown criterion,†_International Journal of Rock Mechanics and
+# Mining Sciences_, vol. 43, no. 6, pp. 920–937, 2006
+#
+# [13] J. Clausen and L. Damkilde, “An exact implementation of the Hoek–Brown criterion for elasto-plastic finite element calculations,†_International Journal of Rock Mechanics and Mining Sciences_, vol. 45, no. 6, pp. 831–847, 2008.
+#
+# [14] E. Hoek, "Practical Rock Engineering. 2007." _Online. ed. Rocscience_ (<https://www.rocscience.com/assets/resources/learning/hoek/Practical-Rock-Engineering-Chapter-11-Rock-Mass-Properties.pdf>), 2007.
+#
+# [15] C. Carranza-Torres, C. Fairhurst, The elasto-plastic response of underground excavations in rock
+# masses that satisfy the Hoek-Brown failure criterion, _International Journal of Rock Mechanics and Mining Sciences_, vol. 36, pp. 77-809, 1999
diff --git a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.py
similarity index 51%
rename from Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md
rename to Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.py
index cdb9b98da38b75e33f2c5fdbd8b729f6dd7b4baa..34eb0fdc62eae642392b38496a3f0c0a4f9461f1 100644
--- a/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.md
+++ b/Tests/Data/Mechanics/Linear/DiscWithHole/Linear_Disc_with_hole.py
@@ -1,52 +1,58 @@
-+++
-title = "Linear elasticity: disc with hole"
-date = "2022-04-27"
-author = "Linda Günther, Sophia Einspänner, Robert Habel, Christoph Lehmann and Thomas Nagel"
-web_subsection = "small-deformations"
-+++
-
-
-|<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
-|---|---|--:|
-
-In this benchmark example we consider a linear elastic small deformation problem. More specifically, a plate with a central hole that is put under tension on its top boundary is simulated. By exploiting symmetries, below we evaluate this problem just for the top right quarter of this disc.
-
-Doing this it is important to define boundary conditions for the regarded section of the plate. For the bottom boundary ($x$-axis, $\theta=-90°$) and the left boundary ($y$-axis, $\theta=0°$) we prescribe Dirichlet boundary conditions that constrain the normal displacement along the edge to be zero. On the top, where a tensile traction is applied to the plate, Neumann boundary conditions are prescribed.
-
-For the description of the plate and its load, some dimensions need to be defined. The quarter disc under consideration can be seen as a square with an edge length of $10\, \text{cm}$. The radius of the circular hole is $a = 2\, \text{cm}$ and the applied tension on the top boundary has a value of $\sigma = 10\, \text{kPa}$.
-
-To fully capture and understand the behaviour of stress and strain distributions around the hole, it is necessary to also define material properties. In case of isotropic linear elasticity, the relevant parameters are the Young's modulus and Poisson's ratio. The following parameters are chosen here:
-
-$ E = 1\,\text{MPa} \qquad \nu = 0.3$
-
-For verification of the numerical implementation, the numerical solution of the problem will be compared to the analytical solution.
-
-## Analytical solution
-
-The overall stress distributions in the plate around the hole can be represented by Kirsch's Solution, which is expressed here in cylindrical coordinates. The following equations are valid for an infinitely extended plate. Since the hole is very small compared to the dimension of the plate, we can consider this condition as fulfilled. The parameter $\sigma$ is the applied tension and $a$ is the radius of the inner hole. $r$ and $\theta$ are the cylindrical coordinates of a point in the plate. $r$ is measured from the origin (i.e., center of the hole), and $\theta$ is measured clockwise from the positive $y$-axis, i.e. the $x$-axis is at $\theta = -90°$.
-
-\begin{align}
- \sigma_{rr} \left( r , \theta \right) &=
- \frac{\sigma}{2}
- \left[ \left( 1 - \frac{a^2}{r^2} \right) + \left( 1 + 3 \frac{a^4}{r^4} - 4 \frac{a^2}{r^2} \right) \cos \left( 2 \theta \right) \right]
-\\
- \sigma_{\theta \theta} \left( r , \theta \right) &=
- \frac{\sigma}{2}
- \left[ \left( 1 + \frac{a^2}{r^2} \right) - \left( 1 + 3 \frac{a^4}{r^4} \right) \cos \left( 2 \theta \right) \right]
-\\
- \sigma_{r \theta} \left( r , \theta \right) &=
- - \frac{\sigma}{2}
- \left[ \left( 1 - 3 \frac{a^4}{r^4} + 2 \frac{a^2}{r^2} \right) \sin \left( 2 \theta \right) \right]
-\end{align}
-
-For visualization of the stress distribution within the plate the stresses $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$ are plotted along the $x$- and $y$-axes as well as along the diagonal ($\theta = -45°$).
-As can be seen below, the hole causes characteristic stress distributions.
+# %% [markdown]
+# +++
+# title = "Linear elasticity: disc with hole"
+# date = "2022-04-27"
+# author = "Linda Günther, Sophia Einspänner, Robert Habel, Christoph Lehmann and Thomas Nagel"
+# web_subsection = "small-deformations"
+# +++
+
+# %% [markdown]
+# |<div style="width:330px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:330px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://upload.wikimedia.org/wikipedia/commons/e/e8/TUBAF_Logo.svg" width="300"/></div>|
+# |---|---|--:|
+#
+# In this benchmark example we consider a linear elastic small deformation problem. More specifically, a plate with a central hole that is put under tension on its top boundary is simulated. By exploiting symmetries, below we evaluate this problem just for the top right quarter of this disc.
+#
+# Doing this it is important to define boundary conditions for the regarded section of the plate. For the bottom boundary ($x$-axis, $\theta=-90°$) and the left boundary ($y$-axis, $\theta=0°$) we prescribe Dirichlet boundary conditions that constrain the normal displacement along the edge to be zero. On the top, where a tensile traction is applied to the plate, Neumann boundary conditions are prescribed.
+#
+# For the description of the plate and its load, some dimensions need to be defined. The quarter disc under consideration can be seen as a square with an edge length of $10\, \text{cm}$. The radius of the circular hole is $a = 2\, \text{cm}$ and the applied tension on the top boundary has a value of $\sigma = 10\, \text{kPa}$.
+#
+# To fully capture and understand the behaviour of stress and strain distributions around the hole, it is necessary to also define material properties. In case of isotropic linear elasticity, the relevant parameters are the Young's modulus and Poisson's ratio. The following parameters are chosen here:
+#
+# $ E = 1\,\text{MPa} \qquad \nu = 0.3$
+#
+# For verification of the numerical implementation, the numerical solution of the problem will be compared to the analytical solution.
+#
+# ## Analytical solution
+#
+# The overall stress distributions in the plate around the hole can be represented by Kirsch's Solution, which is expressed here in cylindrical coordinates. The following equations are valid for an infinitely extended plate. Since the hole is very small compared to the dimension of the plate, we can consider this condition as fulfilled. The parameter $\sigma$ is the applied tension and $a$ is the radius of the inner hole. $r$ and $\theta$ are the cylindrical coordinates of a point in the plate. $r$ is measured from the origin (i.e., center of the hole), and $\theta$ is measured clockwise from the positive $y$-axis, i.e. the $x$-axis is at $\theta = -90°$.
+#
+# \begin{align}
+# \sigma_{rr} \left( r , \theta \right) &=
+# \frac{\sigma}{2}
+# \left[ \left( 1 - \frac{a^2}{r^2} \right) + \left( 1 + 3 \frac{a^4}{r^4} - 4 \frac{a^2}{r^2} \right) \cos \left( 2 \theta \right) \right]
+# \\
+# \sigma_{\theta \theta} \left( r , \theta \right) &=
+# \frac{\sigma}{2}
+# \left[ \left( 1 + \frac{a^2}{r^2} \right) - \left( 1 + 3 \frac{a^4}{r^4} \right) \cos \left( 2 \theta \right) \right]
+# \\
+# \sigma_{r \theta} \left( r , \theta \right) &=
+# - \frac{\sigma}{2}
+# \left[ \left( 1 - 3 \frac{a^4}{r^4} + 2 \frac{a^2}{r^2} \right) \sin \left( 2 \theta \right) \right]
+# \end{align}
+#
+# For visualization of the stress distribution within the plate the stresses $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$ are plotted along the $x$- and $y$-axes as well as along the diagonal ($\theta = -45°$).
+# As can be seen below, the hole causes characteristic stress distributions.
+
+# %% jupyter={"source_hidden": true}
+import os
+from pathlib import Path
-```python jupyter={"source_hidden": true}
-import numpy as np
import matplotlib.pyplot as plt
-from matplotlib import cm
+import numpy as np
+import ogstools as ot
+import pyvista as pv
+# %% jupyter={"source_hidden": true}
# Some plot settings
plt.style.use("seaborn-v0_8-deep")
plt.rcParams["lines.linewidth"] = 2.0
@@ -61,9 +67,9 @@ plt.rcParams["axes.spines.left"] = True
plt.rcParams["axes.spines.bottom"] = True
plt.rcParams["axes.axisbelow"] = True
plt.rcParams["figure.figsize"] = (8, 6)
-```
-```python jupyter={"source_hidden": true}
+
+# %% jupyter={"source_hidden": true}
def kirsch_sig_rr(sig, r, theta, a):
# Änderung in Heaviside-Funktion: Spannung erst ab r=1.99... auf Null setzen, da erster Eintrag von dist_sorted etwas kleiner als 2 ist
return (
@@ -101,15 +107,16 @@ def kirsch_sig_rt(sig, r, theta, a):
)
* np.heaviside(r + 1e-7 - a, 1)
)
-```
-### Stress distribution along the $x$-axis ($\theta= -90°$, orthogonal to the load)
-At larger distances from the hole the stresses are approximately distributed as if the plate had no hole ($\sigma_{rr} \equiv \sigma_{xx} \approx 0$, $\sigma_{r\theta} \equiv \sigma_{xy} = 0$, $\sigma_{\theta\theta} \equiv \sigma_{yy} = 10\,\text{kPa}$).
-As the hole is approached, the tangential stress increases until it reaches its maximum value directly at the contour.
-Interestingly, that value is three times as high as the applied traction. It can be concluded that the hole leads to a **threefold stress concentration**.
+# %% [markdown]
+# ### Stress distribution along the $x$-axis ($\theta= -90°$, orthogonal to the load)
+#
+# At larger distances from the hole the stresses are approximately distributed as if the plate had no hole ($\sigma_{rr} \equiv \sigma_{xx} \approx 0$, $\sigma_{r\theta} \equiv \sigma_{xy} = 0$, $\sigma_{\theta\theta} \equiv \sigma_{yy} = 10\,\text{kPa}$).
+# As the hole is approached, the tangential stress increases until it reaches its maximum value directly at the contour.
+# Interestingly, that value is three times as high as the applied traction. It can be concluded that the hole leads to a **threefold stress concentration**.
-```python
+# %%
r = np.linspace(2, 10, 1000)
plt.plot(r, kirsch_sig_rr(10, r, -90, 2), label=r"$\sigma_{rr} \equiv \sigma_{xx}$")
plt.plot(
@@ -128,11 +135,11 @@ plt.xlabel(r"$r\equiv x$ / cm")
plt.xlim(0, 10)
plt.ylabel(r"$\sigma$ / kPa")
plt.ylim(0, 30)
-```
-### Stress distribution along the diagonal ($\theta= -45°$)
+# %% [markdown]
+# ### Stress distribution along the diagonal ($\theta= -45°$)
-```python
+# %%
r = np.linspace(2, 14.14, 1000)
plt.plot(r, kirsch_sig_rr(10, r, -45, 2), label=r"$\sigma_{rr}$")
plt.plot(r, kirsch_sig_tt(10, r, -45, 2), label=r"$\sigma_{\theta \theta}$")
@@ -145,17 +152,17 @@ plt.xlabel(r"$r$ / cm")
plt.xlim(0, 14.14)
plt.ylabel(r"$\sigma$ / kPa")
plt.ylim(0, 10)
-```
-
-### Stress distribution along the $y$-axis ($\theta=0°$, parallel to the applied tension)
-
-**Attention**: Radius is plotted along the vertical axis. I.e., upwards direction in the plot corresponds to the upwards ($y$) direction of the problem.
-The horizontal stress near the hole amounts to the negative of the applied tension. Far away from the hole ($r\to\infty$) it vanishes.
+# %% [markdown]
+# ### Stress distribution along the $y$-axis ($\theta=0°$, parallel to the applied tension)
+#
+# **Attention**: Radius is plotted along the vertical axis. I.e., upwards direction in the plot corresponds to the upwards ($y$) direction of the problem.
+#
+# The horizontal stress near the hole amounts to the negative of the applied tension. Far away from the hole ($r\to\infty$) it vanishes.
+#
+# The vertical stress is zero at the edge of the hole and exhibits a sign change in the vicinity of the hole.
-The vertical stress is zero at the edge of the hole and exhibits a sign change in the vicinity of the hole.
-
-```python
+# %%
r = np.linspace(2, 10, 1000)
plt.plot(kirsch_sig_rr(10, r, 0, 2), r, label=r"$\sigma_{rr} \equiv \sigma_{yy}$")
plt.plot(
@@ -170,21 +177,22 @@ ax = plt.gca()
plt.ylabel(r"$r\equiv y$ / cm")
plt.xlabel(r"$\sigma$ / kPa")
plt.xlim(-10, 10)
-```
-### 2D color plots
+# %% [markdown]
+# ### 2D color plots
+#
+# For gaining a deeper insight for the planar stress distribution in the plate around the hole the following 2D plots are helpful. The magnitude of the acting stresses is expressed by a chromatic graduated scale in separate plots for $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$.
+
-For gaining a deeper insight for the planar stress distribution in the plate around the hole the following 2D plots are helpful. The magnitude of the acting stresses is expressed by a chromatic graduated scale in separate plots for $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$.
+# %%
+def cart_to_cyl(x, y):
+ return [np.sqrt(x**2 + y**2), np.rad2deg(np.arctan(x / y))]
-```python
-cart_to_cyl = lambda x, y: [np.sqrt(x**2 + y**2), np.rad2deg(np.arctan(x / y))]
-```
-```python
+# %%
X, Y = np.meshgrid(np.linspace(0.1, 10, 1000), np.linspace(0.1, 10, 1000))
-```
-```python
+# %%
# cylindical coordinates from Cartesian grid
r, theta = cart_to_cyl(X, Y)
eps = 1e-2
@@ -206,74 +214,64 @@ for i in range(3):
ax[i].set_aspect("equal")
ax[i].set_xlabel("$x$ / cm")
ax[i].set_ylabel("$y$ / cm")
-ax[0].set_title("$\\sigma_{rr}$ / kPa")
-ax[1].set_title("$\\sigma_{\\theta\\theta}$ / kPa")
-ax[2].set_title("$\\sigma_{r\\theta}$ / kPa")
+ax[0].set_title(r"$\sigma_{rr}$ / kPa")
+ax[1].set_title(r"$\sigma_{\theta\theta}$ / kPa")
+ax[2].set_title(r"$\sigma_{r\theta}$ / kPa")
fig.tight_layout()
-```
-
-## Numerical solution
-
-In the following section the previously analytically considered problem of the plate with a central hole is solved with numerical methods (FEM in OpenGeoSys).
-OpenGeoSys performs computations in Cartesian coordinates.
-In order to control, display and compare the results of the numerical solution with the analytical it is necessary to convert the output values from Cartesian coordinates to polar coordinates.
-The formulae below are used to obtain the coordinates $r$ and $\theta$ of the polar coordinate system.
-
-$$
-\begin{align}
- r &= \sqrt{\left(x^2+y^2\right)}
-\\
- \theta &= \arctan\left(\frac yx \right)
-\\
-\end{align}
-$$
-
-For the simulation of the numerical model, some assumptions have to be made:
-
-* The analytical solution is strictly correct only for an infinitely extended plate.
- This requirement can only be approximated in the numerical method. The dimensions of the considered plate must be much larger than the hole. With a hole diameter of $4\, \text{cm}$ and a size of $20\times20\, \text{cm}$ the plate can be considered as infinitely extended. Indeed, a general rule implied that the range of borders should be about five times the size of the cavity.
-* For numerical simulations, it is necessary to specify the mesh resolution with which the model is calculated. A finer mesh provides more detailed gradients but takes more time for calculation, whereas a coarser mesh yields less accurate values, but in a shorter time.
-* Another assumption that enables us to work with a two-dimensional plate model is the plain strain condition. That means that deformations only occur in one plane.
-
-Simultaneous to the analytical solution the results of the numerical model are plotted along the $x$-, $y$- and diagonal axis. To visualize local differences to the analytically computed stress distribution of $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$ the latter are also added in the plots as dotted lines.<br>
-For a better presentation of these differences and as a basis for the following error analysis, plots for absolute and relative error along the regarded axes are have been made.
-
-```python
-import ogstools as ot
-import os
-from pathlib import Path
+# %% [markdown]
+# ## Numerical solution
+#
+# In the following section the previously analytically considered problem of the plate with a central hole is solved with numerical methods (FEM in OpenGeoSys).
+# OpenGeoSys performs computations in Cartesian coordinates.
+# In order to control, display and compare the results of the numerical solution with the analytical it is necessary to convert the output values from Cartesian coordinates to polar coordinates.
+# The formulae below are used to obtain the coordinates $r$ and $\theta$ of the polar coordinate system.
+#
+# $$
+# \begin{align}
+# r &= \sqrt{\left(x^2+y^2\right)}
+# \\
+# \theta &= \arctan\left(\frac yx \right)
+# \\
+# \end{align}
+# $$
+#
+# For the simulation of the numerical model, some assumptions have to be made:
+#
+# * The analytical solution is strictly correct only for an infinitely extended plate.
+# This requirement can only be approximated in the numerical method. The dimensions of the considered plate must be much larger than the hole. With a hole diameter of $4\, \text{cm}$ and a size of $20\times20\, \text{cm}$ the plate can be considered as infinitely extended. Indeed, a general rule implied that the range of borders should be about five times the size of the cavity.
+# * For numerical simulations, it is necessary to specify the mesh resolution with which the model is calculated. A finer mesh provides more detailed gradients but takes more time for calculation, whereas a coarser mesh yields less accurate values, but in a shorter time.
+# * Another assumption that enables us to work with a two-dimensional plate model is the plain strain condition. That means that deformations only occur in one plane.
+#
+# Simultaneous to the analytical solution the results of the numerical model are plotted along the $x$-, $y$- and diagonal axis. To visualize local differences to the analytically computed stress distribution of $\sigma_{rr}$, $\sigma_{\theta\theta}$ and $\sigma_{r\theta}$ the latter are also added in the plots as dotted lines.<br>
+# For a better presentation of these differences and as a basis for the following error analysis, plots for absolute and relative error along the regarded axes are have been made.
+
+# %%
out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
if not out_dir.exists():
out_dir.mkdir(parents=True)
-```
-```python
+# %%
model = ot.Project(
input_file="../disc_with_hole.prj", output_file="../disc_with_hole.prj"
)
-```
-```python
+# %%
model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir}")
-```
-
-```python
-import pyvista as pv
+# %%
pv.set_plot_theme("document")
pv.set_jupyter_backend("static") # comment out for interactive graphics
-```
-```python
+# %%
reader = pv.get_reader(f"{out_dir}/disc_with_hole.pvd")
reader.set_active_time_value(1.0) # go to 1 s
mesh = reader.read()[0] # nulltes Gitter lesen
-```
-```python
+
+# %%
def vec4_to_mat3x3cart(vec4):
- theta = np.arctan2(ys, xs)
+ # theta = np.arctan2(ys, xs)
m = np.zeros((3, 3))
m[0, 0] = vec4[0]
@@ -299,11 +297,12 @@ def vec4_to_mat3x3polar(
rot[1, 1] = np.cos(theta) # rot = Drehmatrix, Drehung um z-Achse
return rot.T * m_cart * rot
-```
-### Stress distribution along the x-axis ($\theta = -90°$)
-```python
+# %% [markdown]
+# ### Stress distribution along the x-axis ($\theta = -90°$)
+
+# %%
pt1 = (0, 1e-6, 0)
pt2 = (10, 1e-6, 0)
xaxis = pv.Line(pt1, pt2, resolution=2)
@@ -316,13 +315,11 @@ indices_sorted = np.argsort(
dist_from_origin
) # Indizes der Punkte sortieren nach ihrem Abstand vom Ursprung
dist_sorted = dist_from_origin[indices_sorted] # index magic
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig = line_mesh.point_data["sigma"]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
num_points = sig.shape[0] # entspricht len(sig) d.h Anzahl der Punkte, hier 36
sig_rr = np.zeros(num_points) # Liste von 36 Nullen
sig_tt = np.zeros(num_points)
@@ -333,9 +330,8 @@ f_abs_rt = np.zeros(num_points)
f_rel_rr = np.zeros(num_points)
f_rel_tt = np.zeros(num_points)
f_rel_rt = np.zeros(num_points)
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
for pt_idx in range(num_points):
sig_vec = sig[pt_idx, :]
xs = line_mesh.points[pt_idx, 0]
@@ -344,15 +340,13 @@ for pt_idx in range(num_points):
sig_rr[pt_idx] = sig_polar[0, 0]
sig_tt[pt_idx] = sig_polar[1, 1]
sig_rt[pt_idx] = sig_polar[0, 1]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig_rr_sorted = sig_rr[indices_sorted]
sig_tt_sorted = sig_tt[indices_sorted]
sig_rt_sorted = sig_rt[indices_sorted]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
for pt_idx in range(num_points):
f_abs_rr[pt_idx] = sig_rr_sorted[pt_idx] * 1000 - kirsch_sig_rr(
10, dist_sorted[pt_idx], -90, 2
@@ -381,9 +375,8 @@ for pt_idx in range(num_points):
f_rel_rt[pt_idx] = f_abs_rt[pt_idx] / kirsch_sig_rt(
10, dist_sorted[pt_idx], -90, 2
)
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
r = np.linspace(2, 10, 1000)
fig, ax = plt.subplots(ncols=3, figsize=(18, 6))
@@ -398,56 +391,56 @@ ax[0].plot(
kirsch_sig_rr(10, r, -90, 2),
color="deepskyblue",
linestyle=":",
- label="$\sigma_{rr,\mathrm{analytical}}$",
+ label=r"$\sigma_{rr,\mathrm{analytical}}$",
)
ax[0].plot(
r,
kirsch_sig_tt(10, r, -90, 2),
color="yellowgreen",
linestyle=":",
- label="$\sigma_{\\theta\\theta,\mathrm{analytical}}$",
+ label=r"$\sigma_{\theta\theta,\mathrm{analytical}}$",
)
ax[0].plot(
r,
kirsch_sig_rt(10, r, -90, 2),
color="orangered",
linestyle=":",
- label="$\sigma_{r\\theta,\mathrm{analytical}}$",
+ label=r"$\sigma_{r\theta,\mathrm{analytical}}$",
)
-ax[0].plot(dist_sorted, sig_rr_sorted * 1000, label="$\sigma_{rr}$")
-ax[0].plot(dist_sorted, sig_tt_sorted * 1000, label="$\sigma_{\\theta\\theta}$")
-ax[0].plot(dist_sorted, sig_rt_sorted * 1000, label="$\sigma_{r\\theta}$")
+ax[0].plot(dist_sorted, sig_rr_sorted * 1000, label=r"$\sigma_{rr}$")
+ax[0].plot(dist_sorted, sig_tt_sorted * 1000, label=r"$\sigma_{\theta\theta}$")
+ax[0].plot(dist_sorted, sig_rt_sorted * 1000, label=r"$\sigma_{r\theta}$")
-ax[0].set_ylabel("$\\sigma$ / kPa")
+ax[0].set_ylabel(r"$\sigma$ / kPa")
ax[0].set(ylim=(-2, 35.5))
ax[0].legend(loc="upper right")
-ax[1].plot(dist_sorted, f_abs_rr, label="$\Delta\sigma_{rr}$")
-ax[1].plot(dist_sorted, f_abs_tt, label="$\Delta\sigma_{\\theta\\theta}$")
-ax[1].plot(dist_sorted, f_abs_rt, label="$\Delta\sigma_{r\\theta}$")
+ax[1].plot(dist_sorted, f_abs_rr, label=r"$\Delta\sigma_{rr}$")
+ax[1].plot(dist_sorted, f_abs_tt, label=r"$\Delta\sigma_{\theta\theta}$")
+ax[1].plot(dist_sorted, f_abs_rt, label=r"$\Delta\sigma_{r\theta}$")
# ax[1].spines['bottom'].set_position('zero')
ax[1].spines["top"].set_color("none")
ax[1].xaxis.tick_bottom()
-ax[1].set_ylabel("$\Delta\\sigma$ / kPa")
+ax[1].set_ylabel(r"$\Delta\sigma$ / kPa")
ax[1].set(ylim=(-2, 4.5))
ax[1].legend(loc="upper right")
ax[2].plot(
dist_sorted,
f_rel_rr,
- label="$\Delta\sigma_{rr}$ / $\sigma_{rr,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{rr}$ / $\sigma_{rr,\mathrm{analytical}}$",
)
ax[2].plot(
dist_sorted,
f_rel_tt,
- label="$\Delta\sigma_{\\theta\\theta}$ / $\sigma_{\\theta\\theta,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{\theta\theta}$ / $\sigma_{\theta\theta,\mathrm{analytical}}$",
)
-# ax[2].plot(dist_sorted, f_rel_rt, label = "$\Delta\sigma_{r\\theta}$ / $\sigma_{r\\theta,analytical}$")
+# ax[2].plot(dist_sorted, f_rel_rt, label = r"$\Delta\sigma_{r\theta}$ / $\sigma_{r\theta,analytical}$")
ax[2].set_ylim(-0.2, 0.2)
-ax[2].set_ylabel("$\Delta\\sigma$ / $\sigma_{\mathrm{analytical}}$")
+ax[2].set_ylabel(r"$\Delta\sigma$ / $\sigma_{\mathrm{analytical}}$")
ax[2].legend() # loc="upper right")
ax[0].set_title("Stress distribution")
@@ -455,11 +448,11 @@ ax[1].set_title("Absolute error")
ax[2].set_title("Relative error")
fig.tight_layout()
-```
-### Stress distribution along the diagonal ($\theta = 45°$)
+# %% [markdown]
+# ### Stress distribution along the diagonal ($\theta = 45°$)
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
pt1 = (0, 0, 0)
pt2 = (10, 10, 0)
diagonal = pv.Line(
@@ -476,16 +469,14 @@ indices_sorted = np.argsort(
dist_from_origin
) # Indizes der Punkte sortieren nach ihrem Abstand vom Ursprung
dist_sorted = dist_from_origin[indices_sorted] # index magic
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig = line_mesh.point_data[
"sigma"
] # Spannung in Kelvin Darstellung als 4 zeiliger Vector statt 3x3 Matrix
# Achtung (CL): Das ist nicht ganz die Kelvin-Darstellung. Die Kelvin-Darstellung hat noch Faktoren sqrt(2) bei den Nichtdiagonalelementen.
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
num_points = sig.shape[0] # entspricht len(sig) d.h Anzahl der Punkte, hier 36
sig_rr = np.zeros(num_points) # Liste von 36 Nullen
sig_tt = np.zeros(num_points)
@@ -496,9 +487,8 @@ f_abs_rt = np.zeros(num_points)
f_rel_rr = np.zeros(num_points)
f_rel_tt = np.zeros(num_points)
f_rel_rt = np.zeros(num_points)
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
for pt_idx in range(num_points):
sig_vec = sig[pt_idx, :]
xs = line_mesh.points[pt_idx, 0]
@@ -507,15 +497,13 @@ for pt_idx in range(num_points):
sig_rr[pt_idx] = sig_polar[0, 0]
sig_tt[pt_idx] = sig_polar[1, 1]
sig_rt[pt_idx] = sig_polar[0, 1]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig_rr_sorted = sig_rr[indices_sorted]
sig_tt_sorted = sig_tt[indices_sorted]
sig_rt_sorted = sig_rt[indices_sorted]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
for pt_idx in range(num_points):
f_abs_rr[pt_idx] = sig_rr_sorted[pt_idx] * 1000 - kirsch_sig_rr(
10, dist_sorted[pt_idx], -45, 2
@@ -544,9 +532,8 @@ for pt_idx in range(num_points):
f_rel_rt[pt_idx] = f_abs_rt[pt_idx] / kirsch_sig_rt(
10, dist_sorted[pt_idx], -45, 2
)
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
r = np.linspace(2, 14.14, 1000)
fig, ax = plt.subplots(ncols=3, figsize=(18, 6))
@@ -562,59 +549,59 @@ ax[0].plot(
kirsch_sig_rr(10, r, -45, 2),
color="deepskyblue",
linestyle=":",
- label="$\sigma_{rr,\mathrm{analytical}}$",
+ label=r"$\sigma_{rr,\mathrm{analytical}}$",
)
ax[0].plot(
r,
kirsch_sig_tt(10, r, -45, 2),
color="yellowgreen",
linestyle=":",
- label="$\sigma_{\\theta\\theta,\mathrm{analytical}}$",
+ label=r"$\sigma_{\theta\theta,\mathrm{analytical}}$",
)
ax[0].plot(
r,
kirsch_sig_rt(10, r, -45, 2),
color="orangered",
linestyle=":",
- label="$\sigma_{r\\theta,\mathrm{analytical}}$",
+ label=r"$\sigma_{r\theta,\mathrm{analytical}}$",
)
-ax[0].plot(dist_sorted, sig_rr_sorted * 1000, label="$\sigma_{rr}$")
-ax[0].plot(dist_sorted, sig_tt_sorted * 1000, label="$\sigma_{\\theta\\theta}$")
-ax[0].plot(dist_sorted, sig_rt_sorted * 1000, label="$\sigma_{r\\theta}$")
+ax[0].plot(dist_sorted, sig_rr_sorted * 1000, label=r"$\sigma_{rr}$")
+ax[0].plot(dist_sorted, sig_tt_sorted * 1000, label=r"$\sigma_{\theta\theta}$")
+ax[0].plot(dist_sorted, sig_rt_sorted * 1000, label=r"$\sigma_{r\theta}$")
ax[0].legend(loc="upper right")
-ax[0].set_ylabel("$\\sigma$ / kPa")
+ax[0].set_ylabel(r"$\sigma$ / kPa")
ax[0].set(ylim=(0, 12))
-ax[1].plot(dist_sorted, f_abs_rr, label="$\Delta\sigma_{rr}$")
-ax[1].plot(dist_sorted, f_abs_tt, label="$\Delta\sigma_{\\theta\\theta}$")
-ax[1].plot(dist_sorted, f_abs_rt, label="$\Delta\sigma_{r\\theta}$")
+ax[1].plot(dist_sorted, f_abs_rr, label=r"$\Delta\sigma_{rr}$")
+ax[1].plot(dist_sorted, f_abs_tt, label=r"$\Delta\sigma_{\theta\theta}$")
+ax[1].plot(dist_sorted, f_abs_rt, label=r"$\Delta\sigma_{r\theta}$")
ax[1].spines["top"].set_color("none")
ax[1].xaxis.tick_bottom()
-ax[1].set_ylabel("$\Delta\\sigma$ / kPa")
+ax[1].set_ylabel(r"$\Delta\sigma$ / kPa")
ax[1].set(ylim=(-0.5, 1.3))
ax[1].legend(loc="upper right")
ax[2].plot(
dist_sorted,
f_rel_rr,
- label="$\Delta\sigma_{rr}$ / $\sigma_{rr,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{rr}$ / $\sigma_{rr,\mathrm{analytical}}$",
)
ax[2].plot(
dist_sorted,
f_rel_tt,
- label="$\Delta\sigma_{\\theta\\theta}$ / $\sigma_{\\theta\\theta,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{\theta\theta}$ / $\sigma_{\theta\theta,\mathrm{analytical}}$",
)
ax[2].plot(
dist_sorted,
f_rel_rt,
- label="$\Delta\sigma_{r\\theta}$ / $\sigma_{r\\theta,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{r\theta}$ / $\sigma_{r\theta,\mathrm{analytical}}$",
)
ax[2].set(ylim=(-0.1, 0.3))
-ax[2].set_ylabel("$\Delta\\sigma$ / $\sigma_{\mathrm{analytical}}$")
+ax[2].set_ylabel(r"$\Delta\sigma$ / $\sigma_{\mathrm{analytical}}$")
ax[2].legend(loc="upper right")
ax[0].set_title("Stress distribution")
@@ -622,11 +609,11 @@ ax[1].set_title("Absolute error")
ax[2].set_title("Relative error")
fig.tight_layout()
-```
-### Stress distribution along the $y$-axis ($\theta=0°$)
+# %% [markdown]
+# ### Stress distribution along the $y$-axis ($\theta=0°$)
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
pt1 = (0, 0, 0)
pt2 = (0, 10, 0)
yaxis = pv.Line(pt1, pt2, resolution=2)
@@ -637,13 +624,11 @@ ys = line_mesh.points[:, 1]
dist_from_origin = np.hypot(xs, ys)
indices_sorted = np.argsort(dist_from_origin)
dist_sorted = dist_from_origin[indices_sorted]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig = line_mesh.point_data["sigma"]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
num_points = sig.shape[0]
sig_rr = np.zeros(num_points)
sig_tt = np.zeros(num_points)
@@ -654,9 +639,8 @@ f_abs_rt = np.zeros(num_points)
f_rel_rr = np.zeros(num_points)
f_rel_tt = np.zeros(num_points)
f_rel_rt = np.zeros(num_points)
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
for pt_idx in range(num_points):
sig_vec = sig[pt_idx, :]
xs = line_mesh.points[pt_idx, 0]
@@ -665,15 +649,13 @@ for pt_idx in range(num_points):
sig_rr[pt_idx] = sig_polar[0, 0]
sig_tt[pt_idx] = sig_polar[1, 1]
sig_rt[pt_idx] = sig_polar[0, 1]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig_rr_sorted = sig_rr[indices_sorted]
sig_tt_sorted = sig_tt[indices_sorted]
sig_rt_sorted = sig_rt[indices_sorted]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
for pt_idx in range(num_points):
f_abs_rr[pt_idx] = sig_rr_sorted[pt_idx] * 1000 - kirsch_sig_rr(
10, dist_sorted[pt_idx], 0, 2
@@ -702,9 +684,8 @@ for pt_idx in range(num_points):
f_rel_rt[pt_idx] = f_abs_rt[pt_idx] / kirsch_sig_rt(
10, dist_sorted[pt_idx], 0, 2
)
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
r = np.linspace(2, 10, 1000)
fig, ax = plt.subplots(ncols=3, figsize=(18, 6))
@@ -719,28 +700,28 @@ ax[0].plot(
kirsch_sig_rr(10, r, 0, 2),
color="deepskyblue",
linestyle=":",
- label="$\sigma_{rr,\mathrm{analytical}}$",
+ label=r"$\sigma_{rr,\mathrm{analytical}}$",
)
ax[0].plot(
r,
kirsch_sig_tt(10, r, 0, 2),
color="yellowgreen",
linestyle=":",
- label="$\sigma_{\\theta\\theta,\mathrm{analytical}}$",
+ label=r"$\sigma_{\theta\theta,\mathrm{analytical}}$",
)
ax[0].plot(
r,
kirsch_sig_rt(10, r, 0, 2),
color="orangered",
linestyle=":",
- label="$\sigma_{r\\theta,\mathrm{analytical}}$",
+ label=r"$\sigma_{r\theta,\mathrm{analytical}}$",
)
-ax[0].plot(dist_sorted, sig_rr_sorted * 1000, label="$\sigma_{rr}$")
-ax[0].plot(dist_sorted, sig_tt_sorted * 1000, label="$\sigma_{\\theta\\theta}$")
-ax[0].plot(dist_sorted, sig_rt_sorted * 1000, label="$\sigma_{r\\theta}$")
+ax[0].plot(dist_sorted, sig_rr_sorted * 1000, label=r"$\sigma_{rr}$")
+ax[0].plot(dist_sorted, sig_tt_sorted * 1000, label=r"$\sigma_{\theta\theta}$")
+ax[0].plot(dist_sorted, sig_rt_sorted * 1000, label=r"$\sigma_{r\theta}$")
-ax[0].set_ylabel("$\\sigma$ / kPa")
+ax[0].set_ylabel(r"$\sigma$ / kPa")
ax[0].set(ylim=(-14, 10.5))
ax[0].legend(loc="lower right")
ax[0].spines["right"].set_color("none")
@@ -749,29 +730,29 @@ ax[0].xaxis.set_ticks_position("bottom")
ax[0].yaxis.set_ticks_position("left")
ax[0].spines["left"].set_position(("data", 0))
-ax[1].plot(dist_sorted, f_abs_rr, label="$\Delta\sigma_{rr}$")
-ax[1].plot(dist_sorted, f_abs_tt, label="$\Delta\sigma_{\\theta\\theta}$")
-ax[1].plot(dist_sorted, f_abs_rt, label="$\Delta\sigma_{r\\theta}$")
+ax[1].plot(dist_sorted, f_abs_rr, label=r"$\Delta\sigma_{rr}$")
+ax[1].plot(dist_sorted, f_abs_tt, label=r"$\Delta\sigma_{\theta\theta}$")
+ax[1].plot(dist_sorted, f_abs_rt, label=r"$\Delta\sigma_{r\theta}$")
ax[1].spines["top"].set_color("none")
ax[1].xaxis.tick_bottom()
-ax[1].set_ylabel("$\Delta\\sigma$ / kPa")
+ax[1].set_ylabel(r"$\Delta\sigma$ / kPa")
ax[1].set(ylim=(-3.5, 2.5))
ax[1].legend(loc="lower right")
ax[2].plot(
dist_sorted,
f_rel_rr,
- label="$\Delta\sigma_{rr}$ / $\sigma_{rr,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{rr}$ / $\sigma_{rr,\mathrm{analytical}}$",
)
ax[2].plot(
dist_sorted,
f_rel_tt,
- label="$\Delta\sigma_{\\theta\\theta}$ / $\sigma_{\\theta\\theta,\mathrm{analytical}}$",
+ label=r"$\Delta\sigma_{\theta\theta}$ / $\sigma_{\theta\theta,\mathrm{analytical}}$",
)
-# ax[2].plot(dist_sorted, f_rel_rt, label = "$\Delta\sigma_{r\\theta}$ / $\sigma_{r\\theta,analytical}$")
+# ax[2].plot(dist_sorted, f_rel_rt, label = r"$\Delta\sigma_{r\theta}$ / $\sigma_{r\theta,analytical}$")
ax[2].set(ylim=(-1, 1))
-ax[2].set_ylabel("$\Delta\\sigma$ / $\sigma_{\mathrm{analytical}}$")
+ax[2].set_ylabel(r"$\Delta\sigma$ / $\sigma_{\mathrm{analytical}}$")
ax[2].legend(loc="lower right")
ax[0].set_title("Stress distribution")
@@ -779,18 +760,18 @@ ax[1].set_title("Absolute error")
ax[2].set_title("Relative error")
fig.tight_layout()
-```
-
-### Error analysis
-
-In this following section, the differences of the numerical to the analytical solution are discussed. Doing this, one studies various aspects of the numerical solution.
-
-Whenever a problem is solved by finite element method (FEM), discretization errors occur. They arise from the creation of the mesh, since the result is calculated only at a finite number of nodes. Discretization errors can be minimized by refining the mesh. However, this increases computing time, so a compromise must be found.<br>
-Accordingly, this error smaller or larger and has to be acceptable for the field of application.
-
-The following image illustrates the chosen mesh for this benchmark example. In the area around the cavity, the stress gradients are rather high, as can also be seen in the diagrams above, and thus call for a finer mesh resolution.
-```python jupyter={"source_hidden": true}
+# %% [markdown]
+# ### Error analysis
+#
+# In this following section, the differences of the numerical to the analytical solution are discussed. Doing this, one studies various aspects of the numerical solution.
+#
+# Whenever a problem is solved by finite element method (FEM), discretization errors occur. They arise from the creation of the mesh, since the result is calculated only at a finite number of nodes. Discretization errors can be minimized by refining the mesh. However, this increases computing time, so a compromise must be found.<br>
+# Accordingly, this error smaller or larger and has to be acceptable for the field of application.
+#
+# The following image illustrates the chosen mesh for this benchmark example. In the area around the cavity, the stress gradients are rather high, as can also be seen in the diagrams above, and thus call for a finer mesh resolution.
+
+# %% jupyter={"source_hidden": true}
plotter = pv.Plotter(window_size=[1024, 500])
plotter.add_mesh(mesh, show_edges=True, show_scalar_bar=False, color=None, scalars=None)
@@ -798,32 +779,30 @@ plotter.show_bounds(ticks="outside", xlabel="x / cm", ylabel="y / cm")
# plotter.add_axes()
plotter.view_xy()
plotter.show()
-```
-All three stress components of the numerical solution show a more or less pronounced buckle most obviously along the diagonal- and y-axes in the middle of the observed radius range. The reason for this is the fact that the **mesh is refined** in the area around the hole. This increases the accuracy and resolution, but at the same time, the cells become distorted rectangles instead of squares.
-
-Another source of error is the **limitation of the plate's geometry**. In a continuous, infinitely extended plate the stress distribution approaches the primary stress state with distance from the cavity. The example of the stress distribution along the y-axis shows that the numerical solution for $\sigma_{rr}$ reaches the Neumann boundary condition set with a value of $10\, \text{kPa}$ for $\text{r} = 10\, \text{cm}$. With Kirsch's Solution for the analytical calculation, $\sigma_{rr}$ converges only at infinity towards the value of $10\, \text{kPa}$ and amounts approximately $9\, \text{kPa}$ for $r = 10\, \text{cm}$.<br>
-A similar error becomes visible along the $x$-axis. The Neumann boundary condition constrains the normal stress along the right edge to be zero . As a consequence $\sigma_{rr}$ must reach a value of $0\, \text{MPa}$ for $r = 10\, \text{cm}$ to satisfy the boundary condition. In the analytical solution however, $\sigma_{rr}$ approaches zero at infinity and thus is larger than in the numerical solution. Errors of that type can be minimized by increasing the plates dimensions or decreasing the diameter of the cavity, so that the boundary conditions do not interfere with the cavity-disturbed stress distribution.<br>
-As a rule of thumb, it can be assumed in this case that the boundaries of the problem can be drawn after at least five-hole radii.
-
-Further artefacts, called **extrapolation errors** occur at all boundaries of the numerical model. It can be recognised by a discontinuity in the graph for the stress distribution and absolute error next to the cavity ($r = 2\, \text{cm}$) and the plate boundaries ($r = 10\, \text{cm}$ resp. $r = 14.14\, \text{cm}$). The stress is derived from strain and this in turn from the nodal displacement. Therefore, there are no calculated stress magnitudes on the nodes, but they are formed as the average value of the neighboring cells. However, if the node is located on the edge directly next to the cavity or plate boundary, no average value can be formed, but the stress is extrapolated. In contrast to most of the other errors, this one is not based on an defect in the calculation of the numerical model, but rather in the output.
-
-### 2D Color Plots
-
-As can see below the 2D color plots of the numerical solution resemble those of the analytical solution. The plot of the tangential stress $\sigma_{\theta\theta}$ also shows a threefold stress concentration. From the plot of the radial stress $\sigma_{rr}$ it becomes clear that the Neumann boundary condition that respects the applied tensile stress of $10 \text{kPa}$ at the top boundary is satisfied.
-
-```python jupyter={"source_hidden": true}
+# %% [markdown]
+# All three stress components of the numerical solution show a more or less pronounced buckle most obviously along the diagonal- and y-axes in the middle of the observed radius range. The reason for this is the fact that the **mesh is refined** in the area around the hole. This increases the accuracy and resolution, but at the same time, the cells become distorted rectangles instead of squares.
+#
+# Another source of error is the **limitation of the plate's geometry**. In a continuous, infinitely extended plate the stress distribution approaches the primary stress state with distance from the cavity. The example of the stress distribution along the y-axis shows that the numerical solution for $\sigma_{rr}$ reaches the Neumann boundary condition set with a value of $10\, \text{kPa}$ for $\text{r} = 10\, \text{cm}$. With Kirsch's Solution for the analytical calculation, $\sigma_{rr}$ converges only at infinity towards the value of $10\, \text{kPa}$ and amounts approximately $9\, \text{kPa}$ for $r = 10\, \text{cm}$.<br>
+# A similar error becomes visible along the $x$-axis. The Neumann boundary condition constrains the normal stress along the right edge to be zero . As a consequence $\sigma_{rr}$ must reach a value of $0\, \text{MPa}$ for $r = 10\, \text{cm}$ to satisfy the boundary condition. In the analytical solution however, $\sigma_{rr}$ approaches zero at infinity and thus is larger than in the numerical solution. Errors of that type can be minimized by increasing the plates dimensions or decreasing the diameter of the cavity, so that the boundary conditions do not interfere with the cavity-disturbed stress distribution.<br>
+# As a rule of thumb, it can be assumed in this case that the boundaries of the problem can be drawn after at least five-hole radii.
+#
+# Further artefacts, called **extrapolation errors** occur at all boundaries of the numerical model. It can be recognised by a discontinuity in the graph for the stress distribution and absolute error next to the cavity ($r = 2\, \text{cm}$) and the plate boundaries ($r = 10\, \text{cm}$ resp. $r = 14.14\, \text{cm}$). The stress is derived from strain and this in turn from the nodal displacement. Therefore, there are no calculated stress magnitudes on the nodes, but they are formed as the average value of the neighboring cells. However, if the node is located on the edge directly next to the cavity or plate boundary, no average value can be formed, but the stress is extrapolated. In contrast to most of the other errors, this one is not based on an defect in the calculation of the numerical model, but rather in the output.
+#
+# ### 2D Color Plots
+#
+# As can see below the 2D color plots of the numerical solution resemble those of the analytical solution. The plot of the tangential stress $\sigma_{\theta\theta}$ also shows a threefold stress concentration. From the plot of the radial stress $\sigma_{rr}$ it becomes clear that the Neumann boundary condition that respects the applied tensile stress of $10 \text{kPa}$ at the top boundary is satisfied.
+
+# %% jupyter={"source_hidden": true}
points = mesh.point_data["sigma"].shape[0]
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
sig_polar = np.zeros([points, 3, 3])
sig_rr = np.zeros([points, 1])
sig_tt = np.zeros([points, 1])
sig_rt = np.zeros([points, 1])
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
xs = mesh.points[:, 0]
ys = mesh.points[:, 1]
sig_vec = mesh.point_data["sigma"]
@@ -832,27 +811,21 @@ for i in range(xs.shape[0]):
sig_rr[i] = sig_polar[i, 0, 0]
sig_tt[i] = sig_polar[i, 1, 1]
sig_rt[i] = sig_polar[i, 0, 1]
-```
-
-```python jupyter={"source_hidden": true}
-if False:
- print(np.min(sig_rr), np.max(sig_rr))
- print(np.min(sig_rt), np.max(sig_rt))
- print(np.min(sig_tt), np.max(sig_tt))
+# %% jupyter={"source_hidden": true}
mesh2 = mesh.copy(deep=False) # workaround some pyvista bug with multiple plots
mesh3 = mesh.copy(deep=False)
-sargs = dict(
- title="sigma_rr / kPa",
- title_font_size=25,
- label_font_size=15,
- n_labels=4,
- position_x=0.05,
- position_y=0.85,
- fmt="%.1f",
- width=0.9,
-)
+sargs = {
+ "title": "sigma_rr / kPa",
+ "title_font_size": 25,
+ "label_font_size": 15,
+ "n_labels": 4,
+ "position_x": 0.05,
+ "position_y": 0.85,
+ "fmt": "%.1f",
+ "width": 0.9,
+}
clim = [-10, 10]
p = pv.Plotter(shape=(1, 3), border=False)
@@ -871,16 +844,16 @@ p.add_mesh(
p.view_xy()
p.camera.zoom(2)
-sargs1 = dict(
- title="sigma_tt / kPa",
- title_font_size=25,
- label_font_size=15,
- n_labels=4,
- position_x=0.05,
- position_y=0.85,
- fmt="%.1f",
- width=0.9,
-)
+sargs1 = {
+ "title": "sigma_tt / kPa",
+ "title_font_size": 25,
+ "label_font_size": 15,
+ "n_labels": 4,
+ "position_x": 0.05,
+ "position_y": 0.85,
+ "fmt": "%.1f",
+ "width": 0.9,
+}
clim = [-35, 35]
p.subplot(0, 1)
@@ -898,16 +871,16 @@ p.add_mesh(
p.view_xy()
p.camera.zoom(2)
-sargs2 = dict(
- title="sigma_rt / kPa",
- title_font_size=25,
- label_font_size=15,
- n_labels=4,
- position_x=0.05,
- position_y=0.85,
- fmt="%.1f",
- width=0.9,
-)
+sargs2 = {
+ "title": "sigma_rt / kPa",
+ "title_font_size": 25,
+ "label_font_size": 15,
+ "n_labels": 4,
+ "position_x": 0.05,
+ "position_y": 0.85,
+ "fmt": "%.1f",
+ "width": 0.9,
+}
clim = [-10, 10]
p.subplot(0, 2)
@@ -927,21 +900,20 @@ p.camera.zoom(2)
p.window_size = [1000, 400]
p.show()
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
if True:
plotter = pv.Plotter()
- sargs = dict(
- title="sigma_rr",
- title_font_size=20,
- label_font_size=16,
- n_labels=4,
- position_x=0.2,
- position_y=0.9,
- fmt="%.1f",
- )
+ sargs = {
+ "title": "sigma_rr",
+ "title_font_size": 20,
+ "label_font_size": 16,
+ "n_labels": 4,
+ "position_x": 0.2,
+ "position_y": 0.9,
+ "fmt": "%.1f",
+ }
plotter.add_mesh(
mesh,
@@ -954,21 +926,20 @@ if True:
# plotter.add_axes()
plotter.view_xy()
plotter.show()
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
if False:
plotter = pv.Plotter()
- sargs = dict(
- title="sigma_tt",
- title_font_size=20,
- label_font_size=16,
- n_labels=4,
- position_x=0.2,
- position_y=0.9,
- fmt="%.1f",
- )
+ sargs = {
+ "title": "sigma_tt",
+ "title_font_size": 20,
+ "label_font_size": 16,
+ "n_labels": 4,
+ "position_x": 0.2,
+ "position_y": 0.9,
+ "fmt": "%.1f",
+ }
plotter.add_mesh(
mesh,
@@ -981,21 +952,20 @@ if False:
# plotter.add_axes()
plotter.view_xy()
plotter.show()
-```
-```python jupyter={"source_hidden": true}
+# %% jupyter={"source_hidden": true}
if False:
plotter = pv.Plotter()
- sargs = dict(
- title="sigma_rt",
- title_font_size=20,
- label_font_size=16,
- n_labels=4,
- position_x=0.2,
- position_y=0.9,
- fmt="%.1f",
- )
+ sargs = {
+ "title": "sigma_rt",
+ "title_font_size": 20,
+ "label_font_size": 16,
+ "n_labels": 4,
+ "position_x": 0.2,
+ "position_y": 0.9,
+ "fmt": "%.1f",
+ }
plotter.add_mesh(
mesh,
@@ -1008,4 +978,3 @@ if False:
plotter.add_axes()
plotter.view_xy()
plotter.show()
-```
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.md b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.md
deleted file mode 100644
index b5976b34739c68b3281dc92ffa63f4a8d836b592..0000000000000000000000000000000000000000
--- a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.md
+++ /dev/null
@@ -1,396 +0,0 @@
-+++
-title = "OGS-PHREEQC-Benchmark: CO2 injection into Opalinus Clay"
-date = "2023-09-07"
-author = "Shuang Chen, Vinay Kumar, Jobst Maßmann"
-web_subsection = "reactive-transport"
-image = "figures/results_comparision_sce1.png"
-+++
-
-
-$$\require{mhchem}$$
-
-## Overview
-
-The present numerical work is motivated by the CO2LPIE project (shortened as CL-experiment) [1], which is an in-situ experiment that is being conducted at the Mont Terri
-rock laboratory.
-In the experiment, carbon dioxide ($\ce{CO2}$) is injected into the Opalinus Clay leading to changes in its hydraulic, mechanical and chemical properties.
-In general, these processes are of great interest in the evaluation of barrier integrity.
-
-Two scenarios are considered in this benchmark.
-In the first scenario, the pure $\ce{CO2}$ gas induced Calcite dissolution is simulated by OGS-6-IPhreeqc and the results are verified with the related PHREEQC example presented by Appelo and Postma [2].
-A comprehensive information regarding the computational procedure of OGS-6-IPhreeqc can be found in Lu et al. [3].
-In the second scenario, the simulation considers a more accurate representation of chemical environments based on the CL Project.
-This includes the incorporation of primary minerals typically found in Opalinus Clay and the relevant composition of species present in the pore water.
-
-## Scenario #1: Kinetics of $\ce{CO2}$ induced calcite dissolution
-
-### Problem description
-
-In this case, carbon dioxide ($\ce{CO2}$) gas with a partial pressure $\mathrm{10^{-1.5}}$ $\mathrm{bar}$ is injected into a fluid containing calcite.
-The temperature is maintained at 10 $\mathrm{°C}$, the fluid’s pH value is set at 6 and the initial concentration of calcite is 1 $\mathrm{mol\cdot kg^{-1}\cdot water}$ (shortened as $\mathrm{mol/kgw}$).
-The dissolution of the $\ce{CO2}$ gas in water together with the calcite dissolution pathways can be described with the following reactions (Plummer et al. [4], Appelo and Postma [2]).
-
-$$ \ce{CO2_{(g)} + H2O -> H2CO3^{\ast}} $$
-$$ \ce{CaCO3 + H+ <=> Ca^2+ + HCO3-} $$
-$$ \ce{CaCO3 + H2CO3^{\ast} <=> Ca^2+ + 2HCO3-} $$
-$$ \ce{CaCO3 + H2O <=> Ca^2+ + HCO3- + OH-} $$
-
-Two kinetic rates are adopted to describe the calcite dissolution.
-The first approach involves a simplified rate calculation derived directly from the current concentration of calcium,
-employing the following formula: $\text{rate} = 10^{-6.91} - 10^{-1.52} * (c_{\ce{Ca}})^2$, with $c_{\ce{Ca}}$ in $\mathrm{mol/kgw}$.
-And another one is the well known PWP rate, which is proposed by Plummer, Wigley and Parkhurst (denoted further as “PWPâ€) in 1978 for calcite dissolution based on three dissolution reactions [4].
-The formulae of the PWP approach is available to be found in the PHREEQC database e.g. `phreeqc.dat`.
-In the numerical experiment, the total simulation time is 30 000 $\mathrm{s}$.
-Details about the case study is described in the example 5.9 from Appelo and Postma [2].
-The related PHREEQC script is available online to be found [5].
-
-### Model and results
-
-A simple 1D line-element model with one element is constructed in this work.
-The coupled Hydraulic-Chemical (HC) Process is adopted for the OGS-6-IPhreeqc simulation.
-To match the PHREEQC example, the advection and diffusion have been set to zero in the modelling.
-The initial conditions and chemical parameters are provided as outlined in the associated example.
-The Fig. 1 depicts the comparison between the computed results obtained from the OGS-6-IPhreeqc model and the results derived from PHREEQC, which shows a very
-good agreement with each other.
-Preparation of the simulation and results evaluation are presented in what follows.
-
-### 1) Solve the numerical model
-
-```python
-import os
-import vtuIO
-import numpy as np
-import matplotlib.pyplot as plt
-import pandas as pd
-import time
-```
-
-```python
-from pathlib import Path
-
-out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
-if not out_dir.exists():
- out_dir.mkdir(parents=True)
-
-# simple
-prj_name = "calcite_simple.prj"
-
-t0 = time.time()
-print(">>> OGS started execution ... <<<")
-
-!ogs {prj_name} -o {out_dir} > {out_dir}/out_simple.txt
-
-tf = time.time()
-print(">>> OGS terminated execution <<< Elapsed time: ", round(tf - t0, 2), " s.")
-
-# check if the simulation runs successfully
-check_file = os.path.isfile(f"{out_dir}/calcite_simple_ts_43_t_30000.000000.vtu")
-if check_file:
- print("OGS simulation for the scenario simple case runs successfully")
-else:
- msg = "OGS simulation failed."
- raise Exception(msg)
-
-# pwp
-prj_name = "calcite_pwp.prj"
-
-t0 = time.time()
-print(">>> OGS started execution ... <<<")
-
-!ogs {prj_name} -o {out_dir} > {out_dir}/out_pwp.txt
-
-tf = time.time()
-print(">>> OGS terminated execution <<< Elapsed time: ", round(tf - t0, 2), " s.")
-
-# check if the simulation runs successfully
-check_file = os.path.isfile(f"{out_dir}/calcite_pwp_ts_43_t_30000.000000.vtu")
-if check_file:
- print("OGS simulation for the scenario PWP case runs successfully")
-else:
- msg = "OGS simulation failed."
- raise Exception(msg)
-```
-
-```python
-### Read OGS-6 simulation results
-
-# simple
-pvdfile_simple = vtuIO.PVDIO(f"{out_dir}/calcite_simple.pvd", dim=1)
-
-# pwp
-pvdfile_pwp = vtuIO.PVDIO(f"{out_dir}/calcite_pwp.pvd", dim=1)
-
-# Read PHREEQC results
-
-# simple
-phreeqc_result_df_simple = pd.read_csv(
- "./PHREEQC_results_simple.txt", sep="\s+", header=0, skipinitialspace=True
-)
-
-# pwp
-phreeqc_result_df_pwp = pd.read_csv(
- "./PHREEQC_results_PWP.txt", sep="\s+", header=0, skipinitialspace=True
-)
-
-### visulalization
-# results at point [0, 0, 0] are selected for display
-point = {"pt0": (0.0, 0.0, 0.0)}
-
-# plot
-fig1 = plt.figure(figsize=(8, 6))
-ax = plt.subplot(111)
-
-## results of 'Ca'
-# OGS-6
-ax.plot(
- pvdfile_simple.timesteps,
- pvdfile_simple.read_time_series("Ca", point)["pt0"] * 1e3,
- c="k",
- linestyle="--",
- label="Simple " + "OGS-6",
-)
-ax.plot(
- pvdfile_pwp.timesteps,
- pvdfile_pwp.read_time_series("Ca", point)["pt0"] * 1e3,
- c="k",
- label="PWP " + f"OGS-6",
-)
-# PHREEQC
-ax.scatter(
- phreeqc_result_df_simple.loc[:, "x"] * 1e3,
- phreeqc_result_df_simple.loc[:, "Ca"],
- c="r",
- s=80,
- marker="x",
- label="Simple " + f"PHREEQC",
-)
-ax.scatter(
- phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
- phreeqc_result_df_pwp.loc[:, "Ca"],
- c="b",
- s=80,
- marker="x",
- label="PWP " + f"PHREEQC",
-)
-ax.set_xlim([0, 3e4])
-ax.set_ylim([0, 4])
-ax.set_ylabel("Ca+ [mmol/kgw]")
-ax.set_xlabel("Time [s]")
-ann = ax.annotate(
- "Ca+",
- xy=(5000, 0.7),
- xycoords="data",
- va="center",
- ha="center",
- bbox=dict(boxstyle="round", fc="w"),
- fontsize=16,
-)
-ax.legend(loc="upper left", fontsize=10)
-
-## results of 'pH'
-ax2 = ax.twinx()
-# OGS-6
-ax2.plot(
- pvdfile_simple.timesteps,
- -np.log10(pvdfile_simple.read_time_series("H", point)["pt0"]),
- c="k",
- linestyle="--",
- label="Simple " + "OGS-6",
-)
-ax2.plot(
- pvdfile_pwp.timesteps,
- -np.log10(pvdfile_pwp.read_time_series("H", point)["pt0"]),
- c="k",
- label="PWP " + f"OGS-6",
-)
-# PHREEQC
-ax2.scatter(
- phreeqc_result_df_simple.loc[:, "x"] * 1e3,
- phreeqc_result_df_simple.loc[:, "pH"],
- c="r",
- s=80,
- marker="x",
- label="simple " + f"PHREEQC",
-)
-ax2.scatter(
- phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
- phreeqc_result_df_pwp.loc[:, "pH"],
- c="b",
- s=80,
- marker="x",
- label="PWP " + f"PHREEQC",
-)
-ax2.set_xlim([0, 3e4])
-ax2.set_ylim([4, 7.5])
-ax2.set_ylabel("pH")
-ax2.set_xlabel("Time [s]")
-ann1 = ax2.annotate(
- "pH",
- xy=(5000, 6.4),
- xycoords="data",
- va="center",
- ha="center",
- bbox=dict(boxstyle="round", fc="w"),
- fontsize=16,
-)
-```
-
-Fig. 1: pH and calcium increase with kinetic dissolution of Calcite.
-
-```python
-## L2 difference
-
-# OGS-6 results for comparision
-OGS6_pwp_ca = (
- np.delete((pvdfile_pwp.read_time_series("Ca", point)["pt0"]), np.arange(0, 24, 1))
- * 1e3
-)
-OGS6_pwp_H = np.delete(
- pvdfile_pwp.read_time_series("H", point)["pt0"], np.arange(0, 24, 1)
-)
-
-# plot
-fig1 = plt.figure(figsize=(8, 6))
-ax = plt.subplot(111)
-ax.plot(
- phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
- np.log10(
- ((OGS6_pwp_ca - phreeqc_result_df_pwp.loc[1:, "Ca"]) ** 2) ** 0.5
- / phreeqc_result_df_pwp.loc[:, "Ca"]
- ),
- c="g",
- marker="o",
- markersize=4,
- label="Ca+ ",
-)
-ax.plot(
- phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
- np.log10(
- ((OGS6_pwp_H - 10 ** (-phreeqc_result_df_pwp.loc[1:, "pH"])) ** 2) ** 0.5
- / 10 ** (-phreeqc_result_df_pwp.loc[:, "pH"])
- ),
- c="k",
- marker="o",
- markersize=4,
- label="H+ ",
-)
-ax.set_xlim([1500, 3e4])
-ax.set_xticks(np.arange(0, 30001, 5000))
-ax.set_ylim([-5, 0])
-ax.set_xlabel("Time [s]")
-ax.set_ylabel(
- r"Log($\frac{|| \mathbf{c}^\mathrm{OGS-6} - \mathbf{c}^{\mathrm{PHREEQC}}||_{2}}{\mathbf{c}^{\mathrm{PHREEQC}}}$)"
-)
-ax.grid(color="gray", linestyle="dashed")
-ax.legend(loc="best", fontsize=11)
-```
-
-Fig. 2: L2 relative difference norm between the obtained results from OGS-6 and PHREEQC for the scenario with the PWP rate.
-
-## Scenario #2: Modelling of the $\ce{CO2}$ injection into Opalinus Clay
-
-### Model description
-
-The identical 1D line-element model as in scenario #1 is used in the OGS-6-IPhreeqc simulation.
-In the model, the initial porosity is set to 0.15. Similarly to the scenario #1, hydraulic advection and diffusion are excluded from the simulation.
-In this case, chemical environments based on the CL-experiment are considered.
-For the porous medium in the OGS-6 model, in terms of the solid component, the following mineral composition is assumed in this work: 36% illite, 24% kaolinite, 7.5% calcite and 2.5% dolomite-dis (namely sedimentary (disordered) dolomite).
-The remaining 30% of the mineral is mostly quartz and is not considered in this simulation scenario.
-A more detailed description of the mineral composition of Opalinus Clay can be found in the work of Thury [6].
-In OGS-6-IPhreeqc, the molar amount of each reactive solid component per kilogram of water ${b_m}$ [$\mathrm{mol/kgw}$] is calculated by
-
-$$ {b_m} = \frac{\phi_m}{\rho^{l}\phi V_m},$$
-
-where $\mathrm{\phi_m}$ is the volume fraction of the corresponding mineral ${m}$.
-The $\mathrm{\rho^{l}}$ is the density of the fluid in $\mathrm{kg/m^3}$.
-The $\mathrm{\phi}$ denotes the porosity of the porous medium.
-And the ${V_m}$ is the molar volume of the corresponding mineral ${m}$ in $\mathrm{m^3/mol}$.
-Details of the equation description can be found in Lu et al. [3].
-In the simulation, the used parameters of the Opalinus Clay minerals are listed in table 1.
-
-During the simulation, a constant carbon dioxide ($\ce{CO2}$) gas partial pressure of $10^{1.5}$ $\mathrm{bar}$ is applied to the model domain, same as in the equilibrium phase.
-Notice in this scenario, the $\ce{CO2}$ partial pressure is much higher than it in scenario #1.
-Consequently, the dissolved $\ce{CO2}$ results in an alteration in the acidity of the pore water and leads to the different chemical reactions of each mineral.
-The related reaction formulae are described in the PHREEQC database `llnl.dat`.
-The adopted kinetic dissolution rate of each mineral are referenced from the transition state theory-based reaction mechanism following the work by Palandri and Kharaka [7]. The general equation formula reads
-
-$${{\text{Rate}_\mathrm{mineral}}} = [ k_{\mathrm{acid_{mineral}}}^\mathrm{298.15K}{a_\mathrm{H^{+}}} + k_{\mathrm{neutral_{mineral}}}^\mathrm{298.15K} + k_{\mathrm{base_{mineral}}}^\mathrm{298.15K}{a_\mathrm{{H_{2}CO_{3}}^{*}}} ](1 - \mathrm{SR_{mineral}}),
-\label{eq:transition_state}$$ with
-$$\mathrm{SR_{mineral}} = \frac{\mathrm{IAP_{mineral}}}{K_\mathrm{eq,mineral}},$$
-
-where $a$ denotes the activity of the ion, and $\mathrm{SR}$ is the abbreviation for the saturation ratio of a phase, which describes the ion activity product $\mathrm{IAP}$ divided by equilibrium constant $K_\mathrm{eq}$.
-In the simulation, the main species composition within the Opalinus Clay pore water are considered and the corresponding values are listed in table 2, following the
-work of Wersin et al. [8].
-
-| Minerals | Volume fraction [-] | Molar volume in [$\mathrm{m^3/mol}$] | Molality in [$\mathrm{mol/kgw}$] |
-|:--------------:|:-------------------------:|:--------------------------------------:|:---------------------------------------:|
-| Calcite | 0.06375 | 3.6934e-5 | 11.507 |
-| Dolomite-dis | 0.02125 | 6.439e-5 | 2.2 |
-| Illite | 0.306 | 1.4148e-4 | 14.419 |
-| Kaolinite | 0.204 | 9.935e-5 | 13.689 |
-
-Table 1: Parameter of the Opalinus Clay minerals used in the modeling
-
-| Species | Value in [mol/kgw] |
-|:---------:|:-------------------------:|
-| C(4) | $3.89 \times 10^{-3}$ |
-| Ca | $1.89 \times 10^{-2}$ |
-| Mg | $2.197 \times 10^{-2}$ |
-| Cl | $3.2667 \times 10^{-1}$ |
-| K | $1.92 \times 10^{-3}$ |
-| Na | $2.8067 \times 10^{-1}$ |
-| S(6) | $1.679 \times 10^{-2}$ |
-| Al | $2.295 \times 10^{-10}$ |
-| Si | $1.7 \times 10^{-5}$ |
-| Sr | $4.6 \times 10^{-5}$ |
-
-Table 2: concentration of the species in the Opalinus Clay pore water
-
-For verification purposes, a corresponding PHREEQC model was constructed with identical parameter settings.
-
-### Results
-
-On the Fig. 3 a) the evolution of calcium and magnesium molality during the $\ce{CO2}$ injection computed from the OGS-6 and PHREEQC model is shown.
-And the Fig. 3 b) displays the L2 relative difference norm between the results computed by the two software.
-The results of the two software programs match perfectly.
-On the Fig. 4 the evolution of the dissoluted molality of each minerals during the simulation is illustrated.
-Only calcite and dolomite have been partially dissolved due to $\ce{CO2}$ injection.
-In contrast, clay minerals, the illite and kaolinite were hardly affected.
-The kinetic dissolution rate of calcite and dolomite is mostly controlled by their saturation ratio states.
-When the SR value of calcite and dolomite reaches 1, the dissolution process of the correspondences is stopped.
-
-**Note:**
-The scenario #2 is not directly calculated in the notebook due to the lack of the large input file `llnl.dat` database in the commit. Also the post-processing of the scenario #2 is similar to which is showed in the scenario #1.
-
-
-Fig. 3: a) evolution of the calcium and magnesium molality during the $\ce{CO2}$ injection; b) L2 relative difference norm between the results computed by OGS-6 and PHREEQC.
-
-
-Fig. 4: Evolution of the dissoluted minerals molality and the related saturation ratio of calcite and dolomite over the time.
-
-## Literature
-
-<!-- vale off -->
-
-[1] BGR, CO2LPIE project, 2023. URL: [https://www.bgr.bund.de/DE/Themen/Endlagerung/Projekte/Wirtsgesteine_geotechnische_Barrieren/laufend/Nur-Deutsch/mont_terri_experimente.html?nn=1542156](https://www.bgr.bund.de/DE/Themen/Endlagerung/Projekte/Wirtsgesteine_geotechnische_Barrieren/laufend/Nur-Deutsch/mont_terri_experimente.html?nn=1542156).
-
-[2] Appelo, C. A. J. and Postma, D. (2004). Geochemistry, groundwater and pollution. CRC press.
-
-[3] Lu, R., Nagel, T., Poonoosamy, J., Naumov, D., Fischer, T., Montoya, V., Kolditz, O., and Shao, H.
-(2022). A new operator-splitting finite element scheme for reactive transport modeling in saturated porous
-media. Computers and Geosciences, 163(April 2021):105106. URL: https://doi.org/10.1016/j.cageo.2022.105106.
-
-[4] Plummer, L. N., Wigley, T. M., and Parkhurst, D. L. (1978). KINETICS OF CALCITE DISSOLUTION
-IN CO2-WATER SYSTEMS AT 5 degree TO 60 degree C AND 0. 0 TO 1. 0 ATM CO2.
-
-[5] Appelo, C. A. J. and Postma, D. (2023). Example 5.9: Kinetic dissolution of calcite. URL: https://www.hydrochemistry.eu/a&p/5/ex_5_9.phr.
-
-[6] Thury, M. (2002). The characteristics of the Opalinus Clay investigated in the Mont Terri underground
-rock laboratory in Switzerland. doi:10.1016/S1631-0705(02)01372-5.
-
-[7] Palandri, J. L. and Kharaka, Y. K. (2004). A compilation of rate parameters of water-mineral interaction
-kinetics for application to geochemical modeling.
-
-[8] Wersin, P., Mazurek, M., and Gimmi, T. (2022). Porewater chemistry of Opalinus Clay revisited:
-Findings from 25 years of data collection at the Mont Terri Rock Laboratory. Applied Geochemistry,
-138(November 2021):105234. URL: https://doi.org/10.1016/j.apgeochem.2022.105234.
diff --git a/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.py b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.py
new file mode 100644
index 0000000000000000000000000000000000000000..2d622b4aee63eb741654f79f86f122fc23072e08
--- /dev/null
+++ b/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/CO2Injection/CO2Injection.py
@@ -0,0 +1,395 @@
+# %% [markdown]
+# +++
+# title = "OGS-PHREEQC-Benchmark: CO2 injection into Opalinus Clay"
+# date = "2023-09-07"
+# author = "Shuang Chen, Vinay Kumar, Jobst Maßmann"
+# web_subsection = "reactive-transport"
+# image = "figures/results_comparision_sce1.png"
+# +++
+
+# %% [markdown]
+# $$\require{mhchem}$$
+#
+# ## Overview
+#
+# The present numerical work is motivated by the CO2LPIE project (shortened as CL-experiment) [1], which is an in-situ experiment that is being conducted at the Mont Terri
+# rock laboratory.
+# In the experiment, carbon dioxide ($\ce{CO2}$) is injected into the Opalinus Clay leading to changes in its hydraulic, mechanical and chemical properties.
+# In general, these processes are of great interest in the evaluation of barrier integrity.
+#
+# Two scenarios are considered in this benchmark.
+# In the first scenario, the pure $\ce{CO2}$ gas induced Calcite dissolution is simulated by OGS-6-IPhreeqc and the results are verified with the related PHREEQC example presented by Appelo and Postma [2].
+# A comprehensive information regarding the computational procedure of OGS-6-IPhreeqc can be found in Lu et al. [3].
+# In the second scenario, the simulation considers a more accurate representation of chemical environments based on the CL Project.
+# This includes the incorporation of primary minerals typically found in Opalinus Clay and the relevant composition of species present in the pore water.
+#
+# ## Scenario #1: Kinetics of $\ce{CO2}$ induced calcite dissolution
+#
+# ### Problem description
+#
+# In this case, carbon dioxide ($\ce{CO2}$) gas with a partial pressure $\mathrm{10^{-1.5}}$ $\mathrm{bar}$ is injected into a fluid containing calcite.
+# The temperature is maintained at 10 $\mathrm{°C}$, the fluid`s pH value is set at 6 and the initial concentration of calcite is 1 $\mathrm{mol\cdot kg^{-1}\cdot water}$ (shortened as $\mathrm{mol/kgw}$).
+# The dissolution of the $\ce{CO2}$ gas in water together with the calcite dissolution pathways can be described with the following reactions (Plummer et al. [4], Appelo and Postma [2]).
+#
+# $$ \ce{CO2_{(g)} + H2O -> H2CO3^{\ast}} $$
+# $$ \ce{CaCO3 + H+ <=> Ca^2+ + HCO3-} $$
+# $$ \ce{CaCO3 + H2CO3^{\ast} <=> Ca^2+ + 2HCO3-} $$
+# $$ \ce{CaCO3 + H2O <=> Ca^2+ + HCO3- + OH-} $$
+#
+# Two kinetic rates are adopted to describe the calcite dissolution.
+# The first approach involves a simplified rate calculation derived directly from the current concentration of calcium,
+# employing the following formula: $\text{rate} = 10^{-6.91} - 10^{-1.52} * (c_{\ce{Ca}})^2$, with $c_{\ce{Ca}}$ in $\mathrm{mol/kgw}$.
+# And another one is the well known PWP rate, which is proposed by Plummer, Wigley and Parkhurst (denoted further as “PWPâ€) in 1978 for calcite dissolution based on three dissolution reactions [4].
+# The formulae of the PWP approach is available to be found in the PHREEQC database e.g. `phreeqc.dat`.
+# In the numerical experiment, the total simulation time is 30 000 $\mathrm{s}$.
+# Details about the case study is described in the example 5.9 from Appelo and Postma [2].
+# The related PHREEQC script is available online to be found [5].
+#
+# ### Model and results
+#
+# A simple 1D line-element model with one element is constructed in this work.
+# The coupled Hydraulic-Chemical (HC) Process is adopted for the OGS-6-IPhreeqc simulation.
+# To match the PHREEQC example, the advection and diffusion have been set to zero in the modelling.
+# The initial conditions and chemical parameters are provided as outlined in the associated example.
+# The Fig. 1 depicts the comparison between the computed results obtained from the OGS-6-IPhreeqc model and the results derived from PHREEQC, which shows a very
+# good agreement with each other.
+# Preparation of the simulation and results evaluation are presented in what follows.
+#
+# ### 1) Solve the numerical model
+
+# %%
+import os
+import time
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import vtuIO
+
+# %%
+out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
+if not out_dir.exists():
+ out_dir.mkdir(parents=True)
+
+# simple
+prj_name = "calcite_simple.prj"
+
+t0 = time.time()
+print(">>> OGS started execution ... <<<")
+
+# !ogs {prj_name} -o {out_dir} > {out_dir}/out_simple.txt
+
+tf = time.time()
+print(">>> OGS terminated execution <<< Elapsed time: ", round(tf - t0, 2), " s.")
+
+# check if the simulation runs successfully
+check_file = Path(f"{out_dir}/calcite_simple_ts_43_t_30000.000000.vtu").is_file()
+if check_file:
+ print("OGS simulation for the scenario simple case runs successfully")
+else:
+ msg = "OGS simulation failed."
+ raise Exception(msg)
+
+# pwp
+prj_name = "calcite_pwp.prj"
+
+t0 = time.time()
+print(">>> OGS started execution ... <<<")
+
+# !ogs {prj_name} -o {out_dir} > {out_dir}/out_pwp.txt
+
+tf = time.time()
+print(">>> OGS terminated execution <<< Elapsed time: ", round(tf - t0, 2), " s.")
+
+# check if the simulation runs successfully
+check_file = Path(f"{out_dir}/calcite_pwp_ts_43_t_30000.000000.vtu").is_file()
+if check_file:
+ print("OGS simulation for the scenario PWP case runs successfully")
+else:
+ msg = "OGS simulation failed."
+ raise Exception(msg)
+
+# %%
+### Read OGS-6 simulation results
+
+# simple
+pvdfile_simple = vtuIO.PVDIO(f"{out_dir}/calcite_simple.pvd", dim=1)
+
+# pwp
+pvdfile_pwp = vtuIO.PVDIO(f"{out_dir}/calcite_pwp.pvd", dim=1)
+
+# Read PHREEQC results
+
+# simple
+phreeqc_result_df_simple = pd.read_csv(
+ "./PHREEQC_results_simple.txt", sep=r"\s+", header=0, skipinitialspace=True
+)
+
+# pwp
+phreeqc_result_df_pwp = pd.read_csv(
+ "./PHREEQC_results_PWP.txt", sep=r"\s+", header=0, skipinitialspace=True
+)
+
+### visulalization
+# results at point [0, 0, 0] are selected for display
+point = {"pt0": (0.0, 0.0, 0.0)}
+
+# plot
+fig1 = plt.figure(figsize=(8, 6))
+ax = plt.subplot(111)
+
+## results of 'Ca'
+# OGS-6
+ax.plot(
+ pvdfile_simple.timesteps,
+ pvdfile_simple.read_time_series("Ca", point)["pt0"] * 1e3,
+ c="k",
+ linestyle="--",
+ label="Simple " + "OGS-6",
+)
+ax.plot(
+ pvdfile_pwp.timesteps,
+ pvdfile_pwp.read_time_series("Ca", point)["pt0"] * 1e3,
+ c="k",
+ label="PWP " + "OGS-6",
+)
+# PHREEQC
+ax.scatter(
+ phreeqc_result_df_simple.loc[:, "x"] * 1e3,
+ phreeqc_result_df_simple.loc[:, "Ca"],
+ c="r",
+ s=80,
+ marker="x",
+ label="Simple " + "PHREEQC",
+)
+ax.scatter(
+ phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
+ phreeqc_result_df_pwp.loc[:, "Ca"],
+ c="b",
+ s=80,
+ marker="x",
+ label="PWP " + "PHREEQC",
+)
+ax.set_xlim([0, 3e4])
+ax.set_ylim([0, 4])
+ax.set_ylabel("Ca+ [mmol/kgw]")
+ax.set_xlabel("Time [s]")
+ann = ax.annotate(
+ "Ca+",
+ xy=(5000, 0.7),
+ xycoords="data",
+ va="center",
+ ha="center",
+ bbox={"boxstyle": "round", "fc": "w"},
+ fontsize=16,
+)
+ax.legend(loc="upper left", fontsize=10)
+
+## results of 'pH'
+ax2 = ax.twinx()
+# OGS-6
+ax2.plot(
+ pvdfile_simple.timesteps,
+ -np.log10(pvdfile_simple.read_time_series("H", point)["pt0"]),
+ c="k",
+ linestyle="--",
+ label="Simple " + "OGS-6",
+)
+ax2.plot(
+ pvdfile_pwp.timesteps,
+ -np.log10(pvdfile_pwp.read_time_series("H", point)["pt0"]),
+ c="k",
+ label="PWP " + "OGS-6",
+)
+# PHREEQC
+ax2.scatter(
+ phreeqc_result_df_simple.loc[:, "x"] * 1e3,
+ phreeqc_result_df_simple.loc[:, "pH"],
+ c="r",
+ s=80,
+ marker="x",
+ label="simple " + "PHREEQC",
+)
+ax2.scatter(
+ phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
+ phreeqc_result_df_pwp.loc[:, "pH"],
+ c="b",
+ s=80,
+ marker="x",
+ label="PWP " + "PHREEQC",
+)
+ax2.set_xlim([0, 3e4])
+ax2.set_ylim([4, 7.5])
+ax2.set_ylabel("pH")
+ax2.set_xlabel("Time [s]")
+ann1 = ax2.annotate(
+ "pH",
+ xy=(5000, 6.4),
+ xycoords="data",
+ va="center",
+ ha="center",
+ bbox={"boxstyle": "round", "fc": "w"},
+ fontsize=16,
+)
+
+# %% [markdown]
+# Fig. 1: pH and calcium increase with kinetic dissolution of Calcite.
+
+# %%
+## L2 difference
+
+# OGS-6 results for comparision
+OGS6_pwp_ca = (
+ np.delete((pvdfile_pwp.read_time_series("Ca", point)["pt0"]), np.arange(0, 24, 1))
+ * 1e3
+)
+OGS6_pwp_H = np.delete(
+ pvdfile_pwp.read_time_series("H", point)["pt0"], np.arange(0, 24, 1)
+)
+
+# plot
+fig1 = plt.figure(figsize=(8, 6))
+ax = plt.subplot(111)
+ax.plot(
+ phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
+ np.log10(
+ ((OGS6_pwp_ca - phreeqc_result_df_pwp.loc[1:, "Ca"]) ** 2) ** 0.5
+ / phreeqc_result_df_pwp.loc[:, "Ca"]
+ ),
+ c="g",
+ marker="o",
+ markersize=4,
+ label="Ca+ ",
+)
+ax.plot(
+ phreeqc_result_df_pwp.loc[:, "x"] * 1e3,
+ np.log10(
+ ((OGS6_pwp_H - 10 ** (-phreeqc_result_df_pwp.loc[1:, "pH"])) ** 2) ** 0.5
+ / 10 ** (-phreeqc_result_df_pwp.loc[:, "pH"])
+ ),
+ c="k",
+ marker="o",
+ markersize=4,
+ label="H+ ",
+)
+ax.set_xlim([1500, 3e4])
+ax.set_xticks(np.arange(0, 30001, 5000))
+ax.set_ylim([-5, 0])
+ax.set_xlabel("Time [s]")
+ax.set_ylabel(
+ r"Log($\frac{|| \mathbf{c}^\mathrm{OGS-6} - \mathbf{c}^{\mathrm{PHREEQC}}||_{2}}{\mathbf{c}^{\mathrm{PHREEQC}}}$)"
+)
+ax.grid(color="gray", linestyle="dashed")
+ax.legend(loc="best", fontsize=11)
+
+# %% [markdown]
+# Fig. 2: L2 relative difference norm between the obtained results from OGS-6 and PHREEQC for the scenario with the PWP rate.
+#
+# ## Scenario #2: Modelling of the $\ce{CO2}$ injection into Opalinus Clay
+#
+# ### Model description
+#
+# The identical 1D line-element model as in scenario #1 is used in the OGS-6-IPhreeqc simulation.
+# In the model, the initial porosity is set to 0.15. Similarly to the scenario #1, hydraulic advection and diffusion are excluded from the simulation.
+# In this case, chemical environments based on the CL-experiment are considered.
+# For the porous medium in the OGS-6 model, in terms of the solid component, the following mineral composition is assumed in this work: 36% illite, 24% kaolinite, 7.5% calcite and 2.5% dolomite-dis (namely sedimentary (disordered) dolomite).
+# The remaining 30% of the mineral is mostly quartz and is not considered in this simulation scenario.
+# A more detailed description of the mineral composition of Opalinus Clay can be found in the work of Thury [6].
+# In OGS-6-IPhreeqc, the molar amount of each reactive solid component per kilogram of water ${b_m}$ [$\mathrm{mol/kgw}$] is calculated by
+#
+# $$ {b_m} = \frac{\phi_m}{\rho^{l}\phi V_m},$$
+#
+# where $\mathrm{\phi_m}$ is the volume fraction of the corresponding mineral ${m}$.
+# The $\mathrm{\rho^{l}}$ is the density of the fluid in $\mathrm{kg/m^3}$.
+# The $\mathrm{\phi}$ denotes the porosity of the porous medium.
+# And the ${V_m}$ is the molar volume of the corresponding mineral ${m}$ in $\mathrm{m^3/mol}$.
+# Details of the equation description can be found in Lu et al. [3].
+# In the simulation, the used parameters of the Opalinus Clay minerals are listed in table 1.
+#
+# During the simulation, a constant carbon dioxide ($\ce{CO2}$) gas partial pressure of $10^{1.5}$ $\mathrm{bar}$ is applied to the model domain, same as in the equilibrium phase.
+# Notice in this scenario, the $\ce{CO2}$ partial pressure is much higher than it in scenario #1.
+# Consequently, the dissolved $\ce{CO2}$ results in an alteration in the acidity of the pore water and leads to the different chemical reactions of each mineral.
+# The related reaction formulae are described in the PHREEQC database `llnl.dat`.
+# The adopted kinetic dissolution rate of each mineral are referenced from the transition state theory-based reaction mechanism following the work by Palandri and Kharaka [7]. The general equation formula reads
+#
+# $${{\text{Rate}_\mathrm{mineral}}} = [ k_{\mathrm{acid_{mineral}}}^\mathrm{298.15K}{a_\mathrm{H^{+}}} + k_{\mathrm{neutral_{mineral}}}^\mathrm{298.15K} + k_{\mathrm{base_{mineral}}}^\mathrm{298.15K}{a_\mathrm{{H_{2}CO_{3}}^{*}}} ](1 - \mathrm{SR_{mineral}}),
+# \label{eq:transition_state}$$ with
+# $$\mathrm{SR_{mineral}} = \frac{\mathrm{IAP_{mineral}}}{K_\mathrm{eq,mineral}},$$
+#
+# where $a$ denotes the activity of the ion, and $\mathrm{SR}$ is the abbreviation for the saturation ratio of a phase, which describes the ion activity product $\mathrm{IAP}$ divided by equilibrium constant $K_\mathrm{eq}$.
+# In the simulation, the main species composition within the Opalinus Clay pore water are considered and the corresponding values are listed in table 2, following the
+# work of Wersin et al. [8].
+#
+# | Minerals | Volume fraction [-] | Molar volume in [$\mathrm{m^3/mol}$] | Molality in [$\mathrm{mol/kgw}$] |
+# |:--------------:|:-------------------------:|:--------------------------------------:|:---------------------------------------:|
+# | Calcite | 0.06375 | 3.6934e-5 | 11.507 |
+# | Dolomite-dis | 0.02125 | 6.439e-5 | 2.2 |
+# | Illite | 0.306 | 1.4148e-4 | 14.419 |
+# | Kaolinite | 0.204 | 9.935e-5 | 13.689 |
+#
+# Table 1: Parameter of the Opalinus Clay minerals used in the modeling
+#
+# | Species | Value in [mol/kgw] |
+# |:---------:|:-------------------------:|
+# | C(4) | $3.89 \times 10^{-3}$ |
+# | Ca | $1.89 \times 10^{-2}$ |
+# | Mg | $2.197 \times 10^{-2}$ |
+# | Cl | $3.2667 \times 10^{-1}$ |
+# | K | $1.92 \times 10^{-3}$ |
+# | Na | $2.8067 \times 10^{-1}$ |
+# | S(6) | $1.679 \times 10^{-2}$ |
+# | Al | $2.295 \times 10^{-10}$ |
+# | Si | $1.7 \times 10^{-5}$ |
+# | Sr | $4.6 \times 10^{-5}$ |
+#
+# Table 2: concentration of the species in the Opalinus Clay pore water
+#
+# For verification purposes, a corresponding PHREEQC model was constructed with identical parameter settings.
+#
+# ### Results
+#
+# On the Fig. 3 a) the evolution of calcium and magnesium molality during the $\ce{CO2}$ injection computed from the OGS-6 and PHREEQC model is shown.
+# And the Fig. 3 b) displays the L2 relative difference norm between the results computed by the two software.
+# The results of the two software programs match perfectly.
+# On the Fig. 4 the evolution of the dissoluted molality of each minerals during the simulation is illustrated.
+# Only calcite and dolomite have been partially dissolved due to $\ce{CO2}$ injection.
+# In contrast, clay minerals, the illite and kaolinite were hardly affected.
+# The kinetic dissolution rate of calcite and dolomite is mostly controlled by their saturation ratio states.
+# When the SR value of calcite and dolomite reaches 1, the dissolution process of the correspondences is stopped.
+#
+# **Note:**
+# The scenario #2 is not directly calculated in the notebook due to the lack of the large input file `llnl.dat` database in the commit. Also the post-processing of the scenario #2 is similar to which is showed in the scenario #1.
+#
+# 
+# Fig. 3: a) evolution of the calcium and magnesium molality during the $\ce{CO2}$ injection; b) L2 relative difference norm between the results computed by OGS-6 and PHREEQC.
+#
+# 
+# Fig. 4: Evolution of the dissoluted minerals molality and the related saturation ratio of calcite and dolomite over the time.
+#
+# ## Literature
+#
+# <!-- vale off -->
+#
+# [1] BGR, CO2LPIE project, 2023. URL: [https://www.bgr.bund.de/DE/Themen/Endlagerung/Projekte/Wirtsgesteine_geotechnische_Barrieren/laufend/Nur-Deutsch/mont_terri_experimente.html?nn=1542156](https://www.bgr.bund.de/DE/Themen/Endlagerung/Projekte/Wirtsgesteine_geotechnische_Barrieren/laufend/Nur-Deutsch/mont_terri_experimente.html?nn=1542156).
+#
+# [2] Appelo, C. A. J. and Postma, D. (2004). Geochemistry, groundwater and pollution. CRC press.
+#
+# [3] Lu, R., Nagel, T., Poonoosamy, J., Naumov, D., Fischer, T., Montoya, V., Kolditz, O., and Shao, H.
+# (2022). A new operator-splitting finite element scheme for reactive transport modeling in saturated porous
+# media. Computers and Geosciences, 163(April 2021):105106. URL: https://doi.org/10.1016/j.cageo.2022.105106.
+#
+# [4] Plummer, L. N., Wigley, T. M., and Parkhurst, D. L. (1978). KINETICS OF CALCITE DISSOLUTION
+# IN CO2-WATER SYSTEMS AT 5 degree TO 60 degree C AND 0. 0 TO 1. 0 ATM CO2.
+#
+# [5] Appelo, C. A. J. and Postma, D. (2023). Example 5.9: Kinetic dissolution of calcite. URL: https://www.hydrochemistry.eu/a&p/5/ex_5_9.phr.
+#
+# [6] Thury, M. (2002). The characteristics of the Opalinus Clay investigated in the Mont Terri underground
+# rock laboratory in Switzerland. doi:10.1016/S1631-0705(02)01372-5.
+#
+# [7] Palandri, J. L. and Kharaka, Y. K. (2004). A compilation of rate parameters of water-mineral interaction
+# kinetics for application to geochemical modeling.
+#
+# [8] Wersin, P., Mazurek, M., and Gimmi, T. (2022). Porewater chemistry of Opalinus Clay revisited:
+# Findings from 25 years of data collection at the Mont Terri Rock Laboratory. Applied Geochemistry,
+# 138(November 2021):105234. URL: https://doi.org/10.1016/j.apgeochem.2022.105234.
diff --git a/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.md b/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.md
deleted file mode 100644
index 35e8319099143cb464881126e3d8bf650b9e20e7..0000000000000000000000000000000000000000
--- a/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.md
+++ /dev/null
@@ -1,333 +0,0 @@
-+++
-author = "Thomas Fischer, Frieder Loer"
-date = "2019-10-01T11:54:02+01:00"
-title = "Heat conduction (Line Source Term)"
-weight = 121
-project = ["Parabolic/T/1D_line_source_term_tests/line_source_term.prj"]
-image = "temperature_distribution_line_source_term_in_cylinder.png"
-web_subsection = "heatconduction"
-+++
-
-
-## Equations
-
-We consider the Poisson equation:
-
-$$
-\begin{equation}
-\nabla\cdot(\nabla T) + Q_T = 0 \quad \text{in }\Omega
-\end{equation}$$
-w.r.t Dirichlet-type boundary conditions
-$$
-T(x) = 0 \quad \text{on }\Gamma_D
-$$
-
-where $T$ could be temperature, the subscripts $D$ denotes the Dirichlet-type
-boundary conditions.
-Here, the temperature distribution under the impact of a
-line shaped source term should be studied.
-
-## Problem Specifications and Analytical Solution
-
-In OGS there are several benchmarks for line source terms in 2d and 3d domains
-available.
-Here, some of the 3d benchmarks are described.
-
-### Cylindrical domain
-
-The Poisson equation on cylindrical domain of height $1$ and radius
-$r=1$ is solved.
-In the following figure the geometry, partly semi-transparent,
-is sketched.
-Furthermore, the mesh resolution is shown in the cylindrical domain
-within the first quadrant of the coordinate system.
-In the second quadrant the simulated temperature distribution is depicted.
-
-```python
-# Plot the cylindrical domain with pyvista
-
-import numpy as np
-import pyvista as pv
-pv.set_plot_theme("document")
-pv.set_jupyter_backend("static")
-
-mesh = pv.read("49k_prisms/Cylinder_r_1_h_1_prism_49k.vtu")
-plotter = pv.Plotter()
-
-# Create a dict for colorbar arguments
-sargs = dict(title = "Temperature", height=0.05, width=0.3, position_x=0.6, position_y=0.05)
-
-# Plot transparent part of mesh
-# Create clipping box
-clip_box0 = pv.Box([-1, 1, 0, 1, 0, 1])
-# Apply the clip filter to the mesh
-transparent_mesh = mesh.clip_box(clip_box0)
-# Add mesh to the plotter
-plotter.add_mesh(transparent_mesh, show_edges=False, opacity = 0.5, cmap="coolwarm", scalar_bar_args=sargs)
-
-# Plot Solid part of mesh
-clip_box = pv.Box([-1, 1, -1, 0, 0, 1])
-solid_mesh = mesh.clip_box(clip_box)
-plotter.add_mesh(solid_mesh, show_edges=False, cmap="coolwarm", opacity = 1, scalar_bar_args=sargs)
-
-# Plot grid lines in one quarter of cylinder
-clip_box1 = pv.Box([-1, 1, -1, 0, 0, 1])
-partial_mesh1 = mesh.clip_box(clip_box1)
-clip_box2 = pv.Box([-1, 0, -1, 1, 0, 1])
-partial_mesh2 = partial_mesh1.clip_box(clip_box2)
-plotter.add_mesh(partial_mesh2, show_edges=True, edge_color="mediumblue", cmap="coolwarm", scalar_bar_args=sargs)
-
-# Plot line
-start_point = [-1, 0, 0.5]
-end_point = [1, 0, 0.5]
-line = pv.Line(start_point, end_point)
-plotter.add_mesh(line, show_edges=True, color = "white", line_width = 5)
-
-# Set the camera's field of view
-camera = plotter.camera
-plotter.camera.focal_point = (0, 0, 0.5) # Set the focal point (where the camera is looking)
-plotter.camera.position = (-3, -4.5, 3.5) # Set the camera position (-2, -3, 2.5) *0.5 = (-1, -1.5, 1.25)
-
-plotter.show_grid()
-plotter.add_axes()
-plotter.window_size = [1500,1000]
-plotter.show()
-```
-
-The source term is defined along the line in the center of the cylinder:
-$$
-\begin{equation}
-Q(x) = 1 \quad \text{at } x=0, y=0.
-\end{equation}
-$$
-
-In the above figure the source term is the red vertical line in the origin of the coordinate system.
-
-The analytical solution for a line source in the cylinder is
-
-$$
-\begin{equation}
-T(x) = - \frac{1}{2 \pi} \ln \sqrt{x^2 + y^2}.
-\end{equation}
-$$
-
-```python
-# Define analytical solution
-def t_analytical(x, y):
- return -(1 / (2 * np.pi)) * np.log(np.sqrt(x**2 + y**2))
-```
-
-### Analytical solution in ParaView
-
-Since the analytical solution has a singularity at $(x, y) = (0, 0)$ the
-analytical solution in ParaView is generated as follows:
-
-```none
-if (coordsX^2<0.0001 & coordsY^2<0.0001, temperature, -1/(4*asin(1))*ln(sqrt(coordsX^2+coordsY^2)))
-```
-
-## Numerical simulation 286k mesh
-
-The applied mesh has a resolution of 286k cells.
-
-```python
-# Create output path if it doesn't exist yet
-import os
-from pathlib import Path
-
-out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
-if not out_dir.exists():
- out_dir.mkdir(parents=True)
-```
-
-```python
-# Import OGS class
-import ogstools as ot
-
-model = ot.Project(
- input_file="286k_prisms/line_source_term_in_cylinder.prj",
- output_file="286k_prisms/line_source_term_in_cylinder.prj"
- )
-model.run_model(logfile=f"{out_dir}/286k_out.txt", args=f"-o {out_dir} -m 286k_prisms")
-```
-
-## Results and evaluation
-
-The following plot shows the temperature along the white line in the figure above.
-
-```python
-import os
-import vtuIO
-import numpy as np
-import matplotlib.pyplot as plt
-
-# Extract values along a line in the domain
-
-# Load results
-pvdfile_286k = vtuIO.PVDIO(f"{out_dir}/3D_line_source_term_in_cylinder_286k.pvd", dim=3)
-
-# Get point field names
-fields = pvdfile_286k.get_point_field_names()
-
-# Extract values along line
-number_of_subdivisions = 801
-length = np.linspace(-1, 1, number_of_subdivisions)
-
-# Draws a line through the domain for sampling results
-z_axis = [(i, 0, 0.5) for i in length]
-
-# Extract timestep
-timestep = 1
-temp_286k = pvdfile_286k.read_set_data(timestep, "temperature", pointsetarray=z_axis)
-
-# Plot
-fig, ax = plt.subplots(1, 1, figsize=(10, 5), sharey=True)
-ax.plot(length, temp_286k, label="286k prism", color="tab:orange")
-ax.set_title("286k Prism Temperature")
-ax.set_xlabel("x")
-
-plt.show()
-```
-
-### Comparison with analytical solution
-
-The differences of analytical and computed solution are small outside of the center.
-
-```python
-# Combine analytical solution with numerical solution at singularity in (x,y)=(0,0)
-
-# Create reference for error calculation
-
-# Replace 0 with 1 to prevent division by 0, the respective element will be replaced with the numerical solution anyway
-length_replaced = length.copy()
-length_replaced[int((number_of_subdivisions-1)/2)] = 1
-
-# Replace diverging analytical solution in respective interval below a threshold of 0.01
-threshold = 0.01
-below_threshold = np.where(np.abs(length) < threshold)
-analytical286 = t_analytical(length_replaced, 0)
-analytical286[below_threshold[0][0]:(below_threshold[0][-1]+1)] = temp_286k[below_threshold[0][0]:(below_threshold[0][-1]+1)]
-
-# Plot absolute error
-fig, ax = plt.subplots(1, 1, figsize=(10, 5), sharey=True)
-abs_error286 = temp_286k - analytical286
-ax.plot(length, abs_error286)
-ax.grid(True)
-ax.set_title("286k Prism Absolute Error")
-ax.set_xlabel("x")
-
-plt.show()
-```
-
-Due to the numerical evaluation of the relative error of the computed solution the error grows in the vicinity of the boundary and in the center.
-
-```python
-# Relative error
-fig, ax = plt.subplots(1, 1, figsize=(10, 5))
-rel_error286 = (analytical286[1:-1] - temp_286k[1:-1]) / analytical286[1:-1]
-ax.plot(length[1:-1], rel_error286)
-ax.grid(True)
-ax.set_title("286k Prism Relative Error")
-ax.set_xlabel("x")
-
-plt.show()
-```
-
-### Input files
-
-The project file for the described model is [286k.prj](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/286k_prisms/line_source_term_in_cylinder.prj).
-The project file describes the processes to be solved and the related process variables together with their initial and boundary conditions as well as the source terms.
-
-The input mesh is stored in the VTK file format and can be directly visualized in ParaView for example.
-
-### Cylindrical domain - axisymmetric example
-
-The Poisson equation on cylindrical domain of height $1$ and radius $r=1$ is solved.
-The cylindrical domain is defined as axisymmetric.
-
-#### Numerical simulation
-
-```python
-model = ot.Project(
- input_file="../3D_line_source_term_in_cylinder_axisymmetric/line_source_term_in_cylinder.prj",
- output_file="../3D_line_source_term_in_cylinder_axisymmetric/line_source_term_in_cylinder.prj"
-)
-model.run_model(
- logfile=f"{out_dir}/axisym_out.txt",
- args=f"-o {out_dir} -m ../3D_line_source_term_in_cylinder_axisymmetric",
-)
-```
-
-#### Results and evaluation
-
-```python
-# Plot cylinder cross-section
-
-# Load mesh and plot it
-mesh = pv.read("../3D_line_source_term_in_cylinder_axisymmetric/square_1x1_quad_100x100.vtu")
-plotter = pv.Plotter()
-sargs = dict(title = "Temperature", height=0.05, width=0.4, position_x=0.3, position_y=0.05)
-plotter.add_mesh(mesh, show_edges=False, cmap="coolwarm", scalar_bar_args=sargs)
-
-# Plot line
-start_point = [0, 0.5, 0]
-end_point = [1, 0.5, 0]
-line = pv.Line(start_point, end_point)
-plotter.add_mesh(line, show_edges=True, color = "white", line_width = 1)
-
-plotter.add_axes()
-plotter.show_bounds(mesh, ticks = "both")
-plotter.window_size = [1500, 1000]
-plotter.view_xy()
-
-# Add text annotations for axis labels
-plotter.add_text("X", position=(1100, 140, 0), font_size=16, color="k")
-plotter.add_text("Y", position=(400, 850, 0), font_size=16, color="k")
-plotter.show()
-```
-
-The above figure shows the computed temperature distribution.
-
-The following plot shows the temperature along the white line in the figure above.
-
-```python
-pvdfile_ax = vtuIO.PVDIO(f"{out_dir}/3D_line_source_term_in_cylinder.pvd", dim=3)
-
-# Extract values along line
-# Space axis for plotting
-length = np.linspace(0, 1, 101)
-
-# Draws a line through the domain for sampling results
-z_axis = [(i, 0.5, 0) for i in length]
-
-timestep = 1
-temp_ax = pvdfile_ax.read_set_data(timestep, "temperature", pointsetarray=z_axis)
-
-plt.plot(length, temp_ax)
-plt.title("Temperature")
-plt.xlabel("x")
-plt.ylabel("Temperature (°C)")
-plt.show()
-```
-
-The error and relative error shows the same behaviour like in the simulation models above.
-Outside of the center, that has a singularity in the analytical solution, the errors decreases very fast.
-
-```python
-# Create the reference temperature and combine analytical solution with numerical solution at singularity in (x,y)=(0,0)
-t_ref = t_analytical(length[1:], 0)
-t_ref = np.insert(t_ref, 0, temp_ax[0])
-
-plt.plot(length, t_ref - temp_ax, label="Absolute error")
-plt.plot(length[:-1], (t_ref[:-1] - temp_ax[:-1]) / t_ref[:-1], label="Relative error")
-plt.title("Errors")
-plt.ylabel("Error")
-plt.xlabel("x")
-plt.legend()
-plt.show()
-```
-
-### Input files
-
-The project file for the described model is
-[line_source_term_in_cylinder.prj](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder_axisymmetric/line_source_term_in_cylinder.prj).
diff --git a/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.py b/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.py
new file mode 100644
index 0000000000000000000000000000000000000000..5e7a44f810e92b855cdb8f408054cdbec1f4aaaa
--- /dev/null
+++ b/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/heatconduction-line_source_term.py
@@ -0,0 +1,371 @@
+# %% [markdown]
+# +++
+# author = "Thomas Fischer, Frieder Loer"
+# date = "2019-10-01T11:54:02+01:00"
+# title = "Heat conduction (Line Source Term)"
+# weight = 121
+# project = ["Parabolic/T/1D_line_source_term_tests/line_source_term.prj"]
+# image = "temperature_distribution_line_source_term_in_cylinder.png"
+# web_subsection = "heatconduction"
+# +++
+
+# %% [markdown]
+# ## Equations
+#
+# We consider the Poisson equation:
+#
+# $$
+# \begin{equation}
+# \nabla\cdot(\nabla T) + Q_T = 0 \quad \text{in }\Omega
+# \end{equation}$$
+# w.r.t Dirichlet-type boundary conditions
+# $$
+# T(x) = 0 \quad \text{on }\Gamma_D
+# $$
+#
+# where $T$ could be temperature, the subscripts $D$ denotes the Dirichlet-type
+# boundary conditions.
+# Here, the temperature distribution under the impact of a
+# line shaped source term should be studied.
+#
+# ## Problem Specifications and Analytical Solution
+#
+# In OGS there are several benchmarks for line source terms in 2d and 3d domains
+# available.
+# Here, some of the 3d benchmarks are described.
+#
+# ### Cylindrical domain
+#
+# The Poisson equation on cylindrical domain of height $1$ and radius
+# $r=1$ is solved.
+# In the following figure the geometry, partly semi-transparent,
+# is sketched.
+# Furthermore, the mesh resolution is shown in the cylindrical domain
+# within the first quadrant of the coordinate system.
+# In the second quadrant the simulated temperature distribution is depicted.
+
+# %%
+# Plot the cylindrical domain with pyvista
+
+import os
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import ogstools as ot # Import ogstools
+import pyvista as pv
+import vtuIO
+
+# %%
+pv.set_plot_theme("document")
+pv.set_jupyter_backend("static")
+
+mesh = pv.read("49k_prisms/Cylinder_r_1_h_1_prism_49k.vtu")
+plotter = pv.Plotter()
+
+# Create a dict for colorbar arguments
+sargs = {
+ "title": "Temperature",
+ "height": 0.05,
+ "width": 0.3,
+ "position_x": 0.6,
+ "position_y": 0.05,
+}
+
+# Plot transparent part of mesh
+# Create clipping box
+clip_box0 = pv.Box([-1, 1, 0, 1, 0, 1])
+# Apply the clip filter to the mesh
+transparent_mesh = mesh.clip_box(clip_box0)
+# Add mesh to the plotter
+plotter.add_mesh(
+ transparent_mesh,
+ show_edges=False,
+ opacity=0.5,
+ cmap="coolwarm",
+ scalar_bar_args=sargs,
+)
+
+# Plot Solid part of mesh
+clip_box = pv.Box([-1, 1, -1, 0, 0, 1])
+solid_mesh = mesh.clip_box(clip_box)
+plotter.add_mesh(
+ solid_mesh, show_edges=False, cmap="coolwarm", opacity=1, scalar_bar_args=sargs
+)
+
+# Plot grid lines in one quarter of cylinder
+clip_box1 = pv.Box([-1, 1, -1, 0, 0, 1])
+partial_mesh1 = mesh.clip_box(clip_box1)
+clip_box2 = pv.Box([-1, 0, -1, 1, 0, 1])
+partial_mesh2 = partial_mesh1.clip_box(clip_box2)
+plotter.add_mesh(
+ partial_mesh2,
+ show_edges=True,
+ edge_color="mediumblue",
+ cmap="coolwarm",
+ scalar_bar_args=sargs,
+)
+
+# Plot line
+start_point = [-1, 0, 0.5]
+end_point = [1, 0, 0.5]
+line = pv.Line(start_point, end_point)
+plotter.add_mesh(line, show_edges=True, color="white", line_width=5)
+
+# Set the camera's field of view
+camera = plotter.camera
+plotter.camera.focal_point = (
+ 0,
+ 0,
+ 0.5,
+) # Set the focal point (where the camera is looking)
+plotter.camera.position = (
+ -3,
+ -4.5,
+ 3.5,
+) # Set the camera position (-2, -3, 2.5) *0.5 = (-1, -1.5, 1.25)
+
+plotter.show_grid()
+plotter.add_axes()
+plotter.window_size = [1500, 1000]
+plotter.show()
+
+
+# %% [markdown]
+# The source term is defined along the line in the center of the cylinder:
+# $$
+# \begin{equation}
+# Q(x) = 1 \quad \text{at } x=0, y=0.
+# \end{equation}
+# $$
+#
+# In the above figure the source term is the red vertical line in the origin of the coordinate system.
+#
+# The analytical solution for a line source in the cylinder is
+#
+# $$
+# \begin{equation}
+# T(x) = - \frac{1}{2 \pi} \ln \sqrt{x^2 + y^2}.
+# \end{equation}
+# $$
+
+
+# %%
+# Define analytical solution
+def t_analytical(x, y):
+ return -(1 / (2 * np.pi)) * np.log(np.sqrt(x**2 + y**2))
+
+
+# %% [markdown]
+# ### Analytical solution in ParaView
+#
+# Since the analytical solution has a singularity at $(x, y) = (0, 0)$ the
+# analytical solution in ParaView is generated as follows:
+#
+# ```none
+# if (coordsX^2<0.0001 & coordsY^2<0.0001, temperature, -1/(4*asin(1))*ln(sqrt(coordsX^2+coordsY^2)))
+# ```
+#
+# ## Numerical simulation 286k mesh
+#
+# The applied mesh has a resolution of 286k cells.
+
+# %%
+# Create output path if it doesn't exist yet
+out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
+if not out_dir.exists():
+ out_dir.mkdir(parents=True)
+
+# %%
+model = ot.Project(
+ input_file="286k_prisms/line_source_term_in_cylinder.prj",
+ output_file="286k_prisms/line_source_term_in_cylinder.prj",
+)
+model.run_model(logfile=f"{out_dir}/286k_out.txt", args=f"-o {out_dir} -m 286k_prisms")
+
+# %% [markdown]
+# ## Results and evaluation
+#
+# The following plot shows the temperature along the white line in the figure above.
+
+# %%
+# Extract values along a line in the domain
+
+# Load results
+pvdfile_286k = vtuIO.PVDIO(f"{out_dir}/3D_line_source_term_in_cylinder_286k.pvd", dim=3)
+
+# Get point field names
+fields = pvdfile_286k.get_point_field_names()
+
+# Extract values along line
+number_of_subdivisions = 801
+length = np.linspace(-1, 1, number_of_subdivisions)
+
+# Draws a line through the domain for sampling results
+z_axis = [(i, 0, 0.5) for i in length]
+
+# Extract timestep
+timestep = 1
+temp_286k = pvdfile_286k.read_set_data(timestep, "temperature", pointsetarray=z_axis)
+
+# Plot
+fig, ax = plt.subplots(1, 1, figsize=(10, 5), sharey=True)
+ax.plot(length, temp_286k, label="286k prism", color="tab:orange")
+ax.set_title("286k Prism Temperature")
+ax.set_xlabel("x")
+
+plt.show()
+
+# %% [markdown]
+# ### Comparison with analytical solution
+#
+# The differences of analytical and computed solution are small outside of the center.
+
+# %%
+# Combine analytical solution with numerical solution at singularity in (x,y)=(0,0)
+
+# Create reference for error calculation
+
+# Replace 0 with 1 to prevent division by 0, the respective element will be replaced with the numerical solution anyway
+length_replaced = length.copy()
+length_replaced[int((number_of_subdivisions - 1) / 2)] = 1
+
+# Replace diverging analytical solution in respective interval below a threshold of 0.01
+threshold = 0.01
+below_threshold = np.where(np.abs(length) < threshold)
+analytical286 = t_analytical(length_replaced, 0)
+analytical286[below_threshold[0][0] : (below_threshold[0][-1] + 1)] = temp_286k[
+ below_threshold[0][0] : (below_threshold[0][-1] + 1)
+]
+
+# Plot absolute error
+fig, ax = plt.subplots(1, 1, figsize=(10, 5), sharey=True)
+abs_error286 = temp_286k - analytical286
+ax.plot(length, abs_error286)
+ax.grid(True)
+ax.set_title("286k Prism Absolute Error")
+ax.set_xlabel("x")
+
+plt.show()
+
+# %% [markdown]
+# Due to the numerical evaluation of the relative error of the computed solution the error grows in the vicinity of the boundary and in the center.
+
+# %%
+# Relative error
+fig, ax = plt.subplots(1, 1, figsize=(10, 5))
+rel_error286 = (analytical286[1:-1] - temp_286k[1:-1]) / analytical286[1:-1]
+ax.plot(length[1:-1], rel_error286)
+ax.grid(True)
+ax.set_title("286k Prism Relative Error")
+ax.set_xlabel("x")
+
+plt.show()
+
+# %% [markdown]
+# ### Input files
+#
+# The project file for the described model is [286k.prj](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder/286k_prisms/line_source_term_in_cylinder.prj).
+# The project file describes the processes to be solved and the related process variables together with their initial and boundary conditions as well as the source terms.
+#
+# The input mesh is stored in the VTK file format and can be directly visualized in ParaView for example.
+#
+# ### Cylindrical domain - axisymmetric example
+#
+# The Poisson equation on cylindrical domain of height $1$ and radius $r=1$ is solved.
+# The cylindrical domain is defined as axisymmetric.
+#
+# #### Numerical simulation
+
+# %%
+model = ot.Project(
+ input_file="../3D_line_source_term_in_cylinder_axisymmetric/line_source_term_in_cylinder.prj",
+ output_file="../3D_line_source_term_in_cylinder_axisymmetric/line_source_term_in_cylinder.prj",
+)
+model.run_model(
+ logfile=f"{out_dir}/axisym_out.txt",
+ args=f"-o {out_dir} -m ../3D_line_source_term_in_cylinder_axisymmetric",
+)
+
+# %% [markdown]
+# #### Results and evaluation
+
+# %%
+# Plot cylinder cross-section
+
+# Load mesh and plot it
+mesh = pv.read(
+ "../3D_line_source_term_in_cylinder_axisymmetric/square_1x1_quad_100x100.vtu"
+)
+plotter = pv.Plotter()
+sargs = {
+ "title": "Temperature",
+ "height": 0.05,
+ "width": 0.4,
+ "position_x": 0.3,
+ "position_y": 0.05,
+}
+plotter.add_mesh(mesh, show_edges=False, cmap="coolwarm", scalar_bar_args=sargs)
+
+# Plot line
+start_point = [0, 0.5, 0]
+end_point = [1, 0.5, 0]
+line = pv.Line(start_point, end_point)
+plotter.add_mesh(line, show_edges=True, color="white", line_width=1)
+
+plotter.add_axes()
+plotter.show_bounds(mesh, ticks="both")
+plotter.window_size = [1500, 1000]
+plotter.view_xy()
+
+# Add text annotations for axis labels
+plotter.add_text("X", position=(1100, 140, 0), font_size=16, color="k")
+plotter.add_text("Y", position=(400, 850, 0), font_size=16, color="k")
+plotter.show()
+
+# %% [markdown]
+# The above figure shows the computed temperature distribution.
+#
+# The following plot shows the temperature along the white line in the figure above.
+
+# %%
+pvdfile_ax = vtuIO.PVDIO(f"{out_dir}/3D_line_source_term_in_cylinder.pvd", dim=3)
+
+# Extract values along line
+# Space axis for plotting
+length = np.linspace(0, 1, 101)
+
+# Draws a line through the domain for sampling results
+z_axis = [(i, 0.5, 0) for i in length]
+
+timestep = 1
+temp_ax = pvdfile_ax.read_set_data(timestep, "temperature", pointsetarray=z_axis)
+
+plt.plot(length, temp_ax)
+plt.title("Temperature")
+plt.xlabel("x")
+plt.ylabel("Temperature (°C)")
+plt.show()
+
+# %% [markdown]
+# The error and relative error shows the same behaviour like in the simulation models above.
+# Outside of the center, that has a singularity in the analytical solution, the errors decreases very fast.
+
+# %%
+# Create the reference temperature and combine analytical solution with numerical solution at singularity in (x,y)=(0,0)
+t_ref = t_analytical(length[1:], 0)
+t_ref = np.insert(t_ref, 0, temp_ax[0])
+
+plt.plot(length, t_ref - temp_ax, label="Absolute error")
+plt.plot(length[:-1], (t_ref[:-1] - temp_ax[:-1]) / t_ref[:-1], label="Relative error")
+plt.title("Errors")
+plt.ylabel("Error")
+plt.xlabel("x")
+plt.legend()
+plt.show()
+
+# %% [markdown]
+# ### Input files
+#
+# The project file for the described model is
+# [line_source_term_in_cylinder.prj](https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Parabolic/T/3D_line_source_term_tests/3D_line_source_term_in_cylinder_axisymmetric/line_source_term_in_cylinder.prj).
diff --git a/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.md b/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.md
deleted file mode 100644
index 43ec04837f2e68c3fa07081e7b93247d201b5b2a..0000000000000000000000000000000000000000
--- a/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.md
+++ /dev/null
@@ -1,427 +0,0 @@
-+++
-title = "Wellbore Heat Transport - EUBHE"
-date = "2023-08-03"
-author = "Chaofan Chen, Haibing Shao, Frieder Loer"
-image = "figures/pipe_flow_3d_model.png"
-web_subsection = "heat-transport-bhe"
-project = ["Parabolic/T/BHE_1P/BHE_1P.prj"]
-+++
-
-
-## Problem description
-
-Ramey et al. (1962) proposed the analytical solution concerning the wellbore heat transmission, which can be used to quantify the fluid temperature change in the wellbore.
-In order to verify the single pipe flow model in the OGS, the numerical results were compared with Ramey's analytical solution Ramey et al. (1962).
-The detailed calculation of the Ramey's analytical solution is given below.
-A detailed analysis of an enhanced U-tube borehole heat exchanger (EUBHE) can be found in Chen, C. et al. (2021).
-
-## Model Setup
-
-In this benchmark, the length of the wellbore is 30 m as shown in Figure 1 and the cold water is injected into the inlet point of the wellbore with the constant temperature of 20°C.
-The initial temperature of the fluid and grout in the wellbore is 20 $^{\circ}$C, and temperature of the surrounding rock is 55 $^{\circ}$C.
-The wellbore and pipe diameter are 0.28 m and 0.25826 m, respectively.
-And the flow rate is 0.0002 $m^3/s$.
-
-
-
-Figure 1: Single pipe flow model
-
-```python
-# Import required libraries
-
-import numpy as np
-import matplotlib.pyplot as plt
-import pandas as pd
-```
-
-## Ramey's analytical solution
-
-The equations are followed by their implementation in python.
-
-```python
-# Define constants
-
-T_d = 55 # Temperature of undistributed formation
-T_i = 20 # Temperature of inlet fluid
-Z = np.insert(np.linspace(0.3,30,100), 0, 1e-10) # Depth below surface, start with very small value to prevent division by zero
-q = 0.0002 # Flow rate of fluid in wellbore, m3/s
-rho_f = 1000 # Density of fluid, unit kg/m3
-c_p_f = 4190 # Specific heat capacity of fluid in wellbore
-mu_f = 1.14e-3 # Unit kg/m/s;
-lambda_f = 0.59 # Unit W/m/K
-rho_re = 1800 # Density of geothermal reservoir
-c_p_re = 1778 # Specific heat capacity of geothermal reservoir
-lambda_re = 2.78018 # Heat conductivity coefficient of geothermal reservoir
-lambda_g = 0.73 # Heat conductivity coefficient of grout and pipe
-lambda_pi = 1.3
-r_pi = 0.12913 # Inner radius of pipe and wellbore
-r_b = 0.14
-t_pi = 0.00587 # Thickness of the pipe
-t = 86400 * 5 # Operation time
-```
-
-In Ramey's analytical solution Ramey et al. (1962), the outlet temperature of the pipe inside the wellbore can be calculated by
-
-$$
- T_o(t) = T_{s} + (T_i(t) - T_{s})\exp(-\Delta z/X)
-$$
-
-```python
-def temp(T_d, T_i, delta_z, X):
- return T_d + (T_i - T_d) * np.exp(-delta_z / X)
-```
-
-The coefficient $X$ is determined by
-
-$$
- X = \frac{q\rho_fc_{p,f}(\lambda_{s}+r_pUf(t))}{2\pi r_pU \lambda_{s}}
-$$
-
-where $q$ is the flow rate of the fluid in the wellbore.
-
-```python
-def coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t):
- return (q * rho_f * c_p_f) * (lambda_re + r_pi * U * f_t) / (2 * np.pi * r_pi * U * lambda_re)
-```
-
-With dimensionless time
-
-$$
- t_D = \frac{\lambda_{s}t}{(\rho_{s}c_{p,s}r_b^{2})}
-$$
-
-```python
-def dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b):
- return lambda_re * delta_t / (rho_re * c_p_re * r_b * r_b)
-```
-
-the time function $f(t)$ can be calculated by
-
-<!-- markdownlint-disable reference-links-images -->
-
-$$
- f(t) = [0.4063+0.5\ln(t_D)][1+\frac{0.6}{t_D}], t_D > 1.5
-$$
-
-<!-- markdownlint-enable reference-links-images -->
-
-$$
- f(t) = 1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}), t_D \leqslant 1.5
-$$
-
-```python
-def time_function(t_D):
- if t_D > 1.5:
- f_t = (0.4063 + 0.5 * np.log(t_D)) * (1 + 0.6 / t_D)
- else:
- f_t = 1.1281 * np.sqrt(t_D) * (1 - 0.3 * np.sqrt(t_D))
-
- return f_t
-
-```
-
-The Prandtl and Reynolds number can be calculated as follows
-
-$$
- \mathrm{Pr} = \frac{\mu_f c_{p,f}}{\lambda_f}
-$$
-
-$$
- \mathrm{Re} = \frac{\rho_f v d_{pi}}{\mu_f}
-$$
-
-where, $\mu_f, \rho_f$ and $\lambda_f$ is the fluid viscosity, density and thermal conductivity.
-
-```python
-# Calculate the convective film resistance inside pipe
-v = q / (np.pi * r_pi * r_pi)
-
-# Reynolds number
-Pr = mu_f * c_p_f / lambda_f
-Re = rho_f * v * (2 * r_pi) / mu_f
-
-```
-
-The Nusselt number can be determined by the following equation Diersch, (2011):
-
-$$
- \mathrm{Nu} = 4.364,\ \mathrm{Re} < 2300
-$$
-
-$$
- \mathrm{Nu} = \frac{(\xi_{k}/8)\ \mathrm{Re}_{k}\ \mathrm{Pr}}{1+12.7\sqrt{{\xi_k}/8}(\mathrm{Pr}^{2/3}-1)} [ 1+(\frac{{d_k}^{i}}{L})^{2/3}], Re \geq 10^4
-$$
-
-$$
- \mathrm{Nu} = (1-\gamma_{k})\ 4.364 + \gamma_{k} ( \frac{(0.0308/8)10^{4}\mathrm{Pr}}{1+12.7\ \sqrt{0.0308/8}(\mathrm{Pr}^{2/3}-1)} \left[ 1+\left(\frac{d_{k}^{i}}{L}\right)^{2/3} \right] ), 2 300 \leq \mathrm{Re} < 10^{4}
-$$
-
-with
-
-$$
- \xi_{k}=(1.8\ \log_{10}\mathrm{Re}_{k}-1.5)^{-2}
-$$
-
-and
-
-$$
- \gamma_{k} = \frac{\mathrm{Re}_{k}-2300}{10^{4}-2300}
-$$
-
-```python
-# Estimate convective film coefficient
-xi = (1.8 * np.log(Re) - 1.5)**(-2)
-Nu_l = (xi / 8.0 * Re * Pr) / (1.0 + 12.7 * np.sqrt(xi / 8.0) * ((Pr)**(2.0 / 3.0) - 1.0)) * (1.0 + (r_pi / Z)**(2.0 / 3.0))
-Nu_lam = 4.364
-gamma = (Re - 2300) / (10000 - 2300)
-
-Nu_m = (1.0 - gamma) * 4.364 + gamma * ((0.0308 / 8.0 * 1.0e4 * Pr) / (1.0 + 12.7 * np.sqrt(0.0308 / 8.0) * ((Pr)**(2.0 / 3.0) - 1.0)) * (1.0 + (r_pi / Z)**(2.0 / 3.0)))
-
-# Evaluate Nusselt number; based on the value of Re choose appropriate correlation
-if Re < 2300:
- Nu_p = Nu_lam
-elif Re > 10000:
- Nu_p = Nu_l
-else:
- Nu_p = Nu_m
-```
-
-and the overall heat transfer coefficient $U$ is written as follows,
-
-$$
- U = [\frac{r_{pi}+t_{pi}}{r_{pi}h}+(r_{pi}+t_{pi})(\frac{\ln{\frac{r_{pi}+t_{pi}}{r_{pi}}}}{\lambda_{pi}}+\frac{\ln{\frac{r_b}{r_{pi}+t_{pi}}}}{\lambda_{grout}})]^{-1}
-$$
-
-with
-
-$$
- h = \frac{\lambda_f Nu}{2r_{pi}}
-$$
-
-```python
-h = lambda_f * Nu_p / (2 * r_pi) # Unit: W/m2/K
-
-U = 1 / (((r_pi + t_pi) / (r_pi * h)) + (r_pi + t_pi) * (np.log((r_pi + t_pi) / r_pi) / lambda_pi + np.log(r_b / (r_pi + t_pi)) / lambda_g))
-```
-
-The friction factor $f$, is evaluated by Churchill correlation Churchill et al. (1977),
-
-$$
- f = \frac{1}{(\frac{1}{[((\frac{8}{Re})^{10}+(\frac{Re}{36500})^{20})]^{1/2}}+[2.21(\ln{\frac{Re}{7}})]^{10})^{1/5}}
-$$
-
-```python
-# Churchill correlation for friction factor
-f = 1 / ((1 / (np.sqrt((8 / Re) ** 10 + (Re / 36500) ** 20))) + (2.21 * np.log(Re / 7)) ** 10) ** (1 / 5)
-```
-
-With these equations the analytical solution can be computed:
-
-```python
-# Compute solution
-
-time_range = range(0, t+1, 43200)
-
-# Create a dataframe to store results for different z- and t-steps
-row_names = Z
-column_names = time_range
-
-
-df = pd.DataFrame(index=row_names, columns=column_names)
-
-
-for delta_z in Z:
- T_0 = []
- for delta_t in time_range:
-
- # Compute dimensionless time
- t_D = dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b)
-
- # Compute time function
- f_t = time_function(t_D)
-
- # Compute X
- X = coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t)
-
- # Compute outlet temperature
- T_0.append(temp(T_d, T_i, delta_z, X))
-
- df.loc[delta_z] = T_0
-```
-
-## Numerical solution
-
-```python
-# Create output path
-import os
-from pathlib import Path
-
-out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
-if not out_dir.exists():
- out_dir.mkdir(parents=True)
-```
-
-```python
-import ogstools as ot
-
-# Create an OGS-object
-# Pass the project file and set an output file path
-model=ot.Project(input_file="BHE_1P.prj", output_file=f"{out_dir}/BHE_1P.prj")
-
-# set end time (in seconds) and the time increments
-model.replace_text(5 * 24 * 60 * 60, xpath="./time_loop/processes/process/time_stepping/t_end")
-model.replace_text(30, xpath="./time_loop/processes/process/time_stepping/timesteps/pair/repeat")
-model.replace_text(4 * 60 * 60, xpath="./time_loop/processes/process/time_stepping/timesteps/pair/delta_t")
-# Write to the output file
-model.write_input()
-
-# Run OGS
-model.run_model(logfile=os.path.join(out_dir, "log.txt"), args=f"-o {out_dir} -m .")
-```
-
-## Results and discussion
-
-The outlet temperature change over time was compared against the analytical solution and presented in Figure 2.
-After 30 days, the fluid temperature distribution in the wellbore is shown in Figure 3.
-The maximum relative error between the numerical model and Ramey's analytical solution is less than 0.15%.
-
-In numerical model, the outlet temperature at beginning stage is affected by the initial temperature in the pipe inside the wellbore.
-The initial fluid temperature set in the benchmark means there is water with 20°C filled in the wellbore already before injecting water into the wellbore.
-But in the analytical solution, no initial temperature is set and the temperature keeps equilibrium state at every moment.
-The impact of initial temperature condition in numerical model is decreasing with the increasing of the operational time as shown in Figure 2.
-
-```python
-import vtuIO
-import numpy as np
-import os
-
-out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
-
-# Load the result
-pvdfile = vtuIO.PVDIO(f"{out_dir}/BHE_1P.pvd", dim=3)
-
-# Get point field names
-fields = pvdfile.get_point_field_names()
-print(fields)
-
-# Get all written timesteps
-time = pvdfile.timesteps
-
-# Define observation points
-observation_points = {'outlet':(0.0,10.0,-30.0)}
-
-# Read the time series of point field 'temperature_BHE1'
-temperature_BHE1 = pvdfile.read_time_series('temperature_BHE1', observation_points)
-```
-
-```python
-# Compute analytical solution again but at timesteps of numerical solution for error calculation
-for delta_z in Z:
- T_0 = []
- for delta_t in time:
-
- # Compute dimensionless time
- t_D = dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b)
-
- # Compute time function
- f_t = time_function(t_D)
-
- # Compute X
- X = coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t)
-
- # Compute outlet temperature
- T_0.append(temp(T_d, T_i, delta_z, X))
-```
-
-```python
-# Plot numerical vs. analytical result
-fig, ax1 = plt.subplots(figsize=(10, 8))
-
-x_pos = np.arange(0, 432000, 60*60*24)
-x_ticks = np.arange(0,5,1)
-
-
-ax1.plot(time_range, df.iloc[-1, :], 'k.', markersize=10, markerfacecolor='none', label = 'Ramey\'s analytical solution')
-ax1.plot(time, temperature_BHE1['outlet'][:,0], 'r', label = 'OGS')
-ax1.set_ylabel('Outlet temperature (°C)', fontsize=20)
-ax1.set_xlabel('Time (d)', fontsize=20)
-ax1.set_xticks(x_pos, x_ticks)
-ax1.spines['left'].set_linewidth(2)
-ax1.spines['bottom'].set_linewidth(2)
-ax1.tick_params(axis='both', labelsize=16)
-
-color = 'blue'
-ax2 = ax1.twinx()
-ax2.plot(time ,temperature_BHE1['outlet'][:,0]-T_0, color=color, label = 'Absolute error')
-ax2.set_ylabel('Absolute error (°C)', color=color, fontsize = 20)
-ax2.tick_params(axis='y', labelcolor=color)
-ax2.spines['right'].set_color('blue')
-ax2.spines['right'].set_linewidth(2)
-ax2.tick_params(axis='y', labelsize=16)
-
-plt.figlegend(fontsize=18, loc = (0.4, 0.1), frameon=False)
-fig.tight_layout()
-
-plt.show()
-```
-
-Figure 2: Comparison with analytical solution results
-
-```python
-# Extract values along a line in the domain
-
-# Space axis for plotting
-length = np.linspace(0, -30, 101)
-
-# Draws a line through the domain for sampling results
-x_axis=[(0,10,i) for i in length]
-
-# At timestep i
-i = -1
-temp = pvdfile.read_set_data(time[i], "temperature_BHE1", pointsetarray=x_axis)
-
-# Plot
-fig, ax1 = plt.subplots(figsize=(10, 8))
-ax1.set_xlabel('Distance (m)', fontsize = 20)
-ax1.set_ylabel('Temperature (°C)', fontsize = 20)
-ax1.plot(np.linspace(0,30,101), temp[:,0], 'r', linewidth = 2, label = 'OGS')
-ax1.plot(Z,df.iloc[:, -1],'k.', markersize=10, markerfacecolor='none', label = 'Ramey\'s analytical solution')
-ax1.set_ylim(19,25)
-ax1.set_xlim(0,30)
-ax1.spines['left'].set_linewidth(2)
-ax1.spines['bottom'].set_linewidth(2)
-ax1.tick_params(axis='both', labelsize=16)
-
-color = 'blue'
-ax2 = ax1.twinx()
-ax2.set_ylabel('Absolute error (°C)', color=color, fontsize = 20)
-ax2.plot(np.linspace(0,30,101), temp[:,0]-df.iloc[:, -1], linewidth = 2, color=color, label = 'Absolute error')
-ax2.tick_params(axis='y', labelcolor=color)
-ax2.spines['right'].set_color('blue')
-ax2.spines['right'].set_linewidth(2)
-ax2.tick_params(axis='y', labelsize=16)
-
-fig.tight_layout()
-plt.figlegend(fontsize=18, loc = (0.4, 0.1), frameon=False)
-plt.show()
-
-max_error = 0.08
-if np.max(np.abs(temp[:,0]-df.iloc[:, -1])) > max_error:
- raise RuntimeError(f'The maximum temperature error is over {max_error} °C, which is higher than expected')
-```
-
-Figure 3: Distributed temperature of fluid and absolute error.
-
-## References
-
-<!-- vale off -->
-
-[1] Ramey Jr, H. J. (1962). Wellbore heat transmission. Journal of petroleum Technology, 14(04), 427-435.
-
-[2] Diersch, H-JG, et al. (2011). "Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals." Computers & Geosciences 37.8, 1122-1135.
-
-[3] Chaofan Chen, Wanlong Cai, Dmitri Naumov, Kun Tu, Hongwei Zhou, Yuping Zhang, Olaf Kolditz, Haibing Shao (2021). Numerical investigation on the capacity and efficiency of a deep enhanced U-tube borehole heat exchanger system for building heating. Renewable Energy, 169, 557-572.
-
-[4] Churchill, S. W. (1977). Comprehensive correlating equations for heat, mass and momentum transfer in fully developed flow in smooth tubes. Industrial & Engineering Chemistry Fundamentals, 16(1), 109-116.
-
-<!-- vale on -->
diff --git a/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.py b/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.py
new file mode 100644
index 0000000000000000000000000000000000000000..5410ee078467dcabe4f07409034426468efa8578
--- /dev/null
+++ b/Tests/Data/Parabolic/T/BHE_1P/pipe_flow_ebhe.py
@@ -0,0 +1,480 @@
+# %% [markdown]
+# +++
+# title = "Wellbore Heat Transport - EUBHE"
+# date = "2023-08-03"
+# author = "Chaofan Chen, Haibing Shao, Frieder Loer"
+# image = "figures/pipe_flow_3d_model.png"
+# web_subsection = "heat-transport-bhe"
+# project = ["Parabolic/T/BHE_1P/BHE_1P.prj"]
+# +++
+
+# %% [markdown]
+# ## Problem description
+#
+# Ramey et al. (1962) proposed the analytical solution concerning the wellbore heat transmission, which can be used to quantify the fluid temperature change in the wellbore.
+# In order to verify the single pipe flow model in the OGS, the numerical results were compared with Ramey's analytical solution Ramey et al. (1962).
+# The detailed calculation of the Ramey's analytical solution is given below.
+# A detailed analysis of an enhanced U-tube borehole heat exchanger (EUBHE) can be found in Chen, C. et al. (2021).
+#
+# ## Model Setup
+#
+# In this benchmark, the length of the wellbore is 30 m as shown in Figure 1 and the cold water is injected into the inlet point of the wellbore with the constant temperature of 20°C.
+# The initial temperature of the fluid and grout in the wellbore is 20 $^{\circ}$C, and temperature of the surrounding rock is 55 $^{\circ}$C.
+# The wellbore and pipe diameter are 0.28 m and 0.25826 m, respectively.
+# And the flow rate is 0.0002 $m^3/s$.
+#
+# 
+#
+# Figure 1: Single pipe flow model
+
+# %%
+# Import required libraries
+
+import os
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import ogstools as ot
+import pandas as pd
+import vtuIO
+
+# %% [markdown]
+# ## Ramey's analytical solution
+#
+# The equations are followed by their implementation in python.
+
+# %%
+# Define constants
+
+T_d = 55 # Temperature of undistributed formation
+T_i = 20 # Temperature of inlet fluid
+Z = np.insert(
+ np.linspace(0.3, 30, 100), 0, 1e-10
+) # Depth below surface, start with very small value to prevent division by zero
+q = 0.0002 # Flow rate of fluid in wellbore, m3/s
+rho_f = 1000 # Density of fluid, unit kg/m3
+c_p_f = 4190 # Specific heat capacity of fluid in wellbore
+mu_f = 1.14e-3 # Unit kg/m/s;
+lambda_f = 0.59 # Unit W/m/K
+rho_re = 1800 # Density of geothermal reservoir
+c_p_re = 1778 # Specific heat capacity of geothermal reservoir
+lambda_re = 2.78018 # Heat conductivity coefficient of geothermal reservoir
+lambda_g = 0.73 # Heat conductivity coefficient of grout and pipe
+lambda_pi = 1.3
+r_pi = 0.12913 # Inner radius of pipe and wellbore
+r_b = 0.14
+t_pi = 0.00587 # Thickness of the pipe
+t = 86400 * 5 # Operation time
+
+
+# %% [markdown]
+# In Ramey's analytical solution Ramey et al. (1962), the outlet temperature of the pipe inside the wellbore can be calculated by
+#
+# $$
+# T_o(t) = T_{s} + (T_i(t) - T_{s})\exp(-\Delta z/X)
+# $$
+
+
+# %%
+def temp(T_d, T_i, delta_z, X):
+ return T_d + (T_i - T_d) * np.exp(-delta_z / X)
+
+
+# %% [markdown]
+# The coefficient $X$ is determined by
+#
+# $$
+# X = \frac{q\rho_fc_{p,f}(\lambda_{s}+r_pUf(t))}{2\pi r_pU \lambda_{s}}
+# $$
+#
+# where $q$ is the flow rate of the fluid in the wellbore.
+
+
+# %%
+def coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t):
+ return (
+ (q * rho_f * c_p_f)
+ * (lambda_re + r_pi * U * f_t)
+ / (2 * np.pi * r_pi * U * lambda_re)
+ )
+
+
+# %% [markdown]
+# With dimensionless time
+#
+# $$
+# t_D = \frac{\lambda_{s}t}{(\rho_{s}c_{p,s}r_b^{2})}
+# $$
+
+
+# %%
+def dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b):
+ return lambda_re * delta_t / (rho_re * c_p_re * r_b * r_b)
+
+
+# %% [markdown]
+# the time function $f(t)$ can be calculated by
+#
+# <!-- markdownlint-disable reference-links-images -->
+#
+# $$
+# f(t) = [0.4063+0.5\ln(t_D)][1+\frac{0.6}{t_D}], t_D > 1.5
+# $$
+#
+# <!-- markdownlint-enable reference-links-images -->
+#
+# $$
+# f(t) = 1.1281\sqrt{t_D}(1-0.3\sqrt{t_D}), t_D \leqslant 1.5
+# $$
+
+
+# %%
+def time_function(t_D):
+ if t_D > 1.5:
+ f_t = (0.4063 + 0.5 * np.log(t_D)) * (1 + 0.6 / t_D)
+ else:
+ f_t = 1.1281 * np.sqrt(t_D) * (1 - 0.3 * np.sqrt(t_D))
+
+ return f_t
+
+
+# %% [markdown]
+# The Prandtl and Reynolds number can be calculated as follows
+#
+# $$
+# \mathrm{Pr} = \frac{\mu_f c_{p,f}}{\lambda_f}
+# $$
+#
+# $$
+# \mathrm{Re} = \frac{\rho_f v d_{pi}}{\mu_f}
+# $$
+#
+# where, $\mu_f, \rho_f$ and $\lambda_f$ is the fluid viscosity, density and thermal conductivity.
+
+# %%
+# Calculate the convective film resistance inside pipe
+v = q / (np.pi * r_pi * r_pi)
+
+# Reynolds number
+Pr = mu_f * c_p_f / lambda_f
+Re = rho_f * v * (2 * r_pi) / mu_f
+
+
+# %% [markdown]
+# The Nusselt number can be determined by the following equation Diersch, (2011):
+#
+# $$
+# \mathrm{Nu} = 4.364,\ \mathrm{Re} < 2300
+# $$
+#
+# $$
+# \mathrm{Nu} = \frac{(\xi_{k}/8)\ \mathrm{Re}_{k}\ \mathrm{Pr}}{1+12.7\sqrt{{\xi_k}/8}(\mathrm{Pr}^{2/3}-1)} [ 1+(\frac{{d_k}^{i}}{L})^{2/3}], Re \geq 10^4
+# $$
+#
+# $$
+# \mathrm{Nu} = (1-\gamma_{k})\ 4.364 + \gamma_{k} ( \frac{(0.0308/8)10^{4}\mathrm{Pr}}{1+12.7\ \sqrt{0.0308/8}(\mathrm{Pr}^{2/3}-1)} \left[ 1+\left(\frac{d_{k}^{i}}{L}\right)^{2/3} \right] ), 2 300 \leq \mathrm{Re} < 10^{4}
+# $$
+#
+# with
+#
+# $$
+# \xi_{k}=(1.8\ \log_{10}\mathrm{Re}_{k}-1.5)^{-2}
+# $$
+#
+# and
+#
+# $$
+# \gamma_{k} = \frac{\mathrm{Re}_{k}-2300}{10^{4}-2300}
+# $$
+
+# %%
+# Estimate convective film coefficient
+xi = (1.8 * np.log(Re) - 1.5) ** (-2)
+Nu_l = (
+ (xi / 8.0 * Re * Pr)
+ / (1.0 + 12.7 * np.sqrt(xi / 8.0) * ((Pr) ** (2.0 / 3.0) - 1.0))
+ * (1.0 + (r_pi / Z) ** (2.0 / 3.0))
+)
+Nu_lam = 4.364
+gamma = (Re - 2300) / (10000 - 2300)
+
+Nu_m = (1.0 - gamma) * 4.364 + gamma * (
+ (0.0308 / 8.0 * 1.0e4 * Pr)
+ / (1.0 + 12.7 * np.sqrt(0.0308 / 8.0) * ((Pr) ** (2.0 / 3.0) - 1.0))
+ * (1.0 + (r_pi / Z) ** (2.0 / 3.0))
+)
+
+# Evaluate Nusselt number; based on the value of Re choose appropriate correlation
+if Re < 2300:
+ Nu_p = Nu_lam
+elif Re > 10000:
+ Nu_p = Nu_l
+else:
+ Nu_p = Nu_m
+
+# %% [markdown]
+# and the overall heat transfer coefficient $U$ is written as follows,
+#
+# $$
+# U = [\frac{r_{pi}+t_{pi}}{r_{pi}h}+(r_{pi}+t_{pi})(\frac{\ln{\frac{r_{pi}+t_{pi}}{r_{pi}}}}{\lambda_{pi}}+\frac{\ln{\frac{r_b}{r_{pi}+t_{pi}}}}{\lambda_{grout}})]^{-1}
+# $$
+#
+# with
+#
+# $$
+# h = \frac{\lambda_f Nu}{2r_{pi}}
+# $$
+
+# %%
+h = lambda_f * Nu_p / (2 * r_pi) # Unit: W/m2/K
+
+U = 1 / (
+ ((r_pi + t_pi) / (r_pi * h))
+ + (r_pi + t_pi)
+ * (
+ np.log((r_pi + t_pi) / r_pi) / lambda_pi
+ + np.log(r_b / (r_pi + t_pi)) / lambda_g
+ )
+)
+
+# %% [markdown]
+# The friction factor $f$, is evaluated by Churchill correlation Churchill et al. (1977),
+#
+# $$
+# f = \frac{1}{(\frac{1}{[((\frac{8}{Re})^{10}+(\frac{Re}{36500})^{20})]^{1/2}}+[2.21(\ln{\frac{Re}{7}})]^{10})^{1/5}}
+# $$
+
+# %%
+# Churchill correlation for friction factor
+f = 1 / (
+ (1 / (np.sqrt((8 / Re) ** 10 + (Re / 36500) ** 20))) + (2.21 * np.log(Re / 7)) ** 10
+) ** (1 / 5)
+
+# %% [markdown]
+# With these equations the analytical solution can be computed:
+
+# %%
+# Compute solution
+
+time_range = range(0, t + 1, 43200)
+
+# Create a dataframe to store results for different z- and t-steps
+row_names = Z
+column_names = time_range
+
+
+data = pd.DataFrame(index=row_names, columns=column_names)
+
+
+for delta_z in Z:
+ T_0 = []
+ for delta_t in time_range:
+ # Compute dimensionless time
+ t_D = dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b)
+
+ # Compute time function
+ f_t = time_function(t_D)
+
+ # Compute X
+ X = coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t)
+
+ # Compute outlet temperature
+ T_0.append(temp(T_d, T_i, delta_z, X))
+
+ data.loc[delta_z] = T_0
+
+# %% [markdown]
+# ## Numerical solution
+
+# %%
+# Create output path
+out_dir = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out"))
+if not out_dir.exists():
+ out_dir.mkdir(parents=True)
+
+# %%
+# Create an OGS-object
+# Pass the project file and set an output file path
+model = ot.Project(input_file="BHE_1P.prj", output_file=f"{out_dir}/BHE_1P.prj")
+
+# set end time (in seconds) and the time increments
+model.replace_text(
+ 5 * 24 * 60 * 60, xpath="./time_loop/processes/process/time_stepping/t_end"
+)
+model.replace_text(
+ 30, xpath="./time_loop/processes/process/time_stepping/timesteps/pair/repeat"
+)
+model.replace_text(
+ 4 * 60 * 60,
+ xpath="./time_loop/processes/process/time_stepping/timesteps/pair/delta_t",
+)
+# Write to the output file
+model.write_input()
+
+# Run OGS
+model.run_model(logfile=Path(out_dir) / "log.txt", args=f"-o {out_dir} -m .")
+
+# %% [markdown]
+# ## Results and discussion
+#
+# The outlet temperature change over time was compared against the analytical solution and presented in Figure 2.
+# After 30 days, the fluid temperature distribution in the wellbore is shown in Figure 3.
+# The maximum relative error between the numerical model and Ramey's analytical solution is less than 0.15%.
+#
+# In numerical model, the outlet temperature at beginning stage is affected by the initial temperature in the pipe inside the wellbore.
+# The initial fluid temperature set in the benchmark means there is water with 20°C filled in the wellbore already before injecting water into the wellbore.
+# But in the analytical solution, no initial temperature is set and the temperature keeps equilibrium state at every moment.
+# The impact of initial temperature condition in numerical model is decreasing with the increasing of the operational time as shown in Figure 2.
+
+# %%
+out_dir = os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out")
+
+# Load the result
+pvdfile = vtuIO.PVDIO(f"{out_dir}/BHE_1P.pvd", dim=3)
+
+# Get point field names
+fields = pvdfile.get_point_field_names()
+print(fields)
+
+# Get all written timesteps
+time = pvdfile.timesteps
+
+# Define observation points
+observation_points = {"outlet": (0.0, 10.0, -30.0)}
+
+# Read the time series of point field 'temperature_BHE1'
+temperature_BHE1 = pvdfile.read_time_series("temperature_BHE1", observation_points)
+
+# %%
+# Compute analytical solution again but at timesteps of numerical solution for error calculation
+for delta_z in Z:
+ T_0 = []
+ for delta_t in time:
+ # Compute dimensionless time
+ t_D = dimless_t(lambda_re, delta_t, rho_re, c_p_re, r_b)
+
+ # Compute time function
+ f_t = time_function(t_D)
+
+ # Compute X
+ X = coefficient_x(q, rho_f, c_p_f, lambda_re, r_pi, U, f_t)
+
+ # Compute outlet temperature
+ T_0.append(temp(T_d, T_i, delta_z, X))
+
+# %%
+# Plot numerical vs. analytical result
+fig, ax1 = plt.subplots(figsize=(10, 8))
+
+x_pos = np.arange(0, 432000, 60 * 60 * 24)
+x_ticks = np.arange(0, 5, 1)
+
+
+ax1.plot(
+ time_range,
+ data.iloc[-1, :],
+ "k.",
+ markersize=10,
+ markerfacecolor="none",
+ label="Ramey's analytical solution",
+)
+ax1.plot(time, temperature_BHE1["outlet"][:, 0], "r", label="OGS")
+ax1.set_ylabel("Outlet temperature (°C)", fontsize=20)
+ax1.set_xlabel("Time (d)", fontsize=20)
+ax1.set_xticks(x_pos, x_ticks)
+ax1.spines["left"].set_linewidth(2)
+ax1.spines["bottom"].set_linewidth(2)
+ax1.tick_params(axis="both", labelsize=16)
+
+color = "blue"
+ax2 = ax1.twinx()
+ax2.plot(
+ time, temperature_BHE1["outlet"][:, 0] - T_0, color=color, label="Absolute error"
+)
+ax2.set_ylabel("Absolute error (°C)", color=color, fontsize=20)
+ax2.tick_params(axis="y", labelcolor=color)
+ax2.spines["right"].set_color("blue")
+ax2.spines["right"].set_linewidth(2)
+ax2.tick_params(axis="y", labelsize=16)
+
+plt.figlegend(fontsize=18, loc=(0.4, 0.1), frameon=False)
+fig.tight_layout()
+
+plt.show()
+
+# %% [markdown]
+# Figure 2: Comparison with analytical solution results
+
+# %%
+# Extract values along a line in the domain
+
+# Space axis for plotting
+length = np.linspace(0, -30, 101)
+
+# Draws a line through the domain for sampling results
+x_axis = [(0, 10, i) for i in length]
+
+# At timestep i
+i = -1
+temp = pvdfile.read_set_data(time[i], "temperature_BHE1", pointsetarray=x_axis)
+
+# Plot
+fig, ax1 = plt.subplots(figsize=(10, 8))
+ax1.set_xlabel("Distance (m)", fontsize=20)
+ax1.set_ylabel("Temperature (°C)", fontsize=20)
+ax1.plot(np.linspace(0, 30, 101), temp[:, 0], "r", linewidth=2, label="OGS")
+ax1.plot(
+ Z,
+ data.iloc[:, -1],
+ "k.",
+ markersize=10,
+ markerfacecolor="none",
+ label="Ramey's analytical solution",
+)
+ax1.set_ylim(19, 25)
+ax1.set_xlim(0, 30)
+ax1.spines["left"].set_linewidth(2)
+ax1.spines["bottom"].set_linewidth(2)
+ax1.tick_params(axis="both", labelsize=16)
+
+color = "blue"
+ax2 = ax1.twinx()
+ax2.set_ylabel("Absolute error (°C)", color=color, fontsize=20)
+ax2.plot(
+ np.linspace(0, 30, 101),
+ temp[:, 0] - data.iloc[:, -1],
+ linewidth=2,
+ color=color,
+ label="Absolute error",
+)
+ax2.tick_params(axis="y", labelcolor=color)
+ax2.spines["right"].set_color("blue")
+ax2.spines["right"].set_linewidth(2)
+ax2.tick_params(axis="y", labelsize=16)
+
+fig.tight_layout()
+plt.figlegend(fontsize=18, loc=(0.4, 0.1), frameon=False)
+plt.show()
+
+max_error = 0.08
+if np.max(np.abs(temp[:, 0] - data.iloc[:, -1])) > max_error:
+ msg = f"The maximum temperature error is over {max_error} °C, which is higher than expected"
+ raise RuntimeError(msg)
+
+# %% [markdown]
+# Figure 3: Distributed temperature of fluid and absolute error.
+#
+# ## References
+#
+# <!-- vale off -->
+#
+# [1] Ramey Jr, H. J. (1962). Wellbore heat transmission. Journal of petroleum Technology, 14(04), 427-435.
+#
+# [2] Diersch, H-JG, et al. (2011). "Finite element modeling of borehole heat exchanger systems: Part 1. Fundamentals." Computers & Geosciences 37.8, 1122-1135.
+#
+# [3] Chaofan Chen, Wanlong Cai, Dmitri Naumov, Kun Tu, Hongwei Zhou, Yuping Zhang, Olaf Kolditz, Haibing Shao (2021). Numerical investigation on the capacity and efficiency of a deep enhanced U-tube borehole heat exchanger system for building heating. Renewable Energy, 169, 557-572.
+#
+# [4] Churchill, S. W. (1977). Comprehensive correlating equations for heat, mass and momentum transfer in fully developed flow in smooth tubes. Industrial & Engineering Chemistry Fundamentals, 16(1), 109-116.
+#
+# <!-- vale on -->