Coverage for src/gwtransport/recharge.py: 100%
173 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""
2Recharge-Driven Transport for Aquifers with Areal Recharge.
4Concentration at extraction has two sources. 1) Water infiltrates and is
5transported through an aquifer with constant thickness to extraction. 2) During
6transport, rainfall is mixed instantaneously over the height of the aquifer. In
7an unbounded aquifer all extracted water originates as recharge. Transport is
8advective with linear sorption; there is no microdispersion, molecular diffusion,
9or macrodispersion. Only forward modeling is supported. No assumption is made
10about whether the flow is radial or orthogonal. Two conceptual models share one
11entry point:
13- **Unbounded aquifer** (``aquifer_pore_volume=None``): all extracted water
14 originates as recharge. The residence-time distribution is exponential with
15 mean ``retardation_factor * aquifer_pore_depth / N`` — independent of the
16 pumping rate, hydraulic conductivity, capture-zone size, and planform shape
17 (Haitjema, 1995). In the cumulative-recharge clock
18 ``u(t) = ∫ N dt / (retardation_factor * aquifer_pore_depth)`` (pore volumes
19 flushed) the model is the stationary unit filter ``dC/du = cin_recharge - C``,
20 which this module integrates in closed form per bin. No flow rate is needed.
22- **Bounded aquifer** (``aquifer_pore_volume`` set): the aquifer extent is
23 capped at pore volume ``aquifer_pore_volume`` (strip area
24 ``aquifer_pore_volume / aquifer_pore_depth``). Water with concentration
25 ``cin`` enters at the upstream side at rate ``q_b = flow - N * area``
26 whenever extraction exceeds the rainfall on the strip. When rainfall exceeds
27 extraction (``q_b < 0``) the surplus flows out across the upstream boundary
28 and is lost; the outside has no memory, so when extraction later dominates
29 again the inflow carries the current ``cin``. The exact solution is the
30 unbounded exponential kernel acting on ``cin_recharge``, truncated at the
31 boundary-entry time of the extracted water, with the residual tail weight
32 placed as an atom on ``cin`` at the entry time. With zero recharge this
33 reduces exactly to single-pore-volume piston flow
34 (:func:`gwtransport.advection.infiltration_to_extraction`); with the
35 boundary never feeding the well it reduces exactly to the unbounded model.
37Available functions:
39- :func:`recharge_to_extraction` - Compute the extracted concentration from the recharge
40 concentration ``cin_recharge`` and, in the bounded model, the upstream-boundary concentration
41 ``cin`` and the extraction rate ``flow``. Forward only: there is no extraction-to-recharge
42 inverse. The result is one value per ``cout_tedges`` bin, a flow-weighted bin average in the
43 bounded model and a recharge-weighted bin average in the unbounded model, NaN outside the input
44 time range and in bins whose weight is zero. All terms are closed-form, so the output is exact
45 to machine precision for bin-constant inputs.
47References
48----------
49Haitjema, H.M. (1995). On the residence time distribution in idealized
50groundwatersheds. Journal of Hydrology, 172(1-4), 127-146.
51https://doi.org/10.1016/0022-1694(95)02732-5
53This file is part of gwtransport which is released under AGPL-3.0 license.
54See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
55"""
57import numpy as np
58import numpy.typing as npt
59import pandas as pd
61from gwtransport._time import tedges_to_days
62from gwtransport._validation import (
63 _validate_no_nan,
64 _validate_non_negative_array,
65 _validate_positive_scalar,
66 _validate_retardation_factor,
67 _validate_tedges_parity,
68)
70# Kernel weights older than this many pore-volume flushes are below one ulp of the
71# row sum (e^-60 ~ 9e-27); truncating them keeps the gathers banded without
72# changing any double-precision result.
73_KERNEL_CUTOFF = 60.0
76def recharge_to_extraction(
77 *,
78 cin: npt.ArrayLike | None = None,
79 cin_recharge: npt.ArrayLike,
80 flow: npt.ArrayLike | None = None,
81 recharge: npt.ArrayLike,
82 tedges: pd.DatetimeIndex,
83 cout_tedges: pd.DatetimeIndex,
84 aquifer_pore_volume: float | None = None,
85 aquifer_pore_depth: float,
86 retardation_factor: float = 1.0,
87) -> npt.NDArray[np.floating]:
88 """Compute the concentration of extracted water under uniform areal recharge.
90 Unbounded model (``aquifer_pore_volume=None``): exponential residence-time
91 distribution with mean ``retardation_factor * aquifer_pore_depth / N``
92 (Haitjema, 1995), exact for bin-constant inputs. Bounded model
93 (``aquifer_pore_volume`` set, together with ``cin`` and ``flow``): the
94 exponential kernel is truncated at the upstream-boundary entry time and the
95 residual weight is an atom on ``cin``; water pushed out across the boundary
96 during rainfall surplus is lost.
98 Parameters
99 ----------
100 cin : array-like, optional
101 Concentration of the water entering at the upstream side of the
102 bounded aquifer [concentration units]. Required when
103 ``aquifer_pore_volume`` is set; must be None otherwise.
104 cin_recharge : array-like
105 Concentration of the recharge water entering via the surface
106 [concentration units]. Length must equal ``len(tedges) - 1``; constant
107 over each interval ``[tedges[i], tedges[i+1])``.
108 flow : array-like, optional
109 Extraction rate [m3/day]. Required when ``aquifer_pore_volume`` is
110 set; must be None otherwise, because the unbounded model is
111 independent of the pumping rate (see Notes). Must be non-negative and
112 NaN-free.
113 recharge : array-like
114 Areal recharge rate N [m/day; same length unit as
115 ``aquifer_pore_depth``]. Length must equal ``len(tedges) - 1``.
116 Must be non-negative and NaN-free.
117 tedges : pandas.DatetimeIndex
118 Time bin edges for the input series.
119 cout_tedges : pandas.DatetimeIndex
120 Time bin edges for the output series. Bins not fully inside the
121 ``tedges`` range return NaN.
122 aquifer_pore_volume : float, optional
123 Pore volume of the bounded aquifer [m3]. The strip area between the
124 upstream boundary and the well is
125 ``aquifer_pore_volume / aquifer_pore_depth``. Default None (unbounded).
126 aquifer_pore_depth : float
127 Pore volume per unit surface area: porosity times saturated thickness
128 [m]. The only static aquifer parameter of the unbounded model.
129 retardation_factor : float, optional
130 Compound retardation factor (>= 1.0), by default 1.0. Dilates the
131 solute clock; mixing fractions are unaffected.
133 Returns
134 -------
135 numpy.ndarray
136 Extracted concentration per ``cout_tedges`` bin, length
137 ``len(cout_tedges) - 1``. Flow-weighted bin average (bounded model) or
138 recharge-weighted bin average (unbounded model). NaN for bins outside
139 the input time range, for zero-recharge bins (unbounded), and for
140 zero-extraction bins (bounded).
142 Raises
143 ------
144 ValueError
145 If array lengths do not match the bin-edge pattern, inputs contain NaN
146 or negative values, physical parameters are out of range, or only part
147 of the bounded-model triple (``cin``, ``flow``,
148 ``aquifer_pore_volume``) is provided.
150 See Also
151 --------
152 gwtransport.advection.infiltration_to_extraction : Zero-recharge limit of the bounded model.
153 gwtransport.deposition.deposition_to_extraction : Distributed source along the flow path.
154 :ref:`concept-residence-time` : Background on residence times.
155 :ref:`concept-transport-equation` : Flow-weighted averaging approach.
157 Notes
158 -----
159 The unbounded model needs no flow rate because the capture zone
160 self-adjusts: the well always draws exactly its pumping rate from
161 recharge, over a capture area ``flow / N``. Pumping harder widens the
162 capture area proportionally, leaving the age composition of the extracted
163 water -- set by the ratio of pore storage per unit area
164 (``aquifer_pore_depth``) to recharge per unit area (``N``) -- unchanged,
165 so the flow rate cancels exactly (Haitjema, 1995). In the bounded model
166 the area is fixed by ``aquifer_pore_volume`` instead of adjusting to the
167 well, so the flow rate no longer cancels and must be given.
169 Spin-up follows the ``"constant"`` policy: all inputs are treated as
170 constant at their first values before ``tedges[0]``. For the bounded model
171 this is the steady concentration profile
172 ``C(V) = cr0 + (cin0 - cr0) * (V_R - apv) / (V_R - V)`` when the boundary
173 feeds the well (``q_b(0) > 0``, ``V_R = flow[0] * aquifer_pore_depth /
174 recharge[0]``), and the uniform profile ``cin_recharge[0]`` otherwise.
176 Under constant inputs with ``flow > N * area`` the extracted water is the
177 mass-balance mixture ``cin_recharge + (cin - cin_recharge) * q_b / flow``:
178 an exponential residence-time density carrying the recharge fraction plus
179 a piston atom of mass ``q_b / flow`` at the boundary-to-well travel time.
181 The exponential kernel lives on the dimensionless clock ``u`` and is
182 parameter-free; the pumping rate enters the bounded model only through the
183 boundary-entry times. All formulas are closed-form (exp/log of bin-local
184 quantities), exact to machine precision for bin-constant inputs.
186 References
187 ----------
188 Haitjema, H.M. (1995). On the residence time distribution in idealized
189 groundwatersheds. Journal of Hydrology, 172(1-4), 127-146.
190 https://doi.org/10.1016/0022-1694(95)02732-5
192 Examples
193 --------
194 >>> import numpy as np
195 >>> import pandas as pd
196 >>> from gwtransport.recharge import recharge_to_extraction
197 >>> tedges = pd.date_range("2020-01-01", periods=11, freq="D")
198 >>> cout = recharge_to_extraction(
199 ... cin_recharge=np.full(10, 2.5),
200 ... recharge=np.full(10, 0.002),
201 ... tedges=tedges,
202 ... cout_tedges=tedges[3:],
203 ... aquifer_pore_depth=3.0,
204 ... )
205 >>> np.allclose(cout, 2.5)
206 True
207 """
208 tedges, cout_tedges = pd.DatetimeIndex(tedges), pd.DatetimeIndex(cout_tedges)
209 cr = np.asarray(cin_recharge, dtype=float)
210 rech = np.asarray(recharge, dtype=float)
211 bounded = aquifer_pore_volume is not None
213 if (cin is None) != (flow is None) or (cin is None) == bounded:
214 msg = "cin, flow, and aquifer_pore_volume must be provided together (bounded model) or all be None"
215 raise ValueError(msg)
216 _validate_tedges_parity(tedges, cr, tedges_name="tedges", values_name="cin_recharge")
217 _validate_tedges_parity(tedges, rech, tedges_name="tedges", values_name="recharge")
218 _validate_no_nan(cr, name="cin_recharge")
219 _validate_no_nan(rech, name="recharge")
220 _validate_non_negative_array(rech, name="recharge")
221 _validate_positive_scalar(aquifer_pore_depth, name="aquifer_pore_depth")
222 _validate_retardation_factor(retardation_factor)
224 t = tedges_to_days(tedges)
225 tq = tedges_to_days(cout_tedges, ref=tedges[0])
226 dt = np.diff(t)
227 k = rech / (retardation_factor * aquifer_pore_depth)
228 u = np.concatenate([[0.0], np.cumsum(k * dt)])
229 covered = (tq[:-1] >= t[0]) & (tq[1:] <= t[-1])
231 if not bounded:
232 # The unbounded aquifer is the no-boundary special case: an infinite
233 # pore volume puts the boundary beyond reach (pure kernel and spin-up
234 # terms), and a synthetic flow proportional to recharge makes the
235 # flow-weighted bin average the recharge-weighted one.
236 return _bounded_average(t=t, dt=dt, u=u, k=k, q=k, cr=cr, cb=cr, apv=np.inf, tq=tq, covered=covered)
238 cb = np.asarray(cin, dtype=float)
239 q = np.asarray(flow, dtype=float)
240 _validate_tedges_parity(tedges, cb, tedges_name="tedges", values_name="cin")
241 _validate_tedges_parity(tedges, q, tedges_name="tedges", values_name="flow")
242 _validate_no_nan(cb, name="cin")
243 _validate_no_nan(q, name="flow")
244 _validate_non_negative_array(q, name="flow")
245 _validate_positive_scalar(aquifer_pore_volume, name="aquifer_pore_volume")
246 # The solute clock divides flow by the retardation factor; weighting ratios are unaffected.
247 return _bounded_average(
248 t=t,
249 dt=dt,
250 u=u,
251 k=k,
252 q=q / retardation_factor,
253 cr=cr,
254 cb=cb,
255 apv=float(aquifer_pore_volume),
256 tq=tq,
257 covered=covered,
258 )
261def _arrival_times(*, t, dt, k, q, apv):
262 """Forward arrival time at the well of parcels released at ``(t[j], V=apv)``.
264 Trajectories obey ``dV/dt = -(q - k V)`` with the per-bin closed form
265 ``V(s) = V_R + (V - V_R) e^{k s}``. Returns NaN for parcels that are
266 expelled across the boundary (lost) or have not arrived by ``t[-1]``.
267 These arrivals are exactly the output times where the boundary-entry bin
268 of the extracted water changes, including the pre-record transition. The
269 bin loop over the live front is O(n^2) in the number of input bins.
271 Returns
272 -------
273 ndarray
274 Arrival time in days per release edge ``t[0..n-1]``; NaN if the parcel
275 never reaches the well within the record.
276 """
277 n = len(dt)
278 arrivals = np.full(n, np.nan)
279 pos = np.full(n, np.nan)
280 with np.errstate(over="ignore", invalid="ignore"):
281 for i in range(n):
282 pos[i] = apv
283 idx = np.nonzero(np.isfinite(pos[: i + 1]))[0]
284 vi = pos[idx]
285 if k[i] > 0:
286 v_r = q[i] / k[i]
287 v_end = v_r + (vi - v_r) * np.exp(k[i] * dt[i])
288 hit = v_end <= 0.0
289 s_hit = np.log(v_r / (v_r - vi[hit])) / k[i]
290 else:
291 v_end = vi - q[i] * dt[i]
292 hit = v_end <= 0.0
293 s_hit = vi[hit] / q[i] if q[i] > 0 else vi[hit][:0] # q == 0 cannot hit the well
294 arrivals[idx[hit]] = t[i] + s_hit
295 # Resolved (arrived) parcels and parcels at or beyond the boundary
296 # (expelled, or parked exactly on a stagnant boundary -- the next
297 # edge's release duplicates a parked parcel) drop out of the front.
298 pos[idx] = np.where(hit | (v_end >= apv), np.nan, v_end)
299 return arrivals
302def _backward_entries(*, queries, t, dt, k, q, apv, v_start=0.0, edge_side="right"):
303 """Trace backward characteristics from ``(query, V=v_start)`` to the boundary or to ``t[0]``.
305 Returns ``(s, v0)``: the boundary-entry time ``s`` (NaN if the water
306 predates the record) and the landing position ``v0`` at ``t[0]`` for
307 pre-record water (NaN otherwise). Uses the numerically local per-bin form
308 ``V_a = V_R + (V_b - V_R) e^{-k seg}``; differencing global accumulators
309 instead loses precision catastrophically once ``u >> 1``.
311 ``v_start=apv`` with ``edge_side="left"`` traces the grazing trajectory
312 that touches the boundary exactly at an (on-edge) query time backward into
313 the preceding bin, yielding the left-branch limit of the entry-time map at
314 that arrival (its release time); for queries not preceded by outflow the
315 walk exits immediately and returns the query time itself.
317 Returns
318 -------
319 tuple of ndarray
320 ``(s, v0)`` per query: boundary-entry time in days (NaN for pre-record
321 water) and landing position at ``t[0]`` (NaN for entered water).
322 """
323 n = len(dt)
324 nq = len(queries)
325 m = np.clip(np.searchsorted(t, queries, side=edge_side) - 1, 0, n - 1)
326 s_out = np.full(nq, np.nan)
327 v0 = np.full(nq, np.nan)
328 pos = np.full(nq, float(v_start))
329 open_ = np.ones(nq, dtype=bool)
330 with np.errstate(invalid="ignore", divide="ignore"): # stagnant boundary (vi == v_r == apv)
331 for i in range(n - 1, -1, -1):
332 sel = np.nonzero(open_ & (m >= i))[0]
333 if sel.size == 0:
334 continue
335 starts_here = m[sel] == i
336 seg = np.where(starts_here, queries[sel] - t[i], dt[i])
337 t_hi = np.where(starts_here, queries[sel], t[i + 1])
338 vi = pos[sel]
339 if k[i] > 0:
340 v_r = q[i] / k[i]
341 va = v_r + (vi - v_r) * np.exp(-k[i] * seg)
342 ent = va >= apv
343 back = -np.log((apv - v_r) / (vi[ent] - v_r)) / k[i]
344 back = np.where(np.isfinite(back), back, 0.0)
345 else:
346 va = vi + q[i] * seg
347 # The q > 0 guard keeps a parcel parked exactly on the boundary
348 # through a fully stagnant bin (no flow, no recharge) walking into
349 # earlier bins: nothing moves and nothing is lost there.
350 ent = (va >= apv) & (q[i] > 0.0)
351 back = (apv - vi[ent]) / q[i] if q[i] > 0 else vi[ent][:0]
352 s_ent = np.clip(t_hi[ent] - back, t[i], t_hi[ent])
353 s_out[sel[ent]] = s_ent
354 open_[sel[ent]] = False
355 pos[sel[~ent]] = va[~ent]
356 v0[open_] = pos[open_]
357 return s_out, v0
360def _bounded_average(*, t, dt, u, k, q, cr, cb, apv, tq, covered):
361 """Flow-weighted bin averages of the bounded model, exact via piece integration.
363 The output interval is split at: input edges, output edges, and the
364 arrival times of boundary parcels released at the input edges. Between
365 consecutive breakpoints the entry time stays within one source bin and
366 every term integrates in closed form. The boundary atom integrates to
367 ``cin_js * q_b,js * (s2 - s1)`` via the change of variables
368 ``q e^{-u(t)} dt = e^{-u(s)} q_b(s) ds`` along the entry map; the kernel
369 terms reduce to bin-local exponentials times
370 ``I2 = ∫ q e^{-(u(t)-u(t2))} dt``; pre-record water carries the steady
371 spin-up profile evaluated at the landing position ``v0 = G(t)``, whose
372 flow-weighted integral is ``cr0``-linear plus a closed-form logarithm.
374 Returns
375 -------
376 ndarray
377 Flow-weighted average per output bin; NaN where undefined.
378 """
379 n = len(dt)
380 n_out = len(tq) - 1
381 qb = q - k * np.where(k > 0, apv, 0.0) # the where avoids 0 * inf in the unbounded (apv = inf) routing
383 arrivals = _arrival_times(t=t, dt=dt, k=k, q=q, apv=apv)
384 bp = np.unique(np.concatenate([t, tq[(tq >= t[0]) & (tq <= t[-1])], arrivals[np.isfinite(arrivals)]]))
385 mids = 0.5 * (bp[:-1] + bp[1:])
386 s_bp, v0_bp = _backward_entries(queries=bp, t=t, dt=dt, k=k, q=q, apv=apv)
387 s_mid, _ = _backward_entries(queries=mids, t=t, dt=dt, k=k, q=q, apv=apv)
389 # The entry-time map s*(t) is discontinuous at the arrival of a parcel
390 # released at an edge where boundary inflow resumes after an expulsion or
391 # stagnation episode (the skipped span is the lost window), and the walk at
392 # exactly that arrival resolves to one branch by floating-point luck. Both
393 # one-sided limits are computed robustly instead: the right limit is the
394 # release edge t[g] itself; the left limit is the release time of the
395 # grazing trajectory, traced backward from (t[g], V=apv).
396 g_idx = np.nonzero(np.isfinite(arrivals))[0]
397 s_left = np.full(len(arrivals), np.nan)
398 v0_left = np.full(len(arrivals), np.nan)
399 g_pos = g_idx[g_idx >= 1] # the g = 0 arrival grazes the boundary exactly at t[0] (landing position apv)
400 if g_pos.size:
401 s_left[g_pos], v0_left[g_pos] = _backward_entries(
402 queries=t[g_pos], t=t, dt=dt, k=k, q=q, apv=apv, v_start=apv, edge_side="left"
403 )
404 av = arrivals[g_idx]
405 order = np.argsort(av)
406 av, ae = av[order], g_idx[order]
408 def arrival_edge(x):
409 if av.size == 0:
410 return np.full(len(x), -1)
411 pos = np.minimum(np.searchsorted(av, x), av.size - 1)
412 return np.where(av[pos] == x, ae[pos], -1)
414 t1, t2 = bp[:-1], bp[1:]
415 span = t2 - t1
416 m = np.clip(np.searchsorted(t, mids, side="right") - 1, 0, n - 1)
417 kc = np.searchsorted(tq, mids, side="right") - 1
418 in_out = (kc >= 0) & (kc < n_out)
419 kc = np.clip(kc, 0, n_out - 1)
421 pre = np.isnan(s_mid)
422 js = np.clip(np.searchsorted(t, np.where(pre, t[0], s_mid), side="right") - 1, 0, n - 1)
423 # Entry times at piece endpoints: one-sided limits at arrival breakpoints,
424 # the walked values elsewhere, NaN (grazing/pre-record endpoints) falling
425 # back to the entry-bin edge; the final clip into the piece's own entry bin
426 # [t[js], t[js+1]] is a roundoff guard.
427 e1, e2 = arrival_edge(t1), arrival_edge(t2)
428 s1_raw = np.where(e1 >= 0, t[np.maximum(e1, 0)], s_bp[:-1])
429 s2_raw = np.where(e2 >= 0, s_left[np.maximum(e2, 0)], s_bp[1:])
430 s_lo, s_hi = t[js], t[js + 1]
431 s1 = np.clip(np.where(np.isnan(s1_raw), s_lo, s1_raw), s_lo, s_hi)
432 s2 = np.clip(np.where(np.isnan(s2_raw), s_hi, s2_raw), s_lo, s_hi)
433 ds = np.maximum(s2 - s1, 0.0)
435 u1t = u[m] + k[m] * (t1 - t[m])
436 u2t = u[m] + k[m] * (t2 - t[m])
437 kpos = k[m] > 0
438 vol = q[m] * span
439 fac = np.where(kpos, q[m] / np.where(kpos, k[m], 1.0), vol)
441 def piece_integral(x):
442 """Integrate ``q e^{x - u(t)}`` over each piece in closed form (callers keep x <= u(t1) bounded).
444 Returns
445 -------
446 ndarray
447 One integral per piece.
448 """
449 return np.where(kpos, fac * (np.exp(x - u1t) - np.exp(x - u2t)), fac * np.exp(x - u2t))
451 # Kernel mass: full-bin sum for j in [lo, m), then the current-bin and
452 # entry-bin partial corrections (telescopes to vol when cr == cb == 1).
453 lo = np.where(pre, 0, js)
454 # The per-bin weight magnitude scales with exp(u[col+1] - u1t) (u1t <= u2t is
455 # the piece's smallest clock value), so the cutoff must key to u1t: keying it
456 # to u2t drops still-significant bins when one input bin flushes >_KERNEL_CUTOFF
457 # pore volumes.
458 lo_eff = np.maximum(lo, np.clip(np.searchsorted(u, u1t - _KERNEL_CUTOFF, side="right") - 1, 0, None))
459 w_max = int((m - lo_eff).max(initial=0))
460 ker = np.zeros(len(mids))
461 if w_max > 0:
462 cols = lo_eff[:, None] + np.arange(w_max)[None, :]
463 valid = cols < m[:, None]
464 colsc = np.clip(cols, 0, n - 1)
465 # Each bin weight is an adjacent difference of exp values at its two u
466 # edges; consecutive bins share an edge, so evaluating exp once per edge
467 # over the extended [lo_eff, m] edge array halves the exp count.
468 edges = np.clip(lo_eff[:, None] + np.arange(w_max + 1)[None, :], 0, n)
469 e_lo = np.exp(u[edges] - u1t[:, None])
470 e_hi = np.exp(u[edges] - u2t[:, None])
471 w_lo = e_lo[:, 1:] - e_lo[:, :-1]
472 w_hi = e_hi[:, 1:] - e_hi[:, :-1]
473 wgt = np.where(kpos[:, None], w_lo - w_hi, w_hi)
474 ker = np.einsum("pw,pw->p", np.where(valid, wgt, 0.0), cr[colsc])
475 mass = cr[m] * vol - cr[m] * piece_integral(u[m]) + ker * fac
477 # Boundary atom (entered) or steady-profile spin-up atom (pre-record). The
478 # where keeps -inf * 0 (unbounded routing, pre-record pieces) out of the
479 # discarded branch.
480 entered_atom = (cb[js] - cr[js]) * np.where(pre, 0.0, qb[js]) * ds + cr[js] * piece_integral(u[js])
481 qb0 = qb[0]
482 if qb0 > 0 and k[0] > 0:
483 # Landing positions v0 = G(t) at the piece endpoints. At the
484 # pre-record transition the one-sided limit is the landing position of
485 # the grazing continuation from (t[g], apv): apv itself only when the
486 # transition parcel is the t[0] release; when earlier releases were
487 # expelled (delayed transition) the grazing path lands INSIDE at t0.
488 v0_arr = np.where(np.isnan(v0_left), apv, v0_left)
489 v0_1 = np.where(e1 >= 0, v0_arr[np.maximum(e1, 0)], v0_bp[:-1])
490 v0_2 = np.where(e2 >= 0, v0_arr[np.maximum(e2, 0)], v0_bp[1:])
491 v0_1 = np.where(pre, np.where(np.isnan(v0_1), apv, v0_1), 0.0)
492 v0_2 = np.where(pre, np.where(np.isnan(v0_2), apv, v0_2), 0.0)
493 v_r0 = q[0] / k[0]
494 ic_log = np.log((v_r0 - v0_1) / (v_r0 - v0_2))
495 ic_atom = cr[0] * piece_integral(0.0) + (cb[0] - cr[0]) * (v_r0 - apv) * ic_log
496 elif qb0 > 0:
497 ic_atom = cb[0] * piece_integral(0.0) # piston pre-record: domain full of boundary water
498 else:
499 ic_atom = cr[0] * piece_integral(0.0) # boundary never fed the domain before t0
500 mass += np.where(pre, ic_atom, entered_atom)
502 take = in_out & (span > 0)
503 masses = np.zeros(n_out)
504 vols = np.zeros(n_out)
505 np.add.at(masses, kc[take], mass[take])
506 np.add.at(vols, kc[take], vol[take])
507 with np.errstate(invalid="ignore", divide="ignore"):
508 return np.where(covered & (vols > 0), masses / vols, np.nan)