Add LeakyBuildingPit element (#70)

* Add LeakyBuildingPit element
- add LeakyIntHeadDiffEquation mixin class
- add LeakyIntHeadDiffLineSink element for leaky pit boundaries with leaky walls
- add LeakyBuildingPit element, inherits from BuildingPit but adds resistance term
- add notebook testing LeakyBuildingPit

allow BuildingPit to be placed in models with 1 layer
isort __init__

* use self.disveclayers for RHS

* fix #71
- resfac fixed in LineDoublet(1d)
- resfac fixed in LineSink(1d)

* update leaky building pit nb

* allow passing kwargs to contour func

* add integration of normal flux along line
- add disvecnorm: calculate disvec normal to angle theta
- intdisvecnorm: quad/legendre integration of disvecnorm
- add quad integrand function for integrating from 0->r

* rename *disvecnorm to *normflux
- improve code and docstring

* fix overwrite of label kwarg

* add intnormflux notebook

* add intnormflux notebook

* add quad_vec import

* improve normal flux nb

* optimize code a bit and reorder coordinates

* update normalflux nb

* docstrings normflux functions

* add circular leaky building pit notebook

* rename to IntLeakyWallEquation

* rename to IntLeakyWallEquation

* delete nb

* remove py3.7 testing, add 3.10 and 3.11

* fix buildingpit nb

* add shapely to ci reqs

Author: dbrakenhoff

Diff for: .github/workflows/ci.yml

@@ -14,7 +14,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.7, 3.8, 3.9]
+        python-version: [3.8, 3.9, '3.10', '3.11']
     steps:
       - uses: actions/checkout@v3
 

Diff for: notebooks/BuildingPit.ipynb

@@ -3,4 +3,5 @@ scipy>=0.19
 numba>=0.39
 matplotlib>=2.0,<3.5.0
 pytest>=4.6
-jupyter>=1.0.0
+jupyter>=1.0.0
+shapely

Diff for: notebooks/circular_buildingpit.ipynb

@@ -12,36 +12,41 @@
 # --version number
 __name__ = "timml"
 __author__ = "Mark Bakker"
-from .version import __version__
+from . import bessel
 
 # Import all classes and functions
-from .model import ModelMaq, Model3D, Model
-from .well import WellBase, Well, HeadWell
+from .circareasink import CircAreaSink
 from .constant import Constant, ConstantStar
-from .linesink import (
-    LineSinkBase,
-    HeadLineSinkZero,
-    HeadLineSink,
-    LineSinkDitch,
-    HeadLineSinkString,
-    LineSinkDitchString,
-    HeadLineSinkContainer,
+from .inhomogeneity import (
+    BuildingPit,
+    LeakyBuildingPit,
+    PolygonInhom3D,
+    PolygonInhomMaq,
 )
+from .inhomogeneity1d import StripInhom3D, StripInhomMaq
 from .linedoublet import (
     ImpLineDoublet,
     ImpLineDoubletString,
     LeakyLineDoublet,
     LeakyLineDoubletString,
 )
-from .circareasink import CircAreaSink
-from .inhomogeneity import PolygonInhomMaq, PolygonInhom3D, BuildingPit
-from .inhomogeneity1d import StripInhomMaq, StripInhom3D
-from .uflow import Uflow
-from .trace import timtraceline, timtracelines
-from .linesink1d import LineSink1D, HeadLineSink1D
 from .linedoublet1d import ImpLineDoublet1D, LeakyLineDoublet1D
+from .linesink import (
+    HeadLineSink,
+    HeadLineSinkContainer,
+    HeadLineSinkString,
+    HeadLineSinkZero,
+    LineSinkBase,
+    LineSinkDitch,
+    LineSinkDitchString,
+)
+from .linesink1d import HeadLineSink1D, LineSink1D
+from .model import Model, Model3D, ModelMaq
 from .stripareasink import StripAreaSink
-from . import bessel
+from .trace import timtraceline, timtracelines
+from .uflow import Uflow
+from .version import __version__
+from .well import HeadWell, Well, WellBase
 
 __all__ = [s for s in dir() if not s.startswith("_")]
 

Diff for: notebooks/normal_flux.ipynb

@@ -321,3 +321,91 @@ def equation(self):
                     rhs[istart:istart + self.nlayers] -= flux
 
         return mat, rhs
+
+
+class IntLeakyWallEquation:
+    def equation(self):
+        mat = np.empty((self.nunknowns, self.model.neq))
+        rhs = np.zeros(self.nunknowns)  # Needs to be initialized to zero
+        for icp in range(self.ncp):
+            istart = icp * self.nlayers
+            ieq = 0
+            for e in self.model.elementlist:
+                if e.nunknowns > 0:
+                    flux = self.intflux(
+                        e.disvecinflayers, 
+                        self.xc[icp],
+                        self.yc[icp],
+                        self.xc[icp + 1],
+                        self.yc[icp + 1],
+                        self.layers,
+                        aq=self.aq,
+                    )
+                    headin = (
+                        self.intpot(
+                            e.potinflayers,
+                            self.xcin[icp],
+                            self.ycin[icp],
+                            self.xcin[icp + 1],
+                            self.ycin[icp + 1],
+                            self.layers,
+                            aq=self.aqin,
+                        )
+                        / self.aqin.Tcol[self.layers]
+                    )
+                    headout = (
+                        self.intpot(
+                            e.potinflayers,
+                            self.xcout[icp],
+                            self.ycout[icp],
+                            self.xcout[icp + 1],
+                            self.ycout[icp + 1],
+                            self.layers,
+                            aq=self.aqout,
+                        )
+                        / self.aqout.Tcol[self.layers]
+                    )
+
+                    mat[
+                        istart : istart + self.nlayers, ieq : ieq + e.nunknowns
+                    ] = flux - self.resfac * (headin - headout)
+                    ieq += e.nunknowns
+                else:
+                    flux = self.intflux(
+                        e.disveclayers, 
+                        self.xc[icp],
+                        self.yc[icp],
+                        self.xc[icp + 1],
+                        self.yc[icp + 1],
+                        self.layers,
+                        aq=self.aq,
+                    )
+                    headin = (
+                        self.intpot(
+                            e.potentiallayers,
+                            self.xcin[icp],
+                            self.ycin[icp],
+                            self.xcin[icp + 1],
+                            self.ycin[icp + 1],
+                            self.layers,
+                            aq=self.aqin,
+                        )
+                        / self.aqin.T[self.layers]
+                    )
+                    headout = (
+                        self.intpot(
+                            e.potentiallayers,
+                            self.xcout[icp],
+                            self.ycout[icp],
+                            self.xcout[icp + 1],
+                            self.ycout[icp + 1],
+                            self.layers,
+                            aq=self.aqout,
+                        )
+                        / self.aqout.T[self.layers]
+                    )
+
+                    rhs[istart : istart + self.nlayers] += (-flux +
+                        self.resfac * (headin - headout) 
+                    )
+        return mat, rhs

Diff for: requirements.ci.txt

@@ -4,8 +4,12 @@
 from .aquifer_parameters import param_maq, param_3d
 from .element import Element
 from .constant import ConstantInside, ConstantStar
-from .intlinesink import IntHeadDiffLineSink, IntFluxDiffLineSink, IntFluxLineSink
-from .controlpoints import controlpoints
+from .intlinesink import (
+    IntHeadDiffLineSink,
+    IntFluxDiffLineSink,
+    IntFluxLineSink,
+    LeakyIntHeadDiffLineSink,
+)
 
 __all__ = ['PolygonInhomMaq']
 
@@ -147,7 +151,8 @@ def __init__(self, model, xy, kaq=1, z=[1, 0], c=[], npor=0.3,
         kaq, c, npor, ltype, = param_maq(kaq, z, c, npor, topboundary)
         PolygonInhom.__init__(self, model, xy, kaq, c, z, npor, ltype,
                               hstar, N, order, ndeg)
-        
+
+
 class PolygonInhom3D(PolygonInhom):
     """
     Model3D Class to create a multi-layer model object consisting of
@@ -332,7 +337,8 @@ def isinside(self, x, y):
                 return rv
             angles = np.log(bigZmin1 / bigZplus1).imag
             angle = np.sum(angles)
-            if angle > np.pi: rv = 1
+            if angle > np.pi: 
+                rv = 1
         return rv
 
     def create_elements(self):
@@ -357,20 +363,87 @@ def create_elements(self):
                                 order=self.order, ndeg=self.ndeg,
                                 label=None, addtomodel=True,
                                 aq=aqout, aqin=aqin, aqout=aqout)
+                if len(self.nonimplayers) > 0:
+                    # use these conditions for layers without impermeable or leaky walls
+                    IntHeadDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
+                                        x2=self.x[i + 1], y2=self.y[i + 1],
+                                        layers=self.nonimplayers,
+                                        order=self.order, ndeg=self.ndeg,
+                                        label=None, addtomodel=True,
+                                        aq=aqin, aqin=aqin, aqout=aqout)
+                    IntFluxDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
+                                        x2=self.x[i + 1], y2=self.y[i + 1],
+                                        layers=self.nonimplayers,
+                                        order=self.order, ndeg=self.ndeg,
+                                        label=None, addtomodel=True,
+                                        aq=aqout, aqin=aqin, aqout=aqout)
+
+        if aqin.ltype[0] == 'a':  # add constant on inside
+            c = ConstantInside(self.model, self.zcin.real, self.zcin.imag)
+            c.inhomelement = True
+        if aqin.ltype[0] == 'l':
+            assert self.hstar is not None, 'Error: hstar needs to be set'
+            c = ConstantStar(self.model, self.hstar, aq=aqin)
+            c.inhomelement = True
 
-                # use these conditions for layers without impermeable walls
-                IntHeadDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
-                                    x2=self.x[i + 1], y2=self.y[i + 1],
-                                    layers=self.nonimplayers,
-                                    order=self.order, ndeg=self.ndeg,
-                                    label=None, addtomodel=True,
-                                    aq=aqin, aqin=aqin, aqout=aqout)
-                IntFluxDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
-                                    x2=self.x[i + 1], y2=self.y[i + 1],
-                                    layers=self.nonimplayers,
-                                    order=self.order, ndeg=self.ndeg,
-                                    label=None, addtomodel=True,
-                                    aq=aqout, aqin=aqin, aqout=aqout)
+
+class LeakyBuildingPit(BuildingPit):
+    def __init__(self, model, xy, kaq=1, z=[1, 0], c=[], npor=0.3,
+                 topboundary="conf", hstar=None, order=3,
+                 ndeg=3, layers=[0], res=np.inf):
+        super().__init__(
+            model=model,
+            xy=xy,
+            kaq=kaq,
+            z=z,
+            c=c,
+            npor=npor,
+            topboundary=topboundary,
+            hstar=hstar,
+            order=order,
+            ndeg=ndeg,
+            layers=layers,
+        )
+        self.res = res
+
+    def create_elements(self):
+        aqin = self.model.aq.find_aquifer_data(self.zcin[0].real,
+                                               self.zcin[0].imag)
+        for i in range(self.Nsides):
+            aqout = self.model.aq.find_aquifer_data(self.zcout[i].real,
+                                                    self.zcout[i].imag)
+            if (aqout == self.model.aq) or (
+                aqout.inhom_number > self.inhom_number):
+
+                # Conditions for layers with impermeable walls
+                LeakyIntHeadDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
+                                         x2=self.x[i + 1], y2=self.y[i + 1],
+                                         res=self.res, layers=self.layers,
+                                         order=self.order, ndeg=self.ndeg,
+                                         label=None, addtomodel=True,
+                                         aq=aqin, aqin=aqin, aqout=aqout)
+                
+                LeakyIntHeadDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
+                                         x2=self.x[i + 1], y2=self.y[i + 1],
+                                         res=self.res, layers=self.layers,
+                                         order=self.order, ndeg=self.ndeg,
+                                         label=None, addtomodel=True,
+                                         aq=aqout, aqin=aqin, aqout=aqout)
+                
+                if len(self.nonimplayers) > 0:
+                    # use these conditions for layers without impermeable or leaky walls
+                    IntHeadDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
+                                        x2=self.x[i + 1], y2=self.y[i + 1],
+                                        layers=self.nonimplayers,
+                                        order=self.order, ndeg=self.ndeg,
+                                        label=None, addtomodel=True,
+                                        aq=aqin, aqin=aqin, aqout=aqout)
+                    IntFluxDiffLineSink(self.model, x1=self.x[i], y1=self.y[i],
+                                        x2=self.x[i + 1], y2=self.y[i + 1],
+                                        layers=self.nonimplayers,
+                                        order=self.order, ndeg=self.ndeg,
+                                        label=None, addtomodel=True,
+                                        aq=aqout, aqin=aqin, aqout=aqout)
 
         if aqin.ltype[0] == 'a':  # add constant on inside
             c = ConstantInside(self.model, self.zcin.real, self.zcin.imag)
@@ -379,7 +452,8 @@ def create_elements(self):
             assert self.hstar is not None, 'Error: hstar needs to be set'
             c = ConstantStar(self.model, self.hstar, aq=aqin)
             c.inhomelement = True
-            
+        
+
 class AreaSinkInhom(Element):
     def __init__(self, model, N, xc, \
                  name='AreaSinkInhom', label=None, aq=None):

Diff for: timml/__init__.py

@@ -1,6 +1,11 @@
 import numpy as np
 from .linesink import LineSinkHoBase
-from .equation import HeadDiffEquation2, DisvecDiffEquation2, IntDisVecEquation
+from .equation import (
+    HeadDiffEquation2, 
+    DisvecDiffEquation2, 
+    IntDisVecEquation, 
+    IntLeakyWallEquation,
+)
 from .controlpoints import controlpoints
 # needed for testing
 # from .equation import DisvecDiffEquation
@@ -178,3 +183,58 @@ def initialize(self):
 
     def setparams(self, sol):
         self.parameters[:, 0] = sol
+
+
+class LeakyIntHeadDiffLineSink(LineSinkHoBase, IntLeakyWallEquation):
+    """Element to set numerically integrated head
+    along linesink to equal to:
+        Qnormal = H * (headin - headout) / res
+    Used in LeakyBuildingPit element
+    """
+
+    def __init__(self, model, x1=-1, y1=0, x2=1, y2=0, res=np.inf, layers=None,
+                 order=0, ndeg=3, label=None, addtomodel=True, aq=None,
+                 aqin=None, aqout=None):
+        if layers is None:
+            layers = list(range(model.aq.naq))
+        LineSinkHoBase.__init__(self, model, x1, y1, x2, y2, Qls=0,
+                                layers=layers, order=order,
+                                name='LeakyIntHeadDiffLineSink', label=label,
+                                addtomodel=addtomodel, aq=aq)
+        self.res = res
+        self.inhomelement = True
+        self.ndeg = ndeg
+        self.Xleg, self.wleg = np.polynomial.legendre.leggauss(self.ndeg)
+        self.nunknowns = self.nparam
+        self.aqin = aqin
+        self.aqout = aqout
+
+    def initialize(self):
+        LineSinkHoBase.initialize(self)
+
+        # recalculated with ncp - 1 points instead of ncp points
+        self.xcin, self.ycin = controlpoints(self.ncp - 1, self.z1, self.z2,
+                                             eps=1e-6, include_ends=True)
+        self.xcout, self.ycout = controlpoints(self.ncp - 1, self.z1, self.z2,
+                                               eps=-1e-6, include_ends=True)
+
+        if self.aqin is None:
+            self.aqin = self.model.aq.find_aquifer_data(self.xcin[0],
+                                                        self.ycin[0])
+        if self.aqout is None:
+            self.aqout = self.model.aq.find_aquifer_data(self.xcout[0],
+                                                         self.ycout[0])
+        # set control points x,y depending on which side of
+        # the boundary the element is for
+        if self.aq == self.aqin:
+            self.xc = self.xcin
+            self.yc = self.ycin
+        if self.aq == self.aqout:
+            self.xc = self.xcout
+            self.yc = self.ycout
+
+        # set resistance factor
+        self.resfac = self.aq.Haq[self.layers] / self.res
+
+    def setparams(self, sol):
+        self.parameters[:, 0] = sol

Diff for: timml/equation.py

@@ -56,7 +56,7 @@ def initialize(self):
                                                    eps=-1e-6)
         if self.aq is None:
             self.aq = self.model.aq.find_aquifer_data(self.xc[0], self.yc[0])
-        self.resfac = self.aq.T[self.layers] / self.res
+        self.resfac = self.aq.Haq[self.layers] / self.res
         if self.addtomodel:
             self.aq.add_element(self)
         self.parameters = np.empty((self.nparam, 1))