{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building Pit\n",
    "\n",
    "This example demonstrates the modeling of a building pit with dewatering wells using the TimML library. The model calculates groundwater flow towards the pit and evaluates the effectiveness of dewatering strategies by computing total discharge and drawdown distances.\n",
    "\n",
    "<div class=\"alert alert-info\">\n",
    "\n",
    "<b>Note</b>\n",
    "\n",
    "It is highly recommended to use the <code>BuildingPit</code> element if you want to\n",
    "implement different layer properties inside the building pit as compared to the rest of\n",
    "the aquifer. Adding barriers, e.g. `LeakyLineDoublets`, around an inhomogeneity will\n",
    "not work.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timml as tml"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters\n",
    "\n",
    "Define some aquifer parameters "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kh = 2.0  # m/day\n",
    "f_ani = 1 / 10  # anisotropy factor\n",
    "kv = f_ani * kh\n",
    "ctop = 800.0  # resistance top leaky layer in days"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define aquifer top and bottom, the depth of the sheetpile wall and the position of the dewatering wells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ztop = 0.0  # surface elevation\n",
    "z_well = -13.0  # end depth of the wellscreen\n",
    "z_dw = -15.0  # bottom elevation of sheetpile wall\n",
    "z_extra = z_dw - 15.0  # extra layer\n",
    "zbot = -60.0  # bottom elevation of the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Size of the building pit, and the required drawdown at the center of the building pit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "length = 40.0  # length building pit in m\n",
    "width = 30.0  # width building pit in m\n",
    "\n",
    "h_bem = -6.21  # m\n",
    "offset = 5.0  # distance groundwater extraction element from sheetpiles in m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy = [\n",
    "    (-length / 2, -width / 2),\n",
    "    (length / 2, -width / 2),\n",
    "    (length / 2, width / 2),\n",
    "    (-length / 2, width / 2),\n",
    "    (-length / 2, -width / 2),\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the building pit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(p2,) = plt.plot(*np.array(xy).T, \"o-\")\n",
    "plt.axis(\"equal\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model\n",
    "Set up a model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "z = np.array([ztop + 1, ztop, z_dw, z_dw, z_extra, z_extra, zbot])\n",
    "dz = z[1::2] - z[2::2]\n",
    "kh_arr = kh * np.ones(dz.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resistances of the top confining layer and aquitards"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = np.r_[np.array([ctop]), dz[:-1] / (2 * kv) + dz[1:] / (2 * kv)]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tml.ModelMaq(kaq=kh_arr, z=z, c=c, topboundary=\"semi\", hstar=0.0)\n",
    "\n",
    "layers = np.arange(np.sum(z_dw <= ml.aq.zaqbot))\n",
    "last_lay_dw = layers[-1]\n",
    "\n",
    "inhom = tml.BuildingPitMaq(\n",
    "    ml,\n",
    "    xy,\n",
    "    kaq=kh_arr,\n",
    "    z=z[1:],\n",
    "    topboundary=\"conf\",\n",
    "    c=c[1:],\n",
    "    order=4,\n",
    "    ndeg=3,\n",
    "    layers=layers,\n",
    ")\n",
    "\n",
    "tml.River(\n",
    "    ml,\n",
    "    x1=-length / 2 + offset,\n",
    "    y1=width / 2 - offset,\n",
    "    x2=length / 2 - offset,\n",
    "    y2=width / 2 - offset,\n",
    "    hls=h_bem,\n",
    "    layers=np.arange(last_lay_dw + 1),\n",
    ")\n",
    "tml.River(\n",
    "    ml,\n",
    "    x1=-length / 2 + offset,\n",
    "    y1=0,\n",
    "    x2=length / 2 - offset,\n",
    "    y2=0,\n",
    "    hls=h_bem,\n",
    "    layers=np.arange(last_lay_dw + 1),\n",
    ")\n",
    "tml.River(\n",
    "    ml,\n",
    "    x1=-length / 2 + offset,\n",
    "    y1=-width / 2 + offset,\n",
    "    x2=length / 2 - offset,\n",
    "    y2=-width / 2 + offset,\n",
    "    hls=h_bem,\n",
    "    layers=np.arange(last_lay_dw + 1),\n",
    ")\n",
    "\n",
    "# ml.solve_mp(nproc=2)\n",
    "ml.solve()\n",
    "\n",
    "Qtot = 0.0\n",
    "\n",
    "for e in ml.elementlist:\n",
    "    if e.name == \"River\":\n",
    "        Qtot += e.discharge()\n",
    "\n",
    "print(\"\\nDischarge =\", np.round(Qtot.sum(), 2), \"m^3/day\")\n",
    "\n",
    "y = np.linspace(-width / 2 - 25, width / 2 + 1100, 201)\n",
    "hl = ml.headalongline(np.zeros(201), y, layers=[0])\n",
    "y_5cm = np.interp(\n",
    "    -0.05, ml.headalongline(np.zeros(201), y, layers=0).squeeze(), y, right=np.nan\n",
    ")\n",
    "print(\"Distance to 5 cm drawdown contour =\", np.round(y_5cm, 2), \"m\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot an overview of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.topview()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visualizations\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0.0, length / 2 + 1100, 201)\n",
    "hl = ml.headalongline(x, np.zeros(201), layers=[last_lay_dw, last_lay_dw + 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n",
    "ax.plot(x, hl[0].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n",
    "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n",
    "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n",
    "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n",
    "ax.set_xlabel(\"x (m)\")\n",
    "ax.set_ylabel(\"head (m)\")\n",
    "ax.legend(loc=\"best\")\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-length / 2 - 25, 0.0, 201)\n",
    "hl = ml.headalongline(x, np.zeros(201), layers=[last_lay_dw, last_lay_dw + 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n",
    "\n",
    "ax.plot(x, hl[0].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n",
    "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n",
    "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n",
    "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n",
    "ax.set_xlabel(\"x (m)\")\n",
    "ax.set_ylabel(\"head (m)\")\n",
    "ax.legend(loc=\"best\")\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot cross-section around the sheetpile wall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xoffset = 50\n",
    "zoffset = 15\n",
    "x1, x2, y1, y2 = [-length / 2 - xoffset, 0.0, 0, 0]\n",
    "nudge = 1e-6\n",
    "n = 301"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot head contour cross-sections\n",
    "h = ml.headalongline(\n",
    "    np.linspace(x1 + nudge, x2 - nudge, n),\n",
    "    np.linspace(y1 + nudge, y2 - nudge, n),\n",
    ")\n",
    "L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)\n",
    "xg = np.linspace(0, L, n) + x1\n",
    "\n",
    "zg = 0.5 * (ml.aq.zaqbot + ml.aq.zaqtop)\n",
    "zg = np.hstack((ml.aq.zaqtop[0], zg, ml.aq.zaqbot[-1]))\n",
    "h = np.vstack((h[0], h, h[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels = np.linspace(h_bem - 0.1, -0.0, 51)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n",
    "ax.set_aspect(\"equal\")\n",
    "ml.plots.topview(win=[x1, x2, y1, y2], orientation=\"ver\", newfig=False)\n",
    "cf = ax.contourf(xg, zg, h, levels)\n",
    "cs = ax.contour(xg, zg, h, levels, colors=\"k\", linewidths=0.5)\n",
    "ax.set_ylim(z_dw - zoffset, z_dw + zoffset)\n",
    "ax.set_ylabel(\"depth (m NAP)\")\n",
    "ax.set_xlabel(\"m\")\n",
    "\n",
    "cb = plt.colorbar(cf, ax=ax, shrink=0.6)\n",
    "cb.set_label(\"head (m)\")\n",
    "cb.set_ticks(np.arange(-6, 0.1, 1))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model 2: Add more layers\n",
    "Add more layers to the model to get a more accurate solution of the flow towards the building pit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 11  # number of layers around bottom of sheetpile wall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate thickness of each sublayer above and below the sheetpile wall\n",
    "dz_i_top = (z_well - z_dw) / np.sum(np.arange(n + 1))\n",
    "dz_i_bot = (z_dw - z_extra) / np.sum(np.arange(2 * n + 1))\n",
    "\n",
    "# Generate cumulative depths for sublayers above and below the wall\n",
    "z_layers_top = np.cumsum(np.arange(n) * dz_i_top)\n",
    "z_layers_bot = np.cumsum(np.arange(2 * n) * dz_i_bot)\n",
    "\n",
    "# Combine sublayer depths into a single array for the model\n",
    "zgr = np.r_[z_dw + z_layers_top[::-1], (z_dw - z_layers_bot)[1:]]\n",
    "\n",
    "# Build full array of layer boundaries for the model\n",
    "z4 = np.r_[\n",
    "    np.array([ztop + 1, ztop, z_well, z_well]),\n",
    "    np.repeat(zgr, 2, 0),\n",
    "    z_extra,\n",
    "    z_extra,\n",
    "    zbot,\n",
    "]\n",
    "\n",
    "# Calculate thicknesses and hydraulic conductivities for each layer\n",
    "dz4 = z4[1:-1:2] - z4[2::2]\n",
    "kh_arr = kh * np.ones(dz4.shape)\n",
    "\n",
    "# Calculate resistance for each layer\n",
    "c = np.r_[np.array([ctop]), dz4[:-1] / (2 * kv) + dz4[1:] / (2 * kv)]\n",
    "\n",
    "# Set hydraulic conductivity of the top layer to a very low value\n",
    "kh_arr2 = kh_arr.copy()\n",
    "kh_arr2[0] = 1e-5"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build model, solve, and calculate total discharge and distance to the 5 cm drawdown contour."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tml.ModelMaq(kaq=kh_arr, z=z4, c=c, topboundary=\"semi\", hstar=0.0)\n",
    "\n",
    "layers = np.arange(np.sum(z_dw <= ml.aq.zaqbot))\n",
    "last_lay_dw = layers[-1]\n",
    "inhom = tml.BuildingPitMaq(\n",
    "    ml,\n",
    "    xy,\n",
    "    kaq=kh_arr2,\n",
    "    z=z4[1:],\n",
    "    topboundary=\"conf\",\n",
    "    c=c[1:],\n",
    "    order=4,\n",
    "    ndeg=3,\n",
    "    layers=layers,\n",
    ")\n",
    "\n",
    "wlayers = np.arange(np.sum(-14 <= ml.aq.zaqbot))\n",
    "wlayers = wlayers[1:]\n",
    "\n",
    "tml.River(\n",
    "    ml,\n",
    "    x1=-length / 2 + offset,\n",
    "    y1=width / 2 - offset,\n",
    "    x2=length / 2 - offset,\n",
    "    y2=width / 2 - offset,\n",
    "    hls=h_bem,\n",
    "    layers=wlayers,\n",
    ")\n",
    "tml.River(\n",
    "    ml,\n",
    "    x1=-length / 2 + offset,\n",
    "    y1=0,\n",
    "    x2=length / 2 - offset,\n",
    "    y2=0,\n",
    "    hls=h_bem,\n",
    "    layers=wlayers,\n",
    "    order=5,\n",
    ")\n",
    "tml.River(\n",
    "    ml,\n",
    "    x1=-length / 2 + offset,\n",
    "    y1=-width / 2 + offset,\n",
    "    x2=length / 2 - offset,\n",
    "    y2=-width / 2 + offset,\n",
    "    hls=h_bem,\n",
    "    layers=wlayers,\n",
    ")\n",
    "\n",
    "# ml.solve_mp(nproc=2)\n",
    "ml.solve()\n",
    "\n",
    "Qtot_ml = 0.0\n",
    "\n",
    "for e in ml.elementlist:\n",
    "    if e.name == \"River\":\n",
    "        Qtot_ml += e.discharge()\n",
    "\n",
    "print(\"\\nDischarge =\", np.round(Qtot_ml.sum(), 2), \"m^3/day\")\n",
    "\n",
    "y = np.linspace(-width / 2 - 25, width / 2 + 1100, 201)\n",
    "hl = ml.headalongline(np.zeros(201), y, layers=[0])\n",
    "y_5cm_ml = np.interp(\n",
    "    -0.05, ml.headalongline(np.zeros(201), y, layers=0).squeeze(), y, right=np.nan\n",
    ")\n",
    "print(\"Distance to 5 cm drawdown contour =\", np.round(y_5cm_ml, 2), \"m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_lay_dw = layers[-1]\n",
    "x = np.linspace(0.0, length / 2 + 1100, 201)\n",
    "hl = ml.headalongline(x, np.zeros(201), layers=[0, last_lay_dw, last_lay_dw + 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n",
    "\n",
    "ax.plot(x, hl[0].squeeze(), label=\"head layer 0\")\n",
    "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n",
    "ax.plot(x, hl[2].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n",
    "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n",
    "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n",
    "ax.set_xlabel(\"x (m)\")\n",
    "ax.set_ylabel(\"head (m)\")\n",
    "ax.legend(loc=\"best\")\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-length / 2 - 25, 0.0, 201)\n",
    "hl = ml.headalongline(x, np.zeros(201), layers=[0, last_lay_dw, last_lay_dw + 1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(1, 1, figsize=(10, 3))\n",
    "\n",
    "ax.plot(x, hl[0].squeeze(), label=\"head layer 0\")\n",
    "ax.plot(x, hl[1].squeeze(), label=\"head layer {}\".format(last_lay_dw))\n",
    "ax.plot(x, hl[2].squeeze(), label=\"head layer {}\".format(last_lay_dw + 1))\n",
    "ax.axhline(-0.05, color=\"r\", linestyle=\"dashed\", lw=0.75, label=\"-0.05 m\")\n",
    "ax.axhline(-0.5, color=\"k\", linestyle=\"dashed\", lw=0.75, label=\"-0.5 m\")\n",
    "ax.set_xlabel(\"x (m)\")\n",
    "ax.set_ylabel(\"head (m)\")\n",
    "ax.legend(loc=\"best\")\n",
    "ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot head contours in a vertical cross-section around the sheetpile wall."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xoffset = 50\n",
    "zoffset = 15\n",
    "x1, x2, y1, y2 = [-length / 2 - xoffset, 0.0, 0, 0]\n",
    "nudge = 1e-6\n",
    "n = 301"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot head contour cross-sections\n",
    "h = ml.headalongline(\n",
    "    np.linspace(x1 + nudge, x2 - nudge, n), np.linspace(y1 + nudge, y2 - nudge, n)\n",
    ")\n",
    "L = np.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)\n",
    "xg = np.linspace(0, L, n) + x1\n",
    "\n",
    "zg = 0.5 * (ml.aq.zaqbot + ml.aq.zaqtop)\n",
    "zg = np.hstack((ml.aq.zaqtop[0], zg, ml.aq.zaqbot[-1]))\n",
    "h = np.vstack((h[0], h, h[-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "levels = np.linspace(h_bem - 0.1, -0.0, 51)\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n",
    "ax.set_aspect(\"equal\")\n",
    "ml.plots.topview(win=[x1, x2, y1, y2], orientation=\"ver\", newfig=False)\n",
    "cf = ax.contourf(xg, zg, h, levels)\n",
    "cs = ax.contour(xg, zg, h, levels, colors=\"k\", linewidths=0.5)\n",
    "ax.set_ylim(z_dw - zoffset, z_dw + zoffset)\n",
    "ax.set_ylabel(\"depth (m NAP)\")\n",
    "ax.set_xlabel(\"m\")\n",
    "\n",
    "cb = plt.colorbar(cf, ax=ax, shrink=0.6)\n",
    "cb.ax.set_ylabel(\"head (m)\")\n",
    "cb.set_ticks(np.arange(-6, 0.1, 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note the difference between the computed discharge between the model with only a few layers and the model with very fine layers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Number of layers                  |  N=3  |  N=35 | unit\")\n",
    "print(\"-\" * 56)\n",
    "print(\n",
    "    f\"Discharge from building pit       | {Qtot.sum():>5.1f} \"\n",
    "    f\"| {Qtot_ml.sum():>5.1f} | m^3/d\"\n",
    ")\n",
    "print(f\"Distance to 5 cm drawdown contour | {y_5cm:>5.1f} | {y_5cm_ml:>5.1f} | m\")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}