{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# LineSink1D and LineDoublet1D\n",
    "\n",
    "\n",
    "This notebook contains the examples of the LineSink1D and LineDoublet1D elements."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy.special import erfc\n",
    "\n",
    "import timflow.steady as tfs\n",
    "import timflow.transient as tft\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (10, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the analytical solutions for 1D confined flow in a semi-infinite field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ierfc(z, n):\n",
    "    if n == -1:\n",
    "        return 2 / np.sqrt(np.pi) * np.exp(-z * z)\n",
    "    elif n == 0:\n",
    "        return erfc(z)\n",
    "    else:\n",
    "        result = -z / n * ierfc(z, n - 1) + 1 / (2 * n) * ierfc(z, n - 2)\n",
    "        return np.clip(result, a_min=0.0, a_max=None)\n",
    "\n",
    "\n",
    "def bruggeman_123_02(x, t, dh, k, H, S):\n",
    "    \"\"\"Solution for sudden rise of the water table in a confined aquifer.\n",
    "\n",
    "    From Bruggeman 123.02\n",
    "    \"\"\"\n",
    "    beta = np.sqrt(S / (k * H))\n",
    "    u = beta * x / (2 * np.sqrt(t))\n",
    "    return dh * erfc(u)\n",
    "\n",
    "\n",
    "def bruggeman_123_03(x, t, a, k, H, S):\n",
    "    \"\"\"Solution for linear rise of the water table in a confined aquifer.\n",
    "\n",
    "    From Bruggeman 123.03\n",
    "    \"\"\"\n",
    "    beta = np.sqrt(S / (k * H))\n",
    "    u = beta * x / (2 * np.sqrt(t))\n",
    "    return a * t * ierfc(u, 2) / ierfc(0, 2)\n",
    "\n",
    "\n",
    "def bruggeman_123_05_q(x, t, b, k, H, S):\n",
    "    \"\"\"Solution for constant infiltration/extraction in a confined aquifer.\n",
    "\n",
    "    From Olsthoorn, Th. 2006. Van Edelman naar Bruggeman. Stromingen 12 (2006) p5-11.\n",
    "    \"\"\"\n",
    "    beta = np.sqrt(S / (k * H))\n",
    "    u = beta * x / (2 * np.sqrt(t))\n",
    "    s = 2 * b * np.sqrt(t) / np.sqrt(k * H * S) * ierfc(u, 1) / (ierfc(0, 0))\n",
    "    return s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from Analyical Groundwater Modeling, ch. 5\n",
    "\n",
    "\n",
    "def h_edelman(x, t, T, S, delh, t0=0.0):\n",
    "    u = np.sqrt(S * x**2 / (4 * T * (t - t0)))\n",
    "    return delh * erfc(u)\n",
    "\n",
    "\n",
    "def Qx_edelman(x, t, T, S, delh, t0=0.0):\n",
    "    u = np.sqrt(S * x**2 / (4 * T * (t - t0)))\n",
    "    return T * delh * 2 * u / (x * np.sqrt(np.pi)) * np.exp(-(u**2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check implementation of the iterated integral complementary error function.\n",
    "\n",
    "Compared to figure 3 in [Olsthoorn, (2006)](https://edepot.wur.nl/13730)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.linspace(0, 3, 51)\n",
    "ierfc_min1 = ierfc(z, -1)\n",
    "ierfc0 = ierfc(z, 0)\n",
    "plt.plot(z, ierfc_min1, ls=\"dashed\", label=\"ierfc(z, -1)\", lw=1.0, color=\"k\")\n",
    "plt.plot(z, ierfc0, ls=\"solid\", label=\"ierfc(z, 0)\", lw=1.0, color=\"k\")\n",
    "for n in range(1, 8):\n",
    "    plt.plot(z, ierfc(z, n) / ierfc(0, n), label=f\"ierfc(z, {n}) / ierfc(0, {n})\", lw=1.0)\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single layer confined aquifer\n",
    "\n",
    "Example problem confined aquifer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 5\n",
    "H = 10\n",
    "Ss = 1e-3 / H\n",
    "Q = 2.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlconf = tft.ModelMaq(kaq=k, z=[0, -H], Saq=Ss, tmin=1e-3, tmax=1e3, topboundary=\"conf\")\n",
    "ls = tft.DischargeLineSink1D(mlconf, tsandq=[(0, Q)], layers=[0])\n",
    "mlconf.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot head along $x$ for different times $t$. Compare to analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 100, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(-1, 0, 5)\n",
    "\n",
    "for i in range(len(t)):\n",
    "    h = mlconf.headalongline(x, y, t[i])\n",
    "    plt.plot(x, h.squeeze(), label=f\"t={t[i]:.2f} d\")\n",
    "    ha = bruggeman_123_05_q(x, t[i], -Q / 2, k, H, Ss * H)  # Q/2 because 2-sided flow?\n",
    "    plt.plot(x, ha, \"k:\")\n",
    "plt.plot([], [], \"k:\", label=\"Bruggeman 123.05\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"drawdown [m]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Single layer semi-confined aquifer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(kaq=k, z=[1, 0, -H], c=[10], tmin=1e-3, tmax=1e3, topboundary=\"semi\")\n",
    "ls = tft.DischargeLineSink1D(ml, tsandq=[(0, Q)], layers=[0])\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare to timflow steady-state solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlss = tfs.ModelMaq(kaq=k, z=[1, 0, -H], c=[10], topboundary=\"semi\", hstar=0.0)\n",
    "ls = tfs.LineSink1D(mlss, 0.0, sigls=Q)\n",
    "mlss.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot head along $x$ for different times $t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 100, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(-2, 0, 5)\n",
    "\n",
    "for i in range(len(t)):\n",
    "    h = ml.headalongline(x, y, t[i])\n",
    "    plt.plot(x, h.squeeze(), label=f\"t={t[i]:.2f} d\")\n",
    "\n",
    "hss = mlss.headalongline(x, y)\n",
    "plt.plot(x, hss.squeeze(), \"k--\", label=\"Steady-state\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"drawdown [m]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sudden change in water level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from Analyical Groundwater Modeling, ch. 5, p. 72\n",
    "k = 10.0\n",
    "H = 10.0\n",
    "S = 0.2\n",
    "delh = 2.0\n",
    "t0 = 0.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlconf = tft.ModelMaq(\n",
    "    kaq=k, z=[0, -H], Saq=S, tmin=1, tmax=1e2, topboundary=\"conf\", phreatictop=True\n",
    ")\n",
    "hls = tft.River1D(mlconf, tsandh=[(0, delh)], layers=[0])\n",
    "mlconf.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute head along $x$ for t = 1, 10 and 100 days. Compare to analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 200, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(0, 2, 3)\n",
    "\n",
    "for i in range(len(t)):\n",
    "    h = mlconf.headalongline(x, y, t[i])\n",
    "    plt.plot(x, h.squeeze(), label=f\"t={t[i]:.0f} d\")\n",
    "    ha = h_edelman(x, t[i], k * H, S, delh, t0)\n",
    "    plt.plot(x, ha, \"k:\")\n",
    "\n",
    "plt.plot([], [], c=\"k\", ls=\"dotted\", label=\"Edelman\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Leaky Wall\n",
    "\n",
    "Use the previous problem but add a leaky wall at x=100 m with resistance 10 days."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from Analyical Groundwater Modeling, ch. 5, p. 72\n",
    "k = 10.0\n",
    "H = 10.0\n",
    "S = 0.2\n",
    "delh = 2.0\n",
    "t0 = 0.0\n",
    "\n",
    "# add leaky wall with resistance\n",
    "res = 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlconf = tft.ModelMaq(\n",
    "    kaq=k, z=[0, -H], Saq=S, tmin=1, tmax=1e2, topboundary=\"conf\", phreatictop=True\n",
    ")\n",
    "hls = tft.River1D(mlconf, xls=0.0, tsandh=[(0, delh)], layers=[0])\n",
    "ld = tft.LeakyWall1D(mlconf, xld=100.0, res=res, layers=[0])\n",
    "mlconf.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 200, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(0, 2, 3)\n",
    "\n",
    "for i in range(len(t)):\n",
    "    h = mlconf.headalongline(x, y, t[i])\n",
    "    plt.plot(x, h.squeeze(), label=f\"t={t[i]:.0f} d\")\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel(\"x [m]\")\n",
    "plt.ylabel(\"head [m]\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check boundary condition along wall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 1e-6\n",
    "t = np.logspace(0, 2, 100)\n",
    "hleft = mlconf.head(100 - dx, 0, t)\n",
    "hright = mlconf.head(100 + dx, 0, t)\n",
    "disxleft, _ = mlconf.disvec(100 - dx, 0, t)\n",
    "disxright, _ = mlconf.disvec(100 + dx, 0, t)\n",
    "disxnum = mlconf.aq.Haq * (hleft - hright) / res\n",
    "plt.semilogx(t, hleft[0], label=\"hleft\")\n",
    "plt.semilogx(t, hright[0], label=\"hright\")\n",
    "plt.semilogx(t, disxleft[0], label=\"Qxleft timflow\")\n",
    "plt.semilogx(t, disxright[0], \"--\", label=\"Qxright timflow\")\n",
    "plt.semilogx(t, disxnum[0], \"k:\", label=\"Qx numerical\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head or Qx\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multi-layer system"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = [0.0, -5.0, -6.0, -15.0]\n",
    "Saq = [0.2, 1e-4]\n",
    "c = [10.0]\n",
    "k = [5.0, 10.0]\n",
    "delh = 1.0\n",
    "res = 10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlconf = tft.ModelMaq(\n",
    "    kaq=k, z=z, Saq=Saq, c=c, tmin=1, tmax=1e2, topboundary=\"conf\", phreatictop=True\n",
    ")\n",
    "hls_left = tft.River1D(mlconf, xls=0.0, tsandh=[(0, delh)], layers=[0, 1])\n",
    "hls_right = tft.River1D(mlconf, xls=200.0, tsandh=[(0, 0.0)], layers=[0, 1])\n",
    "ld = tft.LeakyWall1D(mlconf, xld=100.0, res=res, layers=[0])\n",
    "mlconf.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = mlconf.plots.xsection([(-10, 0), (210, 0)], params=True)\n",
    "for e in mlconf.elementlist:\n",
    "    e.plot(ax=ax)\n",
    "ax.set_xlim(-10, 210);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 200, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(0, 2, 3)\n",
    "\n",
    "fig, (ax0, ax1) = plt.subplots(2, 1, sharex=True, sharey=True, figsize=(10, 5))\n",
    "\n",
    "for i in range(len(t)):\n",
    "    h = mlconf.headalongline(x, y, t[i])\n",
    "    ax0.plot(x, h[0].squeeze(), label=f\"t={t[i]:.0f} d\")\n",
    "    ax1.plot(x, h[1].squeeze(), label=f\"t={t[i]:.0f} d\", ls=\"dashed\")\n",
    "\n",
    "ax0.legend(loc=(0, 1), frameon=False, ncol=3)\n",
    "ax0.set_title(\"Layer 0\")\n",
    "ax1.set_xlabel(\"x [m]\")\n",
    "ax1.set_title(\"Layer 1\")\n",
    "for iax in [ax0, ax1]:\n",
    "    iax.set_ylabel(\"head [m]\")\n",
    "    iax.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dx = 1e-6\n",
    "t = np.logspace(0, 2, 100)\n",
    "hleft = mlconf.head(100 - dx, 0, t)\n",
    "hright = mlconf.head(100 + dx, 0, t)\n",
    "disxleft, _ = mlconf.disvec(100 - dx, 0, t)\n",
    "disxright, _ = mlconf.disvec(100 + dx, 0, t)\n",
    "disxnum = mlconf.aq.Haq[:, np.newaxis] * (hleft - hright) / res\n",
    "plt.figure(figsize=(10, 5))\n",
    "for i in range(2):\n",
    "    plt.subplot(2, 1, i + 1)\n",
    "    plt.semilogx(t, hleft[i], label=\"hleft\")\n",
    "    plt.semilogx(t, hright[i], \"--\", label=\"hright\")\n",
    "    plt.semilogx(t, disxleft[i], label=\"Qxleft timflow\")\n",
    "    plt.semilogx(t, disxright[i], \"--\", label=\"Qxright timflow\")\n",
    "    if i == 0:  # otherwise no leaky wall\n",
    "        plt.semilogx(t, disxnum[i], \"k:\", label=\"Qx numerical\")\n",
    "    plt.xlabel(\"time (d)\")\n",
    "    plt.ylabel(\"head or Qx\")\n",
    "    plt.title(f\"layer {i}\")\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Multiscreen Linesink1D\n",
    "\n",
    "Linesink pumping from multiple layers simultaneously with specified total discharge."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = [0.0, -5.0, -6.0, -11.0]\n",
    "Saq = [1e-3, 1e-4]\n",
    "c = [10.0]\n",
    "k = [5.0, 10.0]\n",
    "qls = 6.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tft.ModelMaq(\n",
    "    kaq=k,\n",
    "    z=z,\n",
    "    Saq=Saq,\n",
    "    c=c,\n",
    "    tmin=1,\n",
    "    tmax=1e2,\n",
    "    topboundary=\"conf\",\n",
    ")\n",
    "\n",
    "ls = tft.LineSink1D(ml, xls=0.0, tsandq=[(0.0, qls)], layers=[0, 1])\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = ml.plots.xsection(xy=[(-10, 0), (10, 0)], params=True)\n",
    "ml.elementlist[0].plot(ax=ax)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 200, 101)\n",
    "y = np.zeros_like(x)\n",
    "t = np.logspace(0, 2, 3)\n",
    "\n",
    "h = ml.headalongline(x, y, t)\n",
    "qx, _ = ml.disvecalongline(x, y, t)\n",
    "\n",
    "fig, (ax0, ax1) = plt.subplots(2, 1, sharex=True, figsize=(10, 5))\n",
    "\n",
    "for i in range(len(t)):\n",
    "    (p0,) = ax0.plot(x, h[0, i].squeeze(), label=f\"t={t[i]:.0f} d\")\n",
    "    ax0.plot(x, h[1, i].squeeze(), label=f\"t={t[i]:.0f} d\", ls=\"dashed\", c=p0.get_color())\n",
    "    (p1,) = ax1.plot(x, qx[0, i].squeeze(), label=f\"t={t[i]:.0f} d\")\n",
    "    ax1.plot(\n",
    "        x, qx[1, i].squeeze(), label=f\"t={t[i]:.0f} d\", ls=\"dashed\", c=p1.get_color()\n",
    "    )\n",
    "\n",
    "ax0.legend(loc=(0, 1), frameon=False, ncol=3)\n",
    "ax0.set_ylabel(\"head [m]\")\n",
    "ax1.set_xlabel(\"x [m]\")\n",
    "ax1.set_ylabel(\"q [m$^2$/d]\")\n",
    "for iax in [ax0, ax1]:\n",
    "    iax.grid()\n",
    "fig.align_ylabels()"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}