diff --git a/timflow/steady/aquifer.py b/timflow/steady/aquifer.py index 5f884ca1..6e468941 100644 --- a/timflow/steady/aquifer.py +++ b/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 diff --git a/timflow/steady/aquifer_parameters.py b/timflow/steady/aquifer_parameters.py index 6a2322ac..3db800b0 100644 --- a/timflow/steady/aquifer_parameters.py +++ b/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 ------- diff --git a/timflow/steady/constant.py b/timflow/steady/constant.py index dde0baba..4e081e91 100644 --- a/timflow/steady/constant.py +++ b/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. diff --git a/timflow/steady/inhomogeneity.py b/timflow/steady/inhomogeneity.py index 9f94158e..beb61008 100644 --- a/timflow/steady/inhomogeneity.py +++ b/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): diff --git a/timflow/steady/inhomogeneity1d.py b/timflow/steady/inhomogeneity1d.py index e15687dc..10d08a6c 100644 --- a/timflow/steady/inhomogeneity1d.py +++ b/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, diff --git a/timflow/steady/model.py b/timflow/steady/model.py index e8e02062..2d2528b3 100644 --- a/timflow/steady/model.py +++ b/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): diff --git a/timflow/transient/aquifer.py b/timflow/transient/aquifer.py index 3bb6fc9b..f7883948 100644 --- a/timflow/transient/aquifer.py +++ b/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"): diff --git a/timflow/transient/element.py b/timflow/transient/element.py index f9f33baa..78270367 100644 --- a/timflow/transient/element.py +++ b/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] diff --git a/timflow/transient/inhom1d.py b/timflow/transient/inhom1d.py index 3942a477..1c9cdf0f 100644 --- a/timflow/transient/inhom1d.py +++ b/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,