Top boundary reference and method

timflow/steady/aquifer.py

@@ -89,6 +89,8 @@ def __init__(self, model, kaq, c, z, npor, ltype, model3d=False):
         else:
             self.nporll = self.npor[self.ltype == "l"]
         self.model3d = model3d
+        # set reference to top boundary for background aquifers and inhoms
+        self.topbc = None  # top boundary condition element, if any
 
     def initialize(self):
         self.elementlist = []  # Elementlist of aquifer

timflow/steady/aquifer_parameters.py

@@ -20,7 +20,7 @@ def param_maq(kaq, z, c, npor, top):
     npor : float or array of floats
         Porosity of the aquifer(s).
     top : string
-        'conf' for confined aquifer on top, 'leak' for leaky layer on top.
+        'conf' for confined aquifer on top, 'semi' for semi-confined aquifer on top.
 
     Returns
     -------

timflow/steady/constant.py

@@ -176,8 +176,6 @@ def setparams(self, sol):
         self.parameters[:, 0] = sol
 
 
-# class ConstantStar(Element, PotentialEquation):
-# I don't think we need the equation
 class ConstantStar(Element):
     """Constant representing the particular solution inside a semi-confined aquifer.
 

timflow/steady/inhomogeneity.py

@@ -139,10 +139,12 @@ def create_elements(self):
             if self.N is not None:
                 a = AreaSinkInhom(self.model, self.N, self.xcenter, aq=aqin)
                 a.inhomelement = True
+                self.topbc = a
         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
+            self.topbc = c
 
 
 class PolygonInhomMaq(PolygonInhom):
@@ -543,6 +545,7 @@ 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
+            self.topbc = c
 
 
 class BuildingPitMaq(BuildingPit):
@@ -915,6 +918,7 @@ 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
+            self.topbc = c
 
 
 class LeakyBuildingPitMaq(LeakyBuildingPit):

timflow/steady/inhomogeneity1d.py

@@ -135,12 +135,15 @@ def create_elements(self):
                 assert aqin.ilap, (
                     "Error: infiltration can only be added if topboundary='conf'"
                 )
-                XsectionAreaSinkInhom(self.model, self.x1, self.x2, self.N, layer=0)
+                self.topbc = XsectionAreaSinkInhom(
+                    self.model, self.x1, self.x2, self.N, layer=0
+                )
 
         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
+            self.topbc = c
 
     def plot(
         self,

timflow/steady/model.py

@@ -80,9 +80,11 @@ class Model:
         array indicating for each layer whether it is
         'a' aquifer layer
         'l' leaky layer
+    hstar : float, optional
+        head above the top leaky layer, only used if top is semi-confined.
     """
 
-    def __init__(self, kaq, z, c, npor, ltype, model3d=False):
+    def __init__(self, kaq, z, c, npor, ltype, model3d=False, hstar=None):
         # All input variables are numpy arrays
         # That should be checked outside this function
         self.elementlist = []
@@ -91,8 +93,10 @@ def __init__(self, kaq, z, c, npor, ltype, model3d=False):
         self.name = "Model"
         self.model_type = "steady"  # Model type for plotting and other purposes
 
-        self.plots = PlotSteady(self)
+        if self.aq.ltype[0] == "l":
+            self.aq.topbc = ConstantStar(self, hstar, aq=self.aq)
 
+        self.plots = PlotSteady(self)
         self.initialized = False
 
     def initialize(self):
@@ -1036,10 +1040,8 @@ def __init__(self, kaq=1, z=None, c=None, npor=0.3, topboundary="conf", hstar=No
         if z is None:
             z = [1, 0]
         kaq, c, npor, ltype = param_maq(kaq, z, c, npor, topboundary)
-        super().__init__(kaq=kaq, z=z, c=c, npor=npor, ltype=ltype)
+        super().__init__(kaq=kaq, z=z, c=c, npor=npor, ltype=ltype, hstar=hstar)
         self.name = "ModelMaq"
-        if self.aq.ltype[0] == "l":
-            ConstantStar(self, hstar, aq=self.aq)
 
 
 class Model3D(Model):
@@ -1119,11 +1121,11 @@ def __init__(
         if topboundary == "semi":
             z = np.hstack((z[0] + topthick, z))
         model3d = True
-        super().__init__(kaq=kaq, z=z, c=c, npor=npor, ltype=ltype, model3d=model3d)
+        super().__init__(
+            kaq=kaq, z=z, c=c, npor=npor, ltype=ltype, model3d=model3d, hstar=hstar
+        )
         self.aq.kzoverkh = kzoverkh  # add kzoverkh to aquifer object
         self.name = "Model3D"
-        if self.aq.ltype[0] == "l":
-            ConstantStar(self, hstar, aq=self.aq)
 
 
 class ModelXsection(Model):

timflow/transient/aquifer.py

@@ -87,6 +87,8 @@ def __init__(
         # self.D = self.T / self.Saq
         self.area = 1e200  # Smaller than default of ml.aq so that inhom is found
         self.name = name
+        # set reference to top boundary element for background aquifers and inhoms
+        self.topbc = None  # top boundary condition element, if any
 
     def __repr__(self):
         if self.topboundary.startswith("con"):

timflow/transient/element.py

@@ -333,3 +333,31 @@ def run_after_solve(self):  # noqa: B027
     @abstractmethod
     def plot(self, ax=None):
         """Plot the element."""
+
+    def get_bc(self, t):
+        """Get the boundary condition at time(s) t.
+
+        Parameters
+        ----------
+        t : scalar or array
+            Time(s) at which to get the boundary condition.
+
+        Returns
+        -------
+        bc : scalar or array
+            Boundary condition(s) at time t.
+        """
+        t = np.atleast_1d(t)
+        bc = np.zeros_like(t, dtype=float)
+
+        if self.ntstart == 1:
+            bc[t >= self.tstart[0]] = self.bcin[0]
+        else:
+            for itime in range(self.ntstart):
+                if itime == self.ntstart - 1:
+                    mask = t >= self.tstart[itime]
+                else:
+                    mask = (t >= self.tstart[itime]) & (t < self.tstart[itime + 1])
+                bc[mask] = self.bcin[itime]
+
+        return bc if len(bc) > 1 else bc[0]

timflow/transient/inhom1d.py

@@ -205,12 +205,16 @@ def create_elements(self):
             assert self.topboundary == "con" or self.topboundary == "phr", Exception(
                 "Infiltration can only be applied to a confined aquifer."
             )
-            AreaSinkXsection(self.model, self.x1, self.x2, tsandN=self.tsandN)
+            self.topbc = AreaSinkXsection(
+                self.model, self.x1, self.x2, tsandN=self.tsandN
+            )
         if self.tsandhstar is not None:
             assert self.topboundary == "sem", Exception(
                 "hstar can only be implemented on top of a semi-confined aquifer."
             )
-            HstarXsection(self.model, self.x1, self.x2, tsandhstar=self.tsandhstar)
+            self.topbc = HstarXsection(
+                self.model, self.x1, self.x2, tsandhstar=self.tsandhstar
+            )
 
     def plot(
         self,