{
 "cells": [
  {
   "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\n",
    "\n",
    "plt.rcParams[\"font.size\"] = 8.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Circular area-sink\n",
    "Circular area-sink with radius 100 m, located at the origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0.001\n",
    "R = 100\n",
    "ml = tft.ModelMaq(kaq=5, z=[10, 0], Saq=2e-4, tmin=1e-3, tmax=1e4)\n",
    "ca = tft.CircAreaSink(ml, 0, 0, 100, tsandN=[(0, 0.001)])\n",
    "ml.solve()\n",
    "ml.plots.head_along_line(-200, 200, 0, 0, t=[0.1, 1, 10], figsize=(8, 3), sstart=-200);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 3))\n",
    "x = np.linspace(-200, 200, 200)\n",
    "qx = np.zeros_like(x)\n",
    "for t in [0.1, 1, 10]:\n",
    "    for i in range(len(x)):\n",
    "        qx[i] = ml.disvec(x[i], 1e-6, t)[0].item()\n",
    "    plt.plot(x, qx, label=\"time is \" + str(t))\n",
    "qxb = N * np.pi * R**2 / (2 * np.pi * R)\n",
    "plt.axhline(qxb, color=\"r\", ls=\"--\")\n",
    "plt.axhline(-qxb, color=\"r\", ls=\"--\")\n",
    "plt.xlabel(\"x (m)\")\n",
    "plt.ylabel(\"Qx (m^2/d)\")\n",
    "plt.legend(loc=\"best\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circular area-sink and well\n",
    "\n",
    "Discharge of well is the same as total infiltration rate of the circular area-sink.\n",
    "Well and center of area-sink area located at equal distances from $y$-axis, so that the head remains\n",
    "zero along the $y$-axis. Solution approaches steady-state solution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0.001\n",
    "R = 100\n",
    "Q = N * np.pi * R**2\n",
    "ml = tft.ModelMaq(kaq=5, z=[10, 0], Saq=2e-4, tmin=1e-3, tmax=1e4, M=10)\n",
    "ca = tft.CircAreaSink(ml, -200, 0, 100, tsandN=[(0, 0.001)])\n",
    "w = tft.Well(ml, 200, 0, rw=0.1, tsandQ=[(0, Q)])\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.head_along_line(\n",
    "    -400, 300, 0, 0, t=[0.1, 1, 10, 100, 1000], figsize=(8, 3), sstart=-400\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = np.logspace(-3, 4, 100)\n",
    "h = ml.head(-200, 0, t)\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.semilogx(t, h[0])\n",
    "plt.xlabel(\"time\")\n",
    "plt.ylabel(\"head\")\n",
    "plt.title(\"head at center of area-sink\")\n",
    "plt.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0.001\n",
    "R = 100\n",
    "Q = N * np.pi * R**2\n",
    "ml = tft.ModelMaq(kaq=5, z=[10, 0], Saq=2e-4, tmin=10, tmax=100, M=10)\n",
    "ca = tft.CircAreaSink(ml, -200, 0, 100, tsandN=[(0, 0.001)])\n",
    "w = tft.Well(ml, 200, 0, rw=0.1, tsandQ=[(0, Q)])\n",
    "ml.solve()\n",
    "ml.plots.contour([-300, 300, -200, 200], ngr=40, t=20, labels=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two layers\n",
    "Discharge of well is the same as total infiltration rate of the circular area-sink. Center of area-sink and well are at the origin. Circular area-sink in layer 0, well in layer 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0.001\n",
    "R = 100\n",
    "Q = N * np.pi * R**2\n",
    "ml = tft.ModelMaq(\n",
    "    kaq=[5, 20],\n",
    "    z=[20, 12, 10, 0],\n",
    "    c=[1000],\n",
    "    Saq=[2e-4, 1e-4],\n",
    "    tmin=1e-3,\n",
    "    tmax=1e4,\n",
    "    M=10,\n",
    ")\n",
    "ca = tft.CircAreaSink(ml, 0, 0, 100, tsandN=[(0, 0.001)])\n",
    "w = tft.Well(ml, 0, 0, rw=0.1, tsandQ=[(0, Q)], layers=1)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.head_along_line(\n",
    "    -200, 200, 0, 0, t=[0.1, 100], layers=[0, 1], sstart=-200, figsize=(8, 3)\n",
    ");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.head_along_line(\n",
    "    -500, 500, 0, 0, t=[0.1, 100, 1000], layers=[0, 1], sstart=-500, figsize=(8, 3)\n",
    ");"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}