{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Recharge in a circular area\n",
    "\n",
    "This notebook demonstrates how to model a circular area with constant recharge or extraction using the `CircAreaSink` element."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import timflow.steady as tfs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "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 = tfs.ModelMaq(kaq=5, z=[10, 0])\n",
    "ca = tfs.CircAreaSink(ml, xc=0, yc=0, R=100, N=0.001)\n",
    "ml.solve()\n",
    "x = np.linspace(-200, 200, 100)\n",
    "h = ml.headalongline(x, 0)\n",
    "plt.plot(x, h[0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qx = np.zeros_like(x)\n",
    "for i in range(len(x)):\n",
    "    qx[i], qy = ml.disvec(x[i], 1e-6)\n",
    "plt.plot(x, qx)\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=\"--\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Circular area-sink and well\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 = tfs.ModelMaq(kaq=5, z=[10, 0])\n",
    "ca = tfs.CircAreaSink(ml, xc=-200, yc=0, R=100, N=0.001)\n",
    "w = tfs.Well(ml, 200, 0, Qw=Q, rw=0.1)\n",
    "ml.solve()\n",
    "x = np.linspace(-400, 300, 100)\n",
    "h = ml.headalongline(x, 0)\n",
    "plt.plot(x, h[0]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 0.001\n",
    "R = 100\n",
    "Q = N * np.pi * R**2\n",
    "ml = tfs.ModelMaq(kaq=5, z=[10, 0])\n",
    "ca = tfs.CircAreaSink(ml, xc=-200, yc=0, R=100, N=0.001)\n",
    "w = tfs.Well(ml, 200, 0, Qw=Q, rw=0.1)\n",
    "ml.solve()\n",
    "ml.plots.contour([-300, 300, -200, 200], ngr=40)"
   ]
  },
  {
   "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 = tfs.ModelMaq(kaq=[5, 20], z=[20, 12, 10, 0], c=[1000])\n",
    "ca = tfs.CircAreaSink(ml, xc=0, yc=0, R=100, N=0.001)\n",
    "w = tfs.Well(ml, 0, 0, Qw=Q, rw=0.1, layers=1)\n",
    "rf = tfs.Constant(ml, 1000, 0, 0)\n",
    "ml.solve()\n",
    "x = np.linspace(-200, 200, 100)\n",
    "h = ml.headalongline(x, 0)\n",
    "plt.plot(x, h[0])\n",
    "plt.plot(x, h[1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-1000, 1000, 101)\n",
    "h = ml.headalongline(x, 0)\n",
    "plt.plot(x, h[0])\n",
    "plt.plot(x, h[1]);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}