{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "editable": true,
    "slideshow": {
     "slide_type": ""
    },
    "tags": []
   },
   "source": [
    "# HowTo: Model fluctuating head boundaries with `timflow`\n",
    "\n",
    "This notebook shows how fluctuating head boundaries, such as a head with a varying\n",
    "river, can be modeled with `timflow`. The example is for a cross-sectional model, but\n",
    "the same holds for a two-dimensional model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.transient as tft"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define parameters for a single unconfined aquifer. The head in the river fluctuates as a cosine with period $\\tau$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# parameters\n",
    "k = 20.0  # hydraulic conductivity, m/d\n",
    "H = 50.0  # thickness of aquifer, m\n",
    "T = k * H  # transmissivity, m^2/d\n",
    "S = 0.1  # storage coefficient, -\n",
    "tau = 0.5  # tidal period, d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytic solution for the head is given by the following function (see Bakker and Post, 2020)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solution for unit amplitude\n",
    "def head_river(t, tau, tp=0):\n",
    "    return np.cos(2 * np.pi * (t - tp) / tau)\n",
    "\n",
    "\n",
    "def analytic_head(x, t, tau, S, T, tp=0):\n",
    "    B = np.exp(-x * np.sqrt(S * np.pi / (T * tau)))\n",
    "    ts = x * np.sqrt(S * tau / (4 * np.pi * T))\n",
    "    return B * np.cos(2 * np.pi * (t - tp - ts) / tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The analytic solution is for a river with varying head that has been varying for ever. \n",
    "In `timflow.transient`, head is simulated for several days to get the model to spin-up. \n",
    "In interpreting the results, we only look at the heads in the final tidal period."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In `timflow.transient`, the specified head in a river is constant for a period. For example, `tsandh=[(0, 1), (1, 2), (4, -2)]` means the head in the river equal to 1 from $t=0$ till $t=1$, equal to 2 from $t=1$ till $t=4$ and equal to -2 thereafter. The head in a river fluctuates continuously, of course, and will be approximated by a step function. \n",
    "\n",
    "Consider the case that the head is specified using regular time intervals $\\Delta t$ (this is not necessary of course, but practical for most applications). \n",
    "There are three different ways to specify a continuously vary head boundary in `timflow`. \n",
    "\n",
    "- The head at time $t$ is used for the interval from $t$ till $t+\\Delta t$. This is referred to as *begin*.\n",
    "- The head at time $t+\\Delta t$ is used for the interval from $t$ till $t+\\Delta t$. This is referred to as *end*.\n",
    "- The average between the heads at $t$ and $t+\\Delta t$ is used for the interval from $t$ till $t+\\Delta t$. This is referred to as *mid*.\n",
    "\n",
    "A graph for these three options is shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmax = 5  # day\n",
    "delt = 0.025  # day\n",
    "t = np.arange(0, tmax, delt)\n",
    "hexact = head_river(t, tau)\n",
    "hbegin = head_river(t, tau)\n",
    "hend = head_river(t + delt, tau)\n",
    "hmid = 0.5 * (head_river(t - delt / 2, tau) + head_river(t + delt / 2, tau))\n",
    "tmid = np.hstack((0, 0.5 * (t[:-1] + t[1:])))\n",
    "# plot\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.plot(t, hexact, \"k\", label=\"observed\")\n",
    "plt.step(t, hbegin, where=\"post\", label=\"begin\")\n",
    "plt.step(t, hend, where=\"post\", label=\"end\")\n",
    "plt.step(tmid, hmid, where=\"post\", label=\"mid\")\n",
    "plt.xlim(0, 1)\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend()\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A model is created for all three options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlbegin = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "xsection = tft.XsectionMaq(\n",
    "    model=mlbegin, x1=-np.inf, x2=np.inf, kaq=k, z=[0, -H], Saq=S, phreatictop=True\n",
    ")\n",
    "lsriver = tft.River1D(model=mlbegin, xls=0, tsandh=list(zip(t, hbegin, strict=True)))\n",
    "mlbegin.solve(silent=True)\n",
    "#\n",
    "mlmid = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "xsection = tft.XsectionMaq(\n",
    "    model=mlmid, x1=-np.inf, x2=np.inf, kaq=k, z=[0, -H], Saq=S, phreatictop=True\n",
    ")\n",
    "lsriver = tft.River1D(model=mlmid, xls=0, tsandh=list(zip(tmid, hmid, strict=True)))\n",
    "mlmid.solve(silent=True)\n",
    "#\n",
    "mlend = tft.ModelXsection(naq=1, tmin=1e-4, tmax=1e2)\n",
    "xsection = tft.XsectionMaq(\n",
    "    model=mlend, x1=-np.inf, x2=np.inf, kaq=k, z=[0, -H], Saq=S, phreatictop=True\n",
    ")\n",
    "lsriver = tft.River1D(model=mlend, xls=0, tsandh=list(zip(t, hend, strict=True)))\n",
    "mlend.solve(silent=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's compare the output for all three options. First we look at the head along\n",
    "$x$ at 4 moments in the tidal period. As we can see the *mid* option best\n",
    "matches the analytical solution. The `timflow` and exact solutions are closer\n",
    "together when the interval $\\Delta t$ is reduced."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 200, 100)\n",
    "tlist = [4, 4.125, 4.25, 4.375]\n",
    "plt.figure(figsize=(10, 4))\n",
    "for tp in tlist:\n",
    "    hex = analytic_head(x, tp, tau, S, T)\n",
    "    plt.plot(x, hex, \"k\")\n",
    "    for i, ml in enumerate([mlbegin, mlmid, mlend]):\n",
    "        h = ml.headalongline(x, 0, tp)\n",
    "        plt.plot(x, h[0, 0], \"C\" + str(i), linestyle=\"--\")\n",
    "plt.text(-2, 1, \"t=4\", ha=\"right\", va=\"top\")\n",
    "plt.text(-2, 0.01, \"t=4.125\", ha=\"right\", va=\"bottom\")\n",
    "plt.text(-2, -0.01, \"t=4.375\", ha=\"right\", va=\"top\")\n",
    "plt.text(-2, -1, \"t=4.25\", ha=\"right\", va=\"bottom\")\n",
    "plt.xlim(-20, 200)\n",
    "plt.legend([\"exact\", \"begin\", \"mid\", \"end\"])\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we compare the head at $x=50$ m for the last tidal period. Here we can\n",
    "clearly see that the *begin* and *end* result in solutions that trail and lead\n",
    "the analytic solution, respectively. The *mid* option matches nicely with the\n",
    "analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 3))\n",
    "tp = np.linspace(4.5, 5, 100)\n",
    "hex = analytic_head(50, tp, tau, S, T)\n",
    "plt.plot(tp, hex, \"k\")\n",
    "for ml in [mlbegin, mlmid, mlend]:\n",
    "    h = ml.head(50, 0, tp)\n",
    "    plt.plot(tp, h[0], \"--\")\n",
    "plt.legend([\"exact\", \"begin\", \"mid\", \"end\"])\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.title(\"head at x=50 m\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In conclusion, when modeling fluctuating boundary conditions in `timflow`, use the *mid* approximation as the input time series for your model."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}