{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test wells with analytical solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.integrate import quad\n",
    "from scipy.special import exp1\n",
    "\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"font.size\"] = 8.0\n",
    "plt.rcParams[\"figure.figsize\"] = (8, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Theis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def theis(r, t, T, S, Q):\n",
    "    u = r**2 * S / (4 * T * t)\n",
    "    h = -Q / (4 * np.pi * T) * exp1(u)\n",
    "    return h\n",
    "\n",
    "\n",
    "def theisQr(r, t, T, S, Q):\n",
    "    u = r**2 * S / (4 * T * t)\n",
    "    return -Q / (2 * np.pi) * np.exp(-u) / r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 500\n",
    "S = 1e-4\n",
    "t = np.logspace(-5, 0, 100)\n",
    "r = 30\n",
    "Q = 788"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "htheis = theis(r, t, T, S, Q)\n",
    "Qrtheis = theisQr(r, t, T, S, Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=S / 20, tmin=1e-5, tmax=1)\n",
    "w = tft.Well(ml, tsandQ=[(0, Q)], rw=1e-5)\n",
    "ml.solve()\n",
    "h = ml.head(r, 0, t)\n",
    "Qx, Qy = ml.disvec(r, 0, t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(121)\n",
    "plt.semilogx(t, htheis, \"b\", label=\"theis\")\n",
    "plt.semilogx(t, h[0], \"r--\", label=\"timflow\")\n",
    "plt.xlabel(\"time (day)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.subplot(122)\n",
    "plt.semilogx(t, Qrtheis, \"b\", label=\"theis\")\n",
    "plt.semilogx(t, Qx[0], \"r--\", label=\"timflow\")\n",
    "plt.xlabel(\"time (day)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test(M):\n",
    "    ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=S / 20, tmin=1e-5, tmax=1, M=M)\n",
    "    tft.Well(ml, tsandQ=[(0, Q)], rw=1e-5)\n",
    "    ml.solve(silent=True)\n",
    "    h = ml.head(r, 0, t)\n",
    "    return htheis - h[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enumba = test(M=10)\n",
    "plt.plot(t, enumba, \"C1\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head difference Theis - timflow\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, Qrtheis - Qx[0])\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"Qx difference Theis - Ttim\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compare(M=10):\n",
    "    ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=S / 20, tmin=1e-5, tmax=1, M=M)\n",
    "    tft.Well(ml, tsandQ=[(0, Q)], rw=1e-5)\n",
    "    ml.solve(silent=True)\n",
    "    h = ml.head(r, 0, t)\n",
    "    rmse = np.sqrt(np.mean((h[0] - htheis) ** 2))\n",
    "    return rmse\n",
    "\n",
    "\n",
    "Mlist = np.arange(1, 21)\n",
    "rmse = np.zeros(len(Mlist))\n",
    "for i, M in enumerate(Mlist):\n",
    "    rmse[i] = compare(M)\n",
    "plt.semilogy(Mlist, rmse)\n",
    "plt.xlabel(\"Number of terms M\")\n",
    "plt.xticks(np.arange(1, 21))\n",
    "plt.ylabel(\"relative error\")\n",
    "plt.title(\n",
    "    \"comparison between timflow solution and Theis \\n solution using numba and M terms\"\n",
    ")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def volume(r, t=1):\n",
    "    return -2 * np.pi * r * ml.head(r, 0, t) * ml.aq.Scoefaq[0]\n",
    "\n",
    "\n",
    "quad(volume, 1e-5, np.inf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def theis2(r, t, T, S, Q, tend):\n",
    "    u1 = r**2 * S / (4 * T * t)\n",
    "    u2 = r**2 * S / (4 * T * (t[t > tend] - tend))\n",
    "    h = -Q / (4 * np.pi * T) * exp1(u1)\n",
    "    h[t > tend] -= -Q / (4 * np.pi * T) * exp1(u2)\n",
    "    return h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml2 = tft.ModelMaq(kaq=25, z=[20, 0], Saq=S / 20, tmin=1e-5, tmax=10)\n",
    "w2 = tft.Well(ml2, tsandQ=[(0, Q), (1, 0)])\n",
    "ml2.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t2 = np.linspace(0.01, 2, 100)\n",
    "htheis2 = theis2(r, t2, T, S, Q, tend=1)\n",
    "h2 = ml2.head(r, 0, t2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t2, htheis2, \"b\", label=\"theis\")\n",
    "plt.plot(t2, h2[0], \"r--\", label=\"timflow\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hantush"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 500\n",
    "S = 1e-4\n",
    "c = 1000\n",
    "t = np.logspace(-5, 0, 100)\n",
    "r = 30\n",
    "Q = 788"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import quad\n",
    "\n",
    "\n",
    "def integrand_hantush(y, r, lab):\n",
    "    return np.exp(-y - r**2 / (4 * lab**2 * y)) / y\n",
    "\n",
    "\n",
    "def hantush(r, t, T, S, c, Q, tstart=0):\n",
    "    lab = np.sqrt(T * c)\n",
    "    u = r**2 * S / (4 * T * (t - tstart))\n",
    "    F = quad(integrand_hantush, u, np.inf, args=(r, lab))[0]\n",
    "    return -Q / (4 * np.pi * T) * F\n",
    "\n",
    "\n",
    "hantushvec = np.vectorize(hantush)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=25, z=[21, 20, 0], c=[1000], Saq=S / 20, topboundary=\"semi\", tmin=1e-5, tmax=1\n",
    ")\n",
    "w = tft.Well(ml, tsandQ=[(0, Q)])\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hhantush = hantushvec(30, t, T, S, c, Q)\n",
    "h = ml.head(r, 0, t)\n",
    "plt.semilogx(t, hhantush, \"b\", label=\"hantush\")\n",
    "plt.semilogx(t, h[0], \"r--\", label=\"timflow\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Well with welbore storage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 500\n",
    "S = 1e-4\n",
    "t = np.logspace(-5, 0, 100)\n",
    "rw = 0.3\n",
    "Q = 788\n",
    "\n",
    "ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=S / 20, tmin=1e-5, tmax=1)\n",
    "w = tft.Well(ml, rw=rw, tsandQ=[(0, Q)])\n",
    "ml.solve()\n",
    "hnostorage = ml.head(rw, 0, t)\n",
    "\n",
    "ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=S / 20, tmin=1e-5, tmax=1)\n",
    "w = tft.Well(ml, rw=rw, tsandQ=[(0, Q)], rc=rw)\n",
    "ml.solve()\n",
    "hstorage = ml.head(rw, 0, t)\n",
    "\n",
    "plt.semilogx(t, hnostorage[0], label=\"no storage\")\n",
    "plt.semilogx(t, hstorage[0], label=\"with storage\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.xticks(\n",
    "    [1 / (24 * 60 * 60), 1 / (24 * 60), 1 / 24, 1], [\"1 sec\", \"1 min\", \"1 hr\", \"1 d\"]\n",
    ")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slug test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 25\n",
    "H = 20\n",
    "S = 1e-4 / H\n",
    "t = np.logspace(-7, -1, 100)\n",
    "rw = 0.2\n",
    "rc = 0.2\n",
    "delh = 1\n",
    "ml = tft.ModelMaq(kaq=k, z=[H, 0], Saq=S, tmin=1e-7, tmax=1)\n",
    "Qslug = np.pi * rc**2 * delh\n",
    "w = tft.Well(ml, tsandQ=[(0, -Qslug)], rw=rw, rc=rc, wbstype=\"slug\")\n",
    "ml.solve()\n",
    "h = w.headinside(t)\n",
    "plt.semilogx(t, h[0])\n",
    "plt.xticks(\n",
    "    [1 / (24 * 60 * 60) / 10, 1 / (24 * 60 * 60), 1 / (24 * 60), 1 / 24],\n",
    "    [\"0.1 sec\", \"1 sec\", \"1 min\", \"1 hr\"],\n",
    ")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slug test in 5-layer aquifer\n",
    "Well in top 2 layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 25\n",
    "H = 20\n",
    "Ss = 1e-4 / H\n",
    "t = np.logspace(-7, -1, 100)\n",
    "rw = 0.2\n",
    "rc = 0.2\n",
    "delh = 1\n",
    "ml = tft.Model3D(kaq=k, z=np.linspace(H, 0, 6), Saq=Ss, tmin=1e-7, tmax=1)\n",
    "Qslug = np.pi * rc**2 * delh\n",
    "w = tft.Well(ml, tsandQ=[(0, -Qslug)], rw=rw, rc=rc, layers=[0, 1], wbstype=\"slug\")\n",
    "ml.solve()\n",
    "hw = w.headinside(t)\n",
    "plt.semilogx(t, hw[0], label=\"inside well\")\n",
    "h = ml.head(0.2 + 1e-8, 0, t)\n",
    "for i in range(2, 5):\n",
    "    plt.semilogx(t, h[i], label=\"layer\" + str(i))\n",
    "plt.legend()\n",
    "plt.xticks(\n",
    "    [1 / (24 * 60 * 60) / 10, 1 / (24 * 60 * 60), 1 / (24 * 60), 1 / 24],\n",
    "    [\"0.1 sec\", \"1 sec\", \"1 min\", \"1 hr\"],\n",
    ")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "20 layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 25\n",
    "H = 20\n",
    "S = 1e-4 / H\n",
    "t = np.logspace(-7, -1, 100)\n",
    "rw = 0.2\n",
    "rc = 0.2\n",
    "delh = 1\n",
    "ml = tft.Model3D(kaq=k, z=np.linspace(H, 0, 21), Saq=S, tmin=1e-7, tmax=1)\n",
    "Qslug = np.pi * rc**2 * delh\n",
    "w = tft.Well(ml, tsandQ=[(0, -Qslug)], rw=rw, rc=rc, layers=np.arange(8), wbstype=\"slug\")\n",
    "ml.solve()\n",
    "hw = w.headinside(t)\n",
    "plt.semilogx(t, hw[0], label=\"inside well\")\n",
    "h = ml.head(0.2 + 1e-8, 0, t)\n",
    "for i in range(8, 20):\n",
    "    plt.semilogx(t, h[i], label=\"layer\" + str(i))\n",
    "plt.legend()\n",
    "plt.xticks(\n",
    "    [1 / (24 * 60 * 60) / 10, 1 / (24 * 60 * 60), 1 / (24 * 60), 1 / 24],\n",
    "    [\"0.1 sec\", \"1 sec\", \"1 min\", \"1 hr\"],\n",
    ")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Head Well"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=1e-5, tmin=1e-3, tmax=1000)\n",
    "w = tft.HeadWell(ml, tsandh=[(0, -1)], rw=0.2)\n",
    "ml.solve()\n",
    "ax0 = plt.subplot(1, 2, 1)\n",
    "ml.plots.head_along_line(0.2, 100, 0, 0, 100, t=[0.1, 1, 10], sstart=0.2, ax=ax0)\n",
    "t = np.logspace(-3, 3, 100)\n",
    "dis = w.discharge(t)\n",
    "ax1 = plt.subplot(1, 2, 2)\n",
    "plt.semilogx(t, dis[0], label=\"$r_w$=0.2\")\n",
    "ml = tft.ModelMaq(kaq=25, z=[20, 0], Saq=1e-5, tmin=1e-3, tmax=1000)\n",
    "w = tft.HeadWell(ml, tsandh=[(0, -1)], rw=0.3)\n",
    "ml.solve()\n",
    "dis = w.discharge(t)\n",
    "plt.semilogx(t, dis[0], label=\"$r_w$=0.3\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"discharge (m3/d)\")\n",
    "ax0.legend()\n",
    "ax1.legend()\n",
    "ax1.grid()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}