{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inhomogeneities with ModelMaq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a two-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. Aquifer properties are given in Table 1 and the inhomogeneity data is given in table 2. There is a uniform gradient of 0.002 in Southeastern\n",
    "direction.\n",
    "\n",
    "**Table 1: Aquifer data**\n",
    "\n",
    "|Layer        | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) |\n",
    "|------------:|-----------|-----------|-------|------------|\n",
    "|Aquifer 0    |    10     |     0     |   20  |            |\n",
    "|Leaky Layer 1|           |    -10    |   0   |   4000     |\n",
    "|Aquifer 1    |    20     |    -30    |   10  |            |\n",
    "\n",
    "**Table 2: Inhomogeneity 1 data**\n",
    "\n",
    "\n",
    "|Layer        | $k$ (m/d) | $z_b$ (m) | $z_t$ | $c$ (days) |\n",
    "|:-----------:|----------:|----------:|------:|-----------:|\n",
    "|Aquifer 0    |    2      |     0     |   20  |    -       |\n",
    "|Leaky Layer 1|    -      |    -10    |   0   |   500      |\n",
    "|Aquifer 1    |    80     |    -30    |   -10 |   -        |\n",
    "\n",
    "\n",
    "A layout of the nodes of the inhomogeneity are shown in Fig. 1 (inhomogeneity 2 will be added\n",
    "later on). A well is located in the top aquifer inside inhomogeneity 1 (the black dot).\n",
    "\n",
    "**Figure 1: Layout of elements for exercise 3. A well is located inside inhomogeneity 1. Inhomogeneity 2 is added in the second part of the exercise.**\n",
    "\n",
    "![Layout of elements for exercise 3. A well is located inside inhomogeneity 1. Inhomogeneity 2 is added in the second part of the exercise.](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",
    "figsize = (6, 6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])\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.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy1,\n",
    "    kaq=[2, 80],\n",
    "    z=[20, 0, -10, -30],\n",
    "    c=[500],\n",
    "    topboundary=\"conf\",\n",
    "    order=3,\n",
    "    ndeg=2,\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",
    "\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],\n",
    "    ngr=50,\n",
    "    layers=[0, 1],\n",
    "    levels=50,\n",
    "    labels=1,\n",
    "    decimals=2,\n",
    "    figsize=figsize,\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",
    "    figsize=figsize,\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], figsize=figsize)\n",
    "ml.plots.plotcapzone(\n",
    "    well=w, hstepmax=50, nt=20, zstart=10, tmax=20 * 365.25, orientation=\"both\", ax=axes\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3d\n",
    "\n",
    "Change the elevation of the bottom of aquifer 1 from -30 to -20 inside inhomogeneity 1 (you need to recreate the model). Make a\n",
    "contour plot with a cross-section below it and start some  pathlines from $x=-200$, $y=700$. (note that the cross-section\n",
    "shows the elevation layers in the background aquifer, not the inhomogeneity). Note that the pathlines jump when they enter and\n",
    "exit inhomogeneity 1. This is caused by the jump in the base. It meets all continuity conditions and is an\n",
    "approximation of the smooth change in elevation that occurs over a distance of approximately one aquifer\n",
    "thickness from the boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])\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.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy1,\n",
    "    kaq=[2, 80],\n",
    "    z=[20, 0, -10, -40],\n",
    "    c=[500],\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": [
    "# TODO: ADD AQUIFER TOP/BOTTOM JUMP\n",
    "axes = ml.plots.topview_and_xsection(win=[-200, 800, 0, 800], figsize=figsize)\n",
    "ml.plots.tracelines(\n",
    "    -200 * np.ones(2),\n",
    "    700 * np.ones(2),\n",
    "    [-25, -15],\n",
    "    hstepmax=25,\n",
    "    orientation=\"both\",\n",
    "    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.ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])\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.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy1,\n",
    "    kaq=[2, 80],\n",
    "    z=[20, 0, -10, -30],\n",
    "    c=[500],\n",
    "    topboundary=\"conf\",\n",
    "    order=4,\n",
    "    ndeg=2,\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.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy2,\n",
    "    kaq=[2, 8],\n",
    "    z=[20, 0, -10, -30],\n",
    "    c=[50],\n",
    "    topboundary=\"conf\",\n",
    "    order=4,\n",
    "    ndeg=2,\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(\n",
    "    win=[-200, 1000, 0, 1200], ngr=50, layers=[0, 1], levels=20, figsize=figsize\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test for reverse order"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xy1[::-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ml = tfs.ModelMaq(kaq=[10, 20], z=[20, 0, -10, -30], c=[4000])\n",
    "p1 = tfs.PolygonInhomMaq(\n",
    "    ml,\n",
    "    xy=xy1[::-1],\n",
    "    kaq=[2, 80],\n",
    "    z=[20, 0, -10, -30],\n",
    "    c=[500],\n",
    "    topboundary=\"conf\",\n",
    "    order=3,\n",
    "    ndeg=2,\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(\n",
    "    win=[-200, 800, 0, 800],\n",
    "    ngr=50,\n",
    "    layers=[0, 1],\n",
    "    levels=50,\n",
    "    labels=1,\n",
    "    decimals=2,\n",
    "    figsize=figsize,\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}