ogs-jupyter-lab-h-2d-theis.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "id": "ac287b2f",
   "metadata": {},
   "source": [
    "title = \"H process: Theis solution (Pumping well)\"\n",
    "date = \"2022-08-12\"\n",
    "author = \"Wenqing Wang, Olaf Kolditz\"\n",
    "notebook = \"Parabolic/LiquidFlow/AxiSymTheis/axisym_theis.ipynb\"\n",
    "[menu]\n",
    "  [menu.benchmarks]\n",
    "    parent = \"Liquid Flow\"\n",
    "<!--eofm-->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3848074",
   "metadata": {},
   "outputs": [],
   "source": [
    "#modules\n",
    "from IPython.display import Image\n",
    "import os\n",
    "import pyvista as pv\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import exp1\n",
    "import vtk\n",
    "from vtk.util.numpy_support import vtk_to_numpy\n",
    "import matplotlib.tri as tri\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffd2a84d",
   "metadata": {},
   "outputs": [],
   "source": [
    "#settings\n",
    "path='./'\n",
    "fig_dir = \"./figures/\"\n",
    "data_dir = os.environ.get('OGS_DATA_DIR', '../../../Data')\n",
    "out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')\n",
    "exe_dir = \"/home/ok/ogs/release/bin/\" #relative path needed\n",
    "prj_name = \"axisym_theis\"\n",
    "prj_file = f\"{prj_name}.prj\"\n",
    "pvd_name = \"liquid_pcs\"\n",
    "vtu_name = \"axisym_theis.vtu\"\n",
    "title = \"H process: Theis solution (Pumping well)\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "688250b6",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Image(filename = fig_dir + \"ogs-jupyter-lab.png\", width=150, height=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2ff19ce",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "Image(filename = fig_dir + \"h-tet-1.png\", width=150, height=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fae4b91",
   "metadata": {
    "tags": []
   },
   "source": [
    "## H process: Theis solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d6228a95",
   "metadata": {},
   "source": [
    "**Problem description**\n",
    "\n",
    "Theis’ problem examines the transient lowering of the water table induced by a pumping well. \n",
    "The assumptions required by the Theis solution are:\n",
    "\n",
    "The aquifer \n",
    "- is homogeneous, isotropic, confined, infinite in radial extent,\n",
    "- has uniform thickness, horizontal piezometric surface.\n",
    "\n",
    "The well\n",
    "- is fully penetrating the entire aquifer thickness,\n",
    "- has a constant pumping rate,\n",
    "- well storage effects can be neglected,\n",
    "- no other wells or long term changes in regional water levels."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3cb1306f",
   "metadata": {},
   "source": [
    "**Analytical solution**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "296db500",
   "metadata": {},
   "source": [
    "The analytical solution of the drawdown as a function of time and distance is expressed by\n",
    "$$\n",
    "s(r,t) = h_0 - h(r,t) = \\frac{Q}{4\\pi T}W(u), \\quad \\mathrm{where}\\quad u = \\frac{r^2S}{4Tt}.\n",
    "$$\n",
    "\n",
    "where\n",
    "- $s$ [$L$] is the _drawdown_ or change in hydraulic head,\n",
    "- $h_0$ is the constant initial hydraulic head,\n",
    "- $h$ is the hydrauic head at distance $r$ at time $t$\n",
    "- $Q$ [$L^3T^{-1}$] is the constant pumping (discharge) rate\n",
    "- $S$ [$-$] is the aquifer storage coefficient (volume of water released per unit decrease in $h$ per unit area)\n",
    "- $T$ [$L^2T^{-1}$] is the transmissivity (a measure of how much water is transported horizontally per unit time).\n",
    "\n",
    "The _Well Function_, $W(u)$ is the exponential integral, $E_1(u).$\n",
    "$W(u)$ is defined by an infinite series:\n",
    "$$\n",
    "W(u) = - \\gamma - \\ln u + \\sum_{k=1}^\\infty \\frac{(-1)^{k+1} u^k}{k \\cdot k!}\n",
    "$$\n",
    "where\n",
    "- $\\gamma=0.577215664$ is the Euler-Mascheroni constant\n",
    "\n",
    "For practical applications an approximation to the exponential integral is used often:\n",
    "$$W(u) \\approx -\\gamma - \\ln u$$\n",
    "\n",
    "This results in an expression for $s(r,t)$ known as the Jacob equation:\n",
    "$$\n",
    "s(r,t) = -\\frac{Q}{4\\pi T}\\left(\\gamma + \\ln u \\right).\n",
    "$$\n",
    "For more details we refer to Srivastava and Guzman (19??)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f6f31590",
   "metadata": {},
   "outputs": [],
   "source": [
    "#source: https://scipython.com/blog/linear-and-non-linear-fitting-of-the-theis-equation/\n",
    "\n",
    "def calc_u(r, S, T, t):\n",
    "    \"\"\"Calculate and return the dimensionless time parameter, u.\"\"\"\n",
    "\n",
    "    return r**2 * S / 4 / T / t\n",
    "\n",
    "def theis_drawdown(t, S, T, Q, r):\n",
    "    \"\"\"Calculate and return the drawdown s(r,t) for parameters S, T.\n",
    "\n",
    "    This version uses the Theis equation, s(r,t) = Q * W(u) / (4.pi.T),\n",
    "    where W(u) is the Well function for u = Sr^2 / (4Tt).\n",
    "    S is the aquifer storage coefficient,\n",
    "    T is the transmissivity (m2/day),\n",
    "    r is the distance from the well (m), and\n",
    "    Q is the pumping rate (m3/day).\n",
    "\n",
    "    \"\"\"\n",
    "\n",
    "    u = calc_u(r, S, T, t)\n",
    "    s_theis = Q/4/np.pi/T * exp1(u)\n",
    "    return s_theis\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a24764b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q = 2000        # Pumping rate from well (m3/day)\n",
    "r = 10          # Distance from well (m)\n",
    "\n",
    "# Time grid, days.\n",
    "t = np.array([1, 2, 4, 8, 12, 16, 20, 30, 40, 50, 60, 70, 80, 90, 100])\n",
    "\n",
    "# Calculate some synthetic data to fit.\n",
    "S, T = 0.0003, 1000\n",
    "s = theis_drawdown(t, S, T, Q, r)\n",
    "\n",
    "# Plot the data\n",
    "titlestring = \"Theis: Analytical solution\"\n",
    "plt.title(titlestring)\n",
    "plt.plot(t, s, label='r = '+str(r)+' m')\n",
    "plt.xlabel(r'$t\\;/\\mathrm{days}$')\n",
    "plt.ylabel(r'$s\\;/\\mathrm{m}$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f748dd1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recalculation from days in sec\n",
    "Q = 0.016      # Pumping rate from well (m3/s)\n",
    "t = 864000     # Time in s.\n",
    "\n",
    "# Distance from well (m)\n",
    "r = np.array([0.5, 1, 2, 4, 8, 12, 16, 20, 25, 30, 35, 40])\n",
    "\n",
    "# Calculate some synthetic data to fit.\n",
    "S = 0.001\n",
    "T = 9.2903e-4\n",
    "u = calc_u(r, S, T, t)\n",
    "s = theis_drawdown(t, S, T, Q, r)\n",
    "s = s-5\n",
    "# Add some noise.\n",
    "#noise_sd = 4.e-3\n",
    "#s += s * noise_sd * np.random.randn(s.shape[0])\n",
    "\n",
    "# Plot the data\n",
    "titlestring = \"Theis: Analytical solution\"\n",
    "plt.title(titlestring)\n",
    "plt.plot(r, s, label='t = '+str(t)+' days')\n",
    "plt.xlabel(r'$r\\;/\\mathrm{m}$')\n",
    "plt.ylabel(r'$s\\;/\\mathrm{m}$')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "271ac1cc",
   "metadata": {},
   "source": [
    "**Numerical solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f66a6aec",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "mesh = pv.read(vtu_name)\n",
    "print(\"inspecting vtu-file\")\n",
    "mesh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e5ae5ee",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"inspecting mesh and initial conditions\")\n",
    "#file\n",
    "reader = vtk.vtkXMLUnstructuredGridReader()\n",
    "reader.SetFileName(vtu_name)\n",
    "reader.Update()  # Needed because of GetScalarRange\n",
    "data = reader.GetOutput()\n",
    "pressure = data.GetPointData().GetArray(\"OGS5_pressure\")\n",
    "#points\n",
    "points = data.GetPoints()\n",
    "npts = points.GetNumberOfPoints()\n",
    "x = vtk_to_numpy(points.GetData())\n",
    "triang = tri.Triangulation(x[:,0], x[:,1])\n",
    "#plt.triplot(triang, 'go-', lw=1.0)\n",
    "plt.triplot(triang,lw=0.2)\n",
    "plt.tricontour(triang, pressure, 16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b327862b",
   "metadata": {},
   "source": [
    "**Running OGS**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66ef056a",
   "metadata": {},
   "outputs": [],
   "source": [
    "#run ogs\n",
    "t0 = time.time()\n",
    "print(\"run ogs\")\n",
    "print(f\"{exe_dir}ogs {prj_file} > log.txt\")\n",
    "#! {exe_dir}ogs {prj_file} > log.txt\n",
    "! ogs {prj_file} > log.txt\n",
    "tf = time.time()\n",
    "print(\"computation time: \", round(tf - t0, 2), \" s.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50bebe96",
   "metadata": {},
   "source": [
    "**Spatial Profiles**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2f0bf2c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import vtuIO\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "#Read simulation results\n",
    "pvdfile = vtuIO.PVDIO(f\"{pvd_name}.pvd\", dim=2)\n",
    "xaxis =  [(i,0,0) for i in np.linspace(start=0.01, stop=40, num=100)]\n",
    "r_x = np.array(xaxis)[:,0]\n",
    "time = [8.64,86.4,1728,24192,172800,604800,864000]\n",
    "\n",
    "pressure_xaxis_t = pvdfile.read_set_data(t, 'OGS5_pressure', data_type=\"point\", pointsetarray=xaxis)\n",
    "plt.plot(r_x, pressure_xaxis_t, 'x')\n",
    "\n",
    "for t in time:\n",
    "    pressure_xaxis_t = pvdfile.read_set_data(t, 'pressure', data_type=\"point\", pointsetarray=xaxis)\n",
    "    plt.plot(r_x, pressure_xaxis_t, label='t = '+str(t)+' sec')\n",
    "titlestring = \"Theis: Numerical solution\"\n",
    "plt.title(titlestring)\n",
    "plt.xlabel(r'$r\\;/\\mathrm{m}$')\n",
    "plt.ylabel('Hydraulic head')\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "#plt.savefig(\"theis.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5366c257",
   "metadata": {},
   "outputs": [],
   "source": [
    "time = [864000]\n",
    "pressure_xaxis_t = pvdfile.read_set_data(t, 'pressure', data_type=\"point\", pointsetarray=xaxis)\n",
    "plt.plot(r_x, pressure_xaxis_t, 'x', label='numerical solution (ogs6)')\n",
    "plt.plot(r, s, label='analytical solution')\n",
    "titlestring = \"Theis: Comparison analytical and numerical solution\"\n",
    "plt.title(titlestring)\n",
    "plt.legend()\n",
    "plt.xlabel(r'$radius\\;/\\;\\mathrm{m}$')\n",
    "plt.ylabel(r'$hydraulic head\\;/\\;\\mathrm{m}$')\n",
    "plt.grid()\n",
    "#plt.savefig(\"theis-ana+num.png\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb6272f0",
   "metadata": {},
   "source": [
    "**Temporal Profiles**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78afcf25",
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "print(time.ctime())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bef8bd50",
   "metadata": {
    "tags": []
   },
   "source": [
    "**OGS links**\n",
    "- description: https://www.opengeosys.org/docs/benchmarks/\n",
    "- project file: https://gitlab.opengeosys.org/ogs/ogs/-/blob/master/Tests/Data/Parabolic/LiquidFlow/AxiSymTheis/axisym_theis.prj"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ac32c64",
   "metadata": {},
   "source": [
    "**Credits**\n",
    "- Christian for the analytical solution in Python, https://scipython.com/blog/linear-and-non-linear-fitting-of-the-theis-equation/\n",
    "- Wenqing Wang for set-up the OGS benchmark, https://www.opengeosys.org/docs/benchmarks/liquid-flow/liquid-flow-theis-problem/\n",
    "- Jörg Buchwald for ogs6py and VTUInterface, https://joss.theoj.org/papers/6ddcac13634c6373e6e808952fdf4abc"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}