Coverage for src/gwtransport/fronttracking/events.py: 0%
100 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 20:54 +0000
1"""Event detection for front tracking in (V, θ) coordinates.
3All intersections are pure line/line geometry in the (V, θ) plane because
4every wave speed dV/dθ is independent of flow. Functions return θ-coordinates
5of intersections; the solver translates to user-facing t at the API boundary.
7Events include:
9- Characteristic-characteristic collisions
10- Shock-shock collisions
11- Shock-characteristic collisions
12- Rarefaction boundary interactions
13- Outlet crossings
15All calculations return exact floating-point results with machine precision.
16The ``theta`` coordinate is cumulative flow [m³]; see
17:mod:`gwtransport.fronttracking.waves` for the (V, θ) convention.
19Available functions:
21- :class:`EventType` - Enumeration of the event kinds the solver dispatches on: characteristic-characteristic,
22 shock-shock, shock-characteristic, rarefaction-characteristic and shock-rarefaction collisions, decaying-shock
23 fan exhaustion, wave merges from the face calculus, and outlet crossings.
25- :class:`Event` - Record of one scheduled event: the θ at which it occurs, its :class:`EventType`, the waves
26 involved, the position ``location`` [m³], the ``'head'``/``'tail'`` ``boundary_type`` for rarefaction
27 collisions, and the colliding face pair for merges. It defines no ordering, because the solver orders
28 ``(theta, counter, ...)`` tuples rather than the objects themselves.
30- :func:`find_characteristic_intersection` - First strictly-future crossing ``(θ, V)`` of two characteristics,
31 or ``None`` when their speeds are equal or they never meet after ``theta_current``. Both lines are evaluated
32 from the shared reference ``max(θ_start_1, θ_start_2, θ_current)``.
34- :func:`find_shock_shock_intersection` - Same line/line crossing for two shocks, using each shock's constant
35 Rankine-Hugoniot speed.
37- :func:`find_shock_characteristic_intersection` - Same for a shock against a characteristic. The operand
38 order is load-bearing: ``V`` is evaluated on the first line, so it fixes the successor wave's ``v_start``
39 bit for bit.
41- :func:`find_rarefaction_boundary_intersections` - Crossings of a rarefaction's head and tail lines (which
42 travel at ``1/R(c_head)`` and ``1/R(c_tail)``) with another wave, returned as a list of
43 ``(θ, V, boundary_type)`` with ``boundary_type`` in ``{'head', 'tail'}``.
45- :func:`find_outlet_crossing` - Cumulative flow at which a wave reaches ``v_outlet``, assuming positive flow:
46 a straight-line inverse for characteristics and shocks, and the cached-trajectory inverse
47 ``outlet_crossing_theta`` for the fan-fed shocks (both the decaying and the doubly-fed one). Returns
48 ``None`` if the wave is inactive, moves backward, has a pinned characteristic speed, or has already passed
49 the outlet within a relative tolerance that stops a just-processed crossing from re-firing. Rarefactions
50 are excluded — their callers split head and tail into separate boundary crossings.
52- :func:`is_outlet_crossing_pinned` - Whether a boundary state sits on the ``c_min`` retardation floor with an
53 inflated ``R(c_min)``, so that its scheduled outlet crossing is a floor artifact rather than physics and the
54 caller should drop it.
55"""
57from dataclasses import dataclass
58from enum import Enum
60from gwtransport.fronttracking.math import characteristic_position, characteristic_speed
61from gwtransport.fronttracking.waves import (
62 CharacteristicWave,
63 DecayingShockWave,
64 DoubleFanShockWave,
65 RarefactionWave,
66 ShockWave,
67)
69# Numerical tolerance constants
70EPSILON_SPEED = 1e-15 # Tolerance for checking if two speeds are equal (machine precision)
71# A boundary state at/below the c_min retardation floor whose floored retardation
72# exceeds this value is "pinned": for the n>1 dry-soil singularity R(c_min) is
73# inflated to ~1e6, so the state advects orders of magnitude slower than any
74# physical wave and its outlet crossing lands at a non-physical θ (~1e8) that only
75# pollutes the diagnostic event record. n<1 clean water (R(c_min) ≈ 1, fast) stays
76# well below this threshold and is NOT pinned, so its outlet crossing is kept.
77OUTLET_PIN_RETARDATION = 1e4
80def is_outlet_crossing_pinned(concentration: float, sorption) -> bool:
81 """Whether a boundary state is pinned by the ``c_min`` retardation floor.
83 A crossing scheduled for such a state is a non-physical artifact (its speed is
84 a floor artifact, not physics); the caller drops it so it does not pollute the
85 solver's event record / ``theta_current``.
87 Parameters
88 ----------
89 concentration : float
90 Boundary-state concentration [mass/volume].
91 sorption : SorptionModel
92 Sorption model (supplies ``c_min`` and ``retardation``).
94 Returns
95 -------
96 bool
97 ``True`` only when ``concentration`` is at/below ``c_min`` AND the floored
98 retardation ``R(c_min)`` is inflated past :data:`OUTLET_PIN_RETARDATION`.
99 """
100 c_min = getattr(sorption, "c_min", 0.0)
101 if concentration > c_min:
102 return False
103 return float(sorption.retardation(c_min)) > OUTLET_PIN_RETARDATION
106class EventType(Enum):
107 """All possible event types in front tracking simulation."""
109 CHAR_CHAR_COLLISION = "characteristic_collision"
110 """Two characteristics intersect (will form shock)."""
111 SHOCK_SHOCK_COLLISION = "shock_collision"
112 """Two shocks collide (will merge)."""
113 SHOCK_CHAR_COLLISION = "shock_characteristic_collision"
114 """Shock catches or is caught by characteristic."""
115 RAREF_CHAR_COLLISION = "rarefaction_characteristic_collision"
116 """Rarefaction boundary intersects with characteristic."""
117 SHOCK_RAREF_COLLISION = "shock_rarefaction_collision"
118 """Shock intersects with rarefaction boundary."""
119 DSW_FAN_EXHAUSTED = "decaying_shock_fan_exhausted"
120 """A decaying shock's fan is exhausted (c_decay reached c_fan_tail)."""
121 WAVE_MERGE = "wave_merge"
122 """Two faces overtake (universal merge): any interaction involving a decaying/doubly-fed
123 shock — fan-entry, doubly-fed formation, same-apex annihilation, side exhaustion (a
124 doubly-fed shock crossing its own fan boundary line), and their compositions."""
125 OUTLET_CROSSING = "outlet_crossing"
126 """Wave crosses outlet boundary."""
129@dataclass
130class Event:
131 """A single event in the simulation, ordered by cumulative flow θ.
133 The solver's priority queue orders ``(theta, counter, ...)`` tuples, not
134 ``Event`` objects, so this dataclass intentionally defines no ordering.
136 Parameters
137 ----------
138 theta : float
139 Cumulative flow at which the event occurs [m³].
140 event_type : EventType
141 Type of event.
142 waves_involved : list
143 List of wave objects involved in this event.
144 location : float
145 Volumetric position at which the event occurs [m³].
146 boundary_type : str or None
147 Which rarefaction boundary collided: ``'head'`` or ``'tail'``.
148 Set for rarefaction collision events.
149 """
151 theta: float
152 event_type: EventType
153 waves_involved: list # List[Wave] - can't type hint due to circular import
154 location: float
155 boundary_type: str | None = None
156 faces: tuple | None = None # (Face, Face) for WAVE_MERGE
158 def __repr__(self): # noqa: D105
159 return (
160 f"Event(θ={self.theta:.3f}, type={self.event_type.value}, "
161 f"location={self.location:.3f}, n_waves={len(self.waves_involved)})"
162 )
165def _line_intersection(
166 theta_start_a: float,
167 v_start_a: float,
168 speed_a: float,
169 theta_start_b: float,
170 v_start_b: float,
171 speed_b: float,
172 theta_current: float,
173) -> tuple[float, float] | None:
174 """First future crossing of two straight (V, θ) fronts, or ``None``.
176 Every characteristic/shock/rarefaction-boundary travels at a flow-free constant
177 speed, so each is a line ``V = v_start + speed·(θ − θ_start)``. Both are evaluated
178 from the shared reference ``θ_both = max(θ_start_a, θ_start_b, θ_current)`` and the
179 crossing is returned only when strictly in the future (``dθ > 0``).
181 ``V_intersect`` is evaluated on line ``a``: the ``dθ`` is invariant under an a↔b swap
182 (exact IEEE negation of numerator and denominator) but ``V_intersect`` is not, so the
183 operand order each public wrapper passes is load-bearing for bit-reproducibility.
184 """
185 if abs(speed_a - speed_b) < EPSILON_SPEED:
186 return None
187 theta_both = max(theta_start_a, theta_start_b, theta_current)
188 v_a = v_start_a + speed_a * (theta_both - theta_start_a)
189 v_b = v_start_b + speed_b * (theta_both - theta_start_b)
190 dtheta = (v_b - v_a) / (speed_a - speed_b)
191 if dtheta <= 0:
192 return None
193 return (theta_both + dtheta, v_a + speed_a * dtheta)
196def find_characteristic_intersection(char1, char2, theta_current: float) -> tuple[float, float] | None:
197 """Find exact analytical intersection of two characteristics in (V, θ).
199 Returns (θ_intersect, V_intersect) if the intersection lies in the future
200 (θ > θ_current) and both characteristics are active there; otherwise None.
201 """
202 return _line_intersection(
203 char1.theta_start, char1.v_start, char1.speed(), char2.theta_start, char2.v_start, char2.speed(), theta_current
204 )
207def find_shock_shock_intersection(shock1, shock2, theta_current: float) -> tuple[float, float] | None:
208 """Find exact analytical intersection of two shocks in (V, θ)."""
209 return _line_intersection(
210 shock1.theta_start,
211 shock1.v_start,
212 shock1.speed,
213 shock2.theta_start,
214 shock2.v_start,
215 shock2.speed,
216 theta_current,
217 )
220def find_shock_characteristic_intersection(shock, char, theta_current: float) -> tuple[float, float] | None:
221 """Find exact analytical intersection of a shock and a characteristic in (V, θ)."""
222 return _line_intersection(
223 shock.theta_start, shock.v_start, shock.speed, char.theta_start, char.v_start, char.speed(), theta_current
224 )
227def find_rarefaction_boundary_intersections(raref, other_wave, theta_current: float) -> list[tuple[float, float, str]]:
228 """Intersections of a rarefaction's head/tail with another wave.
230 Both rarefaction boundaries propagate at characteristic speeds (head at
231 ``1/R(c_head)``, tail at ``1/R(c_tail)``), so each is a straight (V, θ) line
232 fed directly into :func:`_line_intersection` — no temporary ``CharacteristicWave``
233 objects. The operand order is per-branch: ``a`` is the raref boundary against a
234 characteristic, but the SHOCK against a raref boundary, so ``V_intersect`` — the new
235 wave's ``v_start`` — matches :func:`find_shock_characteristic_intersection` bit for bit.
237 Returns
238 -------
239 list of tuple
240 ``(θ_intersect, V_intersect, boundary_type)`` for each intersection,
241 where boundary_type is ``'head'`` or ``'tail'``.
242 """
243 intersections = []
244 raref_boundaries = ((raref.head_speed(), "head"), (raref.tail_speed(), "tail"))
246 if isinstance(other_wave, CharacteristicWave):
247 for s_raref, tag in raref_boundaries:
248 hit = _line_intersection(
249 raref.theta_start,
250 raref.v_start,
251 s_raref,
252 other_wave.theta_start,
253 other_wave.v_start,
254 other_wave.speed(),
255 theta_current,
256 )
257 if hit:
258 intersections.append((hit[0], hit[1], tag))
260 elif isinstance(other_wave, ShockWave):
261 # a=SHOCK, b=raref boundary (matches find_shock_characteristic_intersection's order).
262 for s_raref, tag in raref_boundaries:
263 hit = _line_intersection(
264 other_wave.theta_start,
265 other_wave.v_start,
266 other_wave.speed,
267 raref.theta_start,
268 raref.v_start,
269 s_raref,
270 theta_current,
271 )
272 if hit:
273 intersections.append((hit[0], hit[1], tag))
275 elif isinstance(other_wave, RarefactionWave):
276 other_boundaries = (other_wave.head_speed(), other_wave.tail_speed())
277 for s_raref, tag in raref_boundaries:
278 for s_other in other_boundaries:
279 hit = _line_intersection(
280 raref.theta_start,
281 raref.v_start,
282 s_raref,
283 other_wave.theta_start,
284 other_wave.v_start,
285 s_other,
286 theta_current,
287 )
288 if hit:
289 intersections.append((hit[0], hit[1], tag))
291 return intersections
294def find_outlet_crossing(wave, v_outlet: float, theta_current: float) -> float | None:
295 """Find the cumulative flow θ at which the wave crosses ``v_outlet``.
297 Handles ``CharacteristicWave``, ``ShockWave``, and ``DecayingShockWave``.
298 Rarefaction outlet crossings are handled by the callers directly (the
299 solver and ``output.py`` split them into head/tail boundary crossings), so
300 a ``RarefactionWave`` never reaches this function and returns ``None``.
302 Assumes positive flow (waves always move toward larger V). Returns None if
303 the wave has already passed the outlet, is not active, or moves backward.
304 The "already past" check uses a relative tolerance so that a wave whose
305 crossing event has just been processed (and is at v_outlet ± a few ULPs)
306 does not re-emit a duplicate crossing one ULP later.
307 """
308 if not wave.is_active:
309 return None
311 # Suppress re-emission when v_current is within FP of v_outlet: the
312 # crossing was already recorded on the prior iteration.
313 tol = 1e-12 * max(abs(v_outlet), abs(wave.v_start), 1.0)
315 if isinstance(wave, CharacteristicWave):
316 theta_eval = max(theta_current, wave.theta_start)
317 v_current = characteristic_position(
318 wave.concentration, wave.sorption, wave.theta_start, wave.v_start, theta_eval
319 )
321 if v_current is None or v_current >= v_outlet - tol:
322 return None
324 speed = characteristic_speed(wave.concentration, wave.sorption)
326 # A c_min-floored (pinned) characteristic — R(c_min) inflated for n>1,
327 # c→0 — advects too slowly to cross at any physical θ; suppress the
328 # artifact crossing rather than scheduling it at θ~1e8.
329 if speed <= 0 or is_outlet_crossing_pinned(wave.concentration, wave.sorption):
330 return None
332 dtheta = (v_outlet - v_current) / speed
333 return theta_eval + dtheta
335 if isinstance(wave, ShockWave):
336 theta_eval = max(theta_current, wave.theta_start)
337 v_current = wave.v_start + wave.speed * (theta_eval - wave.theta_start)
339 if v_current >= v_outlet - tol:
340 return None
342 if wave.speed <= 0:
343 return None
345 dtheta = (v_outlet - v_current) / wave.speed
346 return theta_eval + dtheta
348 if isinstance(wave, (DecayingShockWave, DoubleFanShockWave)):
349 # Closed-form / cached-trajectory inverse V_s(theta) = v_outlet.
350 theta_cross = wave.outlet_crossing_theta(v_outlet)
351 if theta_cross is None:
352 return None
353 # Suppress re-emission within FP of the prior crossing (same convention
354 # as the linear-shock branch above).
355 if theta_cross <= theta_current + 1e-15 * max(abs(theta_current), 1.0):
356 return None
357 return theta_cross
359 return None