Coverage for src/gwtransport/fronttracking/solver.py: 86%
272 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1"""Event-driven front-tracking solver in (V, θ) coordinates.
3The simulation runs entirely in cumulative-flow space θ. Every public
4output — wave attributes, ``state.events[i]['theta']``,
5``theta_first_arrival`` — is in θ. Translation to user-facing time t is
6the caller's responsibility via ``state.t_at_theta``.
7Time-varying flow is absorbed into the precomputed ``theta_edges`` array
8at ``__init__``; there is no flow-change event.
10Algorithm:
121. Initialize waves from inlet boundary conditions (one per cin step at θ_edges[i]).
132. Find next event (earliest collision or outlet crossing in θ).
143. Advance θ to event.
154. Handle event (create new waves, deactivate old ones).
165. Repeat until no more events.
18All calculations are exact analytical with machine precision.
19"""
21import logging
22from dataclasses import dataclass
23from itertools import pairwise
24from operator import itemgetter
26import numpy as np
27import numpy.typing as npt
28import pandas as pd
30from gwtransport._time import tedges_to_days
31from gwtransport.fronttracking.events import (
32 Event,
33 EventType,
34 find_characteristic_intersection,
35 find_outlet_crossing,
36 find_rarefaction_boundary_intersections,
37 find_shock_characteristic_intersection,
38 find_shock_shock_intersection,
39 is_outlet_crossing_pinned,
40)
41from gwtransport.fronttracking.handlers import (
42 EPSILON_CONCENTRATION,
43 create_inlet_waves_at_theta,
44 handle_characteristic_collision,
45 handle_outlet_crossing,
46 handle_rarefaction_characteristic_collision,
47 handle_shock_characteristic_collision,
48 handle_shock_collision,
49 handle_shock_rarefaction_collision,
50)
51from gwtransport.fronttracking.math import (
52 SorptionModel,
53 compute_first_front_arrival_theta,
54)
55from gwtransport.fronttracking.output import (
56 FP_CANCELLATION_CLAMP,
57 compute_cumulative_outlet_mass,
58)
59from gwtransport.fronttracking.waves import (
60 CharacteristicWave,
61 DecayingShockWave,
62 RarefactionWave,
63 ShockWave,
64 Wave,
65)
67logger = logging.getLogger(__name__)
70@dataclass
71class FrontTrackerState:
72 """Complete state of the front-tracking simulation in (V, θ).
74 Parameters
75 ----------
76 waves : list of Wave
77 All waves created during simulation (includes inactive waves).
78 events : list of dict
79 Event history. Records use the ``"theta"`` key carrying the
80 cumulative flow at which the event occurred [m³]. Callers translate
81 to user-facing time via ``FrontTrackerState.t_at_theta``.
82 theta_current : float
83 Current simulation cumulative flow [m³].
84 v_outlet : float
85 Outlet position [m³].
86 sorption : SorptionModel
87 Sorption model.
88 cin : numpy.ndarray
89 Inlet concentration time series [mass/volume].
90 flow : numpy.ndarray
91 Flow rate time series [m³/day], one value per bin.
92 tedges : pandas.DatetimeIndex
93 Time bin edges.
94 tedges_days : numpy.ndarray
95 ``tedges`` as days from ``tedges[0]``, length ``len(flow) + 1``.
96 theta_edges : numpy.ndarray
97 Cumulative flow at every bin edge. ``theta_edges[i] = sum_{k<i} flow[k] *
98 (tedges_days[k+1] - tedges_days[k])``. Length ``len(flow) + 1``.
99 """
101 waves: list[Wave]
102 events: list[dict]
103 theta_current: float
104 v_outlet: float
105 sorption: SorptionModel
106 cin: np.ndarray
107 flow: np.ndarray
108 tedges: pd.DatetimeIndex
109 tedges_days: npt.NDArray[np.floating]
110 theta_edges: npt.NDArray[np.floating]
112 def t_at_theta(self, theta: float) -> float:
113 """Translate cumulative flow θ back to user-facing time t [days].
115 Piecewise linear inversion of the (tedges_days → theta_edges) map.
116 Implementation note on the (rare) zero-flow case: when a bin has
117 ``flow[i] == 0``, θ is constant across ``[tedges_days[i], tedges_days[i+1])``;
118 ``np.searchsorted(..., side='right') - 1`` lands on the rightmost
119 such bin, so this function returns ``tedges_days[i]`` for the
120 right-most i sharing that θ. Events scheduled at zero-flow bin
121 boundaries therefore align with the END of the zero-flow interval
122 — pick one convention and call it documented.
123 """
124 if theta <= self.theta_edges[0]:
125 return float(self.tedges_days[0])
126 if theta >= self.theta_edges[-1]:
127 # Extrapolate: pretend the last bin continues at its current flow.
128 last_flow = float(self.flow[-1])
129 if last_flow > 0:
130 return float(self.tedges_days[-1] + (theta - self.theta_edges[-1]) / last_flow)
131 return float(self.tedges_days[-1])
133 # Find bin index i with theta_edges[i] <= theta < theta_edges[i+1].
134 # np.searchsorted with side='right' returns the smallest i+1 such that
135 # theta_edges[i+1] > theta, so subtracting 1 gives i.
136 # The boundary returns above guarantee strictly-interior theta here, so
137 # searchsorted lands in [0, len(flow)-1] without an extra index clamp.
138 i = int(np.searchsorted(self.theta_edges, theta, side="right")) - 1
139 flow_i = float(self.flow[i])
140 if flow_i <= 0:
141 return float(self.tedges_days[i])
142 return float(self.tedges_days[i] + (theta - self.theta_edges[i]) / flow_i)
144 def theta_at_t(self, t: float) -> float:
145 """Translate user-facing time t [days] to cumulative flow θ [m³].
147 Piecewise linear forward map. Outside the input range the boundary
148 flow is extrapolated.
149 """
150 if t <= self.tedges_days[0]:
151 return float(self.theta_edges[0])
152 if t >= self.tedges_days[-1]:
153 return float(self.theta_edges[-1] + (t - self.tedges_days[-1]) * float(self.flow[-1]))
155 # Boundary returns above guarantee strictly-interior t, so searchsorted
156 # lands in [0, len(flow)-1] without an extra index clamp.
157 i = int(np.searchsorted(self.tedges_days, t, side="right")) - 1
158 return float(self.theta_edges[i] + (t - self.tedges_days[i]) * float(self.flow[i]))
160 def theta_at_t_array(self, t: npt.ArrayLike) -> npt.NDArray[np.floating]:
161 """Vectorized ``theta_at_t``: map an array of times t [days] to θ [m³].
163 Element-wise identical to :meth:`theta_at_t`; replaces per-scalar loops
164 in the plotting/output breakthrough routines.
166 Parameters
167 ----------
168 t : array-like
169 User-facing time points [days].
171 Returns
172 -------
173 ndarray
174 Cumulative flow θ at each ``t`` [m³].
175 """
176 t_arr = np.asarray(t, dtype=float)
177 # Interior map: i = searchsorted(tedges_days, t, 'right') - 1, clipped to
178 # a valid bin so the boundary branches can overwrite the extrapolated tails.
179 i = np.clip(np.searchsorted(self.tedges_days, t_arr, side="right") - 1, 0, len(self.flow) - 1)
180 theta = self.theta_edges[i] + (t_arr - self.tedges_days[i]) * self.flow[i]
181 # Left of the first edge: clamp to theta_edges[0].
182 theta = np.where(t_arr <= self.tedges_days[0], self.theta_edges[0], theta)
183 # Right of the last edge: extrapolate at the final-bin flow.
184 return np.where(
185 t_arr >= self.tedges_days[-1],
186 self.theta_edges[-1] + (t_arr - self.tedges_days[-1]) * float(self.flow[-1]),
187 theta,
188 )
191class FrontTracker:
192 """Event-driven front-tracking solver for nonlinear sorption transport.
194 Parameters
195 ----------
196 cin : numpy.ndarray
197 Inlet concentration time series [mass/volume]; length ``n``.
198 flow : numpy.ndarray
199 Flow rate time series [m³/day]; length ``n`` (one value per bin).
200 tedges : pandas.DatetimeIndex
201 Time bin edges (length ``n+1``).
202 aquifer_pore_volume : float
203 Total pore volume [m³] — used as the outlet position.
204 sorption : SorptionModel
205 Sorption model.
207 Attributes
208 ----------
209 state : FrontTrackerState
210 Complete simulation state.
211 theta_first_arrival : float
212 Cumulative flow θ at which the first nonzero-concentration wave reaches
213 the outlet [m³]. Translate to user-facing time via
214 ``state.t_at_theta(theta_first_arrival)``.
216 Notes
217 -----
218 The solver works exclusively in cumulative flow θ; events appended to
219 ``state.events`` carry ``"theta"``. Translation to user-facing time t is
220 the caller's responsibility (use ``state.t_at_theta``).
221 """
223 def __init__(
224 self,
225 cin: npt.ArrayLike,
226 flow: npt.ArrayLike,
227 tedges: pd.DatetimeIndex,
228 aquifer_pore_volume: float,
229 sorption: SorptionModel,
230 ):
231 cin = np.asarray(cin, dtype=float)
232 flow = np.asarray(flow, dtype=float)
233 if len(tedges) != len(cin) + 1:
234 msg = f"tedges must have length len(cin) + 1, got {len(tedges)} vs {len(cin) + 1}"
235 raise ValueError(msg)
236 if len(flow) != len(cin):
237 msg = f"flow must have same length as cin, got {len(flow)} vs {len(cin)}"
238 raise ValueError(msg)
239 if np.any(cin < 0):
240 msg = "cin must be non-negative"
241 raise ValueError(msg)
242 if np.any(flow < 0):
243 msg = "flow must be non-negative (negative flow not supported)"
244 raise ValueError(msg)
245 if aquifer_pore_volume <= 0:
246 msg = "aquifer_pore_volume must be positive"
247 raise ValueError(msg)
249 tedges_days = tedges_to_days(tedges)
250 dt_days = np.diff(tedges_days)
251 bin_volumes = flow * dt_days
252 theta_edges = np.concatenate(([0.0], np.cumsum(bin_volumes)))
254 self.state = FrontTrackerState(
255 waves=[],
256 events=[],
257 theta_current=0.0,
258 v_outlet=aquifer_pore_volume,
259 sorption=sorption,
260 cin=cin,
261 flow=flow,
262 tedges=tedges,
263 tedges_days=tedges_days,
264 theta_edges=theta_edges,
265 )
267 self.theta_first_arrival = compute_first_front_arrival_theta(cin, theta_edges, aquifer_pore_volume, sorption)
269 self._initialize_inlet_waves()
271 def _initialize_inlet_waves(self):
272 """Emit one wave per nonzero inlet step at the corresponding ``theta_edges[i]``."""
273 c_prev = 0.0
274 theta_edges = self.state.theta_edges
276 for i in range(len(self.state.cin)):
277 c_new = float(self.state.cin[i])
279 if abs(c_new - c_prev) > EPSILON_CONCENTRATION:
280 new_waves = create_inlet_waves_at_theta(
281 c_prev=c_prev,
282 c_new=c_new,
283 theta=float(theta_edges[i]),
284 sorption=self.state.sorption,
285 )
286 self.state.waves.extend(new_waves)
288 c_prev = c_new
290 def find_next_event(self) -> Event | None:
291 """Return the next event in θ-order, or ``None`` if none."""
292 # Each call collects every candidate and selects the single earliest via
293 # ``min`` over the (theta, counter, ...) tuples; the unique ``counter``
294 # breaks θ-ties deterministically (and stops comparison before the
295 # non-orderable EventType/Wave fields). A heap would only pay for its
296 # invariant to extract one minimum, so a flat list + ``min`` is leaner.
297 candidates: list[tuple] = []
298 counter = 0 # Unique counter to break θ-ties deterministically
300 active_waves = [w for w in self.state.waves if w.is_active]
301 theta_current = self.state.theta_current
303 # Defense-in-depth loop guard: reject collision candidates at or below
304 # θ_current + tol so a degenerate geometry cannot re-fire the identical
305 # event forever (same FP-tolerance scale as the outlet-crossing guard
306 # below). The structural fix is routing all shock↔rarefaction collisions
307 # through DecayingShockWave; this backstops any future near-coincident
308 # event from looping the solver.
309 collision_tol = 1e-12 * max(abs(theta_current), 1.0)
311 def push_collision(theta, event_type, waves, v, boundary):
312 nonlocal counter
313 if not (0 <= v <= self.state.v_outlet):
314 return
315 if theta <= theta_current + collision_tol:
316 return
317 candidates.append((theta, counter, event_type, waves, v, boundary))
318 counter += 1
320 chars = [w for w in active_waves if isinstance(w, CharacteristicWave)]
321 for i, w1 in enumerate(chars):
322 for w2 in chars[i + 1 :]:
323 result = find_characteristic_intersection(w1, w2, theta_current)
324 if result:
325 theta, v = result
326 push_collision(theta, EventType.CHAR_CHAR_COLLISION, [w1, w2], v, None)
328 shocks = [w for w in active_waves if isinstance(w, ShockWave)]
329 for i, w1 in enumerate(shocks):
330 for w2 in shocks[i + 1 :]:
331 result = find_shock_shock_intersection(w1, w2, theta_current)
332 if result:
333 theta, v = result
334 push_collision(theta, EventType.SHOCK_SHOCK_COLLISION, [w1, w2], v, None)
336 for shock in shocks:
337 for char in chars:
338 result = find_shock_characteristic_intersection(shock, char, theta_current)
339 if result:
340 theta, v = result
341 push_collision(theta, EventType.SHOCK_CHAR_COLLISION, [shock, char], v, None)
343 rarefs = [w for w in active_waves if isinstance(w, RarefactionWave)]
344 for raref in rarefs:
345 for char in chars:
346 intersections = find_rarefaction_boundary_intersections(raref, char, theta_current)
347 for theta, v, boundary in intersections:
348 push_collision(theta, EventType.RAREF_CHAR_COLLISION, [raref, char], v, boundary)
350 for shock in shocks:
351 for raref in rarefs:
352 intersections = find_rarefaction_boundary_intersections(raref, shock, theta_current)
353 for theta, v, boundary in intersections:
354 push_collision(theta, EventType.SHOCK_RAREF_COLLISION, [shock, raref], v, boundary)
356 for i, raref1 in enumerate(rarefs):
357 for raref2 in rarefs[i + 1 :]:
358 intersections = find_rarefaction_boundary_intersections(raref1, raref2, theta_current)
359 for theta, v, boundary in intersections:
360 push_collision(theta, EventType.RAREF_RAREF_COLLISION, [raref1, raref2], v, boundary)
362 # Fan-exhaustion: a DecayingShockWave is valid only while c_decay stays
363 # above c_fan_tail. When c_decay reaches c_fan_tail the fan is spent and
364 # the wave hands off to a regular ShockWave(c_fan_tail, c_fixed).
365 for wave in active_waves:
366 if isinstance(wave, DecayingShockWave):
367 theta_exhaust = wave.theta_at_fan_exhaustion()
368 if theta_exhaust is None or theta_exhaust <= theta_current + collision_tol:
369 continue
370 v_exhaust = wave.position_at_theta(theta_exhaust)
371 if v_exhaust is None or not (0 <= v_exhaust <= self.state.v_outlet):
372 continue
373 candidates.append((theta_exhaust, counter, EventType.DSW_FAN_EXHAUSTED, [wave], v_exhaust, None))
374 counter += 1
376 v_outlet = self.state.v_outlet
377 # Same FP-tolerance discipline as events.find_outlet_crossing — prevents
378 # re-emitting an outlet crossing for a boundary that's at v_outlet ± ULPs.
379 outlet_tol = 1e-12 * max(abs(v_outlet), 1.0)
381 for wave in active_waves:
382 if isinstance(wave, RarefactionWave):
383 theta_eval = max(theta_current, wave.theta_start)
384 for c_boundary, pos_fn, speed_fn in (
385 (wave.c_head, wave.head_position_at_theta, wave.head_speed),
386 (wave.c_tail, wave.tail_position_at_theta, wave.tail_speed),
387 ):
388 v_pos = pos_fn(theta_eval)
389 if v_pos is None or v_pos >= v_outlet - outlet_tol:
390 continue
391 s = speed_fn()
392 # Skip a c_min-floored (pinned) boundary — R(c_min) inflated
393 # for n>1, c→0: its crossing lands at a non-physical θ~1e8 and
394 # only pollutes the event record.
395 if s <= 0 or is_outlet_crossing_pinned(c_boundary, wave.sorption):
396 continue
397 theta_cross = theta_eval + (v_outlet - v_pos) / s
398 if theta_cross > theta_current:
399 candidates.append((theta_cross, counter, EventType.OUTLET_CROSSING, [wave], v_outlet, None))
400 counter += 1
401 else:
402 theta_cross = find_outlet_crossing(wave, self.state.v_outlet, theta_current)
403 if theta_cross and theta_cross > theta_current:
404 candidates.append((
405 theta_cross,
406 counter,
407 EventType.OUTLET_CROSSING,
408 [wave],
409 self.state.v_outlet,
410 None,
411 ))
412 counter += 1
414 if not candidates:
415 return None
417 theta_event, _, event_type, waves, v, extra = min(candidates, key=itemgetter(slice(2)))
419 raref_types = {
420 EventType.RAREF_CHAR_COLLISION,
421 EventType.SHOCK_RAREF_COLLISION,
422 EventType.RAREF_RAREF_COLLISION,
423 }
424 boundary_type = extra if event_type in raref_types else None
426 return Event(
427 theta=theta_event,
428 event_type=event_type,
429 waves_involved=waves,
430 location=v,
431 boundary_type=boundary_type,
432 )
434 def handle_event(self, event: Event):
435 """Dispatch an event to its handler and record it (with t translated from θ)."""
436 new_waves: list = []
438 if event.event_type == EventType.CHAR_CHAR_COLLISION:
439 new_waves = handle_characteristic_collision(
440 event.waves_involved[0], event.waves_involved[1], event.theta, event.location
441 )
443 elif event.event_type == EventType.SHOCK_SHOCK_COLLISION:
444 new_waves = handle_shock_collision(
445 event.waves_involved[0], event.waves_involved[1], event.theta, event.location
446 )
448 elif event.event_type == EventType.SHOCK_CHAR_COLLISION:
449 new_waves = handle_shock_characteristic_collision(
450 event.waves_involved[0], event.waves_involved[1], event.theta, event.location
451 )
453 elif event.event_type == EventType.RAREF_CHAR_COLLISION:
454 new_waves = handle_rarefaction_characteristic_collision(
455 event.waves_involved[0],
456 event.waves_involved[1],
457 event.theta,
458 event.location,
459 boundary_type=event.boundary_type,
460 )
462 elif event.event_type == EventType.SHOCK_RAREF_COLLISION:
463 new_waves = handle_shock_rarefaction_collision(
464 event.waves_involved[0],
465 event.waves_involved[1],
466 event.theta,
467 event.location,
468 boundary_type=event.boundary_type,
469 )
471 elif event.event_type == EventType.RAREF_RAREF_COLLISION:
472 # Conservative: rarefaction-rarefaction collision records the event
473 # but makes no topology change. Both rarefactions remain active.
474 new_waves = []
476 elif event.event_type == EventType.DSW_FAN_EXHAUSTED:
477 new_waves = self._handle_fan_exhaustion(event.waves_involved[0], event.theta, event.location)
479 elif event.event_type == EventType.OUTLET_CROSSING:
480 event_record = handle_outlet_crossing(event.waves_involved[0], event.theta, event.location)
481 self.state.events.append(event_record)
482 return
484 self.state.waves.extend(new_waves)
486 self.state.events.append({
487 "theta": event.theta,
488 "type": event.event_type.value,
489 "location": event.location,
490 "waves_before": event.waves_involved,
491 "waves_after": new_waves,
492 })
494 def _handle_fan_exhaustion(self, dsw: DecayingShockWave, theta_event: float, v_event: float) -> list[Wave]:
495 """Hand a fan-exhausted decaying shock off to a regular shock.
497 When ``c_decay`` reaches ``c_fan_tail`` the decaying side is no longer
498 fed by the fan; the wave continues as a constant-speed
499 ``ShockWave(c_fan_tail, c_fixed)`` (sides assigned per ``decay_side``).
500 The handoff is C1-continuous (``dV_s/dθ → S(c_fan_tail, c_fixed)`` = the
501 spawned shock's speed). The decaying shock is deactivated.
503 Returns
504 -------
505 list of Wave
506 ``[ShockWave]`` for the continuation, or ``[]`` if it fails entropy.
507 """
508 if dsw.decay_side == "left":
509 c_left, c_right = dsw.c_fan_tail, dsw.c_fixed
510 else:
511 c_left, c_right = dsw.c_fixed, dsw.c_fan_tail
513 dsw.deactivate(theta_event)
515 if abs(c_left - c_right) < EPSILON_CONCENTRATION:
516 # Fan decayed onto the fixed state — no discontinuity remains.
517 return []
519 shock = ShockWave(
520 theta_start=theta_event,
521 v_start=v_event,
522 c_left=c_left,
523 c_right=c_right,
524 sorption=self.state.sorption,
525 )
526 if not shock.satisfies_entropy():
527 return []
528 return [shock]
530 def run(self, max_iterations: int = 10000, *, verbose: bool = False):
531 """Process events in θ-order until the queue is empty or ``max_iterations`` is reached."""
532 iteration = 0
534 if verbose:
535 logger.info("Starting simulation at θ=%.3f", self.state.theta_current)
536 logger.info("Initial waves: %d", len(self.state.waves))
537 logger.info("First arrival: θ=%.3f", self.theta_first_arrival)
539 while iteration < max_iterations:
540 event = self.find_next_event()
542 if event is None:
543 if verbose:
544 logger.info("Simulation complete after %d events at θ=%.6f", iteration, self.state.theta_current)
545 break
547 self.state.theta_current = event.theta
549 try:
550 self.handle_event(event)
551 except Exception:
552 logger.exception("Error handling event at θ=%.3f", event.theta)
553 raise
555 if iteration % 100 == 0:
556 self.verify_physics()
558 if verbose and iteration % 10 == 0:
559 active = sum(1 for w in self.state.waves if w.is_active)
560 logger.debug("Iteration %d: θ=%.3f, active_waves=%d", iteration, event.theta, active)
562 iteration += 1
564 if iteration >= max_iterations:
565 logger.warning("Reached max_iterations=%d", max_iterations)
567 if verbose:
568 logger.info("Final statistics:")
569 logger.info(" Total events: %d", len(self.state.events))
570 logger.info(" Total waves created: %d", len(self.state.waves))
571 logger.info(" Active waves: %d", sum(1 for w in self.state.waves if w.is_active))
572 logger.info(" First arrival: θ=%.6f", self.theta_first_arrival)
574 def verify_physics(self):
575 """Verify physical correctness: every active shock satisfies Lax entropy.
577 Mass conservation is intentionally NOT checked here. The closed-form
578 identity ``m_out(θ) = m_in(θ) − m_dom(θ)`` makes any runtime
579 ``m_in_domain + m_out_cumulative == m_in_cumulative`` test tautological
580 (residual identically zero, regardless of any ``compute_domain_mass``
581 bug), so it cannot catch a conservation error. The non-tautological,
582 integral-based conservation check (an independent breakthrough integral
583 compared to the inlet mass) lives in
584 :func:`gwtransport.fronttracking.validation.verify_physics` check 7 and
585 is exercised by ``TestEndToEndConservation`` /
586 ``TestIndependentDomainMass``.
588 Raises
589 ------
590 RuntimeError
591 If an active shock violates the Lax entropy condition.
592 """
593 for wave in self.state.waves:
594 if isinstance(wave, ShockWave) and wave.is_active and not wave.satisfies_entropy():
595 msg = (
596 f"Shock at θ_start={wave.theta_start:.3f} violates entropy! "
597 f"c_left={wave.c_left:.3f}, c_right={wave.c_right:.3f}, "
598 f"speed={wave.speed:.6g}"
599 )
600 raise RuntimeError(msg)
603def find_unresolved_interaction(state: FrontTrackerState) -> str | None:
604 """Locate an unresolved wave–wave interaction inside the transport domain.
606 The event-driven solver resolves shock↔shock, shock↔characteristic and
607 shock↔rarefaction collisions (the last into a :class:`DecayingShockWave`), but it
608 never collides anything *with* a decaying shock, nor composes two fans that come to
609 occupy the same region. Such an unresolved interaction leaves the non-interacting wave
610 objects overlapping, so the exact nonlinear multi-front field is not represented — the
611 public ``cout`` degrades to a spurious linear superposition of a nonlinear operator
612 (mass-fabricating once the reader clamps the resulting negative bins to zero). This
613 detects the first offending interaction so the public API can refuse the input rather
614 than return a wrong, non-conservative answer. Two complementary detectors run over the
615 input θ-window ``(0, theta_edges[-1]]``:
617 1. **Geometric fan overlap.** When two or more fan-bearing waves (rarefactions /
618 decaying shocks) cover a common point inside ``[0, v_outlet]``, their composite
619 field is wrong even when total mass happens to be conserved (a positive-but-wrong
620 ``cout`` with no negative bin — e.g. two decaying shocks with a zero fan tail, or a
621 later pulse's fan sweeping an earlier one). A symptom-only proxy cannot see this
622 class, so the geometric scan is a necessary complement.
624 2. **Conservation symptom.** The cumulative outlet mass ``m_out(θ) = m_in(θ) − m_dom(θ)``
625 must be non-decreasing in θ (mass exits the column, it never re-enters). A decrease
626 beyond the FP-cancellation band means the reader's domain-mass field transiently
627 over-counts stored mass — the fingerprint of a shock overtaking another shock /
628 rarefaction / decaying-shock fan. The two fans need NOT share an in-domain point
629 (they may only cross beyond ``v_outlet``), so the geometric scan misses this dominant
630 multi-pulse class; the monotonicity check catches it.
632 Parameters
633 ----------
634 state : FrontTrackerState
635 Completed simulation state (after :meth:`FrontTracker.run`).
637 Returns
638 -------
639 str or None
640 A short description (position/θ and mechanism) of the first offending interaction,
641 or ``None`` when the solution is a clean single-front / well-separated superposition
642 that the reader represents exactly.
644 Notes
645 -----
646 Both scans stay strictly inside the inlet θ-window, so the benign out-of-window
647 saturation clamp (a single-front run whose output bins extend past the last injected
648 mass) does not trip the symptom check. Well-separated pulses that clear a short column
649 before overtaking one another (their fans never share an in-domain point and their
650 cumulative outflow stays monotone) are correctly accepted.
651 """
652 waves = state.waves
653 v_outlet = state.v_outlet
654 theta_hi = float(state.theta_edges[-1])
655 v_tol = 1e-9 * max(v_outlet, 1.0)
656 thetas = np.linspace(theta_hi / 400.0, theta_hi, 400)
658 # (1) Geometric fan overlap.
659 def fan_covers(wave: Wave, v: float, theta: float) -> bool:
660 """Whether ``wave``'s self-similar fan geometrically spans ``v`` at ``theta``."""
661 if isinstance(wave, DecayingShockWave):
662 v_s = wave.position_at_theta(theta)
663 if v_s is None or v == wave.v_origin:
664 return False
665 return (wave.decay_side == "left" and v < v_s) or (wave.decay_side == "right" and v > v_s)
666 return isinstance(wave, RarefactionWave) and wave.contains_point(v, theta)
668 for theta in thetas:
669 fans = [w for w in waves if isinstance(w, (DecayingShockWave, RarefactionWave)) and w.was_active_at(theta)]
670 if len(fans) <= 1: # a single fan (or none) cannot overlap another
671 continue
672 boundaries = {0.0, v_outlet}
673 for wave in waves:
674 if not wave.was_active_at(theta):
675 continue
676 if isinstance(wave, RarefactionWave):
677 candidates = (wave.head_position_at_theta(theta), wave.tail_position_at_theta(theta))
678 else:
679 candidates = (wave.position_at_theta(theta),)
680 boundaries.update(p for p in candidates if p is not None and 0.0 < p < v_outlet)
681 for v_lo, v_hi in pairwise(sorted(boundaries)):
682 if v_hi - v_lo < v_tol:
683 continue
684 v_mid = 0.5 * (v_lo + v_hi)
685 if len([w for w in fans if fan_covers(w, v_mid, theta)]) > 1: # ≥ 2 fans share this point
686 return f"two fans overlap at V={v_mid:.4g} m³ (θ={theta:.4g} m³)"
688 # (2) Conservation symptom: cumulative outlet mass must be monotone non-decreasing.
689 # The band mirrors ``compute_bin_averaged_concentration_exact``'s FP-cancellation clamp
690 # (same constant, same ``max(scale, 1.0)`` floor) scaled to the total injected mass, so
691 # pre-breakthrough cancellation dust stays silent while a genuine over-count (orders of
692 # magnitude above the band) refuses the input.
693 grid = np.concatenate(([0.0], thetas))
694 m_out = np.array([
695 compute_cumulative_outlet_mass(
696 float(theta), v_outlet, waves, state.sorption, cin=state.cin, theta_edges=state.theta_edges
697 )
698 for theta in grid
699 ])
700 m_in_total = float(np.sum(state.cin * np.diff(state.theta_edges)))
701 band = FP_CANCELLATION_CLAMP * np.finfo(float).eps * max(m_in_total, 1.0)
702 steps = np.diff(m_out)
703 worst = int(np.argmin(steps))
704 if steps[worst] < -band:
705 return (
706 f"cumulative outlet mass drops by {-steps[worst]:.4g} near θ={grid[worst + 1]:.4g} m³ "
707 "(a shock overtakes another shock / rarefaction / decaying-shock fan; the reader "
708 "over-counts stored mass there)"
709 )
710 return None