{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {},
   "source": [
    "# 133.1x: Confined flow with precipitation\n",
    "\n",
    "## Transient solution (133.16)\n",
    "This notebook shows the Bruggeman solution for:\n",
    "\n",
    "One dimensional finite flow. Given head or drawdown at `x=b` and zero flux at `x=0`.\n",
    "\n",
    "Steady state (133.17) shown as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import timflow.transient as tft\n",
    "from numpy import linspace\n",
    "\n",
    "from bruggeman.flow1d import bruggeman_133_16, bruggeman_133_17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "bruggeman_133_16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3",
   "metadata": {},
   "outputs": [],
   "source": [
    "bruggeman_133_16?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4",
   "metadata": {},
   "source": [
    "## Transient graph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# aquifer parameters\n",
    "b = 100.0 / 2\n",
    "S = 0.1\n",
    "k = 10.0\n",
    "D = 5.0\n",
    "p = 1e-3  # constant precipitation flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "t = linspace(0.0, 10.0, 100)\n",
    "hx0 = bruggeman_133_16(0.0, t, b, S, k, D, p)\n",
    "hxb2 = bruggeman_133_16(b / 2, t, b, S, k, D, p)\n",
    "\n",
    "plt.figure(figsize=(8, 3), layout=\"tight\")\n",
    "plt.plot(t, hx0, label=\"x=0\")\n",
    "plt.plot(t, hxb2, label=\"x=L/4\")\n",
    "\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Time (d)\")\n",
    "plt.ylabel(\"Head (m)\")\n",
    "plt.title(\"Bruggeman 133.16\")\n",
    "plt.xlim(t[0], t[-1])\n",
    "plt.ylim(0.0)\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {},
   "source": [
    "## Steady state solution (133.17)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "bruggeman_133_17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9",
   "metadata": {},
   "outputs": [],
   "source": [
    "bruggeman_133_17?"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "10",
   "metadata": {},
   "source": [
    "## Drainage to canals"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {},
   "source": [
    "This equation can be used to model drainage to drains and canals. \n",
    "\n",
    "Results are equal to Krayenhoff van de Leur - Maasland (e.g.: 'Cultuurtechnisch Vademecum', 1988, page 523)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0e10bc8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "mlconf = tft.ModelXsection(naq=1, tmin=1e-3, tmax=1e3)\n",
    "\n",
    "left = tft.XsectionMaq(\n",
    "    mlconf,\n",
    "    -np.inf,\n",
    "    -b,\n",
    "    kaq=k,\n",
    "    z=[0, -D],\n",
    "    Saq=S,\n",
    "    topboundary=\"phreatic\",\n",
    ")\n",
    "inf = tft.XsectionMaq(\n",
    "    mlconf,\n",
    "    -b,\n",
    "    b,\n",
    "    kaq=k,\n",
    "    z=[0, -D],\n",
    "    Saq=S,\n",
    "    # c=c,\n",
    "    topboundary=\"phreatic\",\n",
    "    tsandN=[(0.0, p)],\n",
    ")\n",
    "right = tft.XsectionMaq(\n",
    "    mlconf,\n",
    "    b,\n",
    "    np.inf,\n",
    "    kaq=k,\n",
    "    z=[0, -D],\n",
    "    Saq=S,\n",
    "    # c=c,\n",
    "    topboundary=\"phreatic\",\n",
    ")\n",
    "\n",
    "d = -1e-3\n",
    "hls_left = tft.River1D(mlconf, xls=-b + d, tsandh=[(0, 0.0)], layers=[0])\n",
    "hls_right = tft.River1D(mlconf, xls=b - d, tsandh=[(0, 0.0)], layers=[0])\n",
    "\n",
    "mlconf.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "x = linspace(-b, b, 51)\n",
    "t_steps = [0, 0.5, 1.0, 2, 5, 10]  # time steps in days\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(8, 3), layout=\"tight\")\n",
    "for t in t_steps:\n",
    "    ht = bruggeman_133_16(x, t, b, S, k, D, p)\n",
    "    (p0,) = ax.plot(x, ht, label=f\"t={t} d\")\n",
    "    ax.plot(\n",
    "        x,\n",
    "        mlconf.headalongline(x, np.zeros_like(x), t)[0, 0],\n",
    "        label=\"Timflow\",\n",
    "        marker=\"x\",\n",
    "        linestyle=\"none\",\n",
    "        markersize=4,\n",
    "        color=p0.get_color(),\n",
    "    )\n",
    "\n",
    "h_ss = bruggeman_133_17(x, b, k, D, p)\n",
    "ax.plot(x, h_ss, label=\"Steady state\", color=\"k\", linestyle=\"--\")\n",
    "\n",
    "ax.grid(True)\n",
    "ax.set_xlabel(\"Distance (m)\")\n",
    "ax.set_ylabel(\"Head (m)\")\n",
    "ax.set_title(\"Bruggeman 133.16 (transient) and 133.17 (steady state)\", y=1.45)\n",
    "ax.set_ylim(0.0)\n",
    "ax.legend(loc=(0, 1), ncol=7, fontsize=\"small\", frameon=False);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32baf003",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bruggeman (broken)",
   "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.13.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}