{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cec081d4-f256-4820-a133-97302d7e7a7c",
   "metadata": {},
   "source": [
    "# Radial Collector Wells\n",
    "\n",
    "This notebook contains multi-layer solutions of examples in 'Multilayer Analytic Element Modeling of Radial Collector Wells'\n",
    "\n",
    "[Bakker, M., Kelson, V.A. and Luther, K.H., 2005. Multilayer analytic element modeling of radial collector wells. Groundwater, 43(6), pp.926-934.](https://doi.org/10.1111/j.1745-6584.2005.00116.x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b958765d-a91f-4bbc-ac44-a397767eefc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35500d1a-fa15-4e42-8e86-4e2d2a633c79",
   "metadata": {},
   "source": [
    "## Example 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4414a820-c3a3-4a7b-93f9-b4303b715c30",
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 150  # hydraulic conductivity, m/d\n",
    "z = [24, 16, 11, 7, 5, 4.05, 3.45, 3.15, 2.85, 2.55, 1.95, 1, 0]\n",
    "kzoverkh = 1\n",
    "layerw = 7\n",
    "Qw = 12_000  # discharge of well, m^3/d\n",
    "L = 60  # length of well, m\n",
    "rw = 0.15  # radius of well, m\n",
    "nls = 10  # number of (equally sized) line-sinks\n",
    "xr, yr, hr = 60, 0, 24  # coordinates and head at reference point\n",
    "#\n",
    "xy = np.array(\n",
    "    list(zip(np.linspace(-L / 2, L / 2, nls + 1), np.zeros(nls + 1), strict=False))\n",
    ")\n",
    "xyalt = np.zeros((nls, 4))  # specify begin and end of each line-sink\n",
    "xyalt[:, 0] = xy[:-1, 0]\n",
    "xyalt[:, 1] = xy[:-1, 1]\n",
    "xyalt[:, 2] = xy[1:, 0]\n",
    "xyalt[:, 3] = xy[1:, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4945dc81-3ae4-4126-b786-c8f53d7aa91e",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.Model3D(kaq=k, z=z, kzoverkh=kzoverkh)\n",
    "w = tfs.CollectorWell(model=ml, xy=xy, Qw=Qw, layers=layerw, rw=rw)\n",
    "rf = tfs.Constant(model=ml, xr=xr, yr=yr, hr=hr, layer=0)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50a5f564-d52d-409c-943b-363bc504fe43",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"xcp, ycp, head in control points\")\n",
    "for i in range(nls):\n",
    "    ls = w.lslist[i]\n",
    "    print(ls.xc[0], ls.yc[0], ml.head(ls.xc[0], ls.yc[0], layers=layerw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ae884e4-f655-4d3f-a290-1f9f2a51fbaa",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"head in top layer at center of well {ml.head(0, 0, layers=0):.2f} m\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "120c9567-5621-45a6-a4e3-6cc0e206b0a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(\n",
    "    win=[-60, 60, -60, 60],\n",
    "    ngr=50,\n",
    "    layers=7,\n",
    "    levels=np.arange(20, 24.1, 0.2),\n",
    "    decimals=1,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91379fa0-3285-4392-9cf9-415f10e8ba1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.Model3D(kaq=k, z=z, kzoverkh=kzoverkh)\n",
    "w = tfs.CollectorWell(model=ml, xy=xyalt, Qw=Qw, layers=layerw, rw=rw)\n",
    "rf = tfs.Constant(model=ml, xr=xr, yr=yr, hr=hr, layer=0)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "35d2b0f7-9be5-4404-ab01-7be14aeb38ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(\n",
    "    win=[-60, 60, -60, 60],\n",
    "    ngr=50,\n",
    "    layers=7,\n",
    "    levels=np.arange(20, 24.1, 0.2),\n",
    "    decimals=1,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3bff817e-2e55-4a9b-a146-19ea431d6d3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"xcp, ycp, head in control points\")\n",
    "for i in range(nls):\n",
    "    ls = w.lslist[i]\n",
    "    print(ls.xc[0], ls.yc[0], ml.head(ls.xc[0], ls.yc[0], layers=layerw))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "119eb69c-2fb8-42ca-bc93-bc70b61cb052",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"head in top layer at center of well {ml.head(0, 0, layers=0):.2f} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e23f307d-9247-4e42-a383-6044dc122132",
   "metadata": {},
   "source": [
    "## Example 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59aa197f-230c-495d-9199-fd10b90fcd8e",
   "metadata": {},
   "outputs": [],
   "source": [
    "rcaisson = 3  # radius caisson, m\n",
    "xr, yr, hr = 100, 0, 24  # coordinates and head at reference point, m\n",
    "Qw = 60_000  # discharge of well, m^3/d\n",
    "narms = 5\n",
    "\n",
    "ml = tfs.Model3D(kaq=k, z=z, kzoverkh=kzoverkh)\n",
    "w = tfs.RadialCollectorWell(\n",
    "    model=ml,\n",
    "    x=0,\n",
    "    y=0,\n",
    "    L=60,\n",
    "    rcaisson=rcaisson,\n",
    "    Qw=Qw,\n",
    "    narms=narms,\n",
    "    rw=rw,\n",
    "    nls=nls,\n",
    "    layers=7,\n",
    ")\n",
    "rf = tfs.Constant(model=ml, xr=xr, yr=yr, hr=hr, layer=0)\n",
    "ml.solve()\n",
    "ml.plots.contour(\n",
    "    win=[-100, 100, -100, 100],\n",
    "    ngr=101,\n",
    "    layers=7,\n",
    "    levels=np.arange(21, 30, 0.5),\n",
    "    decimals=1,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "618a5369-8eb6-4c3a-9035-14fb28e09805",
   "metadata": {},
   "outputs": [],
   "source": [
    "w.lslist[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aeb5ced7-a9b6-462c-8163-803de03bcab0",
   "metadata": {},
   "outputs": [],
   "source": [
    "rcaisson = 3  # radius caisson, m\n",
    "xr, yr, hr = 100, 0, 24  # coordinates and head at reference point, m\n",
    "Qw = 60_000  # discharge of well, m^3/d\n",
    "narms = 5\n",
    "\n",
    "ml = tfs.Model3D(kaq=k, z=z, kzoverkh=kzoverkh)\n",
    "w = tfs.RadialCollectorWell(\n",
    "    model=ml,\n",
    "    x=0,\n",
    "    y=0,\n",
    "    L=[40, 60, 40, 60, 40],\n",
    "    angle=30,\n",
    "    rcaisson=rcaisson,\n",
    "    Qw=Qw,\n",
    "    narms=narms,\n",
    "    rw=rw,\n",
    "    nls=nls,\n",
    "    layers=7,\n",
    ")\n",
    "rf = tfs.Constant(model=ml, xr=xr, yr=yr, hr=hr, layer=0)\n",
    "ml.solve()\n",
    "ml.plots.contour(\n",
    "    win=[-100, 100, -100, 100],\n",
    "    ngr=101,\n",
    "    layers=7,\n",
    "    levels=np.arange(21, 30, 0.5),\n",
    "    decimals=1,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc9b306e-e192-4721-8914-770972a2db9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "rcaisson = 3  # radius caisson, m\n",
    "xr, yr, hr = 100, 0, 24  # coordinates and head at reference point, m\n",
    "Qw = 60_000  # discharge of well, m^3/d\n",
    "narms = 5\n",
    "\n",
    "ml = tfs.Model3D(kaq=k, z=z, kzoverkh=kzoverkh)\n",
    "w = tfs.RadialCollectorWell(\n",
    "    model=ml,\n",
    "    x=0,\n",
    "    y=0,\n",
    "    L=60,\n",
    "    angle=[-40, -20, 0, 20, 40],\n",
    "    rcaisson=rcaisson,\n",
    "    Qw=Qw,\n",
    "    narms=narms,\n",
    "    rw=rw,\n",
    "    nls=nls,\n",
    "    layers=7,\n",
    ")\n",
    "rf = tfs.Constant(model=ml, xr=xr, yr=yr, hr=hr, layer=0)\n",
    "ml.solve()\n",
    "ml.plots.contour(\n",
    "    win=[-100, 100, -100, 100],\n",
    "    ngr=101,\n",
    "    layers=7,\n",
    "    levels=np.arange(21, 30, 0.5),\n",
    "    decimals=1,\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "679aabf2-a23e-49fc-92ee-67e55d4af857",
   "metadata": {},
   "outputs": [],
   "source": [
    "rcaisson = 3  # radius caisson, m\n",
    "xr, yr, hr = 100, 0, 24  # coordinates and head at reference point, m\n",
    "Qw = 60_000  # discharge of well, m^3/d\n",
    "narms = 5\n",
    "\n",
    "ml = tfs.Model3D(kaq=k, z=z, kzoverkh=kzoverkh)\n",
    "w = tfs.RadialCollectorWell(\n",
    "    model=ml,\n",
    "    x=0,\n",
    "    y=0,\n",
    "    L=60,\n",
    "    rcaisson=rcaisson,\n",
    "    Qw=Qw,\n",
    "    narms=narms,\n",
    "    rw=rw,\n",
    "    nls=nls,\n",
    "    layers=[7, 0, 7, 0, 7],  # just for illustration purposes\n",
    ")\n",
    "rf = tfs.Constant(model=ml, xr=xr, yr=yr, hr=hr, layer=0)\n",
    "ml.solve()\n",
    "\n",
    "# contour in layer 0\n",
    "ax = ml.plots.contour(\n",
    "    win=[-100, 100, -100, 100],\n",
    "    ngr=101,\n",
    "    layers=0,\n",
    "    levels=np.arange(21, 30, 0.5),\n",
    "    decimals=1,\n",
    ")\n",
    "# plot arms in layer 0 with dashed lines\n",
    "ax = w.plot(ax=ax, layer=7, ls=\"dashed\")\n",
    "\n",
    "# contour in layer 7\n",
    "ax = ml.plots.contour(\n",
    "    win=[-100, 100, -100, 100],\n",
    "    ngr=101,\n",
    "    layers=7,\n",
    "    levels=np.arange(21, 30, 0.5),\n",
    "    decimals=1,\n",
    "    color=\"C1\",\n",
    ")\n",
    "# plot arms in layer 0 with dashed lines\n",
    "w.plot(ax=ax, layer=0, ls=\"dashed\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "110347d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "headinside = []\n",
    "print(\"head in collector well:\", w.headinside())\n",
    "for i in range(w.nls):\n",
    "    ls = w.lslist[i]\n",
    "    headinside.append(ml.head(ls.xc[0], ls.yc[0], layers=ls.layers[0]))\n",
    "headinside = np.array(headinside)\n",
    "assert np.allclose(headinside[1:], headinside[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a92dd06-da89-4bf3-a1db-8e5bf59fa902",
   "metadata": {},
   "outputs": [],
   "source": [
    "Q_arms = w.discharge_per_arm()\n",
    "print(\"discharge of each arm: \", Q_arms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd1697b8-dee5-4314-b604-541c96ee67c9",
   "metadata": {},
   "source": [
    "Note that the discharge of arms 0, 2, 4, and the discharges of arms 1, 3 are not exactly equal, even though you may expect that from symmetry. This is the case because the control points are put a horizontal distance `rw` away along each arm, always on the left side when going from the point closest to the caisson to the end point."
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}