Coverage for src/gwtransport/fronttracking/output.py: 86%
412 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +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``.
7Functions
8---------
9::
11 concentration_at_point(v, theta, waves, sorption)
12 compute_breakthrough_curve(theta_array, v_outlet, waves, sorption)
13 compute_bin_averaged_concentration_exact(theta_bin_edges, v_outlet, waves, sorption, *, cin=None, theta_edges_inlet=None)
14 compute_domain_mass(theta, v_outlet, waves, sorption)
15 compute_cumulative_inlet_mass(theta, cin, theta_edges)
16 compute_cumulative_outlet_mass(theta, v_outlet, waves, sorption, *, cin, theta_edges)
17 compute_total_outlet_mass(v_outlet, sorption, *, cin, theta_edges) -> float
19Outlet-mass functions use the PDE conservation identity
20``m_out(θ) = m_in(θ) − m_dom(θ)`` (Bear & Cheng 2010, Ch. 3: mass
21conservation for transport with sorption). ``m_dom`` honors historical
22wave activity via ``wave.was_active_at(theta)`` so retrospective queries
23at θ before a collision event correctly attribute c at v_outlet.
25This file is part of gwtransport which is released under AGPL-3.0 license.
26See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
27"""
29import warnings
30from collections.abc import Sequence
31from operator import itemgetter
33import numpy as np
34import numpy.typing as npt
36from gwtransport.fronttracking.events import find_outlet_crossing
37from gwtransport.fronttracking.math import (
38 ConstantRetardation,
39 NonlinearSorption,
40 SorptionModel,
41)
42from gwtransport.fronttracking.waves import CharacteristicWave, DecayingShockWave, RarefactionWave, ShockWave, Wave
44# Numerical tolerance constants
45EPSILON_VELOCITY = 1e-15 # Tolerance for checking if velocity is effectively zero
46EPSILON_TIME = 1e-15 # Tolerance for negligible time segments
47EPSILON_VOLUME = 1e-15 # Tolerance for negligible spatial-segment widths Δv [m³]
48EPSILON_POSITION = 1e-15 # Tolerance for shock-face proximity in position v [m³]
49# Multiplier on the eps·max(|m_in|,|m_dom|)/Δθ cancellation scale below which a conservation-form
50# residual is numerical zero. The m_dom fan-integral sum accumulates ~few·10² ULP over a
51# multi-pulse record (measured max ratio ~460); 65536 clears that with wide margin while staying
52# ~7 orders below any real breakthrough concentration (O(cin)).
53FP_CANCELLATION_CLAMP = 65536.0
56def concentration_at_point(
57 v: float,
58 theta: float,
59 waves: Sequence[Wave],
60 sorption: SorptionModel, # noqa: ARG001
61) -> float:
62 """Compute concentration at point (v, θ) with exact analytical value.
64 The function works entirely in (V, θ) coordinates: public callers must
65 translate user-facing time t → θ at the API boundary (e.g., via
66 ``FrontTrackerState.theta_at_t``).
68 Parameters
69 ----------
70 v : float
71 Position [m³].
72 theta : float
73 Cumulative flow [m³].
74 waves : list of Wave
75 All waves in the simulation (active and inactive).
76 sorption : SorptionModel
77 Sorption model (unused — kept for API symmetry; wave methods carry
78 their own sorption reference).
80 Returns
81 -------
82 concentration : float
83 Concentration at point (v, θ) [mass/volume].
85 Notes
86 -----
87 **Wave priority**: decaying shocks first (closed-form analytical), then
88 rarefaction fans (spatial extent), then most recently crossing shock or
89 rarefaction tail, then characteristics. If no active wave controls the
90 point, returns 0.0 (initial condition).
91 """
92 # Multi-DSW dispatch: when several active DSW fans contain v, the newest
93 # (largest ``theta_start``) wins — same rule as ``compute_domain_mass``.
94 # Iterating chronologically and short-circuiting picks the older DSW's fan
95 # value, which is wrong for overlap regions. The fan predicate mirrors the
96 # in-fan check used downstream.
97 dsw_in_fan: list[DecayingShockWave] = []
98 for wave in waves:
99 if isinstance(wave, DecayingShockWave) and wave.was_active_at(theta):
100 v_s = wave.position_at_theta(theta)
101 if v_s is None:
102 continue
103 in_fan = (wave.decay_side == "left" and v < v_s) or (wave.decay_side == "right" and v > v_s)
104 if in_fan and v != wave.v_origin:
105 dsw_in_fan.append(wave)
107 if dsw_in_fan:
108 newest = max(dsw_in_fan, key=lambda w: w.theta_start)
109 c = newest.concentration_at_point(v, theta)
110 if c is not None:
111 return c
113 # Fallback: no DSW fan contains v — handle shock-face (v ≈ V_s) and
114 # fixed-side (v on the c_fixed side) cases. Chronological iteration is
115 # fine here because (a) shock-face is rare and unique by v ≈ V_s, and
116 # (b) all DSWs in a typical multi-pulse share the same c_fixed baseline.
117 for wave in waves:
118 if isinstance(wave, DecayingShockWave) and wave.was_active_at(theta):
119 c = wave.concentration_at_point(v, theta)
120 if c is not None:
121 return c
123 # Multi-rarefaction dispatch: same shape as DSWs. RarefactionWave returns
124 # None outside its fan, so collecting non-None candidates is equivalent to
125 # an in-fan check. Newest wins.
126 raref_candidates: list[tuple[float, RarefactionWave]] = []
127 for wave in waves:
128 if isinstance(wave, RarefactionWave) and wave.was_active_at(theta):
129 c = wave.concentration_at_point(v, theta)
130 if c is not None:
131 raref_candidates.append((c, wave))
133 if raref_candidates:
134 return max(raref_candidates, key=lambda cw: cw[1].theta_start)[0]
136 latest_wave_theta = -np.inf
137 latest_wave_c = None
138 # Stacked-shock geometry: for v right of multiple shocks, c at v is c_right
139 # of the shock CLOSEST to v from the left (largest V_shock with V_shock < v).
140 # ``theta_start`` would mis-rank stacked shocks (youngest = innermost ≠
141 # closest-to-v), so we track V_shock directly.
142 # Obstruction check: a shock's c_R applies at v only when no other active
143 # wave (rarefaction tail/head, DSW V_s, other shock) sits between V_shock
144 # and v. Otherwise c_R is the c in the immediate-right-of-shock zone, not
145 # at v.
146 rightmost_passed_v_shock = -np.inf
147 rightmost_passed_c_right: float | None = None
149 def _intervening_wave_between(v_a: float, v_b: float) -> bool:
150 """Return True if any other active wave has a boundary V in (v_a, v_b)."""
151 for other in waves:
152 if not other.was_active_at(theta):
153 continue
154 if isinstance(other, RarefactionWave):
155 v_t = other.tail_position_at_theta(theta)
156 v_h = other.head_position_at_theta(theta)
157 if v_t is not None and v_a < v_t < v_b:
158 return True
159 if v_h is not None and v_a < v_h < v_b:
160 return True
161 else:
162 v_o = other.position_at_theta(theta)
163 if v_o is not None and v_a < v_o < v_b:
164 return True
165 return False
167 for wave in waves:
168 if isinstance(wave, ShockWave) and wave.was_active_at(theta):
169 v_shock = wave.position_at_theta(theta)
170 if v_shock is not None:
171 if abs(v - v_shock) < EPSILON_POSITION:
172 return 0.5 * (wave.c_left + wave.c_right)
174 if abs(wave.speed) > EPSILON_VELOCITY:
175 theta_cross = wave.theta_start + (v - wave.v_start) / wave.speed
177 if theta_cross <= theta:
178 if theta_cross > latest_wave_theta:
179 latest_wave_theta = theta_cross
180 latest_wave_c = wave.c_left
181 elif v > v_shock > rightmost_passed_v_shock and not _intervening_wave_between(v_shock, v):
182 rightmost_passed_v_shock = v_shock
183 rightmost_passed_c_right = wave.c_right
185 for wave in waves:
186 if isinstance(wave, RarefactionWave) and wave.was_active_at(theta):
187 v_tail = wave.tail_position_at_theta(theta)
188 if v_tail is not None and v_tail > v + EPSILON_POSITION:
189 tail_speed = wave.tail_speed()
190 if tail_speed > EPSILON_VELOCITY:
191 theta_pass = wave.theta_start + (v - wave.v_start) / tail_speed
192 if theta_pass <= theta and theta_pass > latest_wave_theta:
193 latest_wave_theta = theta_pass
194 latest_wave_c = wave.c_tail
196 # Priority: latest_wave_c is set by the IF-shock branch (theta_cross <=
197 # theta -> c_L of closest-passing shock) or by the rarefaction-tail loop
198 # (rarefaction tail downstream of v passed v at theta_pass -> c_tail).
199 # rightmost_passed_c_right is set by the elif-shock branch (v > V_shock,
200 # shock upstream of v, c_R of closest-from-left shock). The two represent
201 # different geometric truths: latest_wave_c is "most recent event AT v",
202 # rightmost_passed_c_right is "c just past the upstream shock". When a
203 # rarefaction tail is downstream of the shock AND right of v, that
204 # rarefaction's c_tail wins (the tail's passage at v is more recent than
205 # the shock's contribution at v, geometrically). Prefer latest_wave_c.
206 if latest_wave_c is not None:
207 return latest_wave_c
208 if rightmost_passed_c_right is not None:
209 return rightmost_passed_c_right
211 latest_c = 0.0
212 latest_theta = -np.inf
214 for wave in waves:
215 if isinstance(wave, CharacteristicWave) and wave.was_active_at(theta):
216 v_char_at_theta = wave.position_at_theta(theta)
218 if v_char_at_theta is not None and v_char_at_theta >= v - EPSILON_POSITION:
219 speed = wave.speed()
221 if speed > EPSILON_VELOCITY:
222 theta_pass = wave.theta_start + (v - wave.v_start) / speed
224 if theta_pass <= theta and theta_pass > latest_theta:
225 latest_theta = theta_pass
226 latest_c = wave.concentration
228 return latest_c
231def compute_breakthrough_curve(
232 theta_array: npt.NDArray[np.floating],
233 v_outlet: float,
234 waves: Sequence[Wave],
235 sorption: SorptionModel,
236) -> npt.NDArray[np.floating]:
237 """Concentration at the outlet evaluated over a θ-array (breakthrough curve).
239 Parameters
240 ----------
241 theta_array : array-like
242 Cumulative-flow points at which to query the outlet concentration [m³].
243 Must be sorted in ascending order. Callers translate from user-facing
244 time via ``FrontTrackerState.theta_at_t`` before passing.
245 v_outlet : float
246 Outlet position [m³].
247 waves : list of Wave
248 All waves in the simulation.
249 sorption : SorptionModel
250 Sorption model.
252 Returns
253 -------
254 c_out : numpy.ndarray
255 Concentration at ``v_outlet`` for each θ in ``theta_array`` [mass/volume].
257 See Also
258 --------
259 concentration_at_point : Point-wise concentration
260 compute_bin_averaged_concentration_exact : Bin-averaged concentrations
262 Examples
263 --------
264 .. disable_try_examples
266 ::
268 theta_array = np.linspace(0.0, tracker.state.theta_edges[-1], 1000)
269 c_out = compute_breakthrough_curve(
270 theta_array, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption
271 )
272 """
273 theta_arr = np.asarray(theta_array, dtype=float)
274 c_out = np.zeros(len(theta_arr))
275 for i, theta in enumerate(theta_arr):
276 c_out[i] = concentration_at_point(v_outlet, float(theta), waves, sorption)
277 return c_out
280def identify_outlet_segments(
281 theta_start: float,
282 theta_end: float,
283 v_outlet: float,
284 waves: Sequence[Wave],
285 sorption: SorptionModel,
286) -> list[dict]:
287 """Identify which waves control outlet concentration in θ-interval [theta_start, theta_end].
289 Finds all wave crossing events at the outlet and constructs segments where
290 concentration is constant or varying (rarefaction). All times are expressed
291 as cumulative flow θ [m³].
293 Parameters
294 ----------
295 theta_start : float
296 Start of cumulative-flow interval [m³].
297 theta_end : float
298 End of cumulative-flow interval [m³].
299 v_outlet : float
300 Outlet position [m³].
301 waves : list of Wave
302 All waves in the simulation.
303 sorption : SorptionModel
304 Sorption model.
306 Returns
307 -------
308 segments : list of dict
309 List of segment dictionaries, each containing:
311 - 'theta_start' : float
312 Segment start θ [m³]
313 - 'theta_end' : float
314 Segment end θ [m³]
315 - 'type' : str
316 ``'constant'``, ``'rarefaction'``, or ``'decaying_fan'``.
317 ``'decaying_fan'`` is owned by a :class:`~gwtransport.fronttracking.waves.DecayingShockWave` after
318 its head crosses ``v_outlet``; c at ``v_outlet`` then follows the
319 wave's self-similar fan profile.
320 - 'concentration' : float
321 For constant segments
322 - 'wave' : Wave
323 For rarefaction and decaying_fan segments
324 - 'c_start' : float
325 Concentration at segment start
326 - 'c_end' : float
327 Concentration at segment end
329 Notes
330 -----
331 Segments are constructed by:
333 1. Finding all wave crossing events at the outlet for θ in [theta_start, theta_end].
334 2. Sorting events by θ.
335 3. Creating constant-concentration segments between events.
336 4. Handling rarefaction and decaying-fan profiles with θ-varying concentration.
338 The segments completely partition the interval [theta_start, theta_end].
339 """
340 # Find all waves that cross outlet in this θ-range
341 outlet_events: list[dict] = []
343 # Track rarefactions / decaying shocks that already contain the outlet at
344 # theta_start (no crossing event in [theta_start, theta_end]).
345 active_rarefactions_at_start: list[RarefactionWave | DecayingShockWave] = []
347 for wave in waves:
348 # Retrospective filter: ``identify_outlet_segments`` is called over
349 # arbitrary [theta_start, theta_end] windows (e.g., plotting after the
350 # simulation ends). ``is_active`` is the wave's *current* (end-of-sim)
351 # state and skips waves that legitimately crossed v_outlet during the
352 # window but were later deactivated by a collision. Skip only if the
353 # wave's lifetime ended before the window started.
354 if wave.theta_deactivation <= theta_start:
355 continue
357 if isinstance(wave, DecayingShockWave):
358 # The wave's outlet crossing arrival behaves like a rarefaction head
359 # arrival: before arrival, v_outlet is downstream (c=c_fixed for
360 # decay_side='left'); after arrival, v_outlet is inside the fan
361 # whose c follows the self-similar profile and asymptotes to the
362 # fan's tail concentration.
363 theta_cross = wave.outlet_crossing_theta(v_outlet)
364 if theta_cross is None:
365 continue
366 if theta_cross <= theta_start:
367 # Outlet already inside the fan at theta_start.
368 active_rarefactions_at_start.append(wave)
369 elif theta_cross <= theta_end:
370 # c_after is the fan c just past arrival (the decay-side c at
371 # the arrival θ). theta_cross > wave.theta_start by construction
372 # (outlet_crossing_theta enforces v_outlet > v_start), so
373 # c_decay_at_theta does not return None.
374 c_after = wave.c_decay_at_theta(theta_cross)
375 outlet_events.append({
376 "theta": theta_cross,
377 "wave": wave,
378 "boundary": "head",
379 "c_after": c_after,
380 })
381 continue
383 # For rarefactions, detect both head and tail crossings
384 if isinstance(wave, RarefactionWave):
385 # Check if outlet is already inside this rarefaction at theta_start
386 if wave.contains_point(v_outlet, theta_start):
387 active_rarefactions_at_start.append(wave)
388 # Detect when the tail crosses during [theta_start, theta_end]
389 tail_speed = wave.tail_speed()
390 if tail_speed > EPSILON_VELOCITY:
391 theta_cross = wave.theta_start + (v_outlet - wave.v_start) / tail_speed
392 if theta_start < theta_cross <= theta_end:
393 outlet_events.append({
394 "theta": theta_cross,
395 "wave": wave,
396 "boundary": "tail",
397 "c_after": wave.c_tail,
398 })
399 continue
401 # Head crossing
402 head_speed = wave.head_speed()
403 if head_speed > EPSILON_VELOCITY and wave.v_start < v_outlet:
404 theta_cross = wave.theta_start + (v_outlet - wave.v_start) / head_speed
405 if theta_start <= theta_cross <= theta_end:
406 outlet_events.append({
407 "theta": theta_cross,
408 "wave": wave,
409 "boundary": "head",
410 "c_after": wave.c_head,
411 })
413 # Tail crossing
414 tail_speed = wave.tail_speed()
415 if tail_speed > EPSILON_VELOCITY and wave.v_start < v_outlet:
416 theta_cross = wave.theta_start + (v_outlet - wave.v_start) / tail_speed
417 if theta_start <= theta_cross <= theta_end:
418 outlet_events.append({
419 "theta": theta_cross,
420 "wave": wave,
421 "boundary": "tail",
422 "c_after": wave.c_tail,
423 })
424 else:
425 # Characteristics and shocks
426 theta_cross = find_outlet_crossing(wave, v_outlet, theta_start)
428 if theta_cross is not None and theta_start <= theta_cross <= theta_end:
429 if isinstance(wave, CharacteristicWave):
430 c_after = wave.concentration
431 elif isinstance(wave, ShockWave):
432 # After shock passes outlet, outlet sees left (upstream) state
433 c_after = wave.c_left
434 else:
435 c_after = 0.0
437 outlet_events.append({"theta": theta_cross, "wave": wave, "boundary": None, "c_after": c_after})
439 # Sort events by θ
440 outlet_events.sort(key=itemgetter("theta"))
442 # Create segments between events
443 segments: list[dict] = []
444 current_theta = theta_start
445 current_c = concentration_at_point(v_outlet, theta_start, waves, sorption)
447 # Handle case where we start inside a rarefaction or decaying-shock fan.
448 # Multi-fan overlap: pick the newest (largest ``theta_start``) — matches
449 # ``concentration_at_point`` and ``compute_domain_mass`` dispatch.
450 if active_rarefactions_at_start:
451 raref = max(active_rarefactions_at_start, key=lambda w: w.theta_start)
453 if isinstance(raref, RarefactionWave):
454 # Find when tail crosses (if it does)
455 tail_cross_theta = None
456 for event in outlet_events:
457 if event["wave"] is raref and event["boundary"] == "tail" and event["theta"] > theta_start:
458 tail_cross_theta = event["theta"]
459 break
461 raref_end = min(tail_cross_theta or theta_end, theta_end)
462 c_end = raref.c_tail if tail_cross_theta and tail_cross_theta <= theta_end else None
464 segments.append({
465 "theta_start": theta_start,
466 "theta_end": raref_end,
467 "type": "rarefaction",
468 "wave": raref,
469 "c_start": current_c,
470 "c_end": c_end,
471 })
472 else:
473 # DecayingShockWave fan extends to θ=+∞ (or asymptotes to c_fixed
474 # for n>1 with c_min); treat the whole [theta_start, theta_end]
475 # as one decaying-fan segment.
476 raref_end = theta_end
477 c_end = concentration_at_point(v_outlet, theta_end, waves, sorption)
479 segments.append({
480 "theta_start": theta_start,
481 "theta_end": raref_end,
482 "type": "decaying_fan",
483 "wave": raref,
484 "c_start": current_c,
485 "c_end": c_end,
486 })
488 current_theta = raref_end
489 current_c = (
490 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption) if raref_end < theta_end else current_c
491 )
493 for event in outlet_events:
494 # Skip events that fall inside an already-emitted (typically rarefaction)
495 # segment. ``concentration_at_point`` lets active rarefactions "win"
496 # over a behind-shock c_left; the segment list must reflect the same
497 # convention to avoid double-counting.
498 if event["theta"] < current_theta:
499 continue
501 if isinstance(event["wave"], RarefactionWave) and event["boundary"] == "head":
502 if event["theta"] > current_theta:
503 segments.append({
504 "theta_start": current_theta,
505 "theta_end": event["theta"],
506 "type": "constant",
507 "concentration": current_c,
508 "c_start": current_c,
509 "c_end": current_c,
510 })
512 raref = event["wave"]
513 tail_cross_theta = None
514 for later_event in outlet_events:
515 if (
516 later_event["wave"] is raref
517 and later_event["boundary"] == "tail"
518 and later_event["theta"] > event["theta"]
519 ):
520 tail_cross_theta = later_event["theta"]
521 break
523 raref_end = min(tail_cross_theta or theta_end, theta_end)
525 segments.append({
526 "theta_start": event["theta"],
527 "theta_end": raref_end,
528 "type": "rarefaction",
529 "wave": raref,
530 "c_start": raref.c_head,
531 "c_end": raref.c_tail if tail_cross_theta and tail_cross_theta <= theta_end else None,
532 })
534 current_theta = raref_end
535 current_c = (
536 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption)
537 if raref_end < theta_end
538 else current_c
539 )
540 elif isinstance(event["wave"], DecayingShockWave) and event["boundary"] == "head":
541 if event["theta"] > current_theta:
542 segments.append({
543 "theta_start": current_theta,
544 "theta_end": event["theta"],
545 "type": "constant",
546 "concentration": current_c,
547 "c_start": current_c,
548 "c_end": current_c,
549 })
551 decaying = event["wave"]
552 # The decaying_fan segment ends at the next outlet-crossing event
553 # (if any falls in the window) or at theta_end. Without this split,
554 # multi-DSW pulses would have the first DSW's fan swallow every
555 # later wave's arrival.
556 seg_end = theta_end
557 for later_event in outlet_events:
558 if later_event["theta"] > event["theta"] and later_event["theta"] <= theta_end:
559 seg_end = later_event["theta"]
560 break
561 c_end_val = concentration_at_point(v_outlet, seg_end, waves, sorption)
563 segments.append({
564 "theta_start": event["theta"],
565 "theta_end": seg_end,
566 "type": "decaying_fan",
567 "wave": decaying,
568 "c_start": event["c_after"],
569 "c_end": c_end_val,
570 })
572 current_theta = seg_end
573 current_c = c_end_val
574 else:
575 if event["theta"] > current_theta:
576 segments.append({
577 "theta_start": current_theta,
578 "theta_end": event["theta"],
579 "type": "constant",
580 "concentration": current_c,
581 "c_start": current_c,
582 "c_end": current_c,
583 })
585 current_theta = event["theta"]
586 current_c = event["c_after"]
588 # Final segment
589 if theta_end > current_theta:
590 segments.append({
591 "theta_start": current_theta,
592 "theta_end": theta_end,
593 "type": "constant",
594 "concentration": current_c,
595 "c_start": current_c,
596 "c_end": current_c,
597 })
599 return segments
602def integrate_rarefaction_exact(
603 raref: RarefactionWave, v_outlet: float, theta_start: float, theta_end: float, sorption: SorptionModel
604) -> float:
605 """Exact θ-integral ``∫ c(θ) dθ`` of a rarefaction at the outlet.
607 Convenience wrapper over :func:`integrate_fan_exact` that pulls the fan
608 apex from ``raref.theta_start, raref.v_start``. Returns the mass-like
609 quantity ``∫ c dθ`` (= ``∫ c·flow dt`` in time coordinates).
611 Parameters
612 ----------
613 raref : RarefactionWave
614 Rarefaction wave controlling the outlet.
615 v_outlet : float
616 Outlet position [m³].
617 theta_start, theta_end : float
618 Integration range in cumulative flow [m³]. Either can be ``±np.inf``.
619 sorption : SorptionModel
620 Sorption model (any NonlinearSorption subclass).
622 Returns
623 -------
624 integral : float
625 ``∫ c(θ) dθ`` [mass — i.e. concentration × volume].
626 """
627 return integrate_fan_exact(
628 raref.theta_start, raref.v_start, v_outlet, theta_start, theta_end, sorption, c_apex=raref.c_tail
629 )
632def integrate_fan_exact(
633 theta_origin: float,
634 v_origin: float,
635 v_outlet: float,
636 theta_start: float,
637 theta_end: float,
638 sorption: SorptionModel,
639 c_apex: float = 0.0,
640) -> float:
641 """Exact θ-integral ``∫ c(θ) dθ`` for any self-similar fan at the outlet.
643 Decoupled from the wave object so the same closed-form math applies to
644 both :class:`~gwtransport.fronttracking.waves.RarefactionWave` (apex = ``theta_start, v_start``) and
645 :class:`~gwtransport.fronttracking.waves.DecayingShockWave` (apex = ``theta_origin, v_origin``).
647 Parameters
648 ----------
649 theta_origin, v_origin : float
650 Cumulative flow and position at the fan's apex [m³].
651 v_outlet : float
652 Outlet position [m³].
653 theta_start, theta_end : float
654 Integration range in cumulative flow [m³]. ``theta_end`` may be
655 ``+np.inf``; ``theta_start`` must be finite.
656 sorption : SorptionModel
657 Sorption model (any NonlinearSorption subclass).
658 c_apex : float, optional
659 Concentration on the constant side at the fan apex. For
660 ``RarefactionWave`` this is ``raref.c_tail``; for
661 ``DecayingShockWave`` (decay_side='left') this is ``wave.c_fixed``.
662 For ``c_apex > 0`` the fan formula extrapolates past the physical
663 fan range; the integration is clamped at ``θ_tail`` (where
664 ``c(θ_tail) = c_apex``) and the constant-c_apex region beyond
665 contributes ``c_apex · (theta_end − θ_tail)``. Default 0.0
666 preserves the c=0 apex behavior for canonical c_R=0 fans.
668 Returns
669 -------
670 float
671 Mass-like quantity ``∫ c(θ) dθ`` [mass — concentration × volume].
673 Raises
674 ------
675 TypeError
676 If the sorption model does not support exact fan integration.
677 """
678 # Every NonlinearSorption uses one universal IBP antiderivative, which evaluates the fan
679 # kernel k = R·c − C_T at the segment endpoints via c_and_total_from_retardation.
680 if isinstance(sorption, NonlinearSorption):
681 return _integrate_fan_exact_universal(
682 theta_origin, v_origin, v_outlet, theta_start, theta_end, sorption, c_apex
683 )
685 msg = f"Exact fan integration not supported for {type(sorption).__name__}"
686 raise TypeError(msg)
689def _integrate_fan_exact_universal(
690 theta_origin: float,
691 v_origin: float,
692 v_outlet: float,
693 theta_start: float,
694 theta_end: float,
695 sorption: NonlinearSorption,
696 c_apex: float = 0.0,
697) -> float:
698 r"""Exact θ-integral ``∫ c(θ) dθ`` via the universal IBP antiderivative.
700 For any ``NonlinearSorption`` (with ``R = dC_T/dC``), integration by
701 parts on the self-similar fan ``R(c(θ)) = (θ − θ_origin)/Δv`` gives the
702 closed-form antiderivative
704 .. math::
705 F(\\theta) = c(\\theta)\\,(\\theta - \\theta_{\\rm origin})
706 - \\Delta v \\cdot C_T(c(\\theta)).
708 The derivation uses ``∫ c\\,d\\theta = c·(\\theta-\\theta_0) − ∫(\\theta-\\theta_0)·dc
709 = c·(\\theta-\\theta_0) − \\Delta v · ∫ R(c)\\,dc = c·(\\theta-\\theta_0) − \\Delta v · C_T(c)``,
710 where the last equality is the definition of ``C_T`` as the antiderivative
711 of ``R`` (``R = dC_T/dC``).
713 This formula is exact for any sorption (Brooks-Corey, van Genuchten-Mualem,
714 Freundlich, Langmuir). The only sorption-specific call is
715 ``sorption.concentration_from_retardation`` at the two endpoints — for
716 Brooks-Corey this is closed form; for van Genuchten-Mualem it is one
717 ``brentq`` call per endpoint. No quadrature, no integration loop.
719 Convergence at θ → ∞ for ``c_apex = 0``: for any monotone sorption with
720 ``R(0) = ∞`` (BC, vG, Freundlich n > 1), ``c(∞) = 0`` and ``c·θ → 0``
721 faster than ``Δv·C_T → 0`` (verified termwise from the closed-form
722 asymptotic ``c ~ R^{-α}`` for some ``α > 1``), so ``F(∞) = 0``.
724 For ``c_apex > 0`` the fan formula extrapolates to ``c < c_apex`` past
725 ``θ_tail = θ_origin + Δv·R(c_apex)``; clamp the fan portion at
726 ``θ_tail`` and add ``c_apex·(θ_end − θ_tail)`` for any ``θ_end > θ_tail``.
727 """
728 delta_v = v_outlet - v_origin
729 if delta_v <= 0 or theta_end <= theta_start:
730 return 0.0
732 if c_apex > 0.0:
733 theta_tail = theta_origin + delta_v * float(sorption.retardation(c_apex))
734 theta_end_fan = min(theta_end, theta_tail)
735 else:
736 theta_tail = float("inf")
737 theta_end_fan = theta_end
739 # For a c_apex=0 fan the upper bound may be +∞. The integral converges only when
740 # c → 0 as R → ∞ (so base·c → 0). Freundlich n<1 (c → ∞ as R → ∞) diverges — reject it
741 # explicitly. (DecayingShockWave.mass_after_outlet_arrival returns 0 for the n<1 mirror
742 # before reaching here, so this is a defensive guard for direct callers.)
743 if theta_end_fan == float("inf") and not sorption.fan_converges_at_infinity():
744 msg = "Fan integral diverges at θ=+∞ for this sorption (e.g. Freundlich n<1); pass a finite theta_end"
745 raise ValueError(msg)
747 def antiderivative(theta: float) -> float:
748 if theta == float("inf"):
749 # F(∞) = 0 for any sorption with c → 0 as R → ∞ (guarded above for divergent cases).
750 return 0.0
751 base = (theta - theta_origin) / delta_v
752 if base <= 0.0:
753 return 0.0
754 c, ct = sorption.c_and_total_from_retardation(base)
755 return c * (theta - theta_origin) - delta_v * ct
757 fan_integral = antiderivative(theta_end_fan) - antiderivative(theta_start)
758 constant_contrib = c_apex * max(theta_end - theta_tail, 0.0) if c_apex > 0.0 else 0.0
759 return fan_integral + constant_contrib
762def compute_bin_averaged_concentration_exact(
763 theta_bin_edges: npt.NDArray[np.floating],
764 v_outlet: float,
765 waves: Sequence[Wave],
766 sorption: SorptionModel,
767 *,
768 cin: npt.ArrayLike | None = None,
769 theta_edges_inlet: npt.NDArray[np.floating] | None = None,
770) -> npt.NDArray[np.floating]:
771 """θ-bin-averaged outlet concentration.
773 For each θ-bin ``[θ_i, θ_{i+1}]``::
775 C_avg = (1 / Δθ) · ∫_{θ_i}^{θ_{i+1}} C(v_outlet, θ) dθ
777 With ``cin`` + ``theta_edges_inlet`` provided (recommended for multi-DSW
778 cases), uses the conservation-law identity
779 ``C_avg = (Δm_in − Δm_dom) / Δθ`` per bin — analytical and explicit, no
780 outlet-side fan dispatch. Otherwise falls back to outlet-segment
781 integration (correct for canonical single-DSW cases; may miscount
782 multi-DSW or n<1 mirror geometries).
784 Parameters
785 ----------
786 theta_bin_edges : array-like
787 Cumulative-flow OUTPUT bin edges [m³] (where C_avg is reported).
788 Length N+1 for N bins. Callers translate t-bin edges with
789 ``state.theta_at_t``.
790 v_outlet : float
791 Outlet position [m³].
792 waves : list of Wave
793 All waves from front tracking simulation.
794 sorption : SorptionModel
795 Sorption model.
796 cin : array-like, optional (kw-only)
797 Inlet concentration per inlet θ-bin. When provided with
798 ``theta_edges_inlet``, the conservation form is used.
799 theta_edges_inlet : ndarray, optional (kw-only)
800 θ bin edges of the INLET (``state.theta_edges``), length
801 ``len(cin) + 1``.
803 Returns
804 -------
805 c_avg : numpy.ndarray
806 Bin-averaged outlet concentrations [mass/volume]. Length N.
808 Raises
809 ------
810 ValueError
811 If any output θ-bin has non-positive width.
813 See Also
814 --------
815 concentration_at_point : Point-wise concentration
816 compute_breakthrough_curve : Breakthrough curve
817 compute_cumulative_outlet_mass : Cumulative outlet mass via conservation
818 """
819 theta_edges_out = np.asarray(theta_bin_edges, dtype=float)
820 dtheta_out = np.diff(theta_edges_out)
822 if np.any(dtheta_out <= 0):
823 bad = int(np.argmin(dtheta_out))
824 msg = (
825 f"Invalid θ-bin: theta_bin_edges[{bad}]={theta_edges_out[bad]} >= "
826 f"theta_bin_edges[{bad + 1}]={theta_edges_out[bad + 1]}"
827 )
828 raise ValueError(msg)
830 if cin is not None and theta_edges_inlet is not None:
831 # Conservation form: c_avg = Δm_out/Δθ where m_out = m_in − m_dom.
832 # m_in(θ) = ∫₀^θ cin dτ is the piecewise-LINEAR (in θ) integral of the piecewise-constant
833 # cin, so evaluate it at every output edge in O(N+M) from the inlet-bin cumulative sums
834 # plus the partial bin containing θ — instead of the dense N_out×M_in clip-and-matmul
835 # (≈245× slower, a 128 MB temporary at N=M=4000). Edges below te_in[0] contribute 0
836 # (nothing is injected before the record starts, even if it starts mid-window); edges at
837 # or past te_in[-1] saturate at the total. This mirrors ``compute_cumulative_inlet_mass``'s
838 # clip exactly. m_dom stays a per-θ spatial geometry loop.
839 te_in = np.asarray(theta_edges_inlet, dtype=float)
840 cin_arr = np.asarray(cin, dtype=float)
841 cum_in = np.concatenate([[0.0], np.cumsum(cin_arr * np.diff(te_in))])
842 idx = np.clip(np.searchsorted(te_in, theta_edges_out, side="right") - 1, 0, len(cin_arr) - 1)
843 m_in_at_edges = cum_in[idx] + cin_arr[idx] * (theta_edges_out - te_in[idx])
844 m_in_at_edges = np.where(theta_edges_out < te_in[0], 0.0, m_in_at_edges)
845 m_in_at_edges = np.where(theta_edges_out >= te_in[-1], cum_in[-1], m_in_at_edges)
846 m_dom_at_edges = np.array([
847 compute_domain_mass(theta=float(theta_e), v_outlet=v_outlet, waves=waves, sorption=sorption)
848 for theta_e in theta_edges_out
849 ])
850 # ``compute_cumulative_outlet_mass`` short-circuits to 0 for θ ≤ 0;
851 # replicate that clamp so non-positive output edges contribute no mass.
852 m_out_at_edges = np.where(theta_edges_out <= 0.0, 0.0, m_in_at_edges - m_dom_at_edges)
853 result = np.diff(m_out_at_edges) / dtheta_out
854 # FP-noise clamp scaled to the m_in − m_dom CANCELLATION magnitude: each cumulative-mass
855 # edge carries ~eps·max(|m_in|,|m_dom|) rounding (the m_dom fan-integral sum accumulates
856 # several hundred ULP over a multi-pulse record), amplified by 1/Δθ. The former
857 # 1e-12·max(cout) band keyed off the OUTPUT concentration, which collapses to ~0 before
858 # breakthrough — far too tight, so ~1e-11 cancellation dust tripped the diagnostic on
859 # fully in-range inputs. Residuals within this band are numerical zero: clamp and stay
860 # silent.
861 mass_scale = np.maximum(np.abs(m_in_at_edges), np.abs(m_dom_at_edges))
862 eps_band = (
863 FP_CANCELLATION_CLAMP
864 * np.finfo(float).eps
865 * np.maximum(np.maximum(mass_scale[:-1], mass_scale[1:]), 1.0)
866 / dtheta_out
867 )
868 result = np.where(np.abs(result) < eps_band, 0.0, result)
869 # A residual MORE negative than the FP band is a genuine conservation-form violation with two
870 # possible drivers, and BOTH can be present at once. Attribute per offending bin, not off the
871 # single last edge: a bin whose right edge passes ``theta_edges_inlet[-1]`` sees m_in saturate
872 # while the wave list keeps evolving (out-of-window); a fully in-window offending bin is a real
873 # wave-model over-count (e.g. overlapping non-interacting waves from an oscillating inlet, see
874 # ``compute_domain_mass`` Notes). Report every cause that actually produced a negative bin
875 # instead of the previously unconditional — and often factually false — "edges exceed range"
876 # message. Clamp to 0 to preserve the ``cout >= 0`` API contract either way.
877 neg_mask = result < -eps_band
878 if np.any(neg_mask):
879 worst = float(np.min(result))
880 te_in_last = float(te_in[-1])
881 out_of_window = neg_mask & (theta_edges_out[1:] > te_in_last)
882 causes = []
883 if np.any(out_of_window):
884 causes.append(
885 f"output θ-bin edges exceeding theta_edges_inlet[-1]={te_in_last:.3f} (m_in "
886 "saturates at the last injected mass while the wave list keeps evolving — extend "
887 "cin with trailing zeros to cover the output range, or restrict output bins to the "
888 "inlet window)"
889 )
890 if np.any(neg_mask & ~out_of_window):
891 causes.append(
892 "the domain-mass integral growing faster than the inlet-mass integral within the "
893 "inlet window (an m_dom over-count — e.g. overlapping non-interacting waves from a "
894 "continuous/oscillating inlet — or a wave-list / cin inconsistency)"
895 )
896 warnings.warn(
897 f"compute_bin_averaged_concentration_exact produced a concentration as negative as "
898 f"{worst:.3e}, beyond the FP-cancellation band; likely cause(s): {'; '.join(causes)}.",
899 UserWarning,
900 stacklevel=2,
901 )
902 return np.maximum(result, 0.0)
904 # Legacy outlet-segment integration (compatible with hand-constructed
905 # wave-list tests; correct for canonical single-DSW cases).
906 n_bins = len(theta_edges_out) - 1
907 c_avg = np.zeros(n_bins)
908 for i in range(n_bins):
909 theta_start = float(theta_edges_out[i])
910 theta_end = float(theta_edges_out[i + 1])
911 dtheta_bin = theta_end - theta_start
912 segments = identify_outlet_segments(theta_start, theta_end, v_outlet, waves, sorption)
913 total = 0.0
914 for seg in segments:
915 seg_a = max(seg["theta_start"], theta_start)
916 seg_b = min(seg["theta_end"], theta_end)
917 d = seg_b - seg_a
918 if d <= EPSILON_TIME:
919 continue
920 if seg["type"] == "constant":
921 total += seg["concentration"] * d
922 elif seg["type"] == "rarefaction":
923 if isinstance(sorption, NonlinearSorption):
924 total += integrate_rarefaction_exact(seg["wave"], v_outlet, seg_a, seg_b, sorption)
925 else:
926 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_a + seg_b), waves, sorption)
927 total += c_mid * d
928 elif seg["type"] == "decaying_fan":
929 w = seg["wave"]
930 total += integrate_fan_exact(
931 w.theta_origin, w.v_origin, v_outlet, seg_a, seg_b, sorption, c_apex=w.c_fixed
932 )
933 c_avg[i] = total / dtheta_bin
934 return c_avg
937def compute_domain_mass(
938 theta: float,
939 v_outlet: float,
940 waves: Sequence[Wave],
941 sorption: SorptionModel,
942) -> float:
943 """
944 Compute total mass in domain [0, v_outlet] at cumulative flow θ.
946 Integrates concentration over space::
948 M(θ) = ∫₀^v_outlet C_total(v, θ) dv
950 Exact analytical formulas for every wave type: constant regions
951 (``C_total · Δv``), RarefactionWave fan interiors and DecayingShockWave fan
952 interiors (closed-form via :func:`integrate_fan_spatial_exact`).
954 Parameters
955 ----------
956 theta : float
957 Cumulative flow at which to compute domain mass [m³].
958 v_outlet : float
959 Outlet position (domain extent) [m³].
960 waves : list of Wave
961 All waves in the simulation.
962 sorption : SorptionModel
963 Sorption model.
965 Returns
966 -------
967 mass : float
968 Total mass in domain [mass]. Closed-form analytical to machine precision.
970 See Also
971 --------
972 compute_cumulative_inlet_mass : Cumulative inlet mass
973 compute_cumulative_outlet_mass : Cumulative outlet mass
974 concentration_at_point : Point-wise concentration
975 integrate_fan_spatial_exact : Closed-form fan spatial integral
977 Examples
978 --------
979 .. disable_try_examples
981 ::
983 mass = compute_domain_mass(
984 theta=2500.0, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption
985 )
986 mass >= 0.0
987 """
988 # Partition spatial domain into segments at active wave positions.
989 wave_positions = []
991 for wave in waves:
992 # Use was_active_at(theta) — not is_active — so historical wave geometries
993 # contribute to retrospective m_dom queries. Without this, a wave deactivated
994 # by a later collision event is skipped here, which propagates as a "cin echo
995 # at the outlet" bug (m_out = m_in − 0 instead of m_in − m_dom_correct).
996 if not wave.was_active_at(theta):
997 continue
999 if isinstance(wave, (CharacteristicWave, ShockWave)):
1000 v_pos = wave.position_at_theta(theta)
1001 if v_pos is not None and 0 <= v_pos <= v_outlet:
1002 wave_positions.append(v_pos)
1004 elif isinstance(wave, RarefactionWave):
1005 v_head = wave.head_position_at_theta(theta)
1006 v_tail = wave.tail_position_at_theta(theta)
1008 if v_head is not None and 0 <= v_head <= v_outlet:
1009 wave_positions.append(v_head)
1010 if v_tail is not None and 0 <= v_tail <= v_outlet:
1011 wave_positions.append(v_tail)
1013 elif isinstance(wave, DecayingShockWave):
1014 v_pos = wave.position_at_theta(theta)
1015 if v_pos is not None and 0 <= v_pos <= v_outlet:
1016 wave_positions.append(v_pos)
1017 if 0 <= wave.v_origin <= v_outlet:
1018 wave_positions.append(wave.v_origin)
1020 # Add domain boundaries
1021 wave_positions.extend([0.0, v_outlet])
1023 # Sort and remove duplicates; all entries are within [0, v_outlet] by
1024 # construction (each append site is guarded by the bounds check).
1025 wave_positions = sorted(set(wave_positions))
1027 # Compute mass in each segment using refined integration
1028 total_mass = 0.0
1030 for i in range(len(wave_positions) - 1):
1031 v_start = wave_positions[i]
1032 v_end = wave_positions[i + 1]
1033 dv = v_end - v_start
1035 if dv < EPSILON_VOLUME:
1036 continue
1038 # Check whether the midpoint is inside any fan-bearing wave; the fan
1039 # spatial integral is closed-form for both RarefactionWave and
1040 # DecayingShockWave (the latter via the parameterised
1041 # ``integrate_fan_spatial_exact``).
1042 #
1043 # Multi-fan dispatch: when multiple DSWs (or rarefactions) nominally
1044 # contain v_mid, the newer one (later ``theta_start``) wins —
1045 # geometrically equivalent to the exclusion rule
1046 # ``v_mid ∈ [V_apex_W₂, V_s_W₁]`` for simulator-produced layouts.
1047 v_mid = 0.5 * (v_start + v_end)
1048 raref_wave: RarefactionWave | None = None
1049 decaying_wave: DecayingShockWave | None = None
1051 dsw_candidates: list[DecayingShockWave] = []
1052 for wave in waves:
1053 if not wave.was_active_at(theta):
1054 continue
1055 if isinstance(wave, DecayingShockWave):
1056 v_s = wave.position_at_theta(theta)
1057 if v_s is None:
1058 continue
1059 # decay_side='left': fan is upstream of V_s (v < V_s);
1060 # decay_side='right': fan is downstream of V_s (v > V_s).
1061 in_fan = (wave.decay_side == "left" and v_mid < v_s) or (wave.decay_side == "right" and v_mid > v_s)
1062 if in_fan and v_mid != wave.v_origin:
1063 dsw_candidates.append(wave)
1065 if dsw_candidates:
1066 decaying_wave = max(dsw_candidates, key=lambda w: w.theta_start)
1068 if decaying_wave is None:
1069 raref_candidates = [
1070 wave
1071 for wave in waves
1072 if wave.was_active_at(theta) and isinstance(wave, RarefactionWave) and wave.contains_point(v_mid, theta)
1073 ]
1074 if raref_candidates:
1075 raref_wave = max(raref_candidates, key=lambda w: w.theta_start)
1077 if raref_wave is not None:
1078 mass_segment = _integrate_rarefaction_spatial_exact(raref_wave, v_start, v_end, theta, sorption)
1079 elif decaying_wave is not None:
1080 # The decaying fan is bounded at its apex by the parent's tail concentration
1081 # ``c_fan_tail`` (the sustained-inlet plateau), NOT the downstream ``c_fixed``.
1082 # ``integrate_fan_spatial_exact`` then clamps the abandoned region at ``c_fan_tail``;
1083 # passing ``c_fixed`` would integrate the unbounded self-similar fan to the apex and
1084 # over-count the tail, breaking domain-mass conservation once a DSW forms.
1085 mass_segment = integrate_fan_spatial_exact(
1086 decaying_wave.theta_origin,
1087 decaying_wave.v_origin,
1088 v_start,
1089 v_end,
1090 theta,
1091 sorption,
1092 c_apex=decaying_wave.c_fan_tail,
1093 )
1094 else:
1095 # Constant region: c at midpoint is exact for the segment.
1096 c = concentration_at_point(v_mid, theta, waves, sorption)
1097 c_total = sorption.total_concentration(c)
1098 mass_segment = c_total * dv
1100 total_mass += mass_segment
1102 return float(total_mass)
1105def _integrate_rarefaction_spatial_exact(
1106 raref: RarefactionWave,
1107 v_start: float,
1108 v_end: float,
1109 theta: float,
1110 sorption: SorptionModel,
1111) -> float:
1112 """Exact spatial integral of a rarefaction's total concentration at fixed θ.
1114 Thin wrapper over :func:`integrate_fan_spatial_exact` that pulls the fan
1115 apex from ``raref.theta_start, raref.v_start``. For ``ConstantRetardation``
1116 the rarefaction degenerates to a single c value (no fan), so we use the
1117 wave's ``concentration_at_point`` directly.
1119 Parameters
1120 ----------
1121 raref : RarefactionWave
1122 Rarefaction wave.
1123 v_start, v_end : float
1124 Integration range [m³].
1125 theta : float
1126 Cumulative flow at which to evaluate [m³].
1127 sorption : SorptionModel
1128 Sorption model.
1130 Returns
1131 -------
1132 float
1133 Mass in the segment ``[v_start, v_end]``.
1134 """
1135 if isinstance(sorption, ConstantRetardation):
1136 v_mid = 0.5 * (v_start + v_end)
1137 c = raref.concentration_at_point(v_mid, theta) or 0.0
1138 c_total = sorption.total_concentration(c)
1139 return c_total * (v_end - v_start)
1141 return integrate_fan_spatial_exact(
1142 raref.theta_start, raref.v_start, v_start, v_end, theta, sorption, c_apex=raref.c_tail
1143 )
1146def integrate_fan_spatial_exact(
1147 theta_origin: float,
1148 v_origin: float,
1149 v_start: float,
1150 v_end: float,
1151 theta: float,
1152 sorption: SorptionModel,
1153 c_apex: float = 0.0,
1154) -> float:
1155 """Exact spatial integral ``∫ C_total(v, θ) dv`` for any self-similar fan.
1157 Decoupled from the wave object so the same closed-form math applies to
1158 :class:`~gwtransport.fronttracking.waves.RarefactionWave` (apex = ``theta_start, v_start``) and
1159 :class:`~gwtransport.fronttracking.waves.DecayingShockWave` (apex = ``theta_origin, v_origin``).
1161 In (V, θ) the self-similar fan satisfies ``R(C) = (θ - θ_origin)/(v - v_origin)``;
1162 define ``kappa = θ - θ_origin`` and ``u = v - v_origin``. The dissolved and
1163 sorbed contributions reduce to power-law forms in ``u`` that admit closed
1164 forms via incomplete beta functions (Freundlich) or elementary sqrt
1165 operations (Langmuir).
1167 Parameters
1168 ----------
1169 theta_origin, v_origin : float
1170 Cumulative flow and position at the fan's apex [m³].
1171 v_start, v_end : float
1172 Integration range in v [m³].
1173 theta : float
1174 Cumulative flow at which to evaluate [m³].
1175 sorption : SorptionModel
1176 Sorption model (any NonlinearSorption subclass).
1177 c_apex : float, optional
1178 Concentration on the constant side at the fan apex (typically the
1179 parent rarefaction's ``c_tail`` or the DSW's ``c_fixed`` for
1180 ``decay_side='left'``). For ``c_apex > 0`` the fan formula is
1181 unphysical for ``u < u_tail = kappa / R(c_apex)``; the integration
1182 is split into a constant-C_total(c_apex) region for
1183 ``u ∈ [u_start, u_tail]`` plus the fan integral for
1184 ``u ∈ [u_tail, u_end]``. Default 0.0 preserves the c=0 apex
1185 behavior for canonical c_R=0 rarefactions.
1187 Returns
1188 -------
1189 float
1190 Mass in the segment ``[v_start, v_end]``.
1192 Raises
1193 ------
1194 TypeError
1195 If the sorption model does not support exact spatial integration.
1196 """
1197 if theta <= theta_origin:
1198 return 0.0
1200 kappa = theta - theta_origin
1201 u_start = v_start - v_origin
1202 u_end = v_end - v_origin
1204 # The fan only exists for v > v_origin; clip u_start at 0 (the apex
1205 # contributes nothing to the integral since c(v_origin)=0 for n>1 and
1206 # the Beta-function form handles the lower-bound singularity for n<1).
1207 # If the whole segment is upstream of the apex, return 0.
1208 if u_end <= 0:
1209 return 0.0
1210 if u_start < 0:
1211 u_start = 0.0
1213 # Split off the constant-c_apex region near the apex for c_apex > 0.
1214 # The fan formula is only valid for u ≥ u_tail = kappa / R(c_apex);
1215 # below u_tail, c is clamped to c_apex (the parent's tail / DSW's fixed
1216 # concentration). Spatial counterpart of the temporal θ_tail clamp in
1217 # _integrate_fan_exact_universal.
1218 constant_contrib = 0.0
1219 if c_apex > 0.0:
1220 u_tail = kappa / float(sorption.retardation(c_apex))
1221 if u_start < u_tail:
1222 c_total_apex = float(sorption.total_concentration(c_apex))
1223 constant_contrib = c_total_apex * (min(u_end, u_tail) - u_start)
1224 u_start = u_tail
1225 if u_end <= u_start:
1226 return constant_contrib
1228 # One universal IBP antiderivative for every NonlinearSorption (see the rationale in
1229 # ``integrate_fan_exact``, the temporal counterpart).
1230 if isinstance(sorption, NonlinearSorption):
1231 return constant_contrib + _integrate_rarefaction_spatial_universal(sorption, kappa, u_start, u_end)
1233 msg = f"Exact spatial fan integration not supported for {type(sorption).__name__}"
1234 raise TypeError(msg)
1237def _integrate_rarefaction_spatial_universal(
1238 sorption: NonlinearSorption,
1239 kappa: float,
1240 u_start: float,
1241 u_end: float,
1242) -> float:
1243 r"""Exact spatial integral ``∫ C_T(c(u)) du`` via the universal IBP antiderivative.
1245 For any ``NonlinearSorption`` with ``R = dC_T/dC``, integration by parts
1246 on the self-similar fan ``R(c(u)) = κ/u`` gives the closed-form
1247 antiderivative
1249 .. math::
1250 G(u) = C_T(c(u))\\cdot u - \\kappa\\cdot c(u).
1252 Derivation: ``∫ C_T\\,du = C_T·u − ∫ u\\,dC_T = C_T·u − ∫ (κ/R)·R\\,dc =
1253 C_T·u − κ·c`` (the second equality uses ``u = κ/R`` and the third uses
1254 ``dC_T = R\\,dc``).
1256 Sorption-specific calls limited to ``concentration_from_retardation`` and
1257 ``total_concentration`` at the two endpoints. For Brooks-Corey both are
1258 closed form; for van Genuchten-Mualem ``concentration_from_retardation`` is
1259 one ``brentq`` per endpoint and ``total_concentration`` is also one
1260 ``brentq`` per endpoint (chained internally). No quadrature.
1262 At the apex (``u → 0``, ``R → ∞``) ``c → 0`` and ``C_T → 0`` so ``G(0) = 0``;
1263 a segment whose lower bound is at (or below, for a bounded-``R`` sorption like
1264 Langmuir where ``c = 0`` for ``u`` below the fan tail) the apex contributes
1265 ``G(u_end) − 0``.
1266 """
1267 if u_end <= 0.0 or u_end <= u_start or kappa <= 0.0:
1268 return 0.0
1269 if u_start <= 0.0:
1270 g_start = 0.0 # G(0) = C_T(0)·0 − κ·0 = 0 (apex)
1271 else:
1272 c_start, ct_start = sorption.c_and_total_from_retardation(kappa / u_start)
1273 g_start = ct_start * u_start - kappa * c_start
1274 c_end, ct_end = sorption.c_and_total_from_retardation(kappa / u_end)
1275 g_end = ct_end * u_end - kappa * c_end
1276 return g_end - g_start
1279def compute_cumulative_inlet_mass(
1280 theta: float,
1281 cin: npt.ArrayLike,
1282 theta_edges: npt.ArrayLike,
1283) -> float:
1284 """Cumulative inlet mass entering the domain from θ=0 to ``theta``.
1286 In cumulative-flow coordinates ``M_in(θ) = ∫₀^θ cin(τ) dτ``; for
1287 piecewise-constant ``cin`` this is exact under summation over θ-bin
1288 widths.
1290 Parameters
1291 ----------
1292 theta : float
1293 Cumulative flow up to which to integrate [m³].
1294 cin : array-like
1295 Inlet concentration per θ-bin [mass/volume].
1296 theta_edges : array-like
1297 θ bin edges [m³], length ``len(cin) + 1``.
1299 Returns
1300 -------
1301 mass_in : float
1302 Cumulative inlet mass [mass].
1304 Examples
1305 --------
1306 .. disable_try_examples
1308 ::
1310 mass_in = compute_cumulative_inlet_mass(
1311 theta=5000.0, cin=cin, theta_edges=theta_edges
1312 )
1313 mass_in >= 0.0
1314 """
1315 te = np.asarray(theta_edges, dtype=float)
1316 widths = np.clip(theta - te[:-1], 0.0, np.diff(te))
1317 return float(np.sum(np.asarray(cin, dtype=float) * widths))
1320def compute_cumulative_outlet_mass(
1321 theta: float,
1322 v_outlet: float,
1323 waves: Sequence[Wave],
1324 sorption: SorptionModel,
1325 *,
1326 cin: npt.ArrayLike,
1327 theta_edges: npt.NDArray[np.floating],
1328) -> float:
1329 """Cumulative mass exiting through the outlet from θ=0 to ``theta``.
1331 Computed analytically via the conservation-law identity::
1333 m_out(θ) = m_in(θ) − m_dom(θ)
1335 derived from integrating the PDE ``∂_θ C_T + ∂_V c = 0`` over the spatial
1336 domain ``[0, v_outlet]`` (Bear & Cheng 2010, Ch. 3: mass conservation
1337 for advection with sorption). This sidesteps the multi-fan dispatch problem
1338 that the outlet-segment integration faces when several DSWs cover
1339 v_outlet simultaneously — every term on the right is purely spatial or a
1340 closed-form inlet sum, no ownership priority needed.
1342 Parameters
1343 ----------
1344 theta : float
1345 Cumulative flow up to which to integrate [m³].
1346 v_outlet : float
1347 Outlet position [m³].
1348 waves : list of Wave
1349 All waves in the simulation.
1350 sorption : SorptionModel
1351 Sorption model.
1352 cin : array-like (kw-only)
1353 Inlet concentration per θ-bin [mass/volume].
1354 theta_edges : ndarray (kw-only)
1355 θ bin edges [m³], length ``len(cin) + 1``.
1357 Returns
1358 -------
1359 mass_out : float
1360 Cumulative outlet mass [mass].
1362 Examples
1363 --------
1364 .. disable_try_examples
1366 ::
1368 mass_out = compute_cumulative_outlet_mass(
1369 theta=5000.0,
1370 v_outlet=500.0,
1371 waves=tracker.state.waves,
1372 sorption=sorption,
1373 cin=cin,
1374 theta_edges=tracker.state.theta_edges,
1375 )
1376 mass_out >= 0.0
1377 """
1378 if theta <= 0.0:
1379 return 0.0
1380 m_in = compute_cumulative_inlet_mass(theta=theta, cin=cin, theta_edges=theta_edges)
1381 m_dom = compute_domain_mass(theta=theta, v_outlet=v_outlet, waves=waves, sorption=sorption)
1382 return m_in - m_dom
1385def compute_total_outlet_mass(
1386 v_outlet: float, # noqa: ARG001
1387 sorption: SorptionModel, # noqa: ARG001
1388 *,
1389 cin: npt.ArrayLike,
1390 theta_edges: npt.NDArray[np.floating],
1391) -> float:
1392 """Total outlet mass over θ → ∞ (finite only for a returning-to-zero pulse).
1394 The final inlet value ``c_∞ = cin[-1]`` is the sustained boundary state as θ → ∞:
1396 - For ``c_∞ = 0`` (canonical c_R=0 pulse): injection ceases, the domain empties, and
1397 every injected mass unit eventually exits — ``m_out_total = m_in_total`` (the finite
1398 record integral ``Σ cin·Δθ``). The wave list is not needed.
1399 - For ``c_∞ > 0`` (sustained ambient): the inlet keeps injecting ``c_∞`` forever, so the
1400 cumulative outlet mass grows without bound — return ``+inf``. The previous formula
1401 ``m_in_total − C_T(c_∞)·v_outlet`` paired the FINITE record integral with the
1402 infinite-time steady-state fill and went **negative** whenever
1403 ``m_in_total < C_T(c_∞)·v_outlet``, which is not a physical outlet mass.
1405 Parameters
1406 ----------
1407 v_outlet : float
1408 Outlet position [m³] (unused for ``c_∞ = 0``; the ``+inf`` branch does not need it).
1409 sorption : SorptionModel
1410 Sorption model (kept for API symmetry; no ``C_T`` evaluation is required).
1411 cin : array-like (kw-only)
1412 Inlet concentration per θ-bin [mass/volume].
1413 theta_edges : ndarray (kw-only)
1414 θ bin edges [m³], length ``len(cin) + 1``.
1416 Returns
1417 -------
1418 float
1419 ``m_in_total`` for ``cin[-1] = 0``; ``+inf`` for ``cin[-1] > 0``.
1421 See Also
1422 --------
1423 compute_cumulative_outlet_mass : Cumulative outlet mass up to a finite θ (use this for
1424 a sustained ``c_∞ > 0`` boundary, where the θ → ∞ total is unbounded).
1425 compute_domain_mass : Spatial integral of C_total in the aquifer
1426 """
1427 cin_arr = np.asarray(cin, dtype=float)
1428 # An empty cin is a malformed (no-bin) input by the cin/theta_edges contract;
1429 # let cin_arr[-1] raise IndexError rather than masking it as c_inf=0.
1430 if float(cin_arr[-1]) > 0.0:
1431 return float("inf")
1432 te = np.asarray(theta_edges, dtype=float)
1433 return float(np.sum(cin_arr * np.diff(te)))