{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inhomogeneities with Model3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a multi-layer aquifer system that contains one inhomogeneity. Inside the inhomogeneity the transmissivity\n",
    "of the top aquifer is much lower and the transmissivity of the bottom aquifer is much higher than outside the\n",
    "inhomogeneity. A well is located in the top aquifer inside the inhomogeneity (the black dot).\n",
    "\n",
    "The layout of the inhomogeneity is shown in the figure below. Aquifer properties are specified in the code.\n",
    "\n",
    "**Figure 1: Layout of elements for this exercise. A well is located inside inhomogeneity 1.**\n",
    "\n",
    "![Layout of elements for this exercise. A well is located inside inhomogeneity 1.](figs/inhomogeneity_exercise3.png)"
   ]
  },
  {
   "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\n",
    "\n",
    "plt.rcParams[\"figure.figsize\"] = (6, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.Model3D(kaq=10, z=[20, 10, 0], kzoverkh=0.02)\n",
    "xy1 = [\n",
    "    (0, 600),\n",
    "    (-100, 400),\n",
    "    (-100, 200),\n",
    "    (100, 100),\n",
    "    (300, 100),\n",
    "    (500, 100),\n",
    "    (700, 300),\n",
    "    (700, 500),\n",
    "    (600, 700),\n",
    "    (400, 700),\n",
    "    (200, 600),\n",
    "]\n",
    "p1 = tfs.PolygonInhom3D(\n",
    "    ml,\n",
    "    xy=xy1,\n",
    "    kaq=2,\n",
    "    z=[20, 10, 0],\n",
    "    kzoverkh=0.002,\n",
    "    topboundary=\"conf\",\n",
    "    order=5,\n",
    "    ndeg=3,\n",
    ")\n",
    "rf = tfs.Constant(ml, xr=1000, yr=0, hr=40)\n",
    "uf = tfs.Uflow(ml, slope=0.002, angle=-45)\n",
    "w = tfs.Well(ml, xw=400, yw=400, Qw=500, rw=0.2, layers=0)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Questions\n",
    "### Exercise 3a\n",
    "What are the leakage factors of the background aquifer and the inhomogeneity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Leakage factor of the background aquifer is:\", ml.aq.lab)\n",
    "print(\"Leakage factor of the inhomogeneity is:\", p1.lab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3b\n",
    "Make a contour plot of both aquifers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(\n",
    "    win=[-200, 800, 0, 800], ngr=50, layers=[0, 1], levels=20, labels=1, decimals=2\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(\n",
    "    win=[-1200, 1800, -1000, 1800],\n",
    "    ngr=50,\n",
    "    layers=[0, 1],\n",
    "    levels=50,\n",
    "    labels=1,\n",
    "    decimals=2,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "htop = ml.headalongline(np.linspace(101, 499, 100), 100 + 0.001 * np.ones(100))\n",
    "hbot = ml.headalongline(np.linspace(101, 499, 100), 100 - 0.001 * np.ones(100))\n",
    "plt.figure()\n",
    "plt.plot(np.linspace(101, 499, 100), htop[0])\n",
    "plt.plot(np.linspace(101, 499, 100), hbot[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qtop = np.zeros(100)\n",
    "qbot = np.zeros(100)\n",
    "layer = 1\n",
    "x = np.linspace(101, 499, 100)\n",
    "for i in range(100):\n",
    "    qx, qy = ml.disvec(x[i], 100 + 0.001)\n",
    "    qtop[i] = qy[layer]\n",
    "    qx, qy = ml.disvec(x[i], 100 - 0.001)\n",
    "    qbot[i] = qy[layer]\n",
    "plt.figure()\n",
    "plt.plot(x, qtop)\n",
    "plt.plot(x, qbot)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3c\n",
    "Create a 20-year capture zone for the well, starting the pathlines halfway the top aquifer. First create a contour plot with a cross-section below it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "axes = ml.plots.topview_and_xsection(win=[-200, 800, 0, 800])\n",
    "ml.plots.plotcapzone(\n",
    "    well=w, hstepmax=50, nt=20, zstart=15, tmax=20 * 365.25, orientation=\"both\", ax=axes\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two inhomogeneities\n",
    "\n",
    "A second inhomogeneity is added, which shares part of its boundary with the first inhomogeneity. The aquifer properties for the inhomogeneity are provided in table 3. Inside this second inhomogeneity, the transmissivity of both the bottom aquifer and the resistance of the leaky layer are reduced. The input is now somewhat complicated. First the data of the two inhomgeneities is entered. Second, analytic elements are placed along the boundary of the inhomogeneity with `MakeInhomPolySide`. \n",
    "This routine places line elements along a string of points, but it requires that the\n",
    "aquifer data is the same on the left and right sides of the line. Hence, for this case we need to break the boundary up\n",
    "in three sections: One section with the background aquifer on one side and `inhom1` on the other, one section\n",
    "with the background aquifer on one side and `inhom2` on the other, and one section with `inhom1` on one side\n",
    "and `inhom2` on the other. The input file is a bit longer\n",
    "\n",
    "**Table 3: Inhomogeneity 2 data**\n",
    "\n",
    "|Layer        | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) |\n",
    "|------------:|----------:|----------:|------:|-----------:|\n",
    "|Aquifer 0    |    2      |     -20   |   0   |   -        |\n",
    "|Leaky Layer 1|    -      |    -40    |   -20 |   50       |\n",
    "|Aquifer 1    |    8      |    -80    |   -40 |   -        |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.Model3D(kaq=10, z=[20, 10, 0], kzoverkh=0.02)\n",
    "xy1 = [\n",
    "    (0, 600),\n",
    "    (-100, 400),\n",
    "    (-100, 200),\n",
    "    (100, 100),\n",
    "    (300, 100),\n",
    "    (500, 100),\n",
    "    (700, 300),\n",
    "    (700, 500),\n",
    "    (600, 700),\n",
    "    (400, 700),\n",
    "    (200, 600),\n",
    "]\n",
    "p1 = tfs.PolygonInhom3D(\n",
    "    ml,\n",
    "    xy=xy1,\n",
    "    kaq=2,\n",
    "    z=[20, 10, 0],\n",
    "    kzoverkh=0.002,\n",
    "    topboundary=\"conf\",\n",
    "    order=5,\n",
    "    ndeg=3,\n",
    ")\n",
    "xy2 = [\n",
    "    (0, 600),\n",
    "    (200, 600),\n",
    "    (400, 700),\n",
    "    (400, 900),\n",
    "    (200, 1100),\n",
    "    (0, 1000),\n",
    "    (-100, 800),\n",
    "]\n",
    "p2 = tfs.PolygonInhom3D(\n",
    "    ml,\n",
    "    xy=xy2,\n",
    "    kaq=20,\n",
    "    z=[20, 10, 0],\n",
    "    kzoverkh=0.05,\n",
    "    topboundary=\"conf\",\n",
    "    order=5,\n",
    "    ndeg=3,\n",
    ")\n",
    "rf = tfs.Constant(ml, xr=1000, yr=0, hr=40)\n",
    "uf = tfs.Uflow(ml, slope=0.002, angle=-45)\n",
    "w = tfs.Well(ml, xw=400, yw=400, Qw=500, rw=0.2, layers=0)\n",
    "ml.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml.plots.contour(win=[-200, 1000, 0, 1200], ngr=50, layers=[0, 1], levels=20);"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}