Coverage for src / gwtransport / fronttracking / output.py: 83%
459 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-27 06:32 +0000
1"""
2Concentration Extraction from Front Tracking Solutions.
4This module computes outlet concentrations from front tracking wave solutions
5using exact analytical integration. All calculations maintain machine precision
6with no numerical dispersion.
8Functions
9---------
10concentration_at_point(v, t, waves, sorption)
11 Compute concentration at any point (v, t) in domain
12compute_breakthrough_curve(t_array, v_outlet, waves, sorption)
13 Compute concentration at outlet over time array
14compute_bin_averaged_concentration_exact(t_edges, v_outlet, waves, sorption)
15 Compute bin-averaged concentrations using exact analytical integration
16compute_domain_mass(t, v_outlet, waves, sorption)
17 Compute total mass in domain [0, v_outlet] at time t using exact analytical integration
18compute_cumulative_inlet_mass(t, cin, flow, tedges_days)
19 Compute cumulative mass entering domain from t=0 to t
20compute_cumulative_outlet_mass(t, v_outlet, waves, sorption, flow, tedges_days)
21 Compute cumulative mass exiting domain from t=0 to t
22find_last_rarefaction_start_time(v_outlet, waves)
23 Find time when last rarefaction head reaches outlet
24integrate_rarefaction_total_mass(raref, v_outlet, t_start, sorption, flow)
25 Compute total mass exiting through rarefaction to infinity
26compute_total_outlet_mass(v_outlet, waves, sorption, flow, tedges_days)
27 Compute total integrated outlet mass until all mass has exited
29This file is part of gwtransport which is released under AGPL-3.0 license.
30See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
31"""
33from collections.abc import Sequence
34from operator import itemgetter
36import mpmath as mp
37import numpy as np
38import numpy.typing as npt
39from scipy.special import beta as beta_func
40from scipy.special import betainc
42from gwtransport.fronttracking.events import find_outlet_crossing
43from gwtransport.fronttracking.math import (
44 ConstantRetardation,
45 FreundlichSorption,
46 LangmuirSorption,
47 NonlinearSorption,
48 SorptionModel,
49)
50from gwtransport.fronttracking.waves import CharacteristicWave, RarefactionWave, ShockWave, Wave
52# Numerical tolerance constants
53EPSILON_VELOCITY = 1e-15 # Tolerance for checking if velocity is effectively zero
54EPSILON_BETA = 1e-15 # Tolerance for checking if beta is effectively zero (linear case)
55EPSILON_TIME = 1e-15 # Tolerance for negligible time segments
56EPSILON_TIME_MATCH = 1e-6 # Tolerance for matching arrival times (for rarefaction identification)
57EPSILON_CONCENTRATION = 1e-10 # Tolerance for checking if concentration is effectively zero
60def concentration_at_point(
61 v: float,
62 t: float,
63 waves: Sequence[Wave],
64 sorption: SorptionModel, # noqa: ARG001
65) -> float:
66 """
67 Compute concentration at point (v, t) with exact analytical value.
69 Searches through all waves to find which wave controls the concentration
70 at the given point in space-time. Uses exact analytical solutions for
71 characteristics, shocks, and rarefaction fans.
73 Parameters
74 ----------
75 v : float
76 Position [m³].
77 t : float
78 Time [days].
79 waves : list of Wave
80 All waves in the simulation (active and inactive).
81 sorption : SorptionModel
82 Sorption model.
84 Returns
85 -------
86 concentration : float
87 Concentration at point (v, t) [mass/volume].
89 Notes
90 -----
91 **Wave Priority**:
92 The algorithm checks waves in this order:
93 1. Rarefaction waves (if point is inside rarefaction fan)
94 2. Shocks (discontinuities)
95 3. Rarefaction tails (if rarefaction tail has passed point)
96 4. Characteristics (smooth regions)
98 If no active wave controls the point, returns 0.0 (initial condition).
100 **Rarefaction Tail Behavior**:
101 After a rarefaction tail passes a query point, that point maintains the
102 tail concentration as a plateau. This ensures proper concentration behavior
103 after rarefaction waves pass through.
105 **Machine Precision**:
106 All position and concentration calculations use exact analytical formulas.
107 Numerical tolerance is only used for equality checks (v == v_shock).
109 Examples
110 --------
111 ::
113 sorption = FreundlichSorption(
114 k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
115 )
116 # After running simulation with waves...
117 c = concentration_at_point(v=250.0, t=5.0, waves=all_waves, sorption=sorption)
118 c >= 0.0
119 """
120 # Check rarefactions first (they have spatial extent and override other waves)
121 for wave in waves:
122 if isinstance(wave, RarefactionWave) and wave.is_active:
123 c = wave.concentration_at_point(v, t)
124 if c is not None:
125 return c
127 # Track the most recent wave to affect position v
128 # We need to compare crossing times of shocks and rarefaction tails
129 latest_wave_time = -np.inf
130 latest_wave_c = None
132 # Check shocks - track which shocks control this position
133 for wave in waves:
134 if isinstance(wave, ShockWave) and wave.is_active:
135 v_shock = wave.position_at_time(t)
136 if v_shock is not None:
137 # Tolerance for exact shock position
138 tol = 1e-15
140 if abs(v - v_shock) < tol:
141 # Exactly at shock position
142 return 0.5 * (wave.c_left + wave.c_right)
144 # Determine if shock has crossed position v and when
145 if wave.velocity is not None and abs(wave.velocity) > EPSILON_VELOCITY:
146 t_cross = wave.t_start + (v - wave.v_start) / wave.velocity
148 if t_cross <= t:
149 # Shock has crossed position v by time t
150 # After crossing, point sees c_left (concentration behind shock)
151 if t_cross > latest_wave_time:
152 latest_wave_time = t_cross
153 latest_wave_c = wave.c_left
154 elif v > v_shock and wave.t_start > latest_wave_time:
155 # Point is ahead of shock (shock hasn't reached it yet)
156 # Check if this is the closest shock ahead of us
157 # In this case, we see c_right unless overridden by another wave
158 # We track this with a negative time (shock formation time) to indicate
159 # it's a "passive" state (not actively changed)
160 latest_wave_time = wave.t_start
161 latest_wave_c = wave.c_right
163 # Check rarefaction tails - they can override previous waves
164 for wave in waves:
165 if isinstance(wave, RarefactionWave) and wave.is_active:
166 v_tail = wave.tail_position_at_time(t)
167 if v_tail is not None and v_tail > v + 1e-15:
168 # Rarefaction tail has passed position v
169 # Find when it passed: v_start + tail_velocity * (t_pass - t_start) = v
170 tail_vel = wave.tail_velocity()
171 if tail_vel > EPSILON_VELOCITY:
172 t_pass = wave.t_start + (v - wave.v_start) / tail_vel
173 if t_pass <= t and t_pass > latest_wave_time:
174 latest_wave_time = t_pass
175 latest_wave_c = wave.c_tail
177 if latest_wave_c is not None:
178 return latest_wave_c
180 # Check characteristics
181 # Find the most recent characteristic that has reached position v by time t
182 latest_c = 0.0
183 latest_time = -np.inf
185 for wave in waves:
186 if isinstance(wave, CharacteristicWave) and wave.is_active:
187 # Check if this characteristic has reached position v by time t
188 v_char_at_t = wave.position_at_time(t)
190 if v_char_at_t is not None and v_char_at_t >= v - 1e-15:
191 # This characteristic has passed through v
192 # Find when it passed: v_start + vel*(t_pass - t_start) = v
193 vel = wave.velocity()
195 if vel > EPSILON_VELOCITY:
196 t_pass = wave.t_start + (v - wave.v_start) / vel
198 if t_pass <= t and t_pass > latest_time:
199 latest_time = t_pass
200 latest_c = wave.concentration
202 return latest_c
205def compute_breakthrough_curve(
206 t_array: npt.NDArray[np.floating],
207 v_outlet: float,
208 waves: Sequence[Wave],
209 sorption: SorptionModel,
210) -> npt.NDArray[np.floating]:
211 """
212 Compute concentration at outlet over time array.
214 This is the breakthrough curve: C(v_outlet, t) for all t in t_array.
216 Parameters
217 ----------
218 t_array : array-like
219 Time points [days]. Must be sorted in ascending order.
220 v_outlet : float
221 Outlet position [m³].
222 waves : list of Wave
223 All waves in the simulation.
224 sorption : SorptionModel
225 Sorption model.
227 Returns
228 -------
229 c_out : numpy.ndarray
230 Array of concentrations matching t_array [mass/volume].
232 See Also
233 --------
234 concentration_at_point : Point-wise concentration
235 compute_bin_averaged_concentration_exact : Bin-averaged concentrations
237 Examples
238 --------
239 ::
241 t_array = np.linspace(0, 100, 1000)
242 c_out = compute_breakthrough_curve(
243 t_array, v_outlet=500.0, waves=all_waves, sorption=sorption
244 )
245 len(c_out) == len(t_array)
246 """
247 t_array = np.asarray(t_array, dtype=float)
248 c_out = np.zeros(len(t_array))
250 for i, t in enumerate(t_array):
251 c_out[i] = concentration_at_point(v_outlet, t, waves, sorption)
253 return c_out
256def identify_outlet_segments(
257 t_start: float,
258 t_end: float,
259 v_outlet: float,
260 waves: Sequence[Wave],
261 sorption: SorptionModel,
262) -> list[dict]:
263 """
264 Identify which waves control outlet concentration in time interval [t_start, t_end].
266 Finds all wave crossing events at the outlet and constructs segments where
267 concentration is constant or varying (rarefaction).
269 Parameters
270 ----------
271 t_start : float
272 Start of time interval [days].
273 t_end : float
274 End of time interval [days].
275 v_outlet : float
276 Outlet position [m³].
277 waves : list of Wave
278 All waves in the simulation.
279 sorption : SorptionModel
280 Sorption model.
282 Returns
283 -------
284 segments : list of dict
285 List of segment dictionaries, each containing:
287 - 't_start' : float
288 Segment start time
289 - 't_end' : float
290 Segment end time
291 - 'type' : str
292 'constant' or 'rarefaction'
293 - 'concentration' : float
294 For constant segments
295 - 'wave' : Wave
296 For rarefaction segments
297 - 'c_start' : float
298 Concentration at segment start
299 - 'c_end' : float
300 Concentration at segment end
302 Notes
303 -----
304 Segments are constructed by:
305 1. Finding all wave crossing events at outlet in [t_start, t_end]
306 2. Sorting events chronologically
307 3. Creating constant-concentration segments between events
308 4. Handling rarefaction fans with time-varying concentration
310 The segments completely partition the time interval [t_start, t_end].
311 """
312 # Find all waves that cross outlet in this time range
313 outlet_events = []
315 # Track rarefactions that already contain the outlet at t_start
316 # These need to be handled separately since they don't generate crossing events
317 active_rarefactions_at_start = []
319 for wave in waves:
320 if not wave.is_active:
321 continue
323 # For rarefactions, detect both head and tail crossings
324 if isinstance(wave, RarefactionWave):
325 # Check if outlet is already inside this rarefaction at t_start
326 if wave.contains_point(v_outlet, t_start):
327 active_rarefactions_at_start.append(wave)
328 # Don't add crossing events for this wave since we're already inside it
329 # But we still need to detect when the tail crosses during [t_start, t_end]
330 v_tail = wave.tail_position_at_time(t_start)
331 if v_tail is not None and v_tail < v_outlet:
332 vel_tail = wave.tail_velocity()
333 if vel_tail > 0:
334 dt = (v_outlet - v_tail) / vel_tail
335 t_cross = t_start + dt
336 if t_start < t_cross <= t_end:
337 outlet_events.append({
338 "time": t_cross,
339 "wave": wave,
340 "boundary": "tail",
341 "c_after": wave.c_tail,
342 })
343 continue
345 # Head crossing
346 v_head = wave.head_position_at_time(t_start)
347 if v_head is not None and v_head < v_outlet:
348 vel_head = wave.head_velocity()
349 if vel_head > 0:
350 dt = (v_outlet - v_head) / vel_head
351 t_cross = t_start + dt
352 if t_start <= t_cross <= t_end:
353 outlet_events.append({
354 "time": t_cross,
355 "wave": wave,
356 "boundary": "head",
357 "c_after": wave.c_head,
358 })
360 # Tail crossing
361 v_tail = wave.tail_position_at_time(t_start)
362 if v_tail is not None and v_tail < v_outlet:
363 vel_tail = wave.tail_velocity()
364 if vel_tail > 0:
365 dt = (v_outlet - v_tail) / vel_tail
366 t_cross = t_start + dt
367 if t_start <= t_cross <= t_end:
368 outlet_events.append({
369 "time": t_cross,
370 "wave": wave,
371 "boundary": "tail",
372 "c_after": wave.c_tail,
373 })
374 else:
375 # Characteristics and shocks
376 t_cross = find_outlet_crossing(wave, v_outlet, t_start)
378 if t_cross is not None and t_start <= t_cross <= t_end:
379 if isinstance(wave, CharacteristicWave):
380 c_after = wave.concentration
381 elif isinstance(wave, ShockWave):
382 # After shock passes outlet, outlet sees left (upstream) state
383 c_after = wave.c_left
384 else:
385 c_after = 0.0
387 outlet_events.append({"time": t_cross, "wave": wave, "boundary": None, "c_after": c_after})
389 # Sort events by time
390 outlet_events.sort(key=itemgetter("time"))
392 # Create segments between events
393 segments = []
394 current_t = t_start
395 current_c = concentration_at_point(v_outlet, t_start, waves, sorption)
397 # Handle case where we start inside a rarefaction
398 if active_rarefactions_at_start:
399 # Should only be one rarefaction containing the outlet at t_start
400 raref = active_rarefactions_at_start[0]
402 # Find when tail crosses (if it does)
403 tail_cross_time = None
404 for event in outlet_events:
405 if event["wave"] is raref and event["boundary"] == "tail" and event["time"] > t_start:
406 tail_cross_time = event["time"]
407 break
409 # Create rarefaction segment from t_start
410 raref_end = min(tail_cross_time or t_end, t_end)
412 c_start = concentration_at_point(v_outlet, t_start, waves, sorption)
413 c_end = None
414 if tail_cross_time and tail_cross_time <= t_end:
415 c_end = raref.c_tail
417 segments.append({
418 "t_start": t_start,
419 "t_end": raref_end,
420 "type": "rarefaction",
421 "wave": raref,
422 "c_start": c_start,
423 "c_end": c_end,
424 })
426 current_t = raref_end
427 current_c = (
428 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption) if raref_end < t_end else current_c
429 )
431 for event in outlet_events:
432 # Check if we're entering a rarefaction fan
433 if isinstance(event["wave"], RarefactionWave) and event["boundary"] == "head":
434 # Before rarefaction head: constant segment
435 if event["time"] > current_t:
436 segments.append({
437 "t_start": current_t,
438 "t_end": event["time"],
439 "type": "constant",
440 "concentration": current_c,
441 "c_start": current_c,
442 "c_end": current_c,
443 })
445 # Find when tail crosses (if it does)
446 raref = event["wave"]
447 tail_cross_time = None
449 for later_event in outlet_events:
450 if (
451 later_event["wave"] is raref
452 and later_event["boundary"] == "tail"
453 and later_event["time"] > event["time"]
454 ):
455 tail_cross_time = later_event["time"]
456 break
458 # Rarefaction segment
459 raref_end = min(tail_cross_time or t_end, t_end)
461 segments.append({
462 "t_start": event["time"],
463 "t_end": raref_end,
464 "type": "rarefaction",
465 "wave": raref,
466 "c_start": raref.c_head,
467 "c_end": raref.c_tail if tail_cross_time and tail_cross_time <= t_end else None,
468 })
470 current_t = raref_end
471 current_c = (
472 concentration_at_point(v_outlet, raref_end + 1e-10, waves, sorption) if raref_end < t_end else current_c
473 )
474 else:
475 # Regular event (characteristic or shock crossing)
476 # Segment before event
477 if event["time"] > current_t:
478 segments.append({
479 "t_start": current_t,
480 "t_end": event["time"],
481 "type": "constant",
482 "concentration": current_c,
483 "c_start": current_c,
484 "c_end": current_c,
485 })
487 current_t = event["time"]
488 current_c = event["c_after"]
490 # Final segment
491 if t_end > current_t:
492 segments.append({
493 "t_start": current_t,
494 "t_end": t_end,
495 "type": "constant",
496 "concentration": current_c,
497 "c_start": current_c,
498 "c_end": current_c,
499 })
501 return segments
504def integrate_rarefaction_exact(
505 raref: RarefactionWave, v_outlet: float, t_start: float, t_end: float, sorption: SorptionModel
506) -> float:
507 """
508 Exact analytical integral of rarefaction concentration over time at fixed position.
510 Computes integral of C(t) dt from t_start to t_end where C(t) is the
511 self-similar rarefaction solution at the outlet. Dispatches to
512 isotherm-specific closed-form formulas.
514 Parameters
515 ----------
516 raref : RarefactionWave
517 Rarefaction wave controlling the outlet.
518 v_outlet : float
519 Outlet position [m3].
520 t_start : float
521 Integration start time [days]. Can be -np.inf.
522 t_end : float
523 Integration end time [days]. Can be np.inf.
524 sorption : SorptionModel
525 Sorption model (FreundlichSorption or LangmuirSorption).
527 Returns
528 -------
529 integral : float
530 Exact integral value [concentration * time].
532 Raises
533 ------
534 TypeError
535 If sorption model does not support exact rarefaction integration.
536 """
537 if isinstance(sorption, FreundlichSorption):
538 return _integrate_rarefaction_exact_freundlich(raref, v_outlet, t_start, t_end, sorption)
539 if isinstance(sorption, LangmuirSorption):
540 return _integrate_rarefaction_exact_langmuir(raref, v_outlet, t_start, t_end, sorption)
541 msg = f"Exact rarefaction integration not supported for {type(sorption).__name__}"
542 raise TypeError(msg)
545def _integrate_rarefaction_exact_freundlich(
546 raref: RarefactionWave, v_outlet: float, t_start: float, t_end: float, sorption: FreundlichSorption
547) -> float:
548 """Exact temporal integral for Freundlich rarefaction.
550 See `integrate_rarefaction_exact` for parameters.
552 Returns
553 -------
554 float
555 Exact integral value [concentration * time].
557 Raises
558 ------
559 ValueError
560 If sorption is linear (n = 1) or integral diverges.
562 Notes
563 -----
564 For Freundlich: C(t) = [(kappa*t + mu - 1)/alpha]^(1/beta) where
565 kappa = flow/(v_outlet - v_origin), mu = -flow*t_origin/(v_outlet - v_origin),
566 alpha = rho_b*k_f/(n_por*n), beta = 1/n - 1.
568 Antiderivative: F(t) = coeff * (kappa*t + mu - 1)^(1/beta + 1)
569 """
570 t_origin = raref.t_start
571 v_origin = raref.v_start
572 flow = raref.flow
574 kappa = flow / (v_outlet - v_origin)
575 mu = -flow * t_origin / (v_outlet - v_origin)
577 alpha = sorption.bulk_density * sorption.k_f / (sorption.porosity * sorption.n)
578 beta = 1.0 / sorption.n - 1.0
580 if abs(beta) < EPSILON_BETA:
581 msg = "integrate_rarefaction_exact requires nonlinear sorption (n != 1)"
582 raise ValueError(msg)
584 exponent = 1.0 / beta + 1.0
585 coeff = 1.0 / (alpha ** (1.0 / beta) * kappa * exponent)
587 def antiderivative(t: float) -> float:
588 if np.isinf(t):
589 if t > 0:
590 if exponent < 0:
591 return 0.0
592 msg = f"Integral diverges at t=+∞ with exponent={exponent} > 0"
593 raise ValueError(msg)
594 return 0.0
596 base = kappa * t + mu - 1.0
597 if base <= 0:
598 return 0.0
599 return coeff * base**exponent
601 return antiderivative(t_end) - antiderivative(t_start)
604def _integrate_rarefaction_exact_langmuir(
605 raref: RarefactionWave, v_outlet: float, t_start: float, t_end: float, sorption: LangmuirSorption
606) -> float:
607 """Exact temporal integral for Langmuir rarefaction.
609 See `integrate_rarefaction_exact` for parameters.
611 Returns
612 -------
613 float
614 Exact integral value [concentration * time].
616 Notes
617 -----
618 For Langmuir: C(t) = sqrt(A / B(t)) - K_L where
619 B(t) = kappa*t + mu - 1,
620 kappa = flow/(v_outlet - v_origin), mu = -flow*t_origin/(v_outlet - v_origin),
621 A = rho_b * s_max * K_L / n_por.
623 Antiderivative: F(t) = (2*sqrt(A)/kappa) * sqrt(B(t)) - K_L * t
624 """
625 t_origin = raref.t_start
626 v_origin = raref.v_start
627 flow = raref.flow
629 kappa = flow / (v_outlet - v_origin)
630 mu = -flow * t_origin / (v_outlet - v_origin)
631 a_coeff = sorption.a_coeff
632 k_l = sorption.k_l
634 coeff_sqrt = 2.0 * np.sqrt(a_coeff) / kappa
636 def antiderivative(t: float) -> float:
637 if np.isinf(t):
638 if t > 0:
639 # For Langmuir, sqrt(B) → ∞ as t → ∞: integral diverges.
640 # Physical Langmuir rarefactions always have finite tail velocity,
641 # so t_end should always be finite.
642 msg = "Langmuir rarefaction integral diverges at t=+∞"
643 raise ValueError(msg)
644 return 0.0
646 base = kappa * t + mu - 1.0
647 if base <= 0:
648 return 0.0
649 return coeff_sqrt * np.sqrt(base) - k_l * t
651 return antiderivative(t_end) - antiderivative(t_start)
654def compute_bin_averaged_concentration_exact(
655 t_edges: npt.NDArray[np.floating],
656 v_outlet: float,
657 waves: Sequence[Wave],
658 sorption: SorptionModel,
659) -> npt.NDArray[np.floating]:
660 """
661 Compute bin-averaged concentration using EXACT analytical integration.
663 For each time bin [t_i, t_{i+1}], computes:
664 C_avg = (1/(t_{i+1} - t_i)) * ∫_{t_i}^{t_{i+1}} C(v_outlet, t) dt
666 This is the critical function for maintaining machine precision in output.
667 All integrations use exact analytical formulas with no numerical quadrature.
669 Parameters
670 ----------
671 t_edges : array-like
672 Time bin edges [days]. Length N+1 for N bins.
673 v_outlet : float
674 Outlet position [m³].
675 waves : list of Wave
676 All waves from front tracking simulation.
677 sorption : SorptionModel
678 Sorption model.
680 Returns
681 -------
682 c_avg : numpy.ndarray
683 Bin-averaged concentrations [mass/volume]. Length N.
685 Raises
686 ------
687 ValueError
688 If any time bin has non-positive width (``t_edges[i] >= t_edges[i+1]``),
689 or if an unknown segment type is encountered during integration.
691 See Also
692 --------
693 concentration_at_point : Point-wise concentration
694 compute_breakthrough_curve : Breakthrough curve
695 integrate_rarefaction_exact : Exact rarefaction integration
697 Notes
698 -----
699 **Algorithm**:
701 1. For each bin [t_i, t_{i+1}]:
703 a. Identify which wave segments control outlet during this period
704 b. For each segment, compute: Constant C gives integral = C * Δt,
705 Rarefaction C(t) uses exact analytical integral formula
706 c. Sum segment integrals and divide by bin width
708 **Machine Precision**:
710 - Constant segments: exact to floating-point precision
711 - Rarefaction segments: uses closed-form antiderivative
712 - No numerical quadrature or interpolation
713 - Maintains mass balance to ~1e-14 relative error
715 **Rarefaction Integration**:
717 For Freundlich sorption, rarefaction concentration at outlet varies as::
719 C(t) = [(kappa*t + mu - 1)/alpha]^(1/beta)
721 The exact integral is::
723 ∫ C dt = (1/(alpha^(1/beta)*kappa*exponent)) * (kappa*t + mu - 1)^exponent
725 where exponent = 1/beta + 1.
727 Examples
728 --------
729 ::
731 # After running front tracking simulation
732 t_edges = np.array([0.0, 10.0, 20.0, 30.0])
733 c_avg = compute_bin_averaged_concentration_exact(
734 t_edges=t_edges,
735 v_outlet=500.0,
736 waves=tracker.state.waves,
737 sorption=sorption,
738 )
739 len(c_avg) == len(t_edges) - 1
740 np.all(c_avg >= 0)
741 """
742 t_edges = np.asarray(t_edges, dtype=float)
743 n_bins = len(t_edges) - 1
744 c_avg = np.zeros(n_bins)
746 for i in range(n_bins):
747 t_start = t_edges[i]
748 t_end = t_edges[i + 1]
749 dt = t_end - t_start
751 if dt <= 0:
752 msg = f"Invalid time bin: t_edges[{i}]={t_start} >= t_edges[{i + 1}]={t_end}"
753 raise ValueError(msg)
755 # Identify wave segments controlling outlet in this time bin
756 segments = identify_outlet_segments(t_start, t_end, v_outlet, waves, sorption)
758 # Integrate each segment
759 total_integral = 0.0
761 for seg in segments:
762 seg_t_start = max(seg["t_start"], t_start)
763 seg_t_end = min(seg["t_end"], t_end)
764 seg_dt = seg_t_end - seg_t_start
766 if seg_dt <= EPSILON_TIME: # Skip negligible segments
767 continue
769 if seg["type"] == "constant":
770 # C is constant over segment - exact integral
771 integral = seg["concentration"] * seg_dt
773 elif seg["type"] == "rarefaction":
774 # C(t) given by self-similar solution - use exact analytical integral
775 if isinstance(sorption, NonlinearSorption):
776 raref = seg["wave"]
777 integral = integrate_rarefaction_exact(raref, v_outlet, seg_t_start, seg_t_end, sorption)
778 else:
779 # ConstantRetardation - rarefactions shouldn't form, use constant approximation
780 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_t_start + seg_t_end), waves, sorption)
781 integral = c_mid * seg_dt
782 else:
783 msg = f"Unknown segment type: {seg['type']}"
784 raise ValueError(msg)
786 total_integral += integral
788 c_avg[i] = total_integral / dt
790 return c_avg
793def compute_domain_mass(
794 t: float,
795 v_outlet: float,
796 waves: Sequence[Wave],
797 sorption: SorptionModel,
798) -> float:
799 """
800 Compute total mass in domain [0, v_outlet] at time t using exact analytical integration.
802 Implements runtime mass balance verification as described in High Priority #3
803 of FRONT_TRACKING_REBUILD_PLAN.md. Integrates concentration over space:
805 M(t) = ∫₀^v_outlet C(v, t) dv
807 using exact analytical formulas for each wave type.
809 Parameters
810 ----------
811 t : float
812 Time at which to compute domain mass [days].
813 v_outlet : float
814 Outlet position (domain extent) [m³].
815 waves : list of Wave
816 All waves in the simulation.
817 sorption : SorptionModel
818 Sorption model.
820 Returns
821 -------
822 mass : float
823 Total mass in domain [mass]. Computed to machine precision (~1e-14).
825 See Also
826 --------
827 compute_cumulative_inlet_mass : Cumulative inlet mass
828 compute_cumulative_outlet_mass : Cumulative outlet mass
829 concentration_at_point : Point-wise concentration
831 Notes
832 -----
833 **Algorithm**:
835 1. Partition spatial domain [0, v_outlet] into segments where concentration
836 is controlled by a single wave or is constant.
837 2. For each segment, compute mass analytically:
838 - Constant C: mass = C * Δv
839 - Rarefaction C(v): use exact spatial integration formula
840 3. Sum all segment masses.
842 **Wave Priority** (same as concentration_at_point):
844 1. Rarefactions (if position is inside rarefaction fan)
845 2. Shocks and rarefaction tails (most recent to pass)
846 3. Characteristics (smooth regions)
848 **Rarefaction Spatial Integration**:
850 For a rarefaction fan at fixed time t, concentration varies with position v
851 according to the self-similar solution:
853 R(C) = flow*(t - t_origin)/(v - v_origin)
855 The spatial integral ∫ C dv is computed analytically using the inverse
856 retardation relation.
858 **Integration Precision**:
860 - Constant concentration regions: Exact analytical (C_total * dv)
861 - Rarefaction regions: High-precision trapezoidal quadrature (10+ points)
862 - Overall accuracy: ~1e-10 to 1e-12 relative error
863 - Sufficient for runtime verification; primary outputs remain exact
865 Examples
866 --------
867 ::
869 # After running simulation to time t=10.0
870 mass = compute_domain_mass(
871 t=10.0, v_outlet=500.0, waves=tracker.state.waves, sorption=sorption
872 )
873 mass >= 0.0
874 """
875 # Partition spatial domain into segments based on wave structure
876 # We'll evaluate concentration at many points and identify constant regions
878 # Collect all wave positions at time t
879 wave_positions = []
881 for wave in waves:
882 if not wave.is_active:
883 continue
885 if isinstance(wave, (CharacteristicWave, ShockWave)):
886 v_pos = wave.position_at_time(t)
887 if v_pos is not None and 0 <= v_pos <= v_outlet:
888 wave_positions.append(v_pos)
890 elif isinstance(wave, RarefactionWave):
891 v_head = wave.head_position_at_time(t)
892 v_tail = wave.tail_position_at_time(t)
894 if v_head is not None and 0 <= v_head <= v_outlet:
895 wave_positions.append(v_head)
896 if v_tail is not None and 0 <= v_tail <= v_outlet:
897 wave_positions.append(v_tail)
899 # Add domain boundaries
900 wave_positions.extend([0.0, v_outlet])
902 # Sort and remove duplicates
903 wave_positions = sorted(set(wave_positions))
905 # Remove positions outside domain
906 wave_positions = [v for v in wave_positions if 0 <= v <= v_outlet]
908 # Compute mass in each segment using refined integration
909 total_mass = 0.0
911 for i in range(len(wave_positions) - 1):
912 v_start = wave_positions[i]
913 v_end = wave_positions[i + 1]
914 dv = v_end - v_start
916 if dv < EPSILON_VELOCITY:
917 continue
919 # Check if this segment is inside a rarefaction fan
920 # Sample at midpoint to check
921 v_mid = 0.5 * (v_start + v_end)
922 inside_rarefaction = False
923 raref_wave = None
925 for wave in waves:
926 if isinstance(wave, RarefactionWave) and wave.is_active and wave.contains_point(v_mid, t):
927 inside_rarefaction = True
928 raref_wave = wave
929 break
931 if inside_rarefaction and raref_wave is not None:
932 # Rarefaction: concentration varies with position
933 # Use EXACT analytical spatial integration
934 mass_segment = _integrate_rarefaction_spatial_exact(raref_wave, v_start, v_end, t, sorption)
935 else:
936 # Constant concentration region - exact integration
937 v_mid = 0.5 * (v_start + v_end)
938 c = concentration_at_point(v_mid, t, waves, sorption)
939 c_total = sorption.total_concentration(c)
940 mass_segment = c_total * dv
942 total_mass += mass_segment
944 return total_mass
947def _integrate_rarefaction_spatial_exact(
948 raref: RarefactionWave,
949 v_start: float,
950 v_end: float,
951 t: float,
952 sorption: SorptionModel,
953) -> float:
954 """
955 Exact analytical spatial integral of rarefaction total concentration at fixed time.
957 Computes integral of C_total(v) dv from v_start to v_end analytically using
958 closed-form antiderivatives. This maintains machine precision for the mass
959 balance diagnostic.
961 Parameters
962 ----------
963 raref : RarefactionWave
964 Rarefaction wave.
965 v_start : float
966 Integration start position [m3].
967 v_end : float
968 Integration end position [m3].
969 t : float
970 Time [days].
971 sorption : SorptionModel
972 Sorption model.
974 Returns
975 -------
976 mass : float
977 Exact mass in segment to machine precision.
979 Raises
980 ------
981 TypeError
982 If sorption model does not support exact spatial integration.
984 Notes
985 -----
986 For rarefaction at time t: R(C) = kappa/(v - v0) where kappa = flow*(t - t0).
988 For Freundlich: R = 1 + alpha*C^beta where alpha = rho_b*k_f/(n_por*n),
989 beta = 1/n - 1.
990 Total concentration: C_total = C + (rho_b/n_por)*k_f*C^(1/n).
992 Both integrals reduce to power-law forms u^p (kappa-u)^q du which can be
993 expressed using the generalized incomplete beta function via mpmath.betainc().
995 For Langmuir, the integral uses only sqrt operations.
996 """
997 if isinstance(sorption, ConstantRetardation):
998 # Constant retardation: no rarefactions
999 v_mid = 0.5 * (v_start + v_end)
1000 c = raref.concentration_at_point(v_mid, t) or 0.0
1001 c_total = sorption.total_concentration(c)
1002 return c_total * (v_end - v_start)
1004 t_origin = raref.t_start
1005 v_origin = raref.v_start
1006 flow = raref.flow
1008 if t <= t_origin:
1009 return 0.0
1011 kappa = flow * (t - t_origin)
1012 u_start = v_start - v_origin
1013 u_end = v_end - v_origin
1015 if u_start <= 0 or u_end <= 0:
1016 return 0.0
1018 if isinstance(sorption, LangmuirSorption):
1019 return _integrate_rarefaction_spatial_langmuir(sorption, kappa, u_start, u_end)
1021 if isinstance(sorption, FreundlichSorption):
1022 return _integrate_rarefaction_spatial_freundlich(sorption, kappa, u_start, u_end)
1024 msg = f"Exact spatial rarefaction integration not supported for {type(sorption).__name__}"
1025 raise TypeError(msg)
1028def _integrate_rarefaction_spatial_freundlich(
1029 sorption: FreundlichSorption, kappa: float, u_start: float, u_end: float
1030) -> float:
1031 """Exact spatial integral for Freundlich rarefaction using beta functions.
1033 Returns
1034 -------
1035 float
1036 Exact mass in segment.
1037 """
1038 alpha = sorption.bulk_density * sorption.k_f / (sorption.porosity * sorption.n)
1039 rho_b = sorption.bulk_density
1040 n_por = sorption.porosity
1041 k_f = sorption.k_f
1042 n = sorption.n
1044 beta = 1 / n - 1
1045 t_start_norm = u_start / kappa
1046 t_end_norm = u_end / kappa
1048 # Dissolved: ∫ C(u) du via incomplete beta function
1049 a_diss = 1 - 1 / beta
1050 b_diss = 1 + 1 / beta
1052 if a_diss > 0 and b_diss > 0:
1053 beta_diss = betainc(a_diss, b_diss, t_end_norm) - betainc(a_diss, b_diss, t_start_norm)
1054 beta_diss *= beta_func(a_diss, b_diss)
1055 else:
1056 beta_diss = float(mp.betainc(a_diss, b_diss, t_start_norm, t_end_norm, regularized=False))
1058 coeff_diss = (1 / alpha) ** (1 / beta)
1059 mass_dissolved = coeff_diss * kappa * beta_diss
1061 # Sorbed: ∫ (rho_b/n_por)*k_f*C^(1/n) du
1062 exponent_sorb = 1 / (1 - n)
1063 a_sorb = 1 - exponent_sorb
1064 b_sorb = 1 + exponent_sorb
1066 if a_sorb > 0 and b_sorb > 0:
1067 beta_sorb = betainc(a_sorb, b_sorb, t_end_norm) - betainc(a_sorb, b_sorb, t_start_norm)
1068 beta_sorb *= beta_func(a_sorb, b_sorb)
1069 else:
1070 beta_sorb = float(mp.betainc(a_sorb, b_sorb, t_start_norm, t_end_norm, regularized=False))
1072 coeff_sorb = (rho_b / n_por) * k_f / (alpha**exponent_sorb)
1073 mass_sorbed = coeff_sorb * kappa * beta_sorb
1075 return mass_dissolved + mass_sorbed
1078def _integrate_rarefaction_spatial_langmuir(
1079 sorption: LangmuirSorption, kappa: float, u_start: float, u_end: float
1080) -> float:
1081 """Exact spatial integral of C_total(v) for Langmuir rarefaction.
1083 Returns
1084 -------
1085 float
1086 Exact mass in segment.
1088 Notes
1089 -----
1090 With u = v - v_origin, kappa = flow*(t - t_origin):
1091 C(u) = sqrt(a_coeff*u/(kappa-u)) - K_L
1092 C_total = C + (rho_b/n_por) * s_max * C/(K_L + C)
1094 Since K_L + C = sqrt(a_coeff*u/(kappa-u)):
1095 C/(K_L+C) = 1 - K_L*sqrt((kappa-u)/(a_coeff*u))
1097 The integral simplifies to:
1098 integral C_total du = -2*sqrt(a_coeff)*[sqrt(u*(kappa-u))]_start^end
1099 + (rho_b*s_max/n_por - K_L)*(u_end - u_start)
1100 """
1101 a_coeff = sorption.a_coeff
1102 k_l = sorption.k_l
1103 sorbed_max = sorption.bulk_density * sorption.s_max / sorption.porosity
1105 term_sqrt_end = np.sqrt(u_end * (kappa - u_end))
1106 term_sqrt_start = np.sqrt(u_start * (kappa - u_start))
1108 return -2.0 * np.sqrt(a_coeff) * (term_sqrt_end - term_sqrt_start) + (sorbed_max - k_l) * (u_end - u_start)
1111def compute_cumulative_inlet_mass(
1112 t: float,
1113 cin: npt.ArrayLike,
1114 flow: npt.ArrayLike,
1115 tedges_days: npt.ArrayLike,
1116) -> float:
1117 """
1118 Compute cumulative mass entering domain from t=0 to t.
1120 Integrates inlet mass flux over time:
1121 M_in(t) = ∫₀^t cin(τ) * flow(τ) dτ
1123 using exact analytical integration of piecewise-constant functions.
1125 Parameters
1126 ----------
1127 t : float
1128 Time up to which to integrate [days].
1129 cin : array-like
1130 Inlet concentration time series [mass/volume].
1131 Piecewise constant within bins defined by tedges_days.
1132 flow : array-like
1133 Flow rate time series [m³/day].
1134 Piecewise constant within bins defined by tedges_days.
1135 tedges_days : numpy.ndarray
1136 Time bin edges [days]. Length len(cin) + 1.
1138 Returns
1139 -------
1140 mass_in : float
1141 Cumulative mass entered [mass].
1143 Notes
1144 -----
1145 For piecewise-constant cin and flow:
1146 M_in = Σ cin[i] * flow[i] * dt[i]
1148 where the sum is over all bins from tedges_days[0] to t.
1149 Partial bins are handled exactly.
1151 Examples
1152 --------
1153 ::
1155 mass_in = compute_cumulative_inlet_mass(
1156 t=50.0, cin=cin, flow=flow, tedges_days=tedges_days
1157 )
1158 mass_in >= 0.0
1159 """
1160 tedges_arr = np.asarray(tedges_days, dtype=float)
1161 cin_arr = np.asarray(cin, dtype=float)
1162 flow_arr = np.asarray(flow, dtype=float)
1164 # Find which bin t falls into
1165 if t <= tedges_arr[0]:
1166 return 0.0
1168 if t >= tedges_arr[-1]:
1169 # Integrate all bins
1170 dt = np.diff(tedges_arr)
1171 return float(np.sum(cin_arr * flow_arr * dt))
1173 # Find bin containing t
1174 bin_idx = np.searchsorted(tedges_arr, t, side="right") - 1
1176 # Mass flux across inlet boundary = Q * C_in (aqueous concentration)
1177 # This is correct for sorbing solutes: only dissolved mass flows with water
1178 # Integrate complete bins before t
1179 if bin_idx > 0:
1180 dt_complete = np.diff(tedges_arr[: bin_idx + 1])
1181 mass_complete = np.sum(cin_arr[:bin_idx] * flow_arr[:bin_idx] * dt_complete)
1182 else:
1183 mass_complete = 0.0
1185 # Add partial bin
1186 if bin_idx >= 0 and bin_idx < len(cin_arr):
1187 dt_partial = t - tedges_arr[bin_idx]
1188 mass_partial = cin_arr[bin_idx] * flow_arr[bin_idx] * dt_partial
1189 else:
1190 mass_partial = 0.0
1192 return float(mass_complete + mass_partial)
1195def find_last_rarefaction_start_time(
1196 v_outlet: float,
1197 waves: Sequence[Wave],
1198) -> float:
1199 """
1200 Find the time when the last rarefaction head reaches the outlet.
1202 For rarefactions, we integrate analytically so we only need to know
1203 when the rarefaction starts at the outlet (head arrival).
1205 Parameters
1206 ----------
1207 v_outlet : float
1208 Outlet position [m³].
1209 waves : list of Wave
1210 All waves in the simulation.
1212 Returns
1213 -------
1214 t_last : float
1215 Time when last rarefaction head reaches outlet [days].
1216 For non-rarefaction waves, uses their arrival time.
1217 Returns 0.0 if no waves reach the outlet.
1219 Notes
1220 -----
1221 This function finds when the last wave structure *starts* at the outlet.
1222 For rarefactions, this is the head arrival time. The tail may arrive
1223 much later (or at infinite time for rarefactions to C=0), but the
1224 total mass in the rarefaction is computed analytically.
1226 Examples
1227 --------
1228 ::
1230 t_last = find_last_rarefaction_start_time(v_outlet=500.0, waves=all_waves)
1231 t_last >= 0.0
1232 """
1233 t_last = 0.0
1235 for wave in waves:
1236 if not wave.is_active:
1237 continue
1239 if isinstance(wave, RarefactionWave):
1240 # For rarefaction, use head arrival (when rarefaction starts)
1241 head_vel = wave.head_velocity()
1242 if head_vel > EPSILON_VELOCITY:
1243 t_arrival = wave.t_start + (v_outlet - wave.v_start) / head_vel
1244 t_last = max(t_last, t_arrival)
1245 elif isinstance(wave, (CharacteristicWave, ShockWave)):
1246 # For characteristics and shocks, compute arrival time
1247 vel = wave.velocity if isinstance(wave, ShockWave) else wave.velocity()
1248 if vel is not None and vel > EPSILON_VELOCITY:
1249 t_arrival = wave.t_start + (v_outlet - wave.v_start) / vel
1250 t_last = max(t_last, t_arrival)
1252 return t_last
1255def compute_cumulative_outlet_mass(
1256 t: float,
1257 v_outlet: float,
1258 waves: Sequence[Wave],
1259 sorption: SorptionModel,
1260 flow: npt.ArrayLike,
1261 tedges_days: npt.ArrayLike,
1262) -> float:
1263 """
1264 Compute cumulative mass exiting domain from t=0 to t.
1266 Integrates outlet mass flux over time:
1267 M_out(t) = ∫₀^t cout(τ) * flow(τ) dτ
1269 using exact analytical integration. Outlet concentration cout(τ) is obtained
1270 from the wave solution, and flow(τ) is piecewise constant.
1272 Parameters
1273 ----------
1274 t : float
1275 Time up to which to integrate [days].
1276 v_outlet : float
1277 Outlet position [m³].
1278 waves : list of Wave
1279 All waves in the simulation.
1280 sorption : SorptionModel
1281 Sorption model.
1282 flow : array-like
1283 Flow rate time series [m³/day].
1284 Piecewise constant within bins defined by tedges_days.
1285 tedges_days : numpy.ndarray
1286 Time bin edges [days]. Length len(flow) + 1.
1288 Returns
1289 -------
1290 mass_out : float
1291 Cumulative mass exited [mass].
1293 Notes
1294 -----
1295 The outlet concentration is obtained from wave solution via
1296 concentration_at_point(v_outlet, τ, waves, sorption).
1298 For each flow bin [t_i, t_{i+1}], the mass flux integral is computed
1299 exactly using identify_outlet_segments and analytical integration
1300 (constant segments and rarefaction segments).
1302 Examples
1303 --------
1304 ::
1306 mass_out = compute_cumulative_outlet_mass(
1307 t=50.0,
1308 v_outlet=500.0,
1309 waves=tracker.state.waves,
1310 sorption=sorption,
1311 flow=flow,
1312 tedges_days=tedges_days,
1313 )
1314 mass_out >= 0.0
1315 """
1316 tedges_arr = np.asarray(tedges_days, dtype=float)
1317 flow_arr = np.asarray(flow, dtype=float)
1319 if t <= tedges_arr[0]:
1320 return 0.0
1322 # Integrate bin by bin through all flow bins, then continue to t if needed
1323 total_mass = 0.0
1325 # Process bins within tedges range
1326 n_flow_bins = len(flow_arr)
1328 for i in range(n_flow_bins):
1329 t_bin_start = tedges_arr[i]
1330 t_bin_end = tedges_arr[i + 1]
1332 # Skip bins entirely before t
1333 if t_bin_end <= tedges_arr[0]:
1334 continue
1336 # Clip to [tedges[0], t]
1337 t_bin_start = max(t_bin_start, tedges_arr[0])
1338 t_bin_end = min(t_bin_end, t)
1340 flow_i = flow_arr[i]
1341 dt_i = t_bin_end - t_bin_start
1343 if dt_i <= 0:
1344 continue
1346 # Compute ∫_{t_bin_start}^{t_bin_end} cout(τ) dτ using exact integration
1347 segments = identify_outlet_segments(t_bin_start, t_bin_end, v_outlet, waves, sorption)
1349 integral_c_dt = 0.0
1351 for seg in segments:
1352 seg_t_start = max(seg["t_start"], t_bin_start)
1353 seg_t_end = min(seg["t_end"], t_bin_end)
1354 seg_dt = seg_t_end - seg_t_start
1356 if seg_dt <= EPSILON_TIME:
1357 continue
1359 if seg["type"] == "constant":
1360 integral_c_dt += seg["concentration"] * seg_dt
1361 elif seg["type"] == "rarefaction":
1362 if isinstance(sorption, NonlinearSorption):
1363 raref = seg["wave"]
1364 integral_c_dt += integrate_rarefaction_exact(raref, v_outlet, seg_t_start, seg_t_end, sorption)
1365 else:
1366 # ConstantRetardation - use midpoint
1367 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_t_start + seg_t_end), waves, sorption)
1368 integral_c_dt += c_mid * seg_dt
1370 # Mass for this bin = flow * ∫ cout dt
1371 total_mass += flow_i * integral_c_dt
1373 return float(total_mass)
1376def integrate_rarefaction_total_mass(
1377 raref: RarefactionWave,
1378 v_outlet: float,
1379 t_start: float,
1380 sorption: SorptionModel,
1381 flow: float,
1382) -> float:
1383 """
1384 Compute total mass exiting through a rarefaction.
1386 For a rarefaction that starts at the outlet at time t_start, compute the
1387 total mass that will exit through the rarefaction. Integration endpoint
1388 depends on the rarefaction tail concentration:
1390 - If c_tail ≈ 0: Integrate to infinity (tail extends infinitely)
1391 - If c_tail > 0: Integrate to finite tail arrival time
1393 Parameters
1394 ----------
1395 raref : RarefactionWave
1396 Rarefaction wave.
1397 v_outlet : float
1398 Outlet position [m³].
1399 t_start : float
1400 Time when rarefaction head reaches outlet [days].
1401 sorption : SorptionModel
1402 Sorption model.
1403 flow : float
1404 Flow rate [m³/day] (assumed constant).
1406 Returns
1407 -------
1408 total_mass : float
1409 Total mass that exits through rarefaction [mass].
1411 Notes
1412 -----
1413 Uses the exact analytical integral:
1414 M_total = ∫_{t_start}^{t_end} Q * C(t) dt
1416 where C(t) is the concentration at the outlet from the rarefaction wave.
1418 For n > 1 (favorable sorption), rarefactions typically decrease to C=0,
1419 so t_end = ∞ and the integral converges.
1421 For n < 1 (unfavorable sorption), rarefactions typically increase from
1422 low C to high C, so c_tail > 0 and the tail arrives at finite time
1423 t_end = t_start + (v_outlet - v_start) / tail_velocity.
1425 Examples
1426 --------
1427 ::
1429 mass = integrate_rarefaction_total_mass(
1430 raref=raref_wave,
1431 v_outlet=500.0,
1432 t_start=40.0,
1433 sorption=sorption,
1434 flow=100.0,
1435 )
1436 mass >= 0.0
1437 """
1438 if isinstance(sorption, ConstantRetardation):
1439 # No rarefactions with constant retardation
1440 return 0.0
1442 # Determine integration endpoint based on c_tail
1443 # For rarefactions with c_tail ≈ 0, the tail extends to infinity
1444 # For rarefactions with c_tail > 0, the tail arrives at finite time
1445 if raref.c_tail < EPSILON_CONCENTRATION:
1446 # Rarefaction tail goes to C≈0, extends to infinite time
1447 # This is typical for n > 1 rarefactions (concentration decreases)
1448 t_end = np.inf
1449 else:
1450 # Rarefaction tail has finite concentration, arrives at finite time
1451 # This is typical for n < 1 rarefactions (concentration increases)
1452 tail_velocity = raref.tail_velocity()
1453 if tail_velocity < EPSILON_VELOCITY:
1454 # Tail velocity is effectively zero, extends to infinity
1455 t_end = np.inf
1456 else:
1457 # Compute finite tail arrival time
1458 t_end = raref.t_start + (v_outlet - raref.v_start) / tail_velocity
1460 # Integrate from t_start to t_end
1461 integral_c_dt = integrate_rarefaction_exact(raref, v_outlet, t_start, t_end, sorption)
1463 return flow * integral_c_dt
1466def compute_total_outlet_mass(
1467 v_outlet: float,
1468 waves: Sequence[Wave],
1469 sorption: SorptionModel,
1470 flow: npt.ArrayLike,
1471 tedges_days: npt.ArrayLike,
1472) -> tuple[float, float]:
1473 """
1474 Compute total integrated outlet mass until all mass has exited.
1476 Automatically determines when the last wave passes the outlet and
1477 integrates the outlet mass flux until that time, regardless of the
1478 provided tedges extent.
1480 Parameters
1481 ----------
1482 v_outlet : float
1483 Outlet position [m³].
1484 waves : list of Wave
1485 All waves in the simulation.
1486 sorption : SorptionModel
1487 Sorption model.
1488 flow : array-like
1489 Flow rate time series [m³/day].
1490 Piecewise constant within bins defined by tedges_days.
1491 tedges_days : numpy.ndarray
1492 Time bin edges [days]. Length len(flow) + 1.
1494 Returns
1495 -------
1496 total_mass_out : float
1497 Total mass that exits through outlet [mass].
1498 t_integration_end : float
1499 Time until which integration was performed [days].
1500 This is the time when the last wave passes the outlet.
1502 See Also
1503 --------
1504 compute_cumulative_outlet_mass : Cumulative outlet mass up to time t
1505 find_last_rarefaction_start_time : Find when last rarefaction starts
1506 integrate_rarefaction_total_mass : Total mass in rarefaction to infinity
1508 Notes
1509 -----
1510 This function:
1511 1. Finds when the last rarefaction *starts* at the outlet (head arrival)
1512 2. Integrates outlet mass flux until that time
1513 3. Adds analytical integral of rarefaction mass from start to infinity
1515 For rarefactions to C=0, the tail has infinite arrival time but the
1516 total mass is finite and computed analytically.
1518 Examples
1519 --------
1520 ::
1522 total_mass, t_end = compute_total_outlet_mass(
1523 v_outlet=500.0,
1524 waves=tracker.state.waves,
1525 sorption=sorption,
1526 flow=flow,
1527 tedges_days=tedges_days,
1528 )
1529 total_mass >= 0.0
1530 t_end >= tedges_days[0]
1531 """
1532 # Find when the last rarefaction starts at the outlet
1533 t_last_raref_start = find_last_rarefaction_start_time(v_outlet, waves)
1535 tedges_arr = np.asarray(tedges_days, dtype=float)
1536 flow_arr = np.asarray(flow, dtype=float)
1538 # Integrate up to when last rarefaction starts
1539 if t_last_raref_start <= tedges_arr[-1]:
1540 # All waves start within provided time range
1541 mass_up_to_raref_start = compute_cumulative_outlet_mass(
1542 t_last_raref_start, v_outlet, waves, sorption, flow_arr, tedges_arr
1543 )
1544 flow_at_raref_start = flow_arr[-1] # Use last flow value
1545 else:
1546 # Need to extend beyond tedges to reach rarefaction start
1547 # First, compute mass up to tedges[-1]
1548 mass_within_tedges = compute_cumulative_outlet_mass(
1549 tedges_arr[-1], v_outlet, waves, sorption, flow_arr, tedges_arr
1550 )
1552 # Then, integrate from tedges[-1] to t_last_raref_start
1553 flow_extended = flow_arr[-1]
1554 t_start_extended = tedges_arr[-1]
1555 t_end_extended = t_last_raref_start
1557 # Get outlet segments for extended period
1558 segments = identify_outlet_segments(t_start_extended, t_end_extended, v_outlet, waves, sorption)
1560 integral_c_dt = 0.0
1562 for seg in segments:
1563 seg_t_start = max(seg["t_start"], t_start_extended)
1564 seg_t_end = min(seg["t_end"], t_end_extended)
1565 seg_dt = seg_t_end - seg_t_start
1567 if seg_dt <= EPSILON_TIME:
1568 continue
1570 if seg["type"] == "constant":
1571 integral_c_dt += seg["concentration"] * seg_dt
1572 elif seg["type"] == "rarefaction":
1573 if isinstance(sorption, NonlinearSorption):
1574 raref = seg["wave"]
1575 integral_c_dt += integrate_rarefaction_exact(raref, v_outlet, seg_t_start, seg_t_end, sorption)
1576 else:
1577 # ConstantRetardation - use midpoint
1578 c_mid = concentration_at_point(v_outlet, 0.5 * (seg_t_start + seg_t_end), waves, sorption)
1579 integral_c_dt += c_mid * seg_dt
1581 mass_up_to_raref_start = mass_within_tedges + flow_extended * integral_c_dt
1582 flow_at_raref_start = flow_extended
1584 # Find rarefactions that are active at the outlet after t_last_raref_start
1585 # and add their total integrated mass
1586 total_raref_mass = 0.0
1588 for wave in waves:
1589 if not wave.is_active:
1590 continue
1592 if isinstance(wave, RarefactionWave):
1593 # Check if this rarefaction reaches the outlet
1594 head_vel = wave.head_velocity()
1595 if head_vel > EPSILON_VELOCITY and wave.v_start < v_outlet:
1596 t_raref_start_at_outlet = wave.t_start + (v_outlet - wave.v_start) / head_vel
1598 # If this rarefaction starts at or after t_last_raref_start, include its total mass
1599 # (with small tolerance for numerical precision)
1600 if abs(t_raref_start_at_outlet - t_last_raref_start) < EPSILON_TIME_MATCH:
1601 # This is the last rarefaction - integrate to infinity
1602 raref_mass = integrate_rarefaction_total_mass(
1603 raref=wave,
1604 v_outlet=v_outlet,
1605 t_start=t_raref_start_at_outlet,
1606 sorption=sorption,
1607 flow=flow_at_raref_start,
1608 )
1609 total_raref_mass += raref_mass
1611 # For rarefactions with finite tails (c_tail > 0, typical for n < 1),
1612 # all mass is already accounted for in the rarefaction integration
1613 # from head to tail. No additional mass needs to be integrated after
1614 # the tail - if there were more waves, they would be in the wave list.
1615 total_mass = mass_up_to_raref_start + total_raref_mass
1617 return float(total_mass), t_last_raref_start