{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test circular building pit\n",
    "\n",
    "Comparing timflow solutions with `LeakyWalls` and `LeakyBuildingPit` to the exact\n",
    "analytical solution for a circular building pit in a semi-confined aquifer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import packages\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import shapely\n",
    "from scipy.special import i0, i1, k0, k1\n",
    "\n",
    "import timflow.steady as tfs"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select a radius for the circular building pit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = 100.0  # radius of building pit"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify a resolution to control the number of leaky segments used in the timflow models to approximate the exact  analytical solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# subdivide circle into segments, increase resolution to increase no. of elements\n",
    "resolution = 5  # resolution=5 results in 20 segments\n",
    "circle = shapely.Point(0, 0).buffer(R, resolution=resolution)\n",
    "print(\"Number of pts      :\", len(circle.exterior.xy[0]))\n",
    "print(\"Number of segments :\", len(circle.exterior.xy[0]) - 1)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show the points used to approximate the circular building pit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shapely.MultiPoint(shapely.points(shapely.get_coordinates(circle)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create our model we create a list of x,y-coordinates, ordered counter-clockwise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# flip so coordinates are ordered counter clockwise\n",
    "xy = list(\n",
    "    zip(\n",
    "        shapely.get_coordinates(circle)[:, 0],\n",
    "        shapely.get_coordinates(circle)[:, 1],\n",
    "        strict=False,\n",
    "    )\n",
    ")[::-1]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define model parameters:\n",
    "- horizontal and vertical hydraulic conductivity in m/d\n",
    "- top resistance of semi-confined aquifer in days \n",
    "- top and bottom elevations of the aquifer\n",
    "- resistance of the leaky wall in days\n",
    "- order of the elements in timflow (which specifies the no. of control points used in the solution)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model parameters\n",
    "\n",
    "kh = 10  # m/day\n",
    "kv = 0.25 * kh  # m/day\n",
    "\n",
    "ctop = 1000.0  # resistance top leaky layer in days\n",
    "\n",
    "ztop = 0.0  # surface elevation\n",
    "zbot = -20.0  # bottom elevation of the model\n",
    "z = np.array([ztop + 1, ztop, zbot])\n",
    "\n",
    "res = 100.0  # resistance of leaky wall, in days\n",
    "\n",
    "rw = 0.3  # well radius, in m\n",
    "Qw = 100.0  # well discharge, in m3/d\n",
    "\n",
    "o = 7  # order"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a timflow model with a circular building pit using the `LeakyWallString` \n",
    "element. The building pit contains a well at location $(0, 0)$ with a discharge of 100\n",
    "$m^3/d$.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Note</b>\n",
    "\n",
    "The <code>LeakyWallString</code> can be used to model the building pit in this\n",
    "example since the aquifer properties are the same inside and outside the building pit.\n",
    "If aquifer properties change inside the building pit, only the `BuildingPit` elements\n",
    "can be used.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml_lld = tfs.ModelMaq(kaq=kh, z=z, c=ctop, topboundary=\"semi\", hstar=0.0)\n",
    "lld = tfs.LeakyWallString(\n",
    "    ml_lld,\n",
    "    xy=xy,\n",
    "    res=res,\n",
    "    layers=[0],\n",
    "    order=o,\n",
    ")\n",
    "well = tfs.Well(ml_lld, 0.0, 0.0, Qw=Qw, rw=rw)\n",
    "ml_lld.solve()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we build the same model, but using the `LeakyBuildingPit` element. This element\n",
    "is technically an inhomogeneity (allowing different aquifer parameters inside the \n",
    "element as compared to the rest of the model), but for comparison purposes we keep \n",
    "the aquifer homogeneous. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=kh, z=z, c=ctop, topboundary=\"semi\", hstar=0.0)\n",
    "bpit = tfs.LeakyBuildingPitMaq(\n",
    "    ml,\n",
    "    xy,\n",
    "    kaq=kh,\n",
    "    z=z,\n",
    "    topboundary=\"semi\",\n",
    "    hstar=0.0,\n",
    "    c=ctop,\n",
    "    layers=[0],\n",
    "    res=res,\n",
    "    order=o,\n",
    ")\n",
    "well = tfs.Well(ml, 0.0, 0.0, Qw=Qw, rw=rw)\n",
    "ml.solve()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the exact analytical solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# translate some of the model parameters defined earlier\n",
    "k = kh  # m/d\n",
    "H = ztop - zbot  # m\n",
    "c = ctop  # d\n",
    "cwall = res  # d\n",
    "Q = Qw  # m^3/d\n",
    "\n",
    "# computed values\n",
    "T = k * H\n",
    "lab = np.sqrt(c * T)\n",
    "C = H * lab / (cwall * T)\n",
    "I0 = i0(R / lab)\n",
    "I1 = i1(R / lab)\n",
    "K0 = k0(R / lab)\n",
    "K1 = k1(R / lab)\n",
    "B = -Q * (K1 * I0 + I1 * K0) / (K0 * I1 + K1 * I0 + I1 * K1 / C)\n",
    "A = -(Q + B) * K1 / I1"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define functions for the head and discharge as a function of the radial distance $r$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def head_nowall(r):  # for comparison\n",
    "    return -Q / (2 * np.pi * T) * k0(r / lab)\n",
    "\n",
    "\n",
    "def head(r):\n",
    "    if r < R:\n",
    "        h = -Q / (2 * np.pi * T) * k0(r / lab) + A / (2 * np.pi * T) * i0(r / lab)\n",
    "    else:\n",
    "        h = B / (2 * np.pi * T) * k0(r / lab)\n",
    "    return h\n",
    "\n",
    "\n",
    "def disr(r):\n",
    "    if r < R:\n",
    "        Qr = -Q / (2 * np.pi * lab) * k1(r / lab) - A / (2 * np.pi * lab) * i1(r / lab)\n",
    "    else:\n",
    "        Qr = B / (2 * np.pi * lab) * k1(r / lab)\n",
    "    return Qr\n",
    "\n",
    "\n",
    "headvec = np.vectorize(head)\n",
    "disrvec = np.vectorize(disr)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the analytical solution for the head (left) and radial discharge (right) as a function of radial distance $r$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.linspace(rw, 1.5 * R, 301)\n",
    "h = headvec(r)\n",
    "Qr = disrvec(r) * 2 * np.pi * r\n",
    "plt.figure(figsize=(10, 3))\n",
    "plt.subplot(121)\n",
    "plt.plot(r, h, label=\"head with wall\")\n",
    "plt.plot(r, head_nowall(r), \"--\", label=\"head no wall\")\n",
    "plt.xlabel(\"r (m)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.subplot(122)\n",
    "plt.plot(r, Qr)\n",
    "plt.xlabel(\"r (m)\")\n",
    "plt.ylabel(\"$Q_r$ (m$^2$/d)\")\n",
    "plt.grid()\n",
    "plt.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare the exact analytical solution to the timflow models from earlier.\n",
    "\n",
    "We specify the angle of the line relative to the positive x-axis along which the heads (and radial discharge) are computed in the timflow models. At the endpoints the radius of the building pit is exactly equal to $R$, whereas in between two endpoints, the radius is slightly smaller. The timflow solutions at the endpoints can show some deviations from the exact solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle = 0.0  # in degrees"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we compare the head. In the top plot, the heads are plotted as a function of radial distance $r$. In the second plot the differences between the timflow models and the exact solutions are shown.\n",
    "\n",
    "The models correspond closely to the exact solution with differences on the order of $10^{-3}$ inside the building pit, and even smaller outside the building pit. The timflow model with the `LeakyBuildingPit` element is slightly more accurate than the `LeakyWall` solution, though differences are small. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.linspace(rw, 1.5 * R, 301)\n",
    "\n",
    "sin = np.sin(np.deg2rad(angle))\n",
    "cos = np.cos(np.deg2rad(angle))\n",
    "xl = np.linspace(rw * cos, cos * 1.5 * R, 301)\n",
    "yl = np.linspace(rw * sin, sin * 1.5 * R, 301)\n",
    "rl = np.sqrt((xl - 0.0) ** 2 + (yl - 0.0) ** 2)\n",
    "\n",
    "h = headvec(r)\n",
    "h_lld = ml_lld.headalongline(xl, yl)\n",
    "h_lp = ml.headalongline(xl, yl)\n",
    "\n",
    "fig, (ax, ax2) = plt.subplots(2, 1, figsize=(10, 5))\n",
    "ax.plot(r, h, color=\"k\", ls=\"dashed\", label=\"exact\")\n",
    "ax.plot(r, h_lp.squeeze(), color=\"C0\", label=\"LeakyBuildingPit\")\n",
    "ax.plot(r, h_lld.squeeze(), color=\"C1\", ls=\"dashed\", label=\"LeakyWall\")\n",
    "ax.legend(loc=(0, 1), ncol=3, frameon=False)\n",
    "ax.set_ylabel(\"head (m)\")\n",
    "ax.grid(True)\n",
    "\n",
    "ax2.axhline(0.0, linestyle=\"dashed\", color=\"k\", lw=0.75)\n",
    "ax2.plot(r, h_lp.squeeze() - h, color=\"C0\", label=\"LeakyBuildingPit\")\n",
    "ax2.plot(r, h_lld.squeeze() - h, color=\"C1\", ls=\"dashed\", label=\"LeakyWall\")\n",
    "ax2.legend(loc=(0, 1), ncol=3, frameon=False)\n",
    "ax2.grid(True)\n",
    "ax2.set_ylabel(\"difference (m)\")\n",
    "# ax2.set_ylim(-0.005, 0.005)\n",
    "\n",
    "fig.align_ylabels()\n",
    "fig.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we compare the computed radial discharge between the timflow models and the exact solution.\n",
    "\n",
    "The top plot shows the radial discharge, and the bottom plot shows the difference between the timflow models and the exact solution. At the wall the computed radial discharge differs from the exact solution. This is caused by the implementation of the elements in timflow, where control points along a line are specified at which (or between which) certain conditions must be satisfied. The solution at a particular point close to the element may not be entirely accurate. For the `LeakyBuildingPit` implementation, this results in the solution at the edges of the segments to become inaccurate. However, integrating the flux along the entire circular building pit should yield relatively accurate results for the radial discharge, as we will see in the next step. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Qr = disrvec(r) * 2 * np.pi * r\n",
    "qx_lld, qy_lld = ml_lld.disvecalongline(xl, yl)\n",
    "qx_lp, qy_lp = ml.disvecalongline(xl, yl)\n",
    "\n",
    "Q_lld = (qx_lld * cos + qy_lld * sin).squeeze() * 2 * np.pi * rl\n",
    "Q_lp = (qx_lp * cos + qy_lp * sin).squeeze() * 2 * np.pi * rl\n",
    "\n",
    "fig, (ax, ax2) = plt.subplots(2, 1, figsize=(10, 5))\n",
    "ax.plot(r, Qr, color=\"k\", ls=\"dashed\", label=\"exact\")\n",
    "ax.plot(r, Q_lp.squeeze(), color=\"C0\", label=\"LeakyBuildingPit\")\n",
    "ax.plot(r, Q_lld.squeeze(), color=\"C1\", ls=\"dashed\", label=\"LeakyWall\")\n",
    "ax.legend(loc=(0, 1), ncol=3, frameon=False)\n",
    "ax.set_ylabel(\"$Q_r$ (m$^2$/d)\")\n",
    "ax.grid(True)\n",
    "\n",
    "ax2.axhline(0.0, linestyle=\"dashed\", color=\"k\", lw=0.75)\n",
    "ax2.plot(r, Q_lp.squeeze() - Qr, color=\"C0\", label=\"LeakyBuildingPit\")\n",
    "ax2.plot(r, Q_lld.squeeze() - Qr, color=\"C1\", ls=\"dashed\", label=\"LeakyWall\")\n",
    "ax2.legend(loc=(0, 1), ncol=3, frameon=False)\n",
    "ax2.grid(True)\n",
    "ax2.set_ylabel(\"difference (m$^2$/d)\")\n",
    "# ax2.set_ylim(-0.25, 0.25)\n",
    "\n",
    "fig.align_ylabels()\n",
    "fig.tight_layout()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check the solution, we integrate the normal flux along each segment of the circular building pit, and compare it to the exact solution. The calculated fluxes with both `timflow` models should be close to the exact solution, at some very small distance inside and outside the building pit.\n",
    "\n",
    "The `Leakywall` solution shows that the discharge is continuous across the element. The `LeakyBuildingPit` on the other hand shows a small jump in discharge. Both solutions differ slightly from the exact analytical solution, though this is also partly caused by the imperfect representation of the circular buildingpit using N line segments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nudge = 1e-3  # some small distance outside and inside the building pit\n",
    "ndeg = 20  # no. of legendre polynomial terms for integration\n",
    "\n",
    "# get x,y coordinates inside and outside of building pit\n",
    "p = shapely.Polygon(xy).exterior.buffer(nudge, join_style=\"mitre\")\n",
    "xyout = np.array(p.exterior.xy).T\n",
    "xyin = np.array(p.interiors[0].xy).T\n",
    "\n",
    "# flip interior points to be ordered the same as the outside points\n",
    "xyin = xyin[::-1]\n",
    "\n",
    "# calculate integrated normal fluxes inside/outside with both quad and legendre\n",
    "Qtot_in_lld = np.sum(ml_lld.intnormflux(xyin, ndeg=ndeg))\n",
    "Qtot_in_lp = np.sum(ml.intnormflux(xyin, ndeg=ndeg))\n",
    "Qtot_out_lld = np.sum(ml_lld.intnormflux(xyout, ndeg=ndeg))\n",
    "Qtot_out_lp = np.sum(ml.intnormflux(xyout, ndeg=ndeg))\n",
    "\n",
    "# shouldn't matter too much but to be exact add/subtract nudge\n",
    "# for the analytical solution as well:\n",
    "Qexact_in = disr(100 - nudge) * 2 * np.pi * (100.0 - nudge)\n",
    "Qexact_out = disr(100 + nudge) * 2 * np.pi * (100.0 + nudge)\n",
    "\n",
    "# print results\n",
    "df = pd.DataFrame(\n",
    "    index=[\"inside\", \"outside\"],\n",
    "    columns=[\"LeakyWall\", \"LeakyBuildingPit\", \"Exact\"],\n",
    ")\n",
    "df.loc[\"inside\"] = Qtot_in_lld, Qtot_in_lp, Qexact_in\n",
    "df.loc[\"outside\"] = Qtot_out_lld, Qtot_out_lp, Qexact_out\n",
    "df.loc[\"difference\"] = df.loc[\"inside\"] - df.loc[\"outside\"]\n",
    "df.index.name = \"Discharge\"\n",
    "df.columns.name = \"Model\"\n",
    "df.style.format(precision=2).set_caption(\n",
    "    \"Discharge along inside and outside of building pit:\"\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}