Skip to content

fix #64 #90

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
dbrakenhoff wants to merge 1 commit into
base: dev
Choose a base branch
from
Open

fix #64 #90

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 19 additions & 18 deletions timflow/steady/aquifer.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@


class AquiferData:
def __init__(self, model, kaq, c, z, npor, ltype):
def __init__(self, model, kaq, c, z, npor, ltype, model3d=False):
"""Initialize aquifer data.

Parameters
Expand All @@ -37,6 +37,8 @@ def __init__(self, model, kaq, c, z, npor, ltype):
Porosity of the aquifer(s).
ltype : string or array of strings
Type of the layers: 'a' for aquifer, 'l' for leaky layer.
model3d : bool, optional
Whether the model is Model3D (default is False).
"""
# All input variables except model should be numpy arrays
# That should be checked outside this function
Expand Down Expand Up @@ -83,6 +85,7 @@ def __init__(self, model, kaq, c, z, npor, ltype):
self.nporll[1:] = self.npor[self.ltype == "l"]
else:
self.nporll = self.npor[self.ltype == "l"]
self.model3d = model3d

def initialize(self):
self.elementlist = [] # Elementlist of aquifer
Expand Down Expand Up @@ -147,11 +150,9 @@ def summary(self):
Summary of aquifer parameters including layer type, thickness,
hydraulic conductivity, and resistance.
"""
if self.nlayers == self.naq:
model3d = True
if self.model3d:
add_cols = ["kzoverkh"]
else:
model3d = False
add_cols = []
summary = pd.DataFrame(
index=range(self.nlayers),
Expand All @@ -162,22 +163,22 @@ def summary(self):
summary["layer_type"] = [layertype[lt] for lt in self.ltype]
maskaq = self.ltype == "a"
if self.ilap == 1: # confined on top
summary.iloc[maskaq, 2] = self.Haq
summary.iloc[maskaq, 3] = self.kaq
if model3d:
summary.iloc[maskaq, 4] = self.c
summary.iloc[0, 4] = np.nan # reset confined resistance to nan
summary.iloc[maskaq, 5] = self.kzoverkh
summary.loc[maskaq, "H"] = self.Haq
summary.loc[maskaq, "k_h"] = self.kaq
if self.model3d:
summary.loc[maskaq, "c"] = self.c
summary.loc[0, "c"] = np.nan # reset confined resistance to nan
summary.loc[maskaq, "kzoverkh"] = self.kzoverkh
else:
summary.iloc[~maskaq, 2] = self.Hll[1:]
summary.iloc[~maskaq, 4] = self.c[1:]
else:
summary.iloc[maskaq, 2] = self.Haq
summary.iloc[maskaq, 3] = self.kaq
if model3d:
summary.iloc[~maskaq, 2] = self.Hll[0]
summary.iloc[maskaq, 4] = self.c
summary.iloc[maskaq, 5] = self.kzoverkh
summary.loc[maskaq, "H"] = self.Haq
summary.loc[maskaq, "k_h"] = self.kaq
if self.model3d:
summary.loc[~maskaq, "H"] = self.Hll[0]
summary.loc[maskaq, "c"] = self.c
summary.loc[maskaq, "kzoverkh"] = self.kzoverkh
else:
summary.iloc[~maskaq, 2] = self.Hll
summary.iloc[~maskaq, 4] = self.c
Expand All @@ -191,8 +192,8 @@ class Aquifer(AquiferData):
Extends AquiferData and supports inhomogeneities within a model.
"""

def __init__(self, model, kaq, c, z, npor, ltype):
super().__init__(model, kaq, c, z, npor, ltype)
def __init__(self, model, kaq, c, z, npor, ltype, model3d=False):
super().__init__(model, kaq, c, z, npor, ltype, model3d)
self.inhomlist = []
self.area = 1e300 # Needed to find smallest inhom

Expand Down
20 changes: 10 additions & 10 deletions timflow/steady/aquifer_parameters.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ def param_maq(kaq, z, c, npor, top):
z = np.atleast_1d(z).astype("d")
c = np.atleast_1d(c).astype("d")
npor = np.atleast_1d(npor).astype("d")
if top == "conf":
if top.startswith("con"):
Naq = int(len(z) / 2)
ltype = np.array(list((Naq - 1) * "al" + "a"))
else: # leaky layer on top
Expand All @@ -49,7 +49,7 @@ def param_maq(kaq, z, c, npor, top):
assert len(kaq) == Naq, "Error: length of kaq needs to be 1 or" + str(Naq)
H = z[:-1] - z[1:]
assert np.all(H >= 0), "Error: Not all layers thicknesses are non-negative" + str(H)
if top == "conf":
if top.startswith("con"):
if len(c) == 1:
c = c * np.ones(Naq - 1)
if len(npor) == 1:
Expand Down Expand Up @@ -104,10 +104,10 @@ def param_3d(kaq, z, kzoverkh, npor, top="conf", topres=0):
z = np.atleast_1d(z).astype("d")
kzoverkh = np.atleast_1d(kzoverkh).astype("d")
npor = np.atleast_1d(npor).astype("d")
if top == "conf":
if top.startswith("conf"):
Naq = len(z) - 1
ltype = np.array(Naq * ["a"])
elif top == "semi":
elif top.startswith("semi"):
Naq = len(z) - 1
ltype = np.hstack(("l", Naq * ["a"]))
if len(kaq) == 1:
Expand All @@ -117,21 +117,21 @@ def param_3d(kaq, z, kzoverkh, npor, top="conf", topres=0):
kzoverkh = kzoverkh * np.ones(Naq)
assert len(kzoverkh) == Naq, "Error: length of kzoverkh needs to be 1 or" + str(Naq)
if len(npor) == 1:
if top == "conf":
if top.startswith("con"):
npor = npor * np.ones(Naq)
elif top == "semi":
elif top.startswith("sem"):
npor = npor * np.ones(Naq + 1)
if top == "conf":
if top.startswith("con"):
assert len(npor) == Naq, "Error: length of npor needs to be 1 or" + str(Naq)
elif top == "semi":
elif top.startswith("sem"):
assert len(npor) == Naq + 1, "Error: length of npor needs to be 1 or" + str(
Naq + 1
)
H = z[:-1] - z[1:]
assert np.all(H >= 0), "Error: Not all layers thicknesses are non-negative" + str(H)
c = 0.5 * H[:-1] / (kzoverkh[:-1] * kaq[:-1]) + 0.5 * H[1:] / (kzoverkh[1:] * kaq[1:])
if top == "conf":
if top.startswith("con"):
c = np.hstack((1e100, c))
elif top == "semi":
elif top.startswith("sem"):
c = np.hstack((topres + 0.5 * H[0] / (kzoverkh[0] * kaq[0]), c))
return kaq, kzoverkh, c, npor, ltype
16 changes: 12 additions & 4 deletions timflow/steady/inhomogeneity1d.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,10 @@ class Xsection(AquiferData):

tiny = 1e-12

def __init__(self, model, x1, x2, kaq, c, z, npor, ltype, hstar, N, name=None):
super().__init__(model, kaq, c, z, npor, ltype)
def __init__(
self, model, x1, x2, kaq, c, z, npor, ltype, model3d, hstar, N, name=None
):
super().__init__(model, kaq, c, z, npor, ltype, model3d)
self.x1 = x1
self.x2 = x2
self.hstar = hstar
Expand Down Expand Up @@ -339,7 +341,10 @@ def __init__(
npor,
ltype,
) = param_maq(kaq, z, c, npor, topboundary)
super().__init__(model, x1, x2, kaq, c, z, npor, ltype, hstar, N, name=name)
model3d = False
super().__init__(
model, x1, x2, kaq, c, z, npor, ltype, model3d, hstar, N, name=name
)


class Xsection3D(Xsection):
Expand Down Expand Up @@ -418,7 +423,10 @@ def __init__(
) = param_3d(kaq, z, kzoverkh, npor, topboundary, topres)
if topboundary == "semi":
z = np.hstack((z[0] + topthick, z))
super().__init__(model, x1, x2, kaq, c, z, npor, ltype, hstar, N, name=name)
model3d = True
super().__init__(
model, x1, x2, kaq, c, z, npor, ltype, model3d, hstar, N, name=name
)
self.kzoverkh = kzoverkh


Expand Down
7 changes: 4 additions & 3 deletions timflow/steady/model.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -53,12 +53,12 @@ class Model:
'l' leaky layer
"""

def __init__(self, kaq, c, z, npor, ltype):
def __init__(self, kaq, c, z, npor, ltype, model3d=False):
# All input variables are numpy arrays
# That should be checked outside this function
self.elementlist = []
self.elementdict = {} # only elements that have a label
self.aq = Aquifer(self, kaq, c, z, npor, ltype)
self.aq = Aquifer(self, kaq, c, z, npor, ltype, model3d)
self.modelname = "ml" # Used for writing out input
self.name = "Model"
self.model_type = "steady" # Model type for plotting and other purposes
Expand Down Expand Up @@ -838,7 +838,8 @@ def __init__(
)
if topboundary == "semi":
z = np.hstack((z[0] + topthick, z))
super().__init__(kaq, c, z, npor, ltype)
model3d = True
super().__init__(kaq, c, z, npor, ltype, model3d)
self.aq.kzoverkh = kzoverkh # add kzoverkh to aquifer object
self.name = "Model3D"
if self.aq.ltype[0] == "l":
Expand Down