12312733157

Committed 13 Dec 2024 09:08AM UTC coverage: 77.784% (+2.4%) from 75.342%

Build # 12312733157

Build Type

Pull #81

github

Committed by

web-flow

Commit Message

Merge fcd788ee4 into f9283bd02

Pull Request Pull Request #81: Add cross-section models

Run Details

503 of 632 new or added lines in 12 files covered. (79.59%)

42 existing lines in 3 files now uncovered.

3011 of 3871 relevant lines covered (77.78%)

2.33 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

78.33

/ttim/aquifer.py

import numpy as np


class AquiferData:
    def __init__(
        self,
        model,
        kaq,
        z,
        Haq,
        Hll,
        c,
        Saq,
        Sll,
        poraq,
        porll,
        ltype,
        topboundary,
        phreatictop,
        kzoverkh=None,
        model3d=False,
        name=None,
    ):
        """Kzoverkh and model3d only need to be specified when model is model3d."""
        self.model = model
        self.kaq = np.atleast_1d(kaq).astype(float)
        self.z = np.atleast_1d(z).astype(float)
        self.naq = len(self.kaq)
        self.nlayers = len(self.z) - 1
        self.Haq = np.atleast_1d(Haq).astype(float)
        self.Hll = np.atleast_1d(Hll).astype(float)
        self.T = self.kaq * self.Haq
        self.Tcol = self.T.reshape(self.naq, 1)
        self.c = np.atleast_1d(c).astype(float)
        self.c[self.c > 1e100] = 1e100
        self.Saq = np.atleast_1d(Saq).astype(float)
        self.Sll = np.atleast_1d(Sll).astype(float)
        self.Sll[self.Sll < 1e-20] = 1e-20  # Cannot be zero
        self.poraq = np.atleast_1d(poraq).astype(float)
        self.porll = np.atleast_1d(porll).astype(float)
        self.ltype = np.atleast_1d(ltype)
        self.zaqtop = self.z[:-1][self.ltype == "a"]
        self.zaqbot = self.z[1:][self.ltype == "a"]
        self.layernumber = np.zeros(self.nlayers, dtype="int")
        self.layernumber[self.ltype == "a"] = np.arange(self.naq)
        self.layernumber[self.ltype == "l"] = np.arange(self.nlayers - self.naq)
        if self.ltype[0] == "a":
            # first leaky layer below first aquifer layer
            self.layernumber[self.ltype == "l"] += 1
        self.topboundary = topboundary[:3]
        self.phreatictop = phreatictop
        self.kzoverkh = kzoverkh
        if self.kzoverkh is not None:
            self.kzoverkh = np.atleast_1d(self.kzoverkh).astype(float)
            if len(self.kzoverkh) == 1:
                self.kzoverkh = self.kzoverkh * np.ones(self.naq)
        self.model3d = model3d
        if self.model3d:
            assert self.kzoverkh is not None, "model3d specified without kzoverkh"
        # self.D = self.T / self.Saq
        self.area = 1e200  # Smaller than default of ml.aq so that inhom is found
        self.name = name

    def __repr__(self):
        if self.topboundary.startswith("con"):
            topbound = "confined"
        elif self.topboundary.startswith("lea"):
            topbound = "leaky"
        elif self.topboundary.startswith("sem"):
            topbound = "semi-confined"
        else:
            topbound = "unknown"  # should not happen
        return f"Inhom Aquifer: {self.naq} aquifer(s) with {topbound} top boundary"

    def initialize(self):
        """Initialize the aquifer data.

        Attributes
        ----------
        eigval[naq, npval]:
            Array with eigenvalues
        lab[naq, npval]:
            Array with lambda values
        lab2[naq, nint, npint]:
            Array with lambda values reorganized per interval
        eigvec[naq, naq, npval]:
            Array with eigenvector matrices
        coef[naq ,naq, npval]:
            Array with coefficients coef[ilayers, :, np] are the coefficients if the
            element is in ilayers belonging to Laplace parameter number np.
        """
        # Recompute T for when kaq is changed
        self.T = self.kaq * self.Haq
        self.Tcol = self.T.reshape(self.naq, 1)
        # Compute Saq and Sll
        self.Scoefaq = self.Saq * self.Haq
        self.Scoefll = self.Sll * self.Hll
        if (self.topboundary == "con") and self.phreatictop:
            self.Scoefaq[0] = self.Scoefaq[0] / self.Haq[0]
        elif (self.topboundary == "lea") and self.phreatictop:
            self.Scoefll[0] = self.Scoefll[0] / self.Hll[0]
        self.D = self.T / self.Scoefaq
        # Compute c if model3d for when kaq is changed
        if self.model3d:
            self.c[1:] = 0.5 * self.Haq[:-1] / (
                self.kzoverkh[:-1] * self.kaq[:-1]
            ) + 0.5 * self.Haq[1:] / (self.kzoverkh[1:] * self.kaq[1:])
        #
        self.eigval = np.zeros((self.naq, self.model.npval), dtype=complex)
        self.lab = np.zeros((self.naq, self.model.npval), dtype=complex)
        self.eigvec = np.zeros((self.naq, self.naq, self.model.npval), dtype=complex)
        self.coef = np.zeros((self.naq, self.naq, self.model.npval), dtype=complex)
        b = np.diag(np.ones(self.naq))
        for i in range(self.model.npval):
            w, v = self.compute_lab_eigvec(self.model.p[i])
            # Eigenvectors are columns of v
            self.eigval[:, i] = w
            self.eigvec[:, :, i] = v
            self.coef[:, :, i] = np.linalg.solve(v, b).T
        self.lab = 1.0 / np.sqrt(self.eigval)
        self.lab2 = self.lab.copy()
        self.lab2.shape = (self.naq, self.model.nint, self.model.npint)
        self.lababs = np.abs(self.lab2[:, :, 0])  # used to check distances
        self.eigvec2 = self.eigvec.copy()
        self.eigvec2.shape = (self.naq, self.naq, self.model.nint, self.model.npint)

    def compute_lab_eigvec(self, p, returnA=False, B=None):
        sqrtpSc = np.sqrt(p * self.Scoefll * self.c)
        a, b = np.zeros_like(sqrtpSc), np.zeros_like(sqrtpSc)
        small = np.abs(sqrtpSc) < 200
        a[small] = sqrtpSc[small] / np.tanh(sqrtpSc[small])
        b[small] = sqrtpSc[small] / np.sinh(sqrtpSc[small])
        a[~small] = sqrtpSc[~small] / (
            (1.0 - np.exp(-2.0 * sqrtpSc[~small]))
            / (1.0 + np.exp(-2.0 * sqrtpSc[~small]))
        )
        b[~small] = (
            sqrtpSc[~small]
            * 2.0
            * np.exp(-sqrtpSc[~small])
            / (1.0 - np.exp(-2.0 * sqrtpSc[~small]))
        )
        if (self.topboundary[:3] == "sem") or (self.topboundary[:3] == "lea"):
            dzero = sqrtpSc[0] * np.tanh(sqrtpSc[0])

        d0 = p / self.D
        if B is not None:
            d0 = d0 * B  # B is vector of load efficiency paramters
        d0[:-1] += a[1:] / (self.c[1:] * self.T[:-1])
        d0[1:] += a[1:] / (self.c[1:] * self.T[1:])
        if self.topboundary[:3] == "lea":
            d0[0] += dzero / (self.c[0] * self.T[0])
        elif self.topboundary[:3] == "sem":
            d0[0] += a[0] / (self.c[0] * self.T[0])

        dm1 = -b[1:] / (self.c[1:] * self.T[:-1])
        dp1 = -b[1:] / (self.c[1:] * self.T[1:])
        A = np.diag(dm1, -1) + np.diag(d0, 0) + np.diag(dp1, 1)
        if returnA:
            return A
        w, v = np.linalg.eig(A)
        # sorting moved here
        index = np.argsort(abs(w))[::-1]
        w = w[index]
        v = v[:, index]
        return w, v

    def head_to_potential(self, h, layers):
        return h * self.Tcol[layers]

    def potential_to_head(self, pot, layers):
        return pot / self.Tcol[layers]

    def is_inside(self, x, y):
        raise NotImplementedError(
            f"Must overload is_inside in {self.__class__.__name__} class."
        )

    def in_which_layer(self, z):
        """Get layer given elevation z.

        Returns -9999 if above top of system, +9999 if below bottom of system, negative
        for in leaky layer.

        leaky layer -n is on top of aquifer n
        """
        if z > self.zt[0]:
            return -9999
        for i in range(self.naq - 1):
            if z >= self.zb[i]:
                return i
            if z > self.zt[i + 1]:
                return -i - 1
        if z >= self.zb[self.naq - 1]:
            return self.naq - 1
        return +9999

    def findlayer(self, z):
        """Returns layer-number, layer-type and model-layer-number."""
        if z > self.z[0]:
            modellayer = -1
            ltype = "above"
            layernumber = None
        elif z < self.z[-1]:
            modellayer = len(self.layernumber)
            ltype = "below"
            layernumber = None
        else:
            modellayer = np.argwhere((z <= self.z[:-1]) & (z >= self.z[1:]))[0, 0]
            layernumber = self.layernumber[modellayer]
            ltype = self.ltype[modellayer]
        return layernumber, ltype, modellayer


class Aquifer(AquiferData):
    def __init__(
        self,
        model,
        kaq,
        z,
        Haq,
        Hll,
        c,
        Saq,
        Sll,
        poraq,
        porll,
        ltype,
        topboundary,
        phreatictop,
        kzoverkh=None,
        model3d=False,
    ):
        super().__init__(
            model,
            kaq,
            z,
            Haq,
            Hll,
            c,
            Saq,
            Sll,
            poraq,
            porll,
            ltype,
            topboundary,
            phreatictop,
            kzoverkh,
            model3d,
        )
        self.inhomdict = {}
        self.area = 1e300  # Needed to find smallest inhomogeneity

    def __repr__(self):
        if self.topboundary.startswith("con"):
            topbound = "confined"
        elif self.topboundary.startswith("lea"):
            topbound = "leaky"
        elif self.topboundary.startswith("sem"):
            topbound = "semi-confined"
        else:
            topbound = "unknown"  # should not happen
        return f"Background Aquifer: {self.naq} aquifer(s) with {topbound} top boundary"

    def initialize(self):
        super().initialize()
        for inhom in self.inhomdict.values():
            inhom.initialize()

    def is_inside(self, x, y):
        return True

    def find_aquifer_data(self, x, y):
        rv = self
        for aq in self.inhomdict.values():
            if aq.is_inside(x, y):
                if aq.area < rv.area:
                    rv = aq
        return rv

    def add_inhom(self, inhom):
        inhom_number = len(self.inhomdict)
        if inhom.name is None:
            inhom.name = f"inhom{inhom_number}"
        if inhom.name in self.inhomdict:
            raise ValueError(f"Inhomogeneity name '{inhom.name}' already exists.")
        self.inhomdict[inhom.name] = inhom
        return inhom_number


class SimpleAquifer(Aquifer):
    def __init__(self, naq):
        self.naq = naq
        self.inhomdict = {}
        self.area = 1e300  # Needed to find smallest inhomogeneity

    def __repr__(self):
        return f"Simple Aquifer: {self.naq} aquifer(s)"

    def initialize(self):
        for inhom in self.inhomdict.values():
            inhom.initialize()

1	import numpy as np	3✔
2
3
4	class AquiferData:	3✔
5	def __init__(	3✔
6	self,
7	model,
8	kaq,
9	z,
10	Haq,
11	Hll,
12	c,
13	Saq,
14	Sll,
15	poraq,
16	porll,
17	ltype,
18	topboundary,
19	phreatictop,
20	kzoverkh=None,
21	model3d=False,
22	name=None,
23	):
24	"""Kzoverkh and model3d only need to be specified when model is model3d."""
25	self.model = model	3✔
26	self.kaq = np.atleast_1d(kaq).astype(float)	3✔
27	self.z = np.atleast_1d(z).astype(float)	3✔
28	self.naq = len(self.kaq)	3✔
29	self.nlayers = len(self.z) - 1	3✔
30	self.Haq = np.atleast_1d(Haq).astype(float)	3✔
31	self.Hll = np.atleast_1d(Hll).astype(float)	3✔
32	self.T = self.kaq * self.Haq	3✔
33	self.Tcol = self.T.reshape(self.naq, 1)	3✔
34	self.c = np.atleast_1d(c).astype(float)	3✔
35	self.c[self.c > 1e100] = 1e100	3✔
36	self.Saq = np.atleast_1d(Saq).astype(float)	3✔
37	self.Sll = np.atleast_1d(Sll).astype(float)	3✔
38	self.Sll[self.Sll < 1e-20] = 1e-20 # Cannot be zero	3✔
39	self.poraq = np.atleast_1d(poraq).astype(float)	3✔
40	self.porll = np.atleast_1d(porll).astype(float)	3✔
41	self.ltype = np.atleast_1d(ltype)	3✔
42	self.zaqtop = self.z[:-1][self.ltype == "a"]	3✔
43	self.zaqbot = self.z[1:][self.ltype == "a"]	3✔
44	self.layernumber = np.zeros(self.nlayers, dtype="int")	3✔
45	self.layernumber[self.ltype == "a"] = np.arange(self.naq)	3✔
46	self.layernumber[self.ltype == "l"] = np.arange(self.nlayers - self.naq)	3✔
47	if self.ltype[0] == "a":	3✔
48	# first leaky layer below first aquifer layer
49	self.layernumber[self.ltype == "l"] += 1	3✔
50	self.topboundary = topboundary[:3]	3✔
51	self.phreatictop = phreatictop	3✔
52	self.kzoverkh = kzoverkh	3✔
53	if self.kzoverkh is not None:	3✔
54	self.kzoverkh = np.atleast_1d(self.kzoverkh).astype(float)	3✔
55	if len(self.kzoverkh) == 1:	3✔
56	self.kzoverkh = self.kzoverkh * np.ones(self.naq)	3✔
57	self.model3d = model3d	3✔
58	if self.model3d:	3✔
59	assert self.kzoverkh is not None, "model3d specified without kzoverkh"	3✔
60	# self.D = self.T / self.Saq
61	self.area = 1e200 # Smaller than default of ml.aq so that inhom is found	3✔
62	self.name = name	3✔
63
64	def __repr__(self):	3✔
NEW 65	if self.topboundary.startswith("con"):	×
NEW 66	topbound = "confined"	×
NEW 67	elif self.topboundary.startswith("lea"):	×
NEW 68	topbound = "leaky"	×
NEW 69	elif self.topboundary.startswith("sem"):	×
NEW 70	topbound = "semi-confined"	×
71	else:
NEW 72	topbound = "unknown" # should not happen	×
NEW 73	return f"Inhom Aquifer: {self.naq} aquifer(s) with {topbound} top boundary"	×
74
75	def initialize(self):	3✔
76	"""Initialize the aquifer data.
77
78	Attributes
79	----------
80	eigval[naq, npval]:
81	Array with eigenvalues
82	lab[naq, npval]:
83	Array with lambda values
84	lab2[naq, nint, npint]:
85	Array with lambda values reorganized per interval
86	eigvec[naq, naq, npval]:
87	Array with eigenvector matrices
88	coef[naq ,naq, npval]:
89	Array with coefficients coef[ilayers, :, np] are the coefficients if the
90	element is in ilayers belonging to Laplace parameter number np.
91	"""
92	# Recompute T for when kaq is changed
93	self.T = self.kaq * self.Haq	3✔
94	self.Tcol = self.T.reshape(self.naq, 1)	3✔
95	# Compute Saq and Sll
96	self.Scoefaq = self.Saq * self.Haq	3✔
97	self.Scoefll = self.Sll * self.Hll	3✔
98	if (self.topboundary == "con") and self.phreatictop:	3✔
99	self.Scoefaq[0] = self.Scoefaq[0] / self.Haq[0]	3✔
100	elif (self.topboundary == "lea") and self.phreatictop:	3✔
101	self.Scoefll[0] = self.Scoefll[0] / self.Hll[0]	×
102	self.D = self.T / self.Scoefaq	3✔
103	# Compute c if model3d for when kaq is changed
104	if self.model3d:	3✔
105	self.c[1:] = 0.5 * self.Haq[:-1] / (	3✔
106	self.kzoverkh[:-1] * self.kaq[:-1]
107	) + 0.5 * self.Haq[1:] / (self.kzoverkh[1:] * self.kaq[1:])
108	#
109	self.eigval = np.zeros((self.naq, self.model.npval), dtype=complex)	3✔
110	self.lab = np.zeros((self.naq, self.model.npval), dtype=complex)	3✔
111	self.eigvec = np.zeros((self.naq, self.naq, self.model.npval), dtype=complex)	3✔
112	self.coef = np.zeros((self.naq, self.naq, self.model.npval), dtype=complex)	3✔
113	b = np.diag(np.ones(self.naq))	3✔
114	for i in range(self.model.npval):	3✔
115	w, v = self.compute_lab_eigvec(self.model.p[i])	3✔
116	# Eigenvectors are columns of v
117	self.eigval[:, i] = w	3✔
118	self.eigvec[:, :, i] = v	3✔
119	self.coef[:, :, i] = np.linalg.solve(v, b).T	3✔
120	self.lab = 1.0 / np.sqrt(self.eigval)	3✔
121	self.lab2 = self.lab.copy()	3✔
122	self.lab2.shape = (self.naq, self.model.nint, self.model.npint)	3✔
123	self.lababs = np.abs(self.lab2[:, :, 0]) # used to check distances	3✔
124	self.eigvec2 = self.eigvec.copy()	3✔
125	self.eigvec2.shape = (self.naq, self.naq, self.model.nint, self.model.npint)	3✔
126
127	def compute_lab_eigvec(self, p, returnA=False, B=None):	3✔
128	sqrtpSc = np.sqrt(p * self.Scoefll * self.c)	3✔
129	a, b = np.zeros_like(sqrtpSc), np.zeros_like(sqrtpSc)	3✔
130	small = np.abs(sqrtpSc) < 200	3✔
131	a[small] = sqrtpSc[small] / np.tanh(sqrtpSc[small])	3✔
132	b[small] = sqrtpSc[small] / np.sinh(sqrtpSc[small])	3✔
133	a[~small] = sqrtpSc[~small] / (	3✔
134	(1.0 - np.exp(-2.0 * sqrtpSc[~small]))
135	/ (1.0 + np.exp(-2.0 * sqrtpSc[~small]))
136	)
137	b[~small] = (	3✔
138	sqrtpSc[~small]
139	* 2.0
140	* np.exp(-sqrtpSc[~small])
141	/ (1.0 - np.exp(-2.0 * sqrtpSc[~small]))
142	)
143	if (self.topboundary[:3] == "sem") or (self.topboundary[:3] == "lea"):	3✔
144	dzero = sqrtpSc[0] * np.tanh(sqrtpSc[0])	3✔
145
146	d0 = p / self.D	3✔
147	if B is not None:	3✔
148	d0 = d0 * B # B is vector of load efficiency paramters	×
149	d0[:-1] += a[1:] / (self.c[1:] * self.T[:-1])	3✔
150	d0[1:] += a[1:] / (self.c[1:] * self.T[1:])	3✔
151	if self.topboundary[:3] == "lea":	3✔
152	d0[0] += dzero / (self.c[0] * self.T[0])	×
153	elif self.topboundary[:3] == "sem":	3✔
154	d0[0] += a[0] / (self.c[0] * self.T[0])	3✔
155
156	dm1 = -b[1:] / (self.c[1:] * self.T[:-1])	3✔
157	dp1 = -b[1:] / (self.c[1:] * self.T[1:])	3✔
158	A = np.diag(dm1, -1) + np.diag(d0, 0) + np.diag(dp1, 1)	3✔
159	if returnA:	3✔
160	return A	×
161	w, v = np.linalg.eig(A)	3✔
162	# sorting moved here
163	index = np.argsort(abs(w))[::-1]	3✔
164	w = w[index]	3✔
165	v = v[:, index]	3✔
166	return w, v	3✔
167
168	def head_to_potential(self, h, layers):	3✔
169	return h * self.Tcol[layers]	×
170
171	def potential_to_head(self, pot, layers):	3✔
172	return pot / self.Tcol[layers]	3✔
173
174	def is_inside(self, x, y):	3✔
NEW 175	raise NotImplementedError(	×
176	f"Must overload is_inside in {self.__class__.__name__} class."
177	)
178
179	def in_which_layer(self, z):	3✔
180	"""Get layer given elevation z.
181
182	Returns -9999 if above top of system, +9999 if below bottom of system, negative
183	for in leaky layer.
184
185	leaky layer -n is on top of aquifer n
186	"""
187	if z > self.zt[0]:	×
188	return -9999	×
189	for i in range(self.naq - 1):	×
190	if z >= self.zb[i]:	×
191	return i	×
192	if z > self.zt[i + 1]:	×
193	return -i - 1	×
194	if z >= self.zb[self.naq - 1]:	×
195	return self.naq - 1	×
196	return +9999	×
197
198	def findlayer(self, z):	3✔
199	"""Returns layer-number, layer-type and model-layer-number."""
200	if z > self.z[0]:	3✔
201	modellayer = -1	3✔
202	ltype = "above"	3✔
203	layernumber = None	3✔
204	elif z < self.z[-1]:	3✔
205	modellayer = len(self.layernumber)	×
206	ltype = "below"	×
207	layernumber = None	×
208	else:
209	modellayer = np.argwhere((z <= self.z[:-1]) & (z >= self.z[1:]))[0, 0]	3✔
210	layernumber = self.layernumber[modellayer]	3✔
211	ltype = self.ltype[modellayer]	3✔
212	return layernumber, ltype, modellayer	3✔
213
214
215	class Aquifer(AquiferData):	3✔
216	def __init__(	3✔
217	self,
218	model,
219	kaq,
220	z,
221	Haq,
222	Hll,
223	c,
224	Saq,
225	Sll,
226	poraq,
227	porll,
228	ltype,
229	topboundary,
230	phreatictop,
231	kzoverkh=None,
232	model3d=False,
233	):
234	super().__init__(	3✔
235	model,
236	kaq,
237	z,
238	Haq,
239	Hll,
240	c,
241	Saq,
242	Sll,
243	poraq,
244	porll,
245	ltype,
246	topboundary,
247	phreatictop,
248	kzoverkh,
249	model3d,
250	)
251	self.inhomdict = {}	3✔
252	self.area = 1e300 # Needed to find smallest inhomogeneity	3✔
253
254	def __repr__(self):	3✔
NEW 255	if self.topboundary.startswith("con"):	×
NEW 256	topbound = "confined"	×
NEW 257	elif self.topboundary.startswith("lea"):	×
NEW 258	topbound = "leaky"	×
NEW 259	elif self.topboundary.startswith("sem"):	×
NEW 260	topbound = "semi-confined"	×
261	else:
NEW 262	topbound = "unknown" # should not happen	×
NEW 263	return f"Background Aquifer: {self.naq} aquifer(s) with {topbound} top boundary"	×
264
265	def initialize(self):	3✔
266	super().initialize()	3✔
267	for inhom in self.inhomdict.values():	3✔
UNCOV 268	inhom.initialize()	×
269
270	def is_inside(self, x, y):	3✔
NEW 271	return True	×
272
273	def find_aquifer_data(self, x, y):	3✔
274	rv = self	3✔
275	for aq in self.inhomdict.values():	3✔
276	if aq.is_inside(x, y):	3✔
277	if aq.area < rv.area:	3✔
278	rv = aq	3✔
279	return rv	3✔
280
281	def add_inhom(self, inhom):	3✔
282	inhom_number = len(self.inhomdict)	3✔
283	if inhom.name is None:	3✔
284	inhom.name = f"inhom{inhom_number}"	3✔
285	if inhom.name in self.inhomdict:	3✔
NEW 286	raise ValueError(f"Inhomogeneity name '{inhom.name}' already exists.")	×
287	self.inhomdict[inhom.name] = inhom	3✔
288	return inhom_number	3✔
289
290
291	class SimpleAquifer(Aquifer):	3✔
292	def __init__(self, naq):	3✔
293	self.naq = naq	3✔
294	self.inhomdict = {}	3✔
295	self.area = 1e300 # Needed to find smallest inhomogeneity	3✔
296
297	def __repr__(self):	3✔
NEW 298	return f"Simple Aquifer: {self.naq} aquifer(s)"	×
299
300	def initialize(self):	3✔
301	for inhom in self.inhomdict.values():	3✔
302	inhom.initialize()	3✔

mbakker7 / ttim / 12312733157

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous