{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test line-sinks with image wells"
   ]
  },
  {
   "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": [
    "Consider a well pumping in a phreatic aquifer with $S_y=0.1$. The hydraulic conductivity of the aquifer is 10 m/d and the saturated thickness may be approximated as constant and equal to 20 m. The well is located at $(x,y)=(0,0)$. The discharge of the well is 1000 m$^3$/d and the radius is 0.3 m. there is a very long river with a fixed head located along the line $x=50$ m. The head is computed at $(x,y)=(20,0)$ for the first 20 days after the well starts pumping. The solution for an image well is compared to the solution using a `RiverString` element of different lengths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tft.ModelMaq(kaq=10, z=[20, 0], Saq=[0.1], phreatictop=True, tmin=0.001, tmax=100)\n",
    "w1 = tft.Well(ml1, 0, 0, rw=0.3, tsandQ=[(0, 1000)])\n",
    "w2 = tft.Well(ml1, 100, 0, rw=0.3, tsandQ=[(0, -1000)])\n",
    "ml1.solve()\n",
    "t = np.linspace(0.1, 20, 100)\n",
    "h1 = ml1.head(20, 0, t)\n",
    "plt.plot(t, h1[0], label=\"river modeled with image well\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, h1[0], label=\"river modeled with image well\")\n",
    "for ystart in [-50, -100]:\n",
    "    ml2 = tft.ModelMaq(\n",
    "        kaq=10, z=[20, 0], Saq=[0.1], phreatictop=True, tmin=0.001, tmax=100\n",
    "    )\n",
    "    w = tft.Well(ml2, 0, 0, rw=0.3, tsandQ=[(0, 1000)])\n",
    "    yls = np.arange(ystart, -ystart + 1, 20)\n",
    "    xls = 50 * np.ones(len(yls))\n",
    "    lss = tft.RiverString(ml2, xy=list(zip(xls, yls, strict=False)), tsandh=\"fixed\")\n",
    "    ml2.solve()\n",
    "    h2 = ml2.head(20, 0, t)\n",
    "    plt.plot(t, h2[0], \"--\", label=f\"line-sink string from {ystart} to {-ystart}\")\n",
    "plt.title(\"head at (x,y)=(20,0)\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution is repeated for the case where there is a long impermeable wall along $x=50$ m rather than a river."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml1 = tft.ModelMaq(kaq=10, z=[20, 0], Saq=[0.1], phreatictop=True, tmin=0.001, tmax=100)\n",
    "w1 = tft.Well(ml1, 0, 0, rw=0.3, tsandQ=[(0, 1000)])\n",
    "w2 = tft.Well(ml1, 100, 0, rw=0.3, tsandQ=[(0, 1000)])\n",
    "ml1.solve()\n",
    "t = np.linspace(0.1, 20, 100)\n",
    "h1 = ml1.head(20, 0, t)\n",
    "plt.plot(t, h1[0], label=\"impermeable wall modeled with image well\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(t, h1[0], label=\"river modeled with image well\")\n",
    "for ystart in [-100, -200, -400]:\n",
    "    ml2 = tft.ModelMaq(\n",
    "        kaq=10, z=[20, 0], Saq=[0.1], phreatictop=True, tmin=0.001, tmax=100\n",
    "    )\n",
    "    w = tft.Well(ml2, 0, 0, rw=0.3, tsandQ=[(0, 1000)])\n",
    "    yls = np.arange(ystart, -ystart + 1, 20)\n",
    "    xls = 50 * np.ones(len(yls))\n",
    "    lss = tft.LeakyWallString(ml2, xy=list(zip(xls, yls, strict=False)), res=\"imp\")\n",
    "    ml2.solve()\n",
    "    h2 = ml2.head(20, 0, t)\n",
    "    plt.plot(t, h2[0], \"--\", label=f\"line-doublet string from {ystart} to {-ystart}\")\n",
    "plt.title(\"head at (x,y)=(20,0)\")\n",
    "plt.xlabel(\"time (d)\")\n",
    "plt.ylabel(\"head (m)\")\n",
    "plt.legend();"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}