Merge branch 'dev' into intnormflux_per_segment2

mbakker7 · Dec 5, 2023 · 5282197 · 5282197
2 parents 9ce61c8 + bd71006
commit 5282197
Show file tree

Hide file tree

Showing 7 changed files with 549 additions and 92 deletions.
diff --git a/docs/01howto/elements/inhoms.rst b/docs/01howto/elements/inhoms.rst
@@ -1,7 +1,9 @@
 Inhomogeneities
 ===============
 
-Polygonal inhomogeneities can be added in the following ways:
+The following polygonal inhomogeneities are available in TimML. Note that there are 
+inhomogeneities for both multi-aquifer style models (:class:`~timml.model.ModelMaq`)
+and single aquifer models (with multiple layers, :class:`~timml.model.Model3D`).
 
 1. :class:`~timml.inhomogeneity.PolygonInhomMaq`, which is an inhomogeneity consisting
    of a regular sequence of aquifer - leaky layer - aquifer - leaky layer, aquifer, etc.
@@ -12,10 +14,18 @@ Polygonal inhomogeneities can be added in the following ways:
    the resistance from the middle of one layer to the middle of the next layer. Vertical
    anisotropy can be specified. The system may be bounded on top by a leaky layer.
 
-3. :class:`~timml.inhomogeneity.BuildingPit`, which is an inhomogeneity similar to
+3. :class:`~timml.inhomogeneity.BuildingPitMaq`, which is an inhomogeneity similar to
    `PolygonInhomMaq` in which impermeable walls can be placed along the edges of the
    imhomogeneity in specified layers.
 
-4. :class:`~timml.inhomogeneity.LeakyBuildingPit`, which is an inhomogeneity like
-   `BuildingPit` but with leaky walls instead of impermeable walls. The resistance of the
+4. :class:`~timml.inhomogeneity.BuildingPit3D`, which is an inhomogeneity similar to
+   `PolygonInhom3D` in which impermeable walls can be placed along the edges of the
+   imhomogeneity in specified layers.
+
+5. :class:`~timml.inhomogeneity.LeakyBuildingPitMaq`, which is an inhomogeneity like
+   `BuildingPitMaq` but with leaky walls instead of impermeable walls. The resistance of the
+   leaky walls can be specified.
+
+6. :class:`~timml.inhomogeneity.LeakyBuildingPit3D`, which is an inhomogeneity like
+   `BuildingPit3D` but with leaky walls instead of impermeable walls. The resistance of the
    leaky walls can be specified.
diff --git a/notebooks/BuildingPit.ipynb b/notebooks/BuildingPit.ipynb
diff --git a/notebooks/circular_buildingpit.ipynb b/notebooks/circular_buildingpit.ipynb
@@ -18,11 +18,13 @@
    "outputs": [],
    "source": [
     "# import packages\n",
-    "import timml as tml\n",
+    "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
+    "import pandas as pd\n",
     "import shapely\n",
-    "import matplotlib.pyplot as plt\n",
-    "from scipy.special import i0, k0, i1, k1"
+    "from scipy.special import i0, i1, k0, k1\n",
+    "\n",
+    "import timml as tml"
    ]
   },
   {
@@ -229,7 +231,7 @@
    ],
    "source": [
     "ml = tml.ModelMaq(kaq=kh, z=z, c=ctop, topboundary=\"semi\", hstar=0.0)\n",
-    "bpit = tml.LeakyBuildingPit(\n",
+    "bpit = tml.LeakyBuildingPitMaq(\n",
     "    ml,\n",
     "    xy,\n",
     "    kaq=kh,\n",
@@ -504,21 +506,59 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 33,
+   "execution_count": 18,
    "metadata": {},
    "outputs": [
     {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Q inside boundary:\n",
-      "LeakyLineDoublet, LeakyBuildingPit, exact\n",
-      "[-76.69 -76.63 -76.46]\n",
-      "\n",
-      "Q outside boundary:\n",
-      "LeakyLineDoublet, LeakyBuildingPit, exact\n",
-      "[-76.69 -76.6  -76.46]\n"
-     ]
+     "data": {
+      "text/html": [
+       "<style type=\"text/css\">\n",
+       "</style>\n",
+       "<table id=\"T_7a010\">\n",
+       "  <caption>Discharge along inside and outside of building pit:</caption>\n",
+       "  <thead>\n",
+       "    <tr>\n",
+       "      <th class=\"index_name level0\" >Model</th>\n",
+       "      <th id=\"T_7a010_level0_col0\" class=\"col_heading level0 col0\" >LeakyLineDoublet</th>\n",
+       "      <th id=\"T_7a010_level0_col1\" class=\"col_heading level0 col1\" >LeakyBuildingPit</th>\n",
+       "      <th id=\"T_7a010_level0_col2\" class=\"col_heading level0 col2\" >Exact</th>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th class=\"index_name level0\" >Discharge</th>\n",
+       "      <th class=\"blank col0\" >&nbsp;</th>\n",
+       "      <th class=\"blank col1\" >&nbsp;</th>\n",
+       "      <th class=\"blank col2\" >&nbsp;</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th id=\"T_7a010_level0_row0\" class=\"row_heading level0 row0\" >inside</th>\n",
+       "      <td id=\"T_7a010_row0_col0\" class=\"data row0 col0\" >-76.69</td>\n",
+       "      <td id=\"T_7a010_row0_col1\" class=\"data row0 col1\" >-76.63</td>\n",
+       "      <td id=\"T_7a010_row0_col2\" class=\"data row0 col2\" >-76.46</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th id=\"T_7a010_level0_row1\" class=\"row_heading level0 row1\" >outside</th>\n",
+       "      <td id=\"T_7a010_row1_col0\" class=\"data row1 col0\" >-76.69</td>\n",
+       "      <td id=\"T_7a010_row1_col1\" class=\"data row1 col1\" >-76.60</td>\n",
+       "      <td id=\"T_7a010_row1_col2\" class=\"data row1 col2\" >-76.46</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th id=\"T_7a010_level0_row2\" class=\"row_heading level0 row2\" >difference</th>\n",
+       "      <td id=\"T_7a010_row2_col0\" class=\"data row2 col0\" >-0.00</td>\n",
+       "      <td id=\"T_7a010_row2_col1\" class=\"data row2 col1\" >-0.03</td>\n",
+       "      <td id=\"T_7a010_row2_col2\" class=\"data row2 col2\" >-0.00</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n"
+      ],
+      "text/plain": [
+       "<pandas.io.formats.style.Styler at 0x7f6408b27250>"
+      ]
+     },
+     "execution_count": 18,
+     "metadata": {},
+     "output_type": "execute_result"
     }
    ],
    "source": [
@@ -545,11 +585,26 @@
     "Qexact_out = disr(100 + nudge) * 2 * np.pi * (100.0 + nudge)\n",
     "\n",
     "# print results\n",
-    "print(\"Q inside boundary:\\nLeakyLineDoublet, LeakyBuildingPit, exact\")\n",
-    "print(np.round(np.array([Qtot_in_lld, Qtot_in_lp, Qexact_in]), 2))\n",
-    "print(\"\\nQ outside boundary:\\nLeakyLineDoublet, LeakyBuildingPit, exact\")\n",
-    "print(np.round(np.array([Qtot_out_lld, Qtot_out_lp, Qexact_out]), 2))"
+    "df = pd.DataFrame(\n",
+    "    index=[\"inside\", \"outside\"],\n",
+    "    columns=[\"LeakyLineDoublet\", \"LeakyBuildingPit\", \"Exact\"],\n",
+    ")\n",
+    "df.loc[\"inside\"] = Qtot_in_lld, Qtot_in_lp, Qexact_in\n",
+    "df.loc[\"outside\"] = Qtot_out_lld, Qtot_out_lp, Qexact_out\n",
+    "df.loc[\"difference\"] = df.loc[\"inside\"] - df.loc[\"outside\"]\n",
+    "df.index.name = \"Discharge\"\n",
+    "df.columns.name = \"Model\"\n",
+    "df.style.format(precision=2).set_caption(\n",
+    "    \"Discharge along inside and outside of building pit:\"\n",
+    ")"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -568,7 +623,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.9"
+   "version": "3.10.12"
   }
  },
  "nbformat": 4,

diff --git a/timml/__init__.py b/timml/__init__.py
@@ -18,8 +18,10 @@
 from .circareasink import CircAreaSink
 from .constant import Constant, ConstantStar
 from .inhomogeneity import (
-    BuildingPit,
-    LeakyBuildingPit,
+    BuildingPit3D,
+    BuildingPitMaq,
+    LeakyBuildingPit3D,
+    LeakyBuildingPitMaq,
     PolygonInhom3D,
     PolygonInhomMaq,
 )

diff --git a/timml/equation.py b/timml/equation.py
@@ -523,9 +523,9 @@ def equation(self):
                         / self.aqout.Tcol[self.layers]
                     )
 
-                    mat[
-                        istart : istart + self.nlayers, ieq : ieq + e.nunknowns
-                    ] = flux - self.resfac * (headin - headout)
+                    mat[istart : istart + self.nlayers, ieq : ieq + e.nunknowns] = (
+                        flux - self.resfac * (headin - headout)
+                    )
                     ieq += e.nunknowns
                 else:
                     flux = self.intflux(
@@ -562,7 +562,7 @@ def equation(self):
                         / self.aqout.T[self.layers]
                     )
 
-                    rhs[istart : istart + self.nlayers] += -flux + self.resfac * (
-                        headin - headout
+                    rhs[istart : istart + self.nlayers] += (
+                        -flux + self.resfac.squeeze() * (headin - headout)
                     )
         return mat, rhs