Coverage for src/gwtransport/fronttracking/output.py: 83%
283 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""Concentration extraction from front-tracking solutions (V, θ coordinates).
3Every public function in this module takes θ (cumulative flow, m³). Callers
4translate user-facing time t → θ at the API boundary via
5``FrontTrackerState.theta_at_t``.
7Outlet-mass functions use the PDE conservation identity
8``m_out(θ) = m_in(θ) − m_dom(θ)`` (Bear & Cheng 2010, Ch. 3: mass
9conservation for transport with sorption). ``m_dom`` honors historical
10wave activity via ``wave.was_active_at(theta)`` so retrospective queries
11at θ before a collision event correctly attribute c at v_outlet.
13Available functions:
15- :func:`concentration_at_point` - Exact concentration at a single point ``(v, θ)``: the upstream state of the
16 nearest wave face at or downstream of ``v`` among the waves active at ``θ``, the downstream state of the
17 outermost face when nothing lies downstream, and ``0`` in a virgin domain. Exactly on a shock or fan face the
18 two sides are averaged; a contact carries its upstream value at its own position. The ``sorption`` argument is
19 unused — each wave carries its own sorption reference.
21- :func:`compute_breakthrough_curve` - Outlet concentration sampled at every θ of a θ-array, i.e.
22 :func:`concentration_at_point` evaluated at ``v_outlet``. Returns one concentration [mass/volume] per θ; these
23 are point samples of the exact trace, not bin averages.
25- :func:`compute_bin_averaged_concentration_exact` - Outlet concentration averaged over each output θ-bin (a
26 flow-weighted average, since θ is cumulative flow), evaluated through the conservation identity
27 ``c_avg = (Δm_in − Δm_dom) / Δθ`` rather than by integrating the outlet trace, so overlapping fans need no
28 ownership dispatch. A residual inside the floating-point cancellation band is set to zero; a bin more negative
29 than that band warns and is clamped to zero to preserve the ``cout >= 0`` contract.
31- :func:`identify_outlet_segments` - Partition of a θ-interval into the segments over which a single wave
32 controls the outlet, as a list of dicts holding the segment θ-range, its type (``'constant'``,
33 ``'rarefaction'`` or ``'decaying_fan'``), the controlling wave and the concentrations at both segment ends.
34 Crossings extrapolated at or past a wave's deactivation are dropped.
36- :func:`compute_domain_mass` - Total mass ``∫₀^{v_outlet} C_total(v, θ) dv`` held in the domain at one θ [mass],
37 dissolved plus sorbed. The domain is partitioned at the active wave faces; a constant region integrates as
38 ``C_total·Δv`` and a fan region in closed form via :func:`integrate_fan_spatial_exact`.
40- :func:`compute_cumulative_inlet_mass` - Mass injected at the inlet from ``θ = 0`` up to ``theta``,
41 ``∫₀^θ cin dτ`` [mass], exact for the piecewise-constant ``cin`` on ``theta_edges``.
43- :func:`compute_cumulative_outlet_mass` - Mass that has left through the outlet from ``θ = 0`` up to ``theta``
44 [mass], as ``m_in(θ) − m_dom(θ)``; ``0`` for ``θ <= 0``.
46- :func:`compute_total_outlet_mass` - Outlet mass over all θ, from the inlet record alone: the full injected mass
47 ``Σ cin·Δθ`` when the record returns to ``cin[-1] = 0``, and ``+inf`` for a sustained ``cin[-1] > 0`` boundary,
48 which keeps injecting forever.
50- :func:`integrate_fan_exact` - Exact ``∫ c(θ) dθ`` [mass] at a fixed position for any self-similar fan, given
51 the fan apex ``(theta_origin, v_origin)``, via the universal integration-by-parts antiderivative
52 ``F(θ) = c·(θ − θ_origin) − Δv·C_T(c)``. ``c_apex > 0`` clamps the fan at its tail θ and adds the constant
53 plateau beyond it; ``theta_end`` may be ``+inf`` where the fan integral converges, and raises otherwise.
55- :func:`integrate_rarefaction_exact` - :func:`integrate_fan_exact` with the apex and ``c_apex = raref.c_tail``
56 read off a :class:`~gwtransport.fronttracking.waves.RarefactionWave`.
58- :func:`integrate_fan_spatial_exact` - Exact ``∫ C_total(v, θ) dv`` [mass] over a v-segment of a self-similar
59 fan at fixed θ, via the spatial counterpart ``G(u) = C_T(c)·u − κ·c`` with ``κ = θ − θ_origin`` and
60 ``u = v − v_origin``. ``c_apex > 0`` splits off the constant-``C_total(c_apex)`` region near the apex.
62This file is part of gwtransport which is released under AGPL-3.0 license.
63See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
64"""
66import warnings
67from collections.abc import Sequence
68from itertools import pairwise
69from operator import itemgetter
71import numpy as np
72import numpy.typing as npt
74from gwtransport.fronttracking.events import find_outlet_crossing
75from gwtransport.fronttracking.interactions import Face, iter_faces
76from gwtransport.fronttracking.math import (
77 NonlinearSorption,
78 SorptionModel,
79)
80from gwtransport.fronttracking.waves import (
81 CharacteristicWave,
82 DecayingShockWave,
83 Feeder,
84 RarefactionWave,
85 ShockWave,
86 Wave,
87)
89# Numerical tolerance constants
90EPSILON_VELOCITY = 1e-15 # Tolerance for checking if velocity is effectively zero
91EPSILON_VOLUME = 1e-15 # Tolerance for negligible spatial-segment widths Δv [m³]
92EPSILON_POSITION = 1e-15 # Tolerance for shock-face proximity in position v [m³]
93# Multiplier on the eps·max(|m_in|,|m_dom|)/Δθ cancellation scale below which a conservation-form
94# residual is numerical zero. The m_dom fan-integral sum accumulates ~few·10² ULP over a
95# multi-pulse record (measured max ratio ~460); 65536 clears that with wide margin while staying
96# ~7 orders below any real breakthrough concentration (O(cin)).
97FP_CANCELLATION_CLAMP = 65536.0
100def _reader_faces(waves: Sequence[Wave], theta: float) -> list[tuple[float, Face]]:
101 """``(position, face)`` for every active wave face at ``theta`` (reader sweep basis).
103 Fan boundary lines are included: they mark where a downstream-opening fan ends and its
104 plateau begins (a ``decay_side='right'`` decaying/doubly-fed shock has its fan on the
105 downstream side, so its head boundary must partition the domain — the feeder clamp alone
106 only reaches the apex-side plateau).
107 """
108 out: list[tuple[float, Face]] = []
109 for wave in waves:
110 if not wave.was_active_at(theta):
111 continue
112 for face in iter_faces(wave, theta):
113 pos = face.position(theta)
114 if pos is not None:
115 out.append((pos, face))
116 return out
119def concentration_at_point(
120 v: float,
121 theta: float,
122 waves: Sequence[Wave],
123 sorption: SorptionModel, # noqa: ARG001
124) -> float:
125 """Compute concentration at point (v, θ) with exact analytical value.
127 The function works entirely in (V, θ) coordinates: public callers must
128 translate user-facing time t → θ at the API boundary (e.g., via
129 ``FrontTrackerState.theta_at_t``).
131 Parameters
132 ----------
133 v : float
134 Position [m³].
135 theta : float
136 Cumulative flow [m³].
137 waves : list of Wave
138 All waves in the simulation (active and inactive).
139 sorption : SorptionModel
140 Sorption model (unused — kept for API symmetry; wave methods carry
141 their own sorption reference).
143 Returns
144 -------
145 concentration : float
146 Concentration at point (v, θ) [mass/volume].
148 Notes
149 -----
150 **Nearest-downstream-face sweep.** ``C(v, θ)`` is the *left* (upstream) state of the
151 nearest face strictly downstream of ``v`` among the waves active at ``θ``; if no face
152 is downstream, the *right* state of the outermost face; ``0`` in a virgin domain. A
153 face's feeder is a constant or a bounded self-similar fan (clamped to its extent, so the
154 plateau beyond a fan reads correctly). Because every face crossing fires a solver event,
155 the front list is interaction-consistent and each region has a single owner — no
156 ownership heuristic is needed. This is exact for a single front (one owner everywhere)
157 and for genuinely interacting multi-front inputs alike.
158 """
159 faces = _reader_faces(waves, theta)
160 # Nearest face at or downstream of v (within FP): its left (upstream) feeder gives the
161 # state just behind it, which is the state at v. A face strictly upstream of v does not
162 # describe v. If nothing is at/downstream, the outermost face's right (downstream) feeder
163 # holds — the state ahead of every front.
164 downstream = [(p, f) for (p, f) in faces if p >= v - EPSILON_POSITION]
165 if downstream:
166 pos, face = min(downstream, key=itemgetter(0))
167 # Exactly on a shock/fan face: return the average of the two sides (the infinitesimally
168 # thin discontinuity convention). A contact (characteristic) carries its value on the
169 # upstream side, so it returns that (the behind value) at its own position.
170 if abs(pos - v) < EPSILON_POSITION and face.role != "contact":
171 return 0.5 * (face.left.value(v, theta) + face.right.value(v, theta))
172 return face.left.value(v, theta)
173 if faces:
174 _p, face = max(faces, key=itemgetter(0))
175 return face.right.value(v, theta)
176 return 0.0
179def compute_breakthrough_curve(
180 theta_array: npt.NDArray[np.floating],
181 v_outlet: float,
182 waves: Sequence[Wave],
183 sorption: SorptionModel,
184) -> npt.NDArray[np.floating]:
185 """Concentration at the outlet evaluated over a θ-array (breakthrough curve).
187 Parameters
188 ----------
189 theta_array : array-like
190 Cumulative-flow points at which to query the outlet concentration [m³].
191 Must be sorted in ascending order. Callers translate from user-facing
192 time via ``FrontTrackerState.theta_at_t`` before passing.
193 v_outlet : float
194 Outlet position [m³].
195 waves : list of Wave
196 All waves in the simulation.
197 sorption : SorptionModel
198 Sorption model.
200 Returns
201 -------
202 c_out : numpy.ndarray
203 Concentration at ``v_outlet`` for each θ in ``theta_array`` [mass/volume].
205 See Also
206 --------
207 concentration_at_point : Point-wise concentration
208 compute_bin_averaged_concentration_exact : Bin-averaged concentrations
210 Examples
211 --------
212 .. disable_try_examples
214 ::
216 theta_array = np.linspace(0.0, tracker.state.theta_edges[-1], 1000)
217 c_out = compute_breakthrough_curve(
218 theta_array, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption
219 )
220 """
221 theta_arr = np.asarray(theta_array, dtype=float)
222 c_out = np.zeros(len(theta_arr))
223 for i, theta in enumerate(theta_arr):
224 c_out[i] = concentration_at_point(v_outlet, float(theta), waves, sorption)
225 return c_out
228def identify_outlet_segments(
229 theta_start: float,
230 theta_end: float,
231 v_outlet: float,
232 waves: Sequence[Wave],
233 sorption: SorptionModel,
234) -> list[dict]:
235 """Identify which waves control outlet concentration in θ-interval [theta_start, theta_end].
237 Finds all wave crossing events at the outlet and constructs segments where
238 concentration is constant or varying (rarefaction). All times are expressed
239 as cumulative flow θ [m³].
241 Parameters
242 ----------
243 theta_start : float
244 Start of cumulative-flow interval [m³].
245 theta_end : float
246 End of cumulative-flow interval [m³].
247 v_outlet : float
248 Outlet position [m³].
249 waves : list of Wave
250 All waves in the simulation.
251 sorption : SorptionModel
252 Sorption model.
254 Returns
255 -------
256 segments : list of dict
257 List of segment dictionaries, each containing:
259 - 'theta_start' : float
260 Segment start θ [m³]
261 - 'theta_end' : float
262 Segment end θ [m³]
263 - 'type' : str
264 ``'constant'``, ``'rarefaction'``, or ``'decaying_fan'``.
265 ``'decaying_fan'`` is owned by a :class:`~gwtransport.fronttracking.waves.DecayingShockWave` after
266 its head crosses ``v_outlet``; c at ``v_outlet`` then follows the
267 wave's self-similar fan profile.
268 - 'concentration' : float
269 For constant segments
270 - 'wave' : Wave
271 For rarefaction and decaying_fan segments
272 - 'c_start' : float
273 Concentration at segment start
274 - 'c_end' : float
275 Concentration at segment end
277 Notes
278 -----
279 Segments are constructed by:
281 1. Finding all wave crossing events at the outlet for θ in [theta_start, theta_end].
282 2. Sorting events by θ.
283 3. Creating constant-concentration segments between events.
284 4. Handling rarefaction and decaying-fan profiles with θ-varying concentration.
286 The segments completely partition the interval [theta_start, theta_end].
288 Every crossing is clamped to ``theta_cross < theta_deactivation``
289 (matching ``was_active_at`` semantics): a crossing extrapolated past a
290 wave's deactivation is an artifact — after a collision the front belongs
291 to the successor wave, whose own crossing covers the outlet.
292 """
293 # Find all waves that cross outlet in this θ-range
294 outlet_events: list[dict] = []
296 # Track rarefactions / decaying shocks that already contain the outlet at
297 # theta_start (no crossing event in [theta_start, theta_end]).
298 active_rarefactions_at_start: list[RarefactionWave | DecayingShockWave] = []
300 for wave in waves:
301 # Lifetime filter, matching ``Wave.was_active_at``: skip waves that were never
302 # activated (``is_active=False`` with no recorded deactivation) or whose lifetime
303 # ended before the window starts. ``is_active`` alone is the end-of-simulation
304 # state and would also skip waves deactivated after an in-window crossing.
305 never_active = not wave.is_active and wave.theta_deactivation == float("inf")
306 if never_active or wave.theta_deactivation <= theta_start:
307 continue
309 if isinstance(wave, DecayingShockWave):
310 # The wave's outlet crossing arrival behaves like a rarefaction head
311 # arrival: before arrival, v_outlet is downstream (c=c_fixed for
312 # decay_side='left'); after arrival, v_outlet is inside the fan
313 # whose c follows the self-similar profile and asymptotes to the
314 # fan's tail concentration.
315 theta_cross = wave.outlet_crossing_theta(v_outlet)
316 if theta_cross is None or theta_cross >= wave.theta_deactivation:
317 # Crossings at or after deactivation are spurious extrapolations.
318 continue
319 if theta_cross <= theta_start:
320 # Outlet already inside the fan at theta_start.
321 active_rarefactions_at_start.append(wave)
322 elif theta_cross <= theta_end:
323 # c_after is the fan c just past arrival (the decay-side c at
324 # the arrival θ). theta_cross > wave.theta_start by construction
325 # (outlet_crossing_theta enforces v_outlet > v_start), so
326 # c_decay_at_theta does not return None.
327 c_after = wave.c_decay_at_theta(theta_cross)
328 outlet_events.append({
329 "theta": theta_cross,
330 "wave": wave,
331 "boundary": "head",
332 "c_after": c_after,
333 })
334 continue
336 # For rarefactions, detect both head and tail crossings
337 if isinstance(wave, RarefactionWave):
338 # Check if outlet is already inside this rarefaction at theta_start
339 if wave.contains_point(v_outlet, theta_start):
340 active_rarefactions_at_start.append(wave)
341 # Detect when the tail crosses during [theta_start, theta_end]
342 tail_speed = wave.tail_speed()
343 if tail_speed > EPSILON_VELOCITY:
344 theta_cross = wave.theta_start + (v_outlet - wave.v_start) / tail_speed
345 if theta_start < theta_cross <= theta_end and theta_cross < wave.theta_deactivation:
346 outlet_events.append({
347 "theta": theta_cross,
348 "wave": wave,
349 "boundary": "tail",
350 "c_after": wave.c_tail,
351 })
352 continue
354 # Head crossing
355 head_speed = wave.head_speed()
356 if head_speed > EPSILON_VELOCITY and wave.v_start < v_outlet:
357 theta_cross = wave.theta_start + (v_outlet - wave.v_start) / head_speed
358 if theta_start <= theta_cross <= theta_end and theta_cross < wave.theta_deactivation:
359 outlet_events.append({
360 "theta": theta_cross,
361 "wave": wave,
362 "boundary": "head",
363 "c_after": wave.c_head,
364 })
366 # Tail crossing
367 tail_speed = wave.tail_speed()
368 if tail_speed > EPSILON_VELOCITY and wave.v_start < v_outlet:
369 theta_cross = wave.theta_start + (v_outlet - wave.v_start) / tail_speed
370 if theta_start <= theta_cross <= theta_end and theta_cross < wave.theta_deactivation:
371 outlet_events.append({
372 "theta": theta_cross,
373 "wave": wave,
374 "boundary": "tail",
375 "c_after": wave.c_tail,
376 })
377 else:
378 # Characteristics and shocks
379 theta_cross = find_outlet_crossing(wave, v_outlet, theta_start)
381 if theta_cross is not None and theta_start <= theta_cross <= theta_end:
382 if isinstance(wave, CharacteristicWave):
383 c_after = wave.concentration
384 elif isinstance(wave, ShockWave):
385 # After shock passes outlet, outlet sees left (upstream) state
386 c_after = wave.c_left
387 else:
388 c_after = 0.0
390 outlet_events.append({"theta": theta_cross, "wave": wave, "boundary": None, "c_after": c_after})
392 # Sort events by θ
393 outlet_events.sort(key=itemgetter("theta"))
395 # Create segments between events
396 segments: list[dict] = []
397 current_theta = theta_start
398 current_c = concentration_at_point(v_outlet, theta_start, waves, sorption)
400 # Handle case where we start inside a rarefaction or decaying-shock fan.
401 # Multi-fan overlap: pick the newest (largest ``theta_start``) — matches
402 # ``concentration_at_point`` and ``compute_domain_mass`` dispatch.
403 if active_rarefactions_at_start:
404 raref = max(active_rarefactions_at_start, key=lambda w: w.theta_start)
406 if isinstance(raref, RarefactionWave):
407 # Find when tail crosses (if it does)
408 tail_cross_theta = None
409 for event in outlet_events:
410 if event["wave"] is raref and event["boundary"] == "tail" and event["theta"] > theta_start:
411 tail_cross_theta = event["theta"]
412 break
414 raref_end = min(tail_cross_theta or theta_end, theta_end)
415 c_end = raref.c_tail if tail_cross_theta and tail_cross_theta <= theta_end else None
417 segments.append({
418 "theta_start": theta_start,
419 "theta_end": raref_end,
420 "type": "rarefaction",
421 "wave": raref,
422 "c_start": current_c,
423 "c_end": c_end,
424 })
425 else:
426 # DecayingShockWave fan extends to θ=+∞ (or asymptotes to c_fixed
427 # for n>1 with c_min); treat the whole [theta_start, theta_end]
428 # as one decaying-fan segment.
429 raref_end = theta_end
430 c_end = concentration_at_point(v_outlet, theta_end, waves, sorption)
432 segments.append({
433 "theta_start": theta_start,
434 "theta_end": raref_end,
435 "type": "decaying_fan",
436 "wave": raref,
437 "c_start": current_c,
438 "c_end": c_end,
439 })
441 current_theta = raref_end
442 current_c = (
443 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption) if raref_end < theta_end else current_c
444 )
446 for event in outlet_events:
447 # Skip events that fall inside an already-emitted (typically rarefaction)
448 # segment. ``concentration_at_point`` lets active rarefactions "win"
449 # over a behind-shock c_left; the segment list must reflect the same
450 # convention to avoid double-counting.
451 if event["theta"] < current_theta:
452 continue
454 if isinstance(event["wave"], RarefactionWave) and event["boundary"] == "head":
455 if event["theta"] > current_theta:
456 segments.append({
457 "theta_start": current_theta,
458 "theta_end": event["theta"],
459 "type": "constant",
460 "concentration": current_c,
461 "c_start": current_c,
462 "c_end": current_c,
463 })
465 raref = event["wave"]
466 tail_cross_theta = None
467 for later_event in outlet_events:
468 if (
469 later_event["wave"] is raref
470 and later_event["boundary"] == "tail"
471 and later_event["theta"] > event["theta"]
472 ):
473 tail_cross_theta = later_event["theta"]
474 break
476 raref_end = min(tail_cross_theta or theta_end, theta_end)
478 segments.append({
479 "theta_start": event["theta"],
480 "theta_end": raref_end,
481 "type": "rarefaction",
482 "wave": raref,
483 "c_start": raref.c_head,
484 "c_end": raref.c_tail if tail_cross_theta and tail_cross_theta <= theta_end else None,
485 })
487 current_theta = raref_end
488 current_c = (
489 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption)
490 if raref_end < theta_end
491 else current_c
492 )
493 elif isinstance(event["wave"], DecayingShockWave) and event["boundary"] == "head":
494 if event["theta"] > current_theta:
495 segments.append({
496 "theta_start": current_theta,
497 "theta_end": event["theta"],
498 "type": "constant",
499 "concentration": current_c,
500 "c_start": current_c,
501 "c_end": current_c,
502 })
504 decaying = event["wave"]
505 # The decaying_fan segment ends at the next outlet-crossing event
506 # (if any falls in the window) or at theta_end. Without this split,
507 # multi-DSW pulses would have the first DSW's fan swallow every
508 # later wave's arrival.
509 seg_end = theta_end
510 for later_event in outlet_events:
511 if later_event["theta"] > event["theta"] and later_event["theta"] <= theta_end:
512 seg_end = later_event["theta"]
513 break
514 c_end_val = concentration_at_point(v_outlet, seg_end, waves, sorption)
516 segments.append({
517 "theta_start": event["theta"],
518 "theta_end": seg_end,
519 "type": "decaying_fan",
520 "wave": decaying,
521 "c_start": event["c_after"],
522 "c_end": c_end_val,
523 })
525 current_theta = seg_end
526 current_c = c_end_val
527 else:
528 if event["theta"] > current_theta:
529 segments.append({
530 "theta_start": current_theta,
531 "theta_end": event["theta"],
532 "type": "constant",
533 "concentration": current_c,
534 "c_start": current_c,
535 "c_end": current_c,
536 })
538 current_theta = event["theta"]
539 current_c = event["c_after"]
541 # Final segment
542 if theta_end > current_theta:
543 segments.append({
544 "theta_start": current_theta,
545 "theta_end": theta_end,
546 "type": "constant",
547 "concentration": current_c,
548 "c_start": current_c,
549 "c_end": current_c,
550 })
552 return segments
555def integrate_rarefaction_exact(
556 raref: RarefactionWave, v_outlet: float, theta_start: float, theta_end: float, sorption: SorptionModel
557) -> float:
558 """Exact θ-integral ``∫ c(θ) dθ`` of a rarefaction at the outlet.
560 Convenience wrapper over :func:`integrate_fan_exact` that pulls the fan
561 apex from ``raref.theta_start, raref.v_start``. Returns the mass-like
562 quantity ``∫ c dθ`` (= ``∫ c·flow dt`` in time coordinates).
564 Parameters
565 ----------
566 raref : RarefactionWave
567 Rarefaction wave controlling the outlet.
568 v_outlet : float
569 Outlet position [m³].
570 theta_start, theta_end : float
571 Integration range in cumulative flow [m³]. Either can be ``±np.inf``.
572 sorption : SorptionModel
573 Sorption model (any NonlinearSorption subclass).
575 Returns
576 -------
577 integral : float
578 ``∫ c(θ) dθ`` [mass — i.e. concentration × volume].
579 """
580 return integrate_fan_exact(
581 raref.theta_start, raref.v_start, v_outlet, theta_start, theta_end, sorption, c_apex=raref.c_tail
582 )
585def integrate_fan_exact(
586 theta_origin: float,
587 v_origin: float,
588 v_outlet: float,
589 theta_start: float,
590 theta_end: float,
591 sorption: SorptionModel,
592 c_apex: float = 0.0,
593) -> float:
594 """Exact θ-integral ``∫ c(θ) dθ`` for any self-similar fan at the outlet.
596 Decoupled from the wave object so the same closed-form math applies to
597 both :class:`~gwtransport.fronttracking.waves.RarefactionWave` (apex = ``theta_start, v_start``) and
598 :class:`~gwtransport.fronttracking.waves.DecayingShockWave` (apex = ``theta_origin, v_origin``).
600 Parameters
601 ----------
602 theta_origin, v_origin : float
603 Cumulative flow and position at the fan's apex [m³].
604 v_outlet : float
605 Outlet position [m³].
606 theta_start, theta_end : float
607 Integration range in cumulative flow [m³]. ``theta_end`` may be
608 ``+np.inf``; ``theta_start`` must be finite.
609 sorption : SorptionModel
610 Sorption model (any NonlinearSorption subclass).
611 c_apex : float, optional
612 Concentration on the constant side at the fan apex. For
613 ``RarefactionWave`` this is ``raref.c_tail``; for
614 ``DecayingShockWave`` (decay_side='left') this is ``wave.c_fan_tail``
615 (the plateau the outlet holds once the fan's far edge sweeps by).
616 For ``c_apex > 0`` the fan formula extrapolates past the physical
617 fan range; the integration is clamped at ``θ_tail`` (where
618 ``c(θ_tail) = c_apex``) and the constant-c_apex region beyond
619 contributes ``c_apex · (theta_end − θ_tail)``. Default 0.0
620 preserves the c=0 apex behavior for canonical c_R=0 fans.
622 Returns
623 -------
624 float
625 Mass-like quantity ``∫ c(θ) dθ`` [mass — concentration × volume].
627 Raises
628 ------
629 TypeError
630 If the sorption model does not support exact fan integration.
631 """
632 # Every NonlinearSorption uses one universal IBP antiderivative, which evaluates the fan
633 # kernel k = R·c − C_T at the segment endpoints via c_and_total_from_retardation.
634 if isinstance(sorption, NonlinearSorption):
635 return _integrate_fan_exact_universal(
636 theta_origin, v_origin, v_outlet, theta_start, theta_end, sorption, c_apex
637 )
639 msg = f"Exact fan integration not supported for {type(sorption).__name__}"
640 raise TypeError(msg)
643def _integrate_fan_exact_universal(
644 theta_origin: float,
645 v_origin: float,
646 v_outlet: float,
647 theta_start: float,
648 theta_end: float,
649 sorption: NonlinearSorption,
650 c_apex: float = 0.0,
651) -> float:
652 r"""Exact θ-integral ``∫ c(θ) dθ`` via the universal IBP antiderivative.
654 For any ``NonlinearSorption`` (with ``R = dC_T/dC``), integration by
655 parts on the self-similar fan ``R(c(θ)) = (θ − θ_origin)/Δv`` gives the
656 closed-form antiderivative
658 .. math::
659 F(\\theta) = c(\\theta)\\,(\\theta - \\theta_{\\rm origin})
660 - \\Delta v \\cdot C_T(c(\\theta)).
662 The derivation uses ``∫ c\\,d\\theta = c·(\\theta-\\theta_0) − ∫(\\theta-\\theta_0)·dc
663 = c·(\\theta-\\theta_0) − \\Delta v · ∫ R(c)\\,dc = c·(\\theta-\\theta_0) − \\Delta v · C_T(c)``,
664 where the last equality is the definition of ``C_T`` as the antiderivative
665 of ``R`` (``R = dC_T/dC``).
667 This formula is exact for any sorption (Brooks-Corey, van Genuchten-Mualem,
668 Freundlich, Langmuir). The only sorption-specific call is
669 ``sorption.concentration_from_retardation`` at the two endpoints — for
670 Brooks-Corey this is closed form; for van Genuchten-Mualem it is one
671 ``brentq`` call per endpoint. No quadrature, no integration loop.
673 Convergence at θ → ∞ for ``c_apex = 0``: for any monotone sorption with
674 ``R(0) = ∞`` (BC, vG, Freundlich n > 1), ``c(∞) = 0`` and ``c·θ → 0``
675 faster than ``Δv·C_T → 0`` (verified termwise from the closed-form
676 asymptotic ``c ~ R^{-α}`` for some ``α > 1``), so ``F(∞) = 0``.
678 For ``c_apex > 0`` the fan formula extrapolates to ``c < c_apex`` past
679 ``θ_tail = θ_origin + Δv·R(c_apex)``; clamp the fan portion at
680 ``θ_tail`` and add ``c_apex·(θ_end − θ_tail)`` for any ``θ_end > θ_tail``.
681 """
682 delta_v = v_outlet - v_origin
683 if delta_v <= 0 or theta_end <= theta_start:
684 return 0.0
686 if c_apex > 0.0:
687 theta_tail = theta_origin + delta_v * float(sorption.retardation(c_apex))
688 theta_end_fan = min(theta_end, theta_tail)
689 else:
690 theta_tail = float("inf")
691 theta_end_fan = theta_end
693 # For a c_apex=0 fan the upper bound may be +∞. The integral converges only when
694 # c → 0 as R → ∞ (so base·c → 0). Freundlich n<1 (c → ∞ as R → ∞) diverges — reject it
695 # explicitly via the fan_converges_at_infinity guard.
696 if theta_end_fan == float("inf") and not sorption.fan_converges_at_infinity():
697 msg = "Fan integral diverges at θ=+∞ for this sorption (e.g. Freundlich n<1); pass a finite theta_end"
698 raise ValueError(msg)
700 def antiderivative(theta: float) -> float:
701 if theta == float("inf"):
702 # F(∞) = 0 for any sorption with c → 0 as R → ∞ (guarded above for divergent cases).
703 return 0.0
704 base = (theta - theta_origin) / delta_v
705 if base <= 0.0:
706 return 0.0
707 c, ct = sorption.c_and_total_from_retardation(base)
708 return c * (theta - theta_origin) - delta_v * ct
710 fan_integral = antiderivative(theta_end_fan) - antiderivative(theta_start)
711 constant_contrib = c_apex * max(theta_end - theta_tail, 0.0) if c_apex > 0.0 else 0.0
712 return fan_integral + constant_contrib
715def compute_bin_averaged_concentration_exact(
716 theta_bin_edges: npt.NDArray[np.floating],
717 v_outlet: float,
718 waves: Sequence[Wave],
719 sorption: SorptionModel,
720 *,
721 cin: npt.ArrayLike,
722 theta_edges_inlet: npt.NDArray[np.floating],
723) -> npt.NDArray[np.floating]:
724 """θ-bin-averaged outlet concentration.
726 For each θ-bin ``[θ_i, θ_{i+1}]``::
728 C_avg = (1 / Δθ) · ∫_{θ_i}^{θ_{i+1}} C(v_outlet, θ) dθ
730 evaluated through the conservation-law identity
731 ``C_avg = (Δm_in − Δm_dom) / Δθ`` per bin — analytical and explicit, no
732 outlet-side fan dispatch, so multi-DSW and n<1 mirror geometries are
733 handled by the same single path.
735 Parameters
736 ----------
737 theta_bin_edges : array-like
738 Cumulative-flow OUTPUT bin edges [m³] (where C_avg is reported).
739 Length N+1 for N bins. Callers translate t-bin edges with
740 ``state.theta_at_t``.
741 v_outlet : float
742 Outlet position [m³].
743 waves : list of Wave
744 All waves from front tracking simulation.
745 sorption : SorptionModel
746 Sorption model.
747 cin : array-like (kw-only)
748 Inlet concentration per inlet θ-bin.
749 theta_edges_inlet : ndarray (kw-only)
750 θ bin edges of the INLET (``state.theta_edges``), length
751 ``len(cin) + 1``.
753 Returns
754 -------
755 c_avg : numpy.ndarray
756 Bin-averaged outlet concentrations [mass/volume]. Length N.
758 Raises
759 ------
760 ValueError
761 If any output θ-bin has non-positive width.
763 See Also
764 --------
765 concentration_at_point : Point-wise concentration
766 compute_breakthrough_curve : Breakthrough curve
767 compute_cumulative_outlet_mass : Cumulative outlet mass via conservation
768 """
769 theta_edges_out = np.asarray(theta_bin_edges, dtype=float)
770 dtheta_out = np.diff(theta_edges_out)
772 if np.any(dtheta_out <= 0):
773 bad = int(np.argmin(dtheta_out))
774 msg = (
775 f"Invalid θ-bin: theta_bin_edges[{bad}]={theta_edges_out[bad]} >= "
776 f"theta_bin_edges[{bad + 1}]={theta_edges_out[bad + 1]}"
777 )
778 raise ValueError(msg)
780 # Conservation form: c_avg = Δm_out/Δθ where m_out = m_in − m_dom.
781 # m_in(θ) = ∫₀^θ cin dτ is the piecewise-LINEAR (in θ) integral of the piecewise-constant
782 # cin, so evaluate it at every output edge in O(N+M) from the inlet-bin cumulative sums
783 # plus the partial bin containing θ — instead of the dense N_out×M_in clip-and-matmul
784 # (≈245× slower, a 128 MB temporary at N=M=4000). Edges below te_in[0] contribute 0
785 # (nothing is injected before the record starts, even if it starts mid-window); edges at
786 # or past te_in[-1] saturate at the total. This mirrors ``compute_cumulative_inlet_mass``'s
787 # clip exactly. m_dom stays a per-θ spatial geometry loop.
788 te_in = np.asarray(theta_edges_inlet, dtype=float)
789 cin_arr = np.asarray(cin, dtype=float)
790 cum_in = np.concatenate([[0.0], np.cumsum(cin_arr * np.diff(te_in))])
791 idx = np.clip(np.searchsorted(te_in, theta_edges_out, side="right") - 1, 0, len(cin_arr) - 1)
792 m_in_at_edges = cum_in[idx] + cin_arr[idx] * (theta_edges_out - te_in[idx])
793 m_in_at_edges = np.where(theta_edges_out < te_in[0], 0.0, m_in_at_edges)
794 m_in_at_edges = np.where(theta_edges_out >= te_in[-1], cum_in[-1], m_in_at_edges)
795 m_dom_at_edges = np.array([
796 compute_domain_mass(theta=float(theta_e), v_outlet=v_outlet, waves=waves, sorption=sorption)
797 for theta_e in theta_edges_out
798 ])
799 # ``compute_cumulative_outlet_mass`` short-circuits to 0 for θ ≤ 0;
800 # replicate that clamp so non-positive output edges contribute no mass.
801 m_out_at_edges = np.where(theta_edges_out <= 0.0, 0.0, m_in_at_edges - m_dom_at_edges)
802 result = np.diff(m_out_at_edges) / dtheta_out
803 # FP-noise clamp scaled to the m_in − m_dom CANCELLATION magnitude: each cumulative-mass
804 # edge carries ~eps·max(|m_in|,|m_dom|) rounding (the m_dom fan-integral sum accumulates
805 # several hundred ULP over a multi-pulse record), amplified by 1/Δθ. Keying the band off
806 # the OUTPUT concentration instead would collapse it to ~0 before breakthrough, where
807 # ~1e-11 cancellation dust is normal. Residuals within this band are numerical zero:
808 # clamp and stay silent.
809 mass_scale = np.maximum(np.abs(m_in_at_edges), np.abs(m_dom_at_edges))
810 # Band = the larger of the FP-cancellation floor and a 1e-6 relative transient floor.
811 # The latter absorbs the benign ``c_min``-floor artifact at a fan-entry (a decaying
812 # side born at ``c ≈ 1e-12`` rather than exactly 0 perturbs the near-apex fan integral
813 # for one θ-sample); a genuine over-count is O(a pulse's mass), far above the band.
814 eps_band = (
815 max(FP_CANCELLATION_CLAMP * np.finfo(float).eps, 1e-6)
816 * np.maximum(np.maximum(mass_scale[:-1], mass_scale[1:]), 1.0)
817 / dtheta_out
818 )
819 result = np.where(np.abs(result) < eps_band, 0.0, result)
820 # A residual MORE negative than the FP band is a genuine conservation-form violation with two
821 # possible drivers, and BOTH can be present at once. Attribute per offending bin, not off the
822 # single last edge: a bin whose right edge passes ``theta_edges_inlet[-1]`` sees m_in saturate
823 # while the wave list keeps evolving (out-of-window); a fully in-window offending bin is a real
824 # wave-model over-count (e.g. overlapping non-interacting waves from an oscillating inlet, see
825 # ``compute_domain_mass`` Notes). Report every cause that actually produced a negative bin.
826 # Clamp to 0 to preserve the ``cout >= 0`` API contract either way.
827 neg_mask = result < -eps_band
828 if np.any(neg_mask):
829 worst = float(np.min(result))
830 te_in_last = float(te_in[-1])
831 out_of_window = neg_mask & (theta_edges_out[1:] > te_in_last)
832 causes = []
833 if np.any(out_of_window):
834 causes.append(
835 f"output θ-bin edges exceeding theta_edges_inlet[-1]={te_in_last:.3f} (m_in "
836 "saturates at the last injected mass while the wave list keeps evolving — extend "
837 "cin with trailing zeros to cover the output range, or restrict output bins to the "
838 "inlet window)"
839 )
840 if np.any(neg_mask & ~out_of_window):
841 causes.append(
842 "the domain-mass integral growing faster than the inlet-mass integral within the "
843 "inlet window (an m_dom over-count — e.g. overlapping non-interacting waves from a "
844 "continuous/oscillating inlet — or a wave-list / cin inconsistency)"
845 )
846 warnings.warn(
847 f"compute_bin_averaged_concentration_exact produced a concentration as negative as "
848 f"{worst:.3e}, beyond the FP-cancellation band; likely cause(s): {'; '.join(causes)}.",
849 UserWarning,
850 stacklevel=2,
851 )
852 return np.maximum(result, 0.0)
855def compute_domain_mass(
856 theta: float,
857 v_outlet: float,
858 waves: Sequence[Wave],
859 sorption: SorptionModel,
860) -> float:
861 """
862 Compute total mass in domain [0, v_outlet] at cumulative flow θ.
864 Integrates concentration over space::
866 M(θ) = ∫₀^v_outlet C_total(v, θ) dv
868 Exact analytical formulas for every wave type: constant regions
869 (``C_total · Δv``), RarefactionWave fan interiors and DecayingShockWave fan
870 interiors (closed-form via :func:`integrate_fan_spatial_exact`).
872 Parameters
873 ----------
874 theta : float
875 Cumulative flow at which to compute domain mass [m³].
876 v_outlet : float
877 Outlet position (domain extent) [m³].
878 waves : list of Wave
879 All waves in the simulation.
880 sorption : SorptionModel
881 Sorption model.
883 Returns
884 -------
885 mass : float
886 Total mass in domain [mass]. Closed-form analytical to machine precision.
888 See Also
889 --------
890 compute_cumulative_inlet_mass : Cumulative inlet mass
891 compute_cumulative_outlet_mass : Cumulative outlet mass
892 concentration_at_point : Point-wise concentration
893 integrate_fan_spatial_exact : Closed-form fan spatial integral
895 Examples
896 --------
897 .. disable_try_examples
899 ::
901 mass = compute_domain_mass(
902 theta=2500.0, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption
903 )
904 mass >= 0.0
905 """
906 # Partition [0, v_outlet] at the active faces; each segment has a single owner (the
907 # left feeder of its nearest downstream face). A constant owner integrates as
908 # C_total·Δv; a fan owner integrates in closed form via integrate_fan_spatial_exact,
909 # clamped at the fan's near-apex (larger-retardation) bound. Because the wave list is
910 # interaction-consistent (every face crossing fires a solver event), this single-owner
911 # sweep is exact for single-front and genuinely interacting multi-front layouts alike.
912 faces = _reader_faces(waves, theta)
913 boundaries = {0.0, v_outlet}
914 boundaries.update(p for (p, _f) in faces if 0.0 < p < v_outlet)
915 positions = sorted(boundaries)
917 total_mass = 0.0
918 for v_start, v_end in pairwise(positions):
919 dv = v_end - v_start
920 if dv < EPSILON_VOLUME:
921 continue
922 v_mid = 0.5 * (v_start + v_end)
923 feeder = _owner_feeder(faces, v_mid)
925 if feeder is None or feeder.is_const:
926 c = feeder.value(v_mid, theta) if feeder is not None else 0.0
927 total_mass += float(sorption.total_concentration(c)) * dv
928 else:
929 total_mass += integrate_fan_spatial_exact(
930 feeder.theta_apex,
931 feeder.v_apex,
932 v_start,
933 v_end,
934 theta,
935 sorption,
936 c_apex=_near_apex_bound(feeder, sorption),
937 )
939 return float(total_mass)
942def _owner_feeder(faces: list[tuple[float, Face]], v: float) -> Feeder | None:
943 """Return the feeder controlling the state at ``v``: nearest-downstream face's left, else outermost right."""
944 downstream = [(p, f) for (p, f) in faces if p > v + EPSILON_POSITION]
945 if downstream:
946 _p, face = min(downstream, key=itemgetter(0))
947 return face.left
948 if faces:
949 _p, face = max(faces, key=itemgetter(0))
950 return face.right
951 return None
954def _near_apex_bound(feeder: Feeder, sorption: SorptionModel) -> float:
955 """Return the fan boundary concentration with the larger retardation (the near-apex / plateau side).
957 This is the ``c_apex`` clamp for ``integrate_fan_spatial_exact``: the value the fan holds
958 at (and beyond) its apex. Monotonicity-agnostic — the low-c bound for R-decreasing
959 isotherms, the high-c bound for the Freundlich ``n<1`` mirror where R increases with c.
960 """
961 r_a = float(sorption.retardation(feeder.c_a))
962 r_b = float(sorption.retardation(feeder.c_b))
963 return feeder.c_a if r_a >= r_b else feeder.c_b
966def integrate_fan_spatial_exact(
967 theta_origin: float,
968 v_origin: float,
969 v_start: float,
970 v_end: float,
971 theta: float,
972 sorption: SorptionModel,
973 c_apex: float = 0.0,
974) -> float:
975 """Exact spatial integral ``∫ C_total(v, θ) dv`` for any self-similar fan.
977 Decoupled from the wave object so the same closed-form math applies to
978 :class:`~gwtransport.fronttracking.waves.RarefactionWave` (apex = ``theta_start, v_start``) and
979 :class:`~gwtransport.fronttracking.waves.DecayingShockWave` (apex = ``theta_origin, v_origin``).
981 In (V, θ) the self-similar fan satisfies ``R(C) = (θ - θ_origin)/(v - v_origin)``;
982 define ``kappa = θ - θ_origin`` and ``u = v - v_origin``. The dissolved and
983 sorbed contributions reduce to power-law forms in ``u`` that admit closed
984 forms via incomplete beta functions (Freundlich) or elementary sqrt
985 operations (Langmuir).
987 Parameters
988 ----------
989 theta_origin, v_origin : float
990 Cumulative flow and position at the fan's apex [m³].
991 v_start, v_end : float
992 Integration range in v [m³].
993 theta : float
994 Cumulative flow at which to evaluate [m³].
995 sorption : SorptionModel
996 Sorption model (any NonlinearSorption subclass).
997 c_apex : float, optional
998 Concentration on the constant side at the fan apex (typically the
999 parent rarefaction's ``c_tail`` or the DSW's ``c_fixed`` for
1000 ``decay_side='left'``). For ``c_apex > 0`` the fan formula is
1001 unphysical for ``u < u_tail = kappa / R(c_apex)``; the integration
1002 is split into a constant-C_total(c_apex) region for
1003 ``u ∈ [u_start, u_tail]`` plus the fan integral for
1004 ``u ∈ [u_tail, u_end]``. Default 0.0 preserves the c=0 apex
1005 behavior for canonical c_R=0 rarefactions.
1007 Returns
1008 -------
1009 float
1010 Mass in the segment ``[v_start, v_end]``.
1012 Raises
1013 ------
1014 TypeError
1015 If the sorption model does not support exact spatial integration.
1016 """
1017 if theta <= theta_origin:
1018 return 0.0
1020 kappa = theta - theta_origin
1021 u_start = v_start - v_origin
1022 u_end = v_end - v_origin
1024 # The fan only exists for v > v_origin; clip u_start at 0 (the apex
1025 # contributes nothing to the integral since c(v_origin)=0 for n>1 and
1026 # the Beta-function form handles the lower-bound singularity for n<1).
1027 # If the whole segment is upstream of the apex, return 0.
1028 if u_end <= 0:
1029 return 0.0
1030 if u_start < 0:
1031 u_start = 0.0
1033 # Split off the constant-c_apex region near the apex for c_apex > 0.
1034 # The fan formula is only valid for u ≥ u_tail = kappa / R(c_apex);
1035 # below u_tail, c is clamped to c_apex (the parent's tail / DSW's fixed
1036 # concentration). Spatial counterpart of the temporal θ_tail clamp in
1037 # _integrate_fan_exact_universal.
1038 constant_contrib = 0.0
1039 if c_apex > 0.0:
1040 u_tail = kappa / float(sorption.retardation(c_apex))
1041 if u_start < u_tail:
1042 c_total_apex = float(sorption.total_concentration(c_apex))
1043 constant_contrib = c_total_apex * (min(u_end, u_tail) - u_start)
1044 u_start = u_tail
1045 if u_end <= u_start:
1046 return constant_contrib
1048 # One universal IBP antiderivative for every NonlinearSorption (see the rationale in
1049 # ``integrate_fan_exact``, the temporal counterpart).
1050 if isinstance(sorption, NonlinearSorption):
1051 return constant_contrib + _integrate_rarefaction_spatial_universal(sorption, kappa, u_start, u_end)
1053 msg = f"Exact spatial fan integration not supported for {type(sorption).__name__}"
1054 raise TypeError(msg)
1057def _integrate_rarefaction_spatial_universal(
1058 sorption: NonlinearSorption,
1059 kappa: float,
1060 u_start: float,
1061 u_end: float,
1062) -> float:
1063 r"""Exact spatial integral ``∫ C_T(c(u)) du`` via the universal IBP antiderivative.
1065 For any ``NonlinearSorption`` with ``R = dC_T/dC``, integration by parts
1066 on the self-similar fan ``R(c(u)) = κ/u`` gives the closed-form
1067 antiderivative
1069 .. math::
1070 G(u) = C_T(c(u))\\cdot u - \\kappa\\cdot c(u).
1072 Derivation: ``∫ C_T\\,du = C_T·u − ∫ u\\,dC_T = C_T·u − ∫ (κ/R)·R\\,dc =
1073 C_T·u − κ·c`` (the second equality uses ``u = κ/R`` and the third uses
1074 ``dC_T = R\\,dc``).
1076 Sorption-specific calls limited to ``concentration_from_retardation`` and
1077 ``total_concentration`` at the two endpoints. For Brooks-Corey both are
1078 closed form; for van Genuchten-Mualem ``concentration_from_retardation`` is
1079 one ``brentq`` per endpoint and ``total_concentration`` is also one
1080 ``brentq`` per endpoint (chained internally). No quadrature.
1082 At the apex (``u → 0``, ``R → ∞``) ``c → 0`` and ``C_T → 0`` so ``G(0) = 0``;
1083 a segment whose lower bound is at (or below, for a bounded-``R`` sorption like
1084 Langmuir where ``c = 0`` for ``u`` below the fan tail) the apex contributes
1085 ``G(u_end) − 0``.
1086 """
1087 if u_end <= 0.0 or u_end <= u_start or kappa <= 0.0:
1088 return 0.0
1089 if u_start <= 0.0:
1090 g_start = 0.0 # G(0) = C_T(0)·0 − κ·0 = 0 (apex)
1091 else:
1092 c_start, ct_start = sorption.c_and_total_from_retardation(kappa / u_start)
1093 g_start = ct_start * u_start - kappa * c_start
1094 c_end, ct_end = sorption.c_and_total_from_retardation(kappa / u_end)
1095 g_end = ct_end * u_end - kappa * c_end
1096 return g_end - g_start
1099def compute_cumulative_inlet_mass(
1100 theta: float,
1101 cin: npt.ArrayLike,
1102 theta_edges: npt.ArrayLike,
1103) -> float:
1104 """Cumulative inlet mass entering the domain from θ=0 to ``theta``.
1106 In cumulative-flow coordinates ``M_in(θ) = ∫₀^θ cin(τ) dτ``; for
1107 piecewise-constant ``cin`` this is exact under summation over θ-bin
1108 widths.
1110 Parameters
1111 ----------
1112 theta : float
1113 Cumulative flow up to which to integrate [m³].
1114 cin : array-like
1115 Inlet concentration per θ-bin [mass/volume].
1116 theta_edges : array-like
1117 θ bin edges [m³], length ``len(cin) + 1``.
1119 Returns
1120 -------
1121 mass_in : float
1122 Cumulative inlet mass [mass].
1124 Examples
1125 --------
1126 .. disable_try_examples
1128 ::
1130 mass_in = compute_cumulative_inlet_mass(
1131 theta=5000.0, cin=cin, theta_edges=theta_edges
1132 )
1133 mass_in >= 0.0
1134 """
1135 te = np.asarray(theta_edges, dtype=float)
1136 widths = np.clip(theta - te[:-1], 0.0, np.diff(te))
1137 return float(np.sum(np.asarray(cin, dtype=float) * widths))
1140def compute_cumulative_outlet_mass(
1141 theta: float,
1142 v_outlet: float,
1143 waves: Sequence[Wave],
1144 sorption: SorptionModel,
1145 *,
1146 cin: npt.ArrayLike,
1147 theta_edges: npt.NDArray[np.floating],
1148) -> float:
1149 """Cumulative mass exiting through the outlet from θ=0 to ``theta``.
1151 Computed analytically via the conservation-law identity::
1153 m_out(θ) = m_in(θ) − m_dom(θ)
1155 derived from integrating the PDE ``∂_θ C_T + ∂_V c = 0`` over the spatial
1156 domain ``[0, v_outlet]`` (Bear & Cheng 2010, Ch. 3: mass conservation
1157 for advection with sorption). This sidesteps the multi-fan dispatch problem
1158 that the outlet-segment integration faces when several DSWs cover
1159 v_outlet simultaneously — every term on the right is purely spatial or a
1160 closed-form inlet sum, no ownership priority needed.
1162 Parameters
1163 ----------
1164 theta : float
1165 Cumulative flow up to which to integrate [m³].
1166 v_outlet : float
1167 Outlet position [m³].
1168 waves : list of Wave
1169 All waves in the simulation.
1170 sorption : SorptionModel
1171 Sorption model.
1172 cin : array-like (kw-only)
1173 Inlet concentration per θ-bin [mass/volume].
1174 theta_edges : ndarray (kw-only)
1175 θ bin edges [m³], length ``len(cin) + 1``.
1177 Returns
1178 -------
1179 mass_out : float
1180 Cumulative outlet mass [mass].
1182 Examples
1183 --------
1184 .. disable_try_examples
1186 ::
1188 mass_out = compute_cumulative_outlet_mass(
1189 theta=5000.0,
1190 v_outlet=500.0,
1191 waves=tracker.state.waves,
1192 sorption=sorption,
1193 cin=cin,
1194 theta_edges=tracker.state.theta_edges,
1195 )
1196 mass_out >= 0.0
1197 """
1198 if theta <= 0.0:
1199 return 0.0
1200 m_in = compute_cumulative_inlet_mass(theta=theta, cin=cin, theta_edges=theta_edges)
1201 m_dom = compute_domain_mass(theta=theta, v_outlet=v_outlet, waves=waves, sorption=sorption)
1202 return m_in - m_dom
1205def compute_total_outlet_mass(
1206 *,
1207 cin: npt.ArrayLike,
1208 theta_edges: npt.NDArray[np.floating],
1209) -> float:
1210 """Total outlet mass over θ → ∞ (finite only for a returning-to-zero pulse).
1212 The final inlet value ``c_∞ = cin[-1]`` is the sustained boundary state as θ → ∞:
1214 - For ``c_∞ = 0`` (canonical c_R=0 pulse): injection ceases, the domain empties, and
1215 every injected mass unit eventually exits — ``m_out_total = m_in_total`` (the finite
1216 record integral ``Σ cin·Δθ``). The wave list is not needed.
1217 - For ``c_∞ > 0`` (sustained ambient): the inlet keeps injecting ``c_∞`` forever, so the
1218 cumulative outlet mass grows without bound — return ``+inf``. Pairing the FINITE record
1219 integral with the infinite-time steady-state fill ``C_T(c_∞)·v_outlet`` is not an option:
1220 that difference goes negative whenever ``m_in_total < C_T(c_∞)·v_outlet``.
1222 Parameters
1223 ----------
1224 cin : array-like (kw-only)
1225 Inlet concentration per θ-bin [mass/volume].
1226 theta_edges : ndarray (kw-only)
1227 θ bin edges [m³], length ``len(cin) + 1``.
1229 Returns
1230 -------
1231 float
1232 ``m_in_total`` for ``cin[-1] = 0``; ``+inf`` for ``cin[-1] > 0``.
1234 See Also
1235 --------
1236 compute_cumulative_outlet_mass : Cumulative outlet mass up to a finite θ (use this for
1237 a sustained ``c_∞ > 0`` boundary, where the θ → ∞ total is unbounded).
1238 compute_domain_mass : Spatial integral of C_total in the aquifer
1239 """
1240 cin_arr = np.asarray(cin, dtype=float)
1241 # An empty cin is a malformed (no-bin) input by the cin/theta_edges contract;
1242 # let cin_arr[-1] raise IndexError rather than masking it as c_inf=0.
1243 if float(cin_arr[-1]) > 0.0:
1244 return float("inf")
1245 te = np.asarray(theta_edges, dtype=float)
1246 return float(np.sum(cin_arr * np.diff(te)))