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

92.0

/ttim/linesink1d.py

import matplotlib.pyplot as plt
import numpy as np

from ttim.element import Element
from ttim.equation import (
    FluxDiffEquation,
    HeadDiffEquation,
    HeadEquation,
    MscreenEquation,
)


class LineSink1DBase(Element):
    """LineSink1D Base Class.

    All LineSink1D elements are derived from this class
    """

    tiny = 1e-8

    def __init__(
        self,
        model,
        xls=0,
        tsandbc=[(0, 1)],
        res=0,
        wh="H",
        layers=0,
        type="",
        name="LineSink1DBase",
        label=None,
        aq=None,
        inhomelement=False,
    ):
        super().__init__(
            model,
            nparam=1,
            nunknowns=0,
            layers=layers,
            tsandbc=tsandbc,
            type=type,
            name=name,
            label=label,
            inhomelement=inhomelement,
        )
        # Defined here and not in Element as other elements can have multiple
        # parameters per layers
        self.nparam = len(self.layers)
        self.xls = float(xls)
        self.res = np.atleast_1d(res).astype(np.float64)
        self.wh = wh
        self.aq = aq
        self.model.addelement(self)

    def __repr__(self):
        return self.name + " at " + str(self.xls)

    def initialize(self):
        # Control point to make sure the point is always the same for
        # all elements
        self.xc = np.array([self.xls])
        self.yc = np.zeros(1)
        self.ncp = 1
        if self.aq is None:
            self.aq = self.model.aq.find_aquifer_data(self.xc[0], self.yc[0])
        self.setbc()
        coef = self.aq.coef[self.layers, :]
        self.setflowcoef()
        # term is shape (self.nparam,self.aq.naq,self.model.npval)
        self.term = self.flowcoef * coef
        self.term2 = self.term.reshape(
            self.nparam, self.aq.naq, self.model.nint, self.model.npint
        )
        self.dischargeinf = self.flowcoef * coef
        self.dischargeinflayers = np.sum(
            self.dischargeinf * self.aq.eigvec[self.layers, :, :], 1
        )
        if isinstance(self.wh, str):
            if self.wh == "H":
                self.wh = self.aq.Haq[self.layers]
            elif self.wh == "2H":
                self.wh = 2.0 * self.aq.Haq[self.layers]
        else:
            self.wh = np.atleast_1d(self.wh) * np.ones(self.nlayers)
        # Q = (h - hls) / resfach
        self.resfach = self.res / (self.wh)
        # Q = (Phi - Phils) / resfacp
        self.resfacp = self.resfach * self.aq.T[self.layers]

    def setflowcoef(self):
        """Separate function so that this can be overloaded for other types."""
        self.flowcoef = 1.0 / self.model.p  # Step function

    def potinf(self, x, _, aq=None):
        """Can be called with only one x value."""
        if aq is None:
            aq = self.model.aq.find_aquifer_data(x, 0.0)
        rv = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
        if aq == self.aq:
            if (x - self.xls) < 0.0:
                pot = -0.5 * aq.lab * np.exp((x - self.xls) / aq.lab)
            elif (x - self.xls) >= 0.0:
                pot = -0.5 * aq.lab * np.exp(-(x - self.xls) / aq.lab)
            else:
                raise ValueError("Something wrong with passed x value.")
            rv[:] = self.term * pot
        return rv

    def disvecinf(self, x, y, aq=None):
        """Can be called with only one x,y value."""
        if aq is None:
            aq = self.model.aq.find_aquifer_data(x, y)
        rvx = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
        rvy = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)
        if aq == self.aq:
            if (x - self.xls) < 0.0:
                qx = 0.5 * np.exp((x - self.xls) / aq.lab)
            elif (x - self.xls) >= 0.0:
                qx = -0.5 * np.exp(-(x - self.xls) / aq.lab)
            else:
                raise ValueError("Something wrong with passed x value.")
            rvx[:] = self.term * qx
        return rvx, rvy

    def changetrace(
        self, xyzt1, xyzt2, aq, layer, ltype, modellayer, direction, hstepmax
    ):
        raise NotImplementedError("changetrace not implemented for this element")

    def plot(self, ax=None):
        if ax is None:
            _, ax = plt.subplots()
        aq = self.model.aq.find_aquifer_data(self.xls, 0.0)
        for ilay in self.layers:
            ax.plot(
                [self.xls, self.xls],
                [aq.zaqtop[ilay], aq.zaqbot[ilay]],
                "k-",
            )


class DischargeLineSink1D(LineSink1DBase):
    r"""Linesink1D with a specified discharge for each layer that the linesink.

    Parameters
    ----------
    model : Model object
        model to which the element is added
    x : float
        x-coordinate of the linesink
    tsandq : list of tuples
        tuples of starting time and specific discharge after starting time
    res : float
        resistance of the linesink
    layers : int, array or list
        layer (int) or layers (list or array) in which linesink is located
    label : string or None (default: None)
        label of the linesink

    Examples
    --------
    Example of an infinitely long linesink that pumps with a specific discharge of
    100 between times 10 and 50, with a specific discharge of 20 between
    times 50 and 200, and zero speficic discharge after time 200.

    >>> DischargeLineSink1D(ml, tsandq=[(10, 100), (50, 20), (200, 0)])
    """

    def __init__(
        self, model, xls=0, tsandq=[(0, 1)], res=0, wh="H", layers=0, label=None
    ):
        super().__init__(
            model,
            xls,
            tsandbc=tsandq,
            res=res,
            wh=wh,
            layers=layers,
            type="g",
            name="DischargeLineSink1D",
            label=label,
        )


class LineSink1D(LineSink1DBase, MscreenEquation):
    r"""Linesink1D with a specified discharge.

    Parameters
    ----------
    model : Model object
        model to which the element is added
    x : float
        x-coordinate of the linesink
    tsandq : list of tuples
        tuples of starting time and specific discharge after starting time
    res : float
        resistance of the linesink
    layers : int, array or list
        layer (int) or layers (list or array) in which linesink is located
    label : string or None (default: None)
        label of the linesink

    Examples
    --------
    Example of an infinitely long linesink that pumps with a specific discharge of
    100 between times 10 and 50, with a specific discharge of 20 between
    times 50 and 200, and zero speficic discharge after time 200.

    >>> LineSink1D(ml, tsandq=[(10, 100), (50, 20), (200, 0)])
    """

    def __init__(
        self,
        model,
        xls=0,
        tsandq=[(0, 1)],
        res=0,
        wh="H",
        vres=0.0,
        wv=1.0,
        layers=0,
        label=None,
    ):
        super().__init__(
            model,
            xls,
            tsandbc=tsandq,
            res=res,
            wh=wh,
            layers=layers,
            type="v",
            name="LineSink1D",
            label=label,
        )
        self.nunknowns = self.nparam
        # Vertical resistance inside line-sink
        self.vres = np.atleast_1d(vres)
        self.wv = wv
        if len(self.vres) == 1:
            self.vres = self.vres[0] * np.ones(self.nlayers - 1)

    def initialize(self):
        super().initialize()
        self.parameters = np.zeros(
            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
        )
        # qv = (hn - hn-1) / vresfac[n - 1]
        self.vresfac = self.vres / self.wv


class HeadLineSink1D(LineSink1DBase, HeadEquation):
    def __init__(
        self,
        model,
        xls=0,
        tsandh=[(0, 1)],
        res=0,
        wh="H",
        layers=0,
        label=None,
    ):
        super().__init__(
            model,
            xls,
            tsandbc=tsandh,
            res=res,
            wh=wh,
            layers=layers,
            type="v",
            name="HeadLineSink1D",
            label=label,
        )
        self.nunknowns = self.nparam

    def initialize(self):
        super().initialize()
        self.parameters = np.zeros(
            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
        )
        # Needed in solving, solve for a unit head
        self.pc = self.aq.T[self.layers]


class HeadDiffLineSink1D(LineSink1DBase, HeadDiffEquation):
    def __init__(
        self,
        model,
        xls=0,
        layers=0,
        label=None,
        aq=None,
    ):
        super().__init__(
            model,
            xls,
            tsandbc=[(0, 0)],
            res=0.0,
            wh="H",
            layers=layers,
            type="z",
            name="HeadDiffLineSink1D",
            label=label,
            aq=aq,
            inhomelement=True,
        )
        self.nunknowns = self.nparam

    def initialize(self):
        super().initialize()
        self.xcout = self.xc + self.tiny
        self.xcin = self.xc - self.tiny
        self.ycout = np.zeros(1)
        self.ycin = np.zeros(1)
        self.cosout = -np.ones(1)
        self.sinout = np.zeros(1)
        self.aqout = self.model.aq.find_aquifer_data(self.xcout[0], 0)
        self.aqin = self.model.aq.find_aquifer_data(self.xcin[0], 0)

        self.parameters = np.zeros(
            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
        )


class FluxDiffLineSink1D(LineSink1DBase, FluxDiffEquation):
    def __init__(self, model, xls=0, layers=0, label=None, aq=None):
        super().__init__(
            model,
            xls,
            tsandbc=[(0, 0)],
            res=0.0,
            wh="H",
            layers=layers,
            type="z",
            name="FluxDiffLineSink1D",
            label=label,
            aq=aq,
            inhomelement=True,
        )
        self.nunknowns = self.nparam

    def initialize(self):
        super().initialize()
        self.xcout = self.xc - self.tiny
        self.xcin = self.xc + self.tiny
        self.ycout = np.zeros(1)
        self.ycin = np.zeros(1)
        self.cosout = -np.ones(1)
        self.sinout = np.zeros(1)
        self.aqout = self.model.aq.find_aquifer_data(self.xcout[0], 0)
        self.aqin = self.model.aq.find_aquifer_data(self.xcin[0], 0)
        self.parameters = np.zeros(
            (self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
        )

1	import matplotlib.pyplot as plt	3✔
2	import numpy as np	3✔
3
4	from ttim.element import Element	3✔
5	from ttim.equation import (	3✔
6	FluxDiffEquation,
7	HeadDiffEquation,
8	HeadEquation,
9	MscreenEquation,
10	)
11
12
13	class LineSink1DBase(Element):	3✔
14	"""LineSink1D Base Class.
15
16	All LineSink1D elements are derived from this class
17	"""
18
19	tiny = 1e-8	3✔
20
21	def __init__(	3✔
22	self,
23	model,
24	xls=0,
25	tsandbc=[(0, 1)],
26	res=0,
27	wh="H",
28	layers=0,
29	type="",
30	name="LineSink1DBase",
31	label=None,
32	aq=None,
33	inhomelement=False,
34	):
35	super().__init__(	3✔
36	model,
37	nparam=1,
38	nunknowns=0,
39	layers=layers,
40	tsandbc=tsandbc,
41	type=type,
42	name=name,
43	label=label,
44	inhomelement=inhomelement,
45	)
46	# Defined here and not in Element as other elements can have multiple
47	# parameters per layers
48	self.nparam = len(self.layers)	3✔
49	self.xls = float(xls)	3✔
50	self.res = np.atleast_1d(res).astype(np.float64)	3✔
51	self.wh = wh	3✔
52	self.aq = aq	3✔
53	self.model.addelement(self)	3✔
54
55	def __repr__(self):	3✔
NEW 56	return self.name + " at " + str(self.xls)	×
57
58	def initialize(self):	3✔
59	# Control point to make sure the point is always the same for
60	# all elements
61	self.xc = np.array([self.xls])	3✔
62	self.yc = np.zeros(1)	3✔
63	self.ncp = 1	3✔
64	if self.aq is None:	3✔
65	self.aq = self.model.aq.find_aquifer_data(self.xc[0], self.yc[0])	3✔
66	self.setbc()	3✔
67	coef = self.aq.coef[self.layers, :]	3✔
68	self.setflowcoef()	3✔
69	# term is shape (self.nparam,self.aq.naq,self.model.npval)
70	self.term = self.flowcoef * coef	3✔
71	self.term2 = self.term.reshape(	3✔
72	self.nparam, self.aq.naq, self.model.nint, self.model.npint
73	)
74	self.dischargeinf = self.flowcoef * coef	3✔
75	self.dischargeinflayers = np.sum(	3✔
76	self.dischargeinf * self.aq.eigvec[self.layers, :, :], 1
77	)
78	if isinstance(self.wh, str):	3✔
79	if self.wh == "H":	3✔
80	self.wh = self.aq.Haq[self.layers]	3✔
NEW 81	elif self.wh == "2H":	×
NEW 82	self.wh = 2.0 * self.aq.Haq[self.layers]	×
83	else:
NEW 84	self.wh = np.atleast_1d(self.wh) * np.ones(self.nlayers)	×
85	# Q = (h - hls) / resfach
86	self.resfach = self.res / (self.wh)	3✔
87	# Q = (Phi - Phils) / resfacp
88	self.resfacp = self.resfach * self.aq.T[self.layers]	3✔
89
90	def setflowcoef(self):	3✔
91	"""Separate function so that this can be overloaded for other types."""
92	self.flowcoef = 1.0 / self.model.p # Step function	3✔
93
94	def potinf(self, x, _, aq=None):	3✔
95	"""Can be called with only one x value."""
96	if aq is None:	3✔
NEW 97	aq = self.model.aq.find_aquifer_data(x, 0.0)	×
98	rv = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)	3✔
99	if aq == self.aq:	3✔
100	if (x - self.xls) < 0.0:	3✔
101	pot = -0.5 * aq.lab * np.exp((x - self.xls) / aq.lab)	3✔
102	elif (x - self.xls) >= 0.0:	3✔
103	pot = -0.5 * aq.lab * np.exp(-(x - self.xls) / aq.lab)	3✔
104	else:
NEW 105	raise ValueError("Something wrong with passed x value.")	×
106	rv[:] = self.term * pot	3✔
107	return rv	3✔
108
109	def disvecinf(self, x, y, aq=None):	3✔
110	"""Can be called with only one x,y value."""
111	if aq is None:	3✔
NEW 112	aq = self.model.aq.find_aquifer_data(x, y)	×
113	rvx = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)	3✔
114	rvy = np.zeros((self.nparam, aq.naq, self.model.npval), dtype=complex)	3✔
115	if aq == self.aq:	3✔
116	if (x - self.xls) < 0.0:	3✔
117	qx = 0.5 * np.exp((x - self.xls) / aq.lab)	3✔
118	elif (x - self.xls) >= 0.0:	3✔
119	qx = -0.5 * np.exp(-(x - self.xls) / aq.lab)	3✔
120	else:
NEW 121	raise ValueError("Something wrong with passed x value.")	×
122	rvx[:] = self.term * qx	3✔
123	return rvx, rvy	3✔
124
125	def changetrace(	3✔
126	self, xyzt1, xyzt2, aq, layer, ltype, modellayer, direction, hstepmax
127	):
NEW 128	raise NotImplementedError("changetrace not implemented for this element")	×
129
130	def plot(self, ax=None):	3✔
131	if ax is None:	3✔
NEW 132	_, ax = plt.subplots()	×
133	aq = self.model.aq.find_aquifer_data(self.xls, 0.0)	3✔
134	for ilay in self.layers:	3✔
135	ax.plot(	3✔
136	[self.xls, self.xls],
137	[aq.zaqtop[ilay], aq.zaqbot[ilay]],
138	"k-",
139	)
140
141
142	class DischargeLineSink1D(LineSink1DBase):	3✔
143	r"""Linesink1D with a specified discharge for each layer that the linesink.
144
145	Parameters
146	----------
147	model : Model object
148	model to which the element is added
149	x : float
150	x-coordinate of the linesink
151	tsandq : list of tuples
152	tuples of starting time and specific discharge after starting time
153	res : float
154	resistance of the linesink
155	layers : int, array or list
156	layer (int) or layers (list or array) in which linesink is located
157	label : string or None (default: None)
158	label of the linesink
159
160	Examples
161	--------
162	Example of an infinitely long linesink that pumps with a specific discharge of
163	100 between times 10 and 50, with a specific discharge of 20 between
164	times 50 and 200, and zero speficic discharge after time 200.
165
166	>>> DischargeLineSink1D(ml, tsandq=[(10, 100), (50, 20), (200, 0)])
167	"""
168
169	def __init__(	3✔
170	self, model, xls=0, tsandq=[(0, 1)], res=0, wh="H", layers=0, label=None
171	):
172	super().__init__(	3✔
173	model,
174	xls,
175	tsandbc=tsandq,
176	res=res,
177	wh=wh,
178	layers=layers,
179	type="g",
180	name="DischargeLineSink1D",
181	label=label,
182	)
183
184
185	class LineSink1D(LineSink1DBase, MscreenEquation):	3✔
186	r"""Linesink1D with a specified discharge.
187
188	Parameters
189	----------
190	model : Model object
191	model to which the element is added
192	x : float
193	x-coordinate of the linesink
194	tsandq : list of tuples
195	tuples of starting time and specific discharge after starting time
196	res : float
197	resistance of the linesink
198	layers : int, array or list
199	layer (int) or layers (list or array) in which linesink is located
200	label : string or None (default: None)
201	label of the linesink
202
203	Examples
204	--------
205	Example of an infinitely long linesink that pumps with a specific discharge of
206	100 between times 10 and 50, with a specific discharge of 20 between
207	times 50 and 200, and zero speficic discharge after time 200.
208
209	>>> LineSink1D(ml, tsandq=[(10, 100), (50, 20), (200, 0)])
210	"""
211
212	def __init__(	3✔
213	self,
214	model,
215	xls=0,
216	tsandq=[(0, 1)],
217	res=0,
218	wh="H",
219	vres=0.0,
220	wv=1.0,
221	layers=0,
222	label=None,
223	):
224	super().__init__(	3✔
225	model,
226	xls,
227	tsandbc=tsandq,
228	res=res,
229	wh=wh,
230	layers=layers,
231	type="v",
232	name="LineSink1D",
233	label=label,
234	)
235	self.nunknowns = self.nparam	3✔
236	# Vertical resistance inside line-sink
237	self.vres = np.atleast_1d(vres)	3✔
238	self.wv = wv	3✔
239	if len(self.vres) == 1:	3✔
240	self.vres = self.vres[0] * np.ones(self.nlayers - 1)	3✔
241
242	def initialize(self):	3✔
243	super().initialize()	3✔
244	self.parameters = np.zeros(	3✔
245	(self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
246	)
247	# qv = (hn - hn-1) / vresfac[n - 1]
248	self.vresfac = self.vres / self.wv	3✔
249
250
251	class HeadLineSink1D(LineSink1DBase, HeadEquation):	3✔
252	def __init__(	3✔
253	self,
254	model,
255	xls=0,
256	tsandh=[(0, 1)],
257	res=0,
258	wh="H",
259	layers=0,
260	label=None,
261	):
262	super().__init__(	3✔
263	model,
264	xls,
265	tsandbc=tsandh,
266	res=res,
267	wh=wh,
268	layers=layers,
269	type="v",
270	name="HeadLineSink1D",
271	label=label,
272	)
273	self.nunknowns = self.nparam	3✔
274
275	def initialize(self):	3✔
276	super().initialize()	3✔
277	self.parameters = np.zeros(	3✔
278	(self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
279	)
280	# Needed in solving, solve for a unit head
281	self.pc = self.aq.T[self.layers]	3✔
282
283
284	class HeadDiffLineSink1D(LineSink1DBase, HeadDiffEquation):	3✔
285	def __init__(	3✔
286	self,
287	model,
288	xls=0,
289	layers=0,
290	label=None,
291	aq=None,
292	):
293	super().__init__(	3✔
294	model,
295	xls,
296	tsandbc=[(0, 0)],
297	res=0.0,
298	wh="H",
299	layers=layers,
300	type="z",
301	name="HeadDiffLineSink1D",
302	label=label,
303	aq=aq,
304	inhomelement=True,
305	)
306	self.nunknowns = self.nparam	3✔
307
308	def initialize(self):	3✔
309	super().initialize()	3✔
310	self.xcout = self.xc + self.tiny	3✔
311	self.xcin = self.xc - self.tiny	3✔
312	self.ycout = np.zeros(1)	3✔
313	self.ycin = np.zeros(1)	3✔
314	self.cosout = -np.ones(1)	3✔
315	self.sinout = np.zeros(1)	3✔
316	self.aqout = self.model.aq.find_aquifer_data(self.xcout[0], 0)	3✔
317	self.aqin = self.model.aq.find_aquifer_data(self.xcin[0], 0)	3✔
318
319	self.parameters = np.zeros(	3✔
320	(self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
321	)
322
323
324	class FluxDiffLineSink1D(LineSink1DBase, FluxDiffEquation):	3✔
325	def __init__(self, model, xls=0, layers=0, label=None, aq=None):	3✔
326	super().__init__(	3✔
327	model,
328	xls,
329	tsandbc=[(0, 0)],
330	res=0.0,
331	wh="H",
332	layers=layers,
333	type="z",
334	name="FluxDiffLineSink1D",
335	label=label,
336	aq=aq,
337	inhomelement=True,
338	)
339	self.nunknowns = self.nparam	3✔
340
341	def initialize(self):	3✔
342	super().initialize()	3✔
343	self.xcout = self.xc - self.tiny	3✔
344	self.xcin = self.xc + self.tiny	3✔
345	self.ycout = np.zeros(1)	3✔
346	self.ycin = np.zeros(1)	3✔
347	self.cosout = -np.ones(1)	3✔
348	self.sinout = np.zeros(1)	3✔
349	self.aqout = self.model.aq.find_aquifer_data(self.xcout[0], 0)	3✔
350	self.aqin = self.model.aq.find_aquifer_data(self.xcin[0], 0)	3✔
351	self.parameters = np.zeros(	3✔
352	(self.model.ngvbc, self.nparam, self.model.npval), dtype=complex
353	)

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