Coverage for src/gwtransport/fronttracking/solver.py: 0%
248 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 20:54 +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.
20Available functions:
22- :class:`FrontTrackerState` - Dataclass carrying the whole simulation state: every wave ever created
23 (active and deactivated), the event log, the current θ, the outlet position ``v_outlet``, the sorption
24 model, the inlet series, and the ``tedges_days``/``theta_edges`` grids. Its ``t_at_theta``, ``theta_at_t``
25 and ``theta_at_t_array`` methods are the piecewise-linear θ↔t maps a caller needs to read any solver output
26 in user-facing days; events at a zero-flow θ plateau map to the end of that plateau.
28- :class:`FrontTracker` - The event-driven solver. Construction validates the inputs (non-negative ``cin``
29 and ``flow``, ``len(tedges) == len(cin) + 1``, positive pore volume), builds ``theta_edges`` by cumulating
30 ``flow · dt``, sets the resolution horizon ``theta_horizon`` (default: the last inlet edge), computes
31 ``theta_first_arrival``, and emits the inlet waves. ``find_next_event`` returns the earliest candidate in
32 θ — characteristic, shock and rarefaction collisions, face merges involving a decaying or doubly-fed shock,
33 fan exhaustion, and outlet crossings — ``handle_event`` dispatches it to the handler that retires the
34 parents and appends the successors, and ``run`` repeats until the candidate set is empty or
35 ``max_iterations`` is reached. ``verify_physics`` asserts that every active shock still satisfies the Lax
36 entropy condition.
38- :func:`find_unresolved_interaction` - Invariant tripwire over a completed state. The cumulative outlet mass
39 ``m_out(θ) = m_in(θ) − m_dom(θ)`` must never decrease; this returns a short description of the first θ
40 where it does beyond the floating-point cancellation band — the fingerprint of an interaction the solver
41 failed to resolve — and ``None`` for a mass-conserving wave field.
42"""
44import logging
45from dataclasses import dataclass
46from operator import itemgetter
48import numpy as np
49import numpy.typing as npt
50import pandas as pd
52from gwtransport._time import tedges_to_days
53from gwtransport.fronttracking.events import (
54 Event,
55 EventType,
56 find_characteristic_intersection,
57 find_outlet_crossing,
58 find_rarefaction_boundary_intersections,
59 find_shock_characteristic_intersection,
60 find_shock_shock_intersection,
61 is_outlet_crossing_pinned,
62)
63from gwtransport.fronttracking.handlers import (
64 EPSILON_CONCENTRATION,
65 create_inlet_waves_at_theta,
66 handle_characteristic_collision,
67 handle_rarefaction_characteristic_collision,
68 handle_shock_characteristic_collision,
69 handle_shock_collision,
70 handle_shock_rarefaction_collision,
71)
72from gwtransport.fronttracking.interactions import (
73 find_face_crossing,
74 iter_faces,
75 resolve_merge,
76)
77from gwtransport.fronttracking.math import (
78 SorptionModel,
79 compute_first_front_arrival_theta,
80)
81from gwtransport.fronttracking.output import (
82 FP_CANCELLATION_CLAMP,
83 compute_cumulative_outlet_mass,
84)
85from gwtransport.fronttracking.waves import (
86 CharacteristicWave,
87 DecayingShockWave,
88 DoubleFanShockWave,
89 RarefactionWave,
90 ShockWave,
91 Wave,
92)
94logger = logging.getLogger(__name__)
97@dataclass
98class FrontTrackerState:
99 """Complete state of the front-tracking simulation in (V, θ).
101 Parameters
102 ----------
103 waves : list of Wave
104 All waves created during simulation (includes inactive waves).
105 events : list of dict
106 Event history. Records use the ``"theta"`` key carrying the
107 cumulative flow at which the event occurred [m³]. Callers translate
108 to user-facing time via ``FrontTrackerState.t_at_theta``.
109 theta_current : float
110 Current simulation cumulative flow [m³].
111 v_outlet : float
112 Outlet position [m³].
113 sorption : SorptionModel
114 Sorption model.
115 cin : numpy.ndarray
116 Inlet concentration time series [mass/volume].
117 flow : numpy.ndarray
118 Flow rate time series [m³/day], one value per bin.
119 tedges : pandas.DatetimeIndex
120 Time bin edges.
121 tedges_days : numpy.ndarray
122 ``tedges`` as days from ``tedges[0]``, length ``len(flow) + 1``.
123 theta_edges : numpy.ndarray
124 Cumulative flow at every bin edge. ``theta_edges[i] = sum_{k<i} flow[k] *
125 (tedges_days[k+1] - tedges_days[k])``. Length ``len(flow) + 1``.
126 """
128 waves: list[Wave]
129 events: list[dict]
130 theta_current: float
131 v_outlet: float
132 sorption: SorptionModel
133 cin: np.ndarray
134 flow: np.ndarray
135 tedges: pd.DatetimeIndex
136 tedges_days: npt.NDArray[np.floating]
137 theta_edges: npt.NDArray[np.floating]
138 theta_horizon: float = float("inf")
140 def t_at_theta(self, theta: float) -> float:
141 """Translate cumulative flow θ back to user-facing time t [days].
143 Piecewise linear inversion of the (tedges_days → theta_edges) map.
144 Implementation note on the (rare) zero-flow case: when a bin has
145 ``flow[i] == 0``, θ is constant across ``[tedges_days[i], tedges_days[i+1])``;
146 ``np.searchsorted(..., side='right') - 1`` lands on the rightmost
147 such bin, so this function returns ``tedges_days[i]`` for the
148 right-most i sharing that θ. Convention: events at a zero-flow θ
149 plateau map to the END of the zero-flow interval.
150 """
151 if theta <= self.theta_edges[0]:
152 return float(self.tedges_days[0])
153 if theta >= self.theta_edges[-1]:
154 # Extrapolate: pretend the last bin continues at its current flow.
155 last_flow = float(self.flow[-1])
156 if last_flow > 0:
157 return float(self.tedges_days[-1] + (theta - self.theta_edges[-1]) / last_flow)
158 return float(self.tedges_days[-1])
160 # Find bin index i with theta_edges[i] <= theta < theta_edges[i+1].
161 # np.searchsorted with side='right' returns the smallest i+1 such that
162 # theta_edges[i+1] > theta, so subtracting 1 gives i.
163 # The boundary returns above guarantee strictly-interior theta here, so
164 # searchsorted lands in [0, len(flow)-1] without an extra index clamp.
165 i = int(np.searchsorted(self.theta_edges, theta, side="right")) - 1
166 flow_i = float(self.flow[i])
167 if flow_i <= 0:
168 return float(self.tedges_days[i])
169 return float(self.tedges_days[i] + (theta - self.theta_edges[i]) / flow_i)
171 def theta_at_t(self, t: float) -> float:
172 """Translate user-facing time t [days] to cumulative flow θ [m³].
174 Scalar form of :meth:`theta_at_t_array`.
175 """
176 return float(self.theta_at_t_array(t))
178 def theta_at_t_array(self, t: npt.ArrayLike) -> npt.NDArray[np.floating]:
179 """Map times t [days] to cumulative flow θ [m³].
181 Piecewise linear forward map. Outside the input range the boundary flow is
182 extrapolated.
184 Parameters
185 ----------
186 t : array-like
187 User-facing time points [days].
189 Returns
190 -------
191 ndarray
192 Cumulative flow θ at each ``t`` [m³].
193 """
194 t_arr = np.asarray(t, dtype=float)
195 # Interior map: i = searchsorted(tedges_days, t, 'right') - 1, clipped to
196 # a valid bin so the boundary branches can overwrite the extrapolated tails.
197 i = np.clip(np.searchsorted(self.tedges_days, t_arr, side="right") - 1, 0, len(self.flow) - 1)
198 theta = self.theta_edges[i] + (t_arr - self.tedges_days[i]) * self.flow[i]
199 # Left of the first edge: clamp to theta_edges[0].
200 theta = np.where(t_arr <= self.tedges_days[0], self.theta_edges[0], theta)
201 # Right of the last edge: extrapolate at the final-bin flow.
202 return np.where(
203 t_arr >= self.tedges_days[-1],
204 self.theta_edges[-1] + (t_arr - self.tedges_days[-1]) * float(self.flow[-1]),
205 theta,
206 )
209class FrontTracker:
210 """Event-driven front-tracking solver for nonlinear sorption transport.
212 Parameters
213 ----------
214 cin : numpy.ndarray
215 Inlet concentration time series [mass/volume]; length ``n``.
216 flow : numpy.ndarray
217 Flow rate time series [m³/day]; length ``n`` (one value per bin).
218 tedges : pandas.DatetimeIndex
219 Time bin edges (length ``n+1``).
220 aquifer_pore_volume : float
221 Total pore volume [m³] — used as the outlet position.
222 sorption : SorptionModel
223 Sorption model.
225 Attributes
226 ----------
227 state : FrontTrackerState
228 Complete simulation state.
229 theta_first_arrival : float
230 Cumulative flow θ at which the first nonzero-concentration wave reaches
231 the outlet [m³]. Translate to user-facing time via
232 ``state.t_at_theta(theta_first_arrival)``.
234 Notes
235 -----
236 The solver works exclusively in cumulative flow θ; events appended to
237 ``state.events`` carry ``"theta"``. Translation to user-facing time t is
238 the caller's responsibility (use ``state.t_at_theta``).
239 """
241 def __init__(
242 self,
243 cin: npt.ArrayLike,
244 flow: npt.ArrayLike,
245 tedges: pd.DatetimeIndex,
246 aquifer_pore_volume: float,
247 sorption: SorptionModel,
248 theta_horizon: float | None = None,
249 ):
250 cin = np.asarray(cin, dtype=float)
251 flow = np.asarray(flow, dtype=float)
252 if len(tedges) != len(cin) + 1:
253 msg = f"tedges must have length len(cin) + 1, got {len(tedges)} vs {len(cin) + 1}"
254 raise ValueError(msg)
255 if len(flow) != len(cin):
256 msg = f"flow must have same length as cin, got {len(flow)} vs {len(cin)}"
257 raise ValueError(msg)
258 if np.any(cin < 0):
259 msg = "cin must be non-negative"
260 raise ValueError(msg)
261 if np.any(flow < 0):
262 msg = "flow must be non-negative (negative flow not supported)"
263 raise ValueError(msg)
264 if aquifer_pore_volume <= 0:
265 msg = "aquifer_pore_volume must be positive"
266 raise ValueError(msg)
268 tedges_days = tedges_to_days(tedges)
269 dt_days = np.diff(tedges_days)
270 bin_volumes = flow * dt_days
271 theta_edges = np.concatenate(([0.0], np.cumsum(bin_volumes)))
273 # Global interaction resolution runs to a θ-horizon: an unresolved crossing at
274 # θ > horizon cannot affect any reader query at θ ≤ horizon (information propagates
275 # strictly downstream, all speeds ≥ 0). Default = the last inlet edge; the public
276 # API passes the last requested output edge so beyond-outlet merges that still feed
277 # an in-window outlet query are resolved.
278 resolved_horizon = float(theta_edges[-1]) if theta_horizon is None else float(theta_horizon)
280 self.state = FrontTrackerState(
281 waves=[],
282 events=[],
283 theta_current=0.0,
284 v_outlet=aquifer_pore_volume,
285 sorption=sorption,
286 cin=cin,
287 flow=flow,
288 tedges=tedges,
289 tedges_days=tedges_days,
290 theta_edges=theta_edges,
291 theta_horizon=resolved_horizon,
292 )
294 self.theta_first_arrival = compute_first_front_arrival_theta(cin, theta_edges, aquifer_pore_volume, sorption)
296 self._initialize_inlet_waves()
298 def _initialize_inlet_waves(self):
299 """Emit one wave per nonzero inlet step at the corresponding ``theta_edges[i]``."""
300 c_prev = 0.0
301 theta_edges = self.state.theta_edges
303 for i in range(len(self.state.cin)):
304 # A zero-flow (pump-off) bin carries no water: its θ-interval has zero
305 # width, so it is transport-invisible. Skipping it carries the inlet
306 # step through the gap (effective transition c_before_gap → c_after_gap)
307 # instead of emitting coincident spurious waves at the degenerate θ.
308 if theta_edges[i + 1] <= theta_edges[i]:
309 continue
311 c_new = float(self.state.cin[i])
313 if abs(c_new - c_prev) > EPSILON_CONCENTRATION:
314 new_waves = create_inlet_waves_at_theta(
315 c_prev=c_prev,
316 c_new=c_new,
317 theta=float(theta_edges[i]),
318 sorption=self.state.sorption,
319 )
320 self.state.waves.extend(new_waves)
322 c_prev = c_new
324 def find_next_event(self) -> Event | None:
325 """Return the next event in θ-order, or ``None`` if none."""
326 # Each call collects every candidate and selects the single earliest via
327 # ``min`` over the (theta, counter, ...) tuples; the unique ``counter``
328 # breaks θ-ties deterministically (and stops comparison before the
329 # non-orderable EventType/Wave fields). A heap would only pay for its
330 # invariant to extract one minimum, so a flat list + ``min`` is leaner.
331 candidates: list[tuple] = []
332 counter = 0 # Unique counter to break θ-ties deterministically
334 active_waves = [w for w in self.state.waves if w.is_active]
335 theta_current = self.state.theta_current
337 # Defense-in-depth loop guard: reject collision candidates at or below
338 # θ_current + tol so a degenerate geometry cannot re-fire the identical
339 # event forever (same FP-tolerance scale as the outlet-crossing guard
340 # below). The structural fix is routing all shock↔rarefaction collisions
341 # through DecayingShockWave; this backstops any future near-coincident
342 # event from looping the solver.
343 collision_tol = 1e-12 * max(abs(theta_current), 1.0)
345 def push_collision(theta, event_type, waves, v, boundary):
346 nonlocal counter
347 # Global resolution: collisions are resolved wherever they occur,
348 # up to θ_horizon — not only within [0, v_outlet]. The single-owner reader needs a
349 # globally interaction-consistent front order, since the nearest downstream face
350 # of an in-domain query may lie beyond the outlet (e.g. a rarefaction head that
351 # catches a shock past v_outlet still bounds the in-domain fan). Information
352 # propagates strictly downstream (all speeds ≥ 0), so an event at θ > θ_horizon
353 # cannot affect any reader query at θ ≤ θ_horizon.
354 if v < -collision_tol:
355 return
356 if theta <= theta_current + collision_tol or theta > self.state.theta_horizon:
357 return
358 candidates.append((theta, counter, event_type, waves, v, boundary, None))
359 counter += 1
361 chars = [w for w in active_waves if isinstance(w, CharacteristicWave)]
362 for i, w1 in enumerate(chars):
363 for w2 in chars[i + 1 :]:
364 result = find_characteristic_intersection(w1, w2, theta_current)
365 if result:
366 theta, v = result
367 push_collision(theta, EventType.CHAR_CHAR_COLLISION, [w1, w2], v, None)
369 shocks = [w for w in active_waves if isinstance(w, ShockWave)]
370 for i, w1 in enumerate(shocks):
371 for w2 in shocks[i + 1 :]:
372 result = find_shock_shock_intersection(w1, w2, theta_current)
373 if result:
374 theta, v = result
375 push_collision(theta, EventType.SHOCK_SHOCK_COLLISION, [w1, w2], v, None)
377 for shock in shocks:
378 for char in chars:
379 result = find_shock_characteristic_intersection(shock, char, theta_current)
380 if result:
381 theta, v = result
382 push_collision(theta, EventType.SHOCK_CHAR_COLLISION, [shock, char], v, None)
384 rarefs = [w for w in active_waves if isinstance(w, RarefactionWave)]
385 for raref in rarefs:
386 for char in chars:
387 intersections = find_rarefaction_boundary_intersections(raref, char, theta_current)
388 for theta, v, boundary in intersections:
389 push_collision(theta, EventType.RAREF_CHAR_COLLISION, [raref, char], v, boundary)
391 for shock in shocks:
392 for raref in rarefs:
393 intersections = find_rarefaction_boundary_intersections(raref, shock, theta_current)
394 for theta, v, boundary in intersections:
395 push_collision(theta, EventType.SHOCK_RAREF_COLLISION, [shock, raref], v, boundary)
397 # Interaction events: every face pair with at least one decaying/doubly-fed shock. The
398 # closed-form loops above cover char/shock/rarefaction pairs; here a face of a
399 # DSW/DFSW (its curved shock face or a free fan boundary line) meets any other wave's face.
400 # Resolution is GLOBAL — allowed beyond v_outlet up to θ_horizon — because the sweep reader
401 # needs a globally interaction-consistent front order for the nearest-downstream lookup.
402 special_types = (DecayingShockWave, DoubleFanShockWave)
403 if any(isinstance(w, special_types) for w in active_waves):
404 all_faces = [face for wave in active_waves for face in iter_faces(wave, theta_current)]
405 theta_horizon = self.state.theta_horizon
406 for i, fa in enumerate(all_faces):
407 for fb in all_faces[i + 1 :]:
408 if not (isinstance(fa.wave, special_types) or isinstance(fb.wave, special_types)):
409 continue
410 same_wave = fa.wave is fb.wave
411 if same_wave and not (
412 isinstance(fa.wave, DoubleFanShockWave) and {fa.role, fb.role} == {"shock", "boundary"}
413 ):
414 # Same-wave crossings other than a doubly-fed shock meeting its own fan
415 # boundary line never occur (a DSW's fan exhaustion has its own exact
416 # detector; two boundary lines of one wave diverge).
417 continue
418 # A doubly-fed shock crossing its OWN fan boundary line is that side's
419 # exhaustion: the fan between them is spent, and the universal merge
420 # (const far-bound feeder | surviving fan feeder) degrades the wave to a
421 # DecayingShockWave while retiring the boundary — uniformly with every
422 # other face merge. Entropy (λ(c_L) ≥ σ ≥ λ(c_R)) makes these crossings
423 # the only finite-θ side endings: the left fan's slow (upstream-far) char
424 # can only catch the shock from behind, the shock can only catch the
425 # right fan's downstream-far char — exactly the two exposed lines.
426 crossing = find_face_crossing(fa, fb, theta_current, theta_horizon)
427 if crossing is None:
428 continue
429 theta, v = crossing
430 if v < -collision_tol or theta > self.state.theta_horizon:
431 continue
432 # A born-coincident merge is reported at theta == theta_current (a wave
433 # created this step is already touching/past a neighbour it must merge
434 # with). Admit it — resolve_merge deactivates both parents, so it cannot
435 # re-fire — but keep the strict-future filter for every ordinary crossing
436 # so a just-processed event is not re-emitted one ULP later.
437 born_coincident = theta <= theta_current + collision_tol and (
438 abs(fa.wave.theta_start - theta_current) <= collision_tol
439 or abs(fb.wave.theta_start - theta_current) <= collision_tol
440 )
441 if theta <= theta_current + collision_tol and not born_coincident:
442 continue
443 candidates.append((theta, counter, EventType.WAVE_MERGE, [fa.wave, fb.wave], v, None, (fa, fb)))
444 counter += 1
446 # Fan-exhaustion: a DecayingShockWave is valid only while c_decay stays
447 # above c_fan_tail. When c_decay reaches c_fan_tail the fan is spent and
448 # the wave hands off to a regular ShockWave(c_fan_tail, c_fixed).
449 for wave in active_waves:
450 if isinstance(wave, DecayingShockWave):
451 theta_exhaust = wave.theta_at_fan_exhaustion()
452 if theta_exhaust is None or theta_exhaust <= theta_current + collision_tol:
453 continue
454 if theta_exhaust > self.state.theta_horizon:
455 continue
456 v_exhaust = wave.position_at_theta(theta_exhaust)
457 if v_exhaust is None or v_exhaust < -collision_tol:
458 continue
459 candidates.append((theta_exhaust, counter, EventType.DSW_FAN_EXHAUSTED, [wave], v_exhaust, None, None))
460 counter += 1
462 v_outlet = self.state.v_outlet
463 # Same FP-tolerance discipline as events.find_outlet_crossing — prevents
464 # re-emitting an outlet crossing for a boundary that's at v_outlet ± ULPs.
465 outlet_tol = 1e-12 * max(abs(v_outlet), 1.0)
467 for wave in active_waves:
468 if isinstance(wave, RarefactionWave):
469 theta_eval = max(theta_current, wave.theta_start)
470 for c_boundary, pos_fn, speed_fn in (
471 (wave.c_head, wave.head_position_at_theta, wave.head_speed),
472 (wave.c_tail, wave.tail_position_at_theta, wave.tail_speed),
473 ):
474 v_pos = pos_fn(theta_eval)
475 if v_pos is None or v_pos >= v_outlet - outlet_tol:
476 continue
477 s = speed_fn()
478 # Skip a c_min-floored (pinned) boundary — R(c_min) inflated
479 # for n>1, c→0: its crossing lands at a non-physical θ~1e8 and
480 # only pollutes the event record.
481 if s <= 0 or is_outlet_crossing_pinned(c_boundary, wave.sorption):
482 continue
483 theta_cross = theta_eval + (v_outlet - v_pos) / s
484 if theta_cross > theta_current:
485 candidates.append((
486 theta_cross,
487 counter,
488 EventType.OUTLET_CROSSING,
489 [wave],
490 v_outlet,
491 None,
492 None,
493 ))
494 counter += 1
495 else:
496 theta_cross = find_outlet_crossing(wave, self.state.v_outlet, theta_current)
497 if theta_cross and theta_cross > theta_current:
498 candidates.append((
499 theta_cross,
500 counter,
501 EventType.OUTLET_CROSSING,
502 [wave],
503 self.state.v_outlet,
504 None,
505 None,
506 ))
507 counter += 1
509 if not candidates:
510 return None
512 theta_event, _, event_type, waves, v, extra, faces = min(candidates, key=itemgetter(slice(2)))
514 raref_types = {EventType.RAREF_CHAR_COLLISION, EventType.SHOCK_RAREF_COLLISION}
515 boundary_type = extra if event_type in raref_types else None
517 return Event(
518 theta=theta_event,
519 event_type=event_type,
520 waves_involved=waves,
521 location=v,
522 boundary_type=boundary_type,
523 faces=faces,
524 )
526 def handle_event(self, event: Event):
527 """Dispatch an event to its handler and record it (with t translated from θ)."""
528 new_waves: list = []
530 if event.event_type == EventType.CHAR_CHAR_COLLISION:
531 new_waves = handle_characteristic_collision(
532 event.waves_involved[0], event.waves_involved[1], event.theta, event.location
533 )
535 elif event.event_type == EventType.SHOCK_SHOCK_COLLISION:
536 new_waves = handle_shock_collision(
537 event.waves_involved[0], event.waves_involved[1], event.theta, event.location
538 )
540 elif event.event_type == EventType.SHOCK_CHAR_COLLISION:
541 new_waves = handle_shock_characteristic_collision(
542 event.waves_involved[0], event.waves_involved[1], event.theta, event.location
543 )
545 elif event.event_type == EventType.RAREF_CHAR_COLLISION:
546 new_waves = handle_rarefaction_characteristic_collision(
547 event.waves_involved[0],
548 event.waves_involved[1],
549 event.theta,
550 event.location,
551 boundary_type=event.boundary_type,
552 )
554 elif event.event_type == EventType.SHOCK_RAREF_COLLISION:
555 new_waves = handle_shock_rarefaction_collision(
556 event.waves_involved[0],
557 event.waves_involved[1],
558 event.theta,
559 event.location,
560 boundary_type=event.boundary_type,
561 )
563 elif event.event_type == EventType.WAVE_MERGE:
564 assert event.faces is not None # noqa: S101 # WAVE_MERGE always carries its two faces
565 face_a, face_b = event.faces
566 new_waves = resolve_merge(face_a, face_b, event.theta, event.location, self.state.sorption)
568 elif event.event_type == EventType.DSW_FAN_EXHAUSTED:
569 new_waves = self._handle_fan_exhaustion(event.waves_involved[0], event.theta, event.location)
571 elif event.event_type == EventType.OUTLET_CROSSING:
572 # The wave is NOT deactivated: it stays queryable for concentrations
573 # between its origin and the outlet.
574 wave = event.waves_involved[0]
575 self.state.events.append({
576 "theta": event.theta,
577 "type": "outlet_crossing",
578 "wave": wave,
579 "location": event.location,
580 "concentration_left": wave.concentration_left(),
581 "concentration_right": wave.concentration_right(),
582 })
583 return
585 self.state.waves.extend(new_waves)
587 self.state.events.append({
588 "theta": event.theta,
589 "type": event.event_type.value,
590 "location": event.location,
591 "waves_before": event.waves_involved,
592 "waves_after": new_waves,
593 })
595 def _handle_fan_exhaustion(self, dsw: DecayingShockWave, theta_event: float, v_event: float) -> list[Wave]:
596 """Hand a fan-exhausted decaying shock off to a regular shock.
598 When ``c_decay`` reaches ``c_fan_tail`` the decaying side is no longer
599 fed by the fan; the wave continues as a constant-speed
600 ``ShockWave(c_fan_tail, c_fixed)`` (sides assigned per ``decay_side``).
601 The handoff is C1-continuous (``dV_s/dθ → S(c_fan_tail, c_fixed)`` = the
602 spawned shock's speed). The decaying shock is deactivated.
604 Returns
605 -------
606 list of Wave
607 ``[ShockWave]`` for the continuation, or ``[]`` if it fails entropy.
608 """
609 if dsw.decay_side == "left":
610 c_left, c_right = dsw.c_fan_tail, dsw.c_fixed
611 else:
612 c_left, c_right = dsw.c_fixed, dsw.c_fan_tail
614 dsw.deactivate(theta_event)
616 if abs(c_left - c_right) < EPSILON_CONCENTRATION:
617 # Fan decayed onto the fixed state — no discontinuity remains.
618 return []
620 shock = ShockWave(
621 theta_start=theta_event,
622 v_start=v_event,
623 c_left=c_left,
624 c_right=c_right,
625 sorption=self.state.sorption,
626 )
627 if not shock.satisfies_entropy():
628 return []
629 return [shock]
631 def run(self, max_iterations: int = 10000):
632 """Process events in θ-order until the queue is empty or ``max_iterations`` is reached."""
633 iteration = 0
635 while iteration < max_iterations:
636 event = self.find_next_event()
638 if event is None:
639 break
641 self.state.theta_current = event.theta
642 self.handle_event(event)
643 iteration += 1
645 if iteration >= max_iterations:
646 logger.warning("Reached max_iterations=%d", max_iterations)
648 def verify_physics(self):
649 """Verify physical correctness: every active shock satisfies Lax entropy.
651 Mass conservation is intentionally NOT checked here. The closed-form
652 identity ``m_out(θ) = m_in(θ) − m_dom(θ)`` makes any
653 ``m_in_domain + m_out_cumulative == m_in_cumulative`` test tautological
654 (residual identically zero, regardless of any ``compute_domain_mass``
655 bug), so it cannot catch a conservation error. The non-tautological,
656 integral-based conservation check (an independent breakthrough integral
657 compared to the inlet mass) lives in
658 :func:`gwtransport.fronttracking.validation.verify_physics` check 7.
660 Every shock the solver creates is entropy-checked at construction and a
661 shock's ``(c_left, c_right)`` never change, so this scan is an assertion
662 for externally assembled wave lists, not a solver self-check.
664 Raises
665 ------
666 RuntimeError
667 If an active shock violates the Lax entropy condition.
668 """
669 for wave in self.state.waves:
670 if isinstance(wave, ShockWave) and wave.is_active and not wave.satisfies_entropy():
671 msg = (
672 f"Shock at θ_start={wave.theta_start:.3f} violates entropy! "
673 f"c_left={wave.c_left:.3f}, c_right={wave.c_right:.3f}, "
674 f"speed={wave.speed:.6g}"
675 )
676 raise RuntimeError(msg)
679def find_unresolved_interaction(state: FrontTrackerState) -> str | None:
680 """Tripwire for a solver-left inconsistency in the resolved wave field.
682 The solver resolves *every* wave interaction (shock↔shock, fan-entry, doubly-fed
683 formation, same-apex annihilation, and their compositions), so the wave list is
684 interaction-consistent and the single-owner sweep reader is exact — overlapping fans are
685 normal and correct, not an error. This is an internal invariant check: the cumulative
686 outlet mass ``m_out(θ) = m_in(θ) − m_dom(θ)`` must be non-decreasing in θ (mass leaves
687 the column, it never re-enters). A decrease beyond the FP-cancellation band means the
688 reader's domain-mass field transiently over-counts stored mass — the fingerprint of an
689 interaction the solver failed to resolve (a bug). The public API turns a non-``None``
690 return into a fail-loud ``RuntimeError`` rather than returning a silently wrong ``cout``.
692 Parameters
693 ----------
694 state : FrontTrackerState
695 Completed simulation state (after :meth:`FrontTracker.run`).
697 Returns
698 -------
699 str or None
700 A short description (θ and mechanism) of the first monotonicity violation, or
701 ``None`` when the resolved field conserves mass (the normal case).
703 Notes
704 -----
705 The scan stays strictly inside the inlet θ-window, so the benign out-of-window saturation
706 clamp (a run whose output bins extend past the last injected mass) does not trip it.
707 """
708 waves = state.waves
709 v_outlet = state.v_outlet
710 theta_hi = float(state.theta_edges[-1])
711 thetas = np.linspace(theta_hi / 400.0, theta_hi, 400)
713 # Conservation symptom: cumulative outlet mass must be monotone non-decreasing. The band
714 # mirrors ``compute_bin_averaged_concentration_exact``'s FP-cancellation clamp (same
715 # constant, same ``max(scale, 1.0)`` floor) scaled to the total injected mass, so
716 # pre-breakthrough cancellation dust stays silent while a genuine over-count (orders of
717 # magnitude above the band) is caught.
718 grid = np.concatenate(([0.0], thetas))
719 m_out = np.array([
720 compute_cumulative_outlet_mass(
721 float(theta), v_outlet, waves, state.sorption, cin=state.cin, theta_edges=state.theta_edges
722 )
723 for theta in grid
724 ])
725 m_in_total = float(np.sum(state.cin * np.diff(state.theta_edges)))
726 # Tolerance band: the larger of the FP-cancellation floor and a 1e-6 relative transient
727 # allowance. The latter absorbs the benign ``c_min``-floor artifact at a fan-entry (the
728 # decaying side is born at ``c ≈ 1e-12`` rather than exactly 0, perturbing the near-apex
729 # fan integral by ~1e-3 for one θ-sample). A genuine unresolved interaction over-counts by
730 # O(a pulse's mass) — orders of magnitude above this band — so the tripwire still fires.
731 band = max(FP_CANCELLATION_CLAMP * np.finfo(float).eps, 1e-6) * max(m_in_total, 1.0)
732 # A real over-count is SUSTAINED — m_out stays below its running maximum across a θ-range.
733 # A single isolated sub-maximum sample is a measure-zero knife-edge at an exact collision θ
734 # (the parents are deactivated and the successor just born); it does not affect the
735 # integral-based bin-averaged cout, so require two consecutive samples below the running
736 # max before flagging. (The actual cout the API returns never sees the knife-edge.)
737 running_max = np.maximum.accumulate(m_out)
738 below = m_out < running_max - band
739 sustained = below[:-1] & below[1:]
740 if np.any(sustained):
741 idx = int(np.argmax(sustained))
742 drop = float(running_max[idx] - m_out[idx + 1])
743 return (
744 f"cumulative outlet mass drops by {drop:.4g} near θ={grid[idx + 1]:.4g} m³ "
745 "(the reader over-counts stored mass there — an interaction the solver failed to resolve)"
746 )
747 return None