Coverage for src/gwtransport/fronttracking/handlers.py: 80%
109 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""Event handlers for front tracking in (V, θ) coordinates.
3Each handler receives the waves involved in an event and returns the new
4waves created by the interaction. In (V, θ) coordinates every wave speed is
5flow-free, so handlers depend only on concentrations and the sorption
6isotherm — flow does not appear.
8All handlers enforce physical correctness:
10- Mass conservation (Rankine-Hugoniot condition)
11- Entropy conditions (Lax condition for shocks)
12- Causality (no backward-traveling information)
14Handlers modify wave states in-place by deactivating parent waves and
15creating new child waves. Positions are volumetric [m³] and ``theta_event``
16is cumulative flow [m³]; see :mod:`gwtransport.fronttracking.waves` for the
17(V, θ) convention.
19Available functions:
21- :func:`handle_characteristic_collision` - Two characteristics meet, the faster catching the slower from
22 behind. This compression always emits a single shock spanning the two concentrations, in every sorption
23 regime; a resulting shock that fails the Lax condition raises ``RuntimeError``.
25- :func:`handle_shock_collision` - Two shocks merge into one connecting the outer states: ``c_left`` from the
26 upstream shock, ``c_right`` from the downstream one, with the speed recomputed from Rankine-Hugoniot. An
27 entropy-violating merge raises ``RuntimeError``.
29- :func:`handle_shock_characteristic_collision` - A shock and a characteristic meet; the characteristic's
30 concentration replaces ``c_right`` when the shock does the catching and ``c_left`` when the characteristic
31 does. The successor is the modified shock if it satisfies entropy, otherwise a rarefaction, so mass balance
32 is preserved either way.
34- :func:`handle_shock_rarefaction_collision` - A shock meets a rarefaction head or tail. The pair is replaced
35 by one :class:`~gwtransport.fronttracking.waves.DecayingShockWave` that subsumes fan and shock together: a
36 head collision decays the left side from ``raref.c_head`` with ``c_fixed = shock.c_right``, a tail collision
37 decays the right side from ``raref.c_tail`` with ``c_fixed = shock.c_left``, and the fan's opposite boundary
38 becomes ``c_fan_tail`` so partial drying is resolved exactly. Degenerate input where the boundary is not
39 faster than the shock deactivates both waves and emits nothing.
41- :func:`handle_rarefaction_characteristic_collision` - A rarefaction boundary meets a characteristic. The
42 characteristic is absorbed when its concentration matches the boundary value within tolerance; otherwise
43 ``RuntimeError`` is raised rather than silently destroying mass.
45- :func:`create_inlet_waves_at_theta` - Emit the wave a step from ``c_prev`` to ``c_new`` produces at the
46 inlet face ``V = 0``: a shock when the new characteristic speed is higher (compression), a rarefaction when
47 it is lower (expansion), and a characteristic when the two speeds tie. Returns an empty list for a
48 negligible step or for a shock that fails the entropy check.
49"""
51from gwtransport.fronttracking.math import (
52 NonlinearSorption,
53 SorptionModel,
54 characteristic_speed,
55)
56from gwtransport.fronttracking.waves import CharacteristicWave, DecayingShockWave, RarefactionWave, ShockWave
58# Numerical tolerance constants
59# Shared single source for the negligible-concentration-change tolerance; the
60# solver imports this rather than redefining it.
61EPSILON_CONCENTRATION = 1e-15
64def handle_characteristic_collision(
65 char1: CharacteristicWave,
66 char2: CharacteristicWave,
67 theta_event: float,
68 v_event: float,
69) -> list[ShockWave]:
70 """Two characteristics collide → emit a shock.
72 The faster characteristic catches the slower one from behind. By the
73 entropy condition this compressive interaction is always a shock,
74 independently of the sorption regime (Freundlich n>1, n<1, or constant
75 retardation).
77 Parameters
78 ----------
79 char1, char2 : CharacteristicWave
80 Colliding characteristics.
81 theta_event : float
82 Cumulative flow at which the collision occurs [m³].
83 v_event : float
84 Position at which the collision occurs [m³].
86 Returns
87 -------
88 list of ShockWave
89 Single shock created at the collision point.
91 Raises
92 ------
93 RuntimeError
94 If the resulting shock fails the Lax entropy condition.
95 """
96 s1 = characteristic_speed(char1.concentration, char1.sorption)
97 s2 = characteristic_speed(char2.concentration, char2.sorption)
99 if s1 > s2:
100 c_left = char1.concentration
101 c_right = char2.concentration
102 else:
103 c_left = char2.concentration
104 c_right = char1.concentration
106 shock = ShockWave(
107 theta_start=theta_event,
108 v_start=v_event,
109 c_left=c_left,
110 c_right=c_right,
111 sorption=char1.sorption,
112 )
114 if not shock.satisfies_entropy():
115 msg = (
116 f"Characteristic collision created non-entropic shock at θ={theta_event:.3f}, "
117 f"V={v_event:.3f}. c_left={c_left:.3f}, c_right={c_right:.3f}, "
118 f"shock_speed={shock.speed:.6g}"
119 )
120 raise RuntimeError(msg)
122 char1.deactivate(theta_event)
123 char2.deactivate(theta_event)
124 return [shock]
127def handle_shock_collision(
128 shock1: ShockWave,
129 shock2: ShockWave,
130 theta_event: float,
131 v_event: float,
132) -> list[ShockWave]:
133 """Two shocks collide → merge into a single shock connecting outer states.
135 The merged shock has ``c_left`` from the faster (upstream) shock,
136 ``c_right`` from the slower (downstream) shock; its speed is recomputed
137 via Rankine-Hugoniot.
139 Parameters
140 ----------
141 shock1, shock2 : ShockWave
142 Colliding shocks.
143 theta_event, v_event : float
144 Cumulative flow [m³] and position [m³] of the collision.
146 Returns
147 -------
148 list of ShockWave
149 Single merged shock.
151 Raises
152 ------
153 RuntimeError
154 If the merged shock violates the entropy condition.
155 """
156 if shock1.speed > shock2.speed:
157 c_left = shock1.c_left
158 c_right = shock2.c_right
159 else:
160 c_left = shock2.c_left
161 c_right = shock1.c_right
163 merged = ShockWave(
164 theta_start=theta_event,
165 v_start=v_event,
166 c_left=c_left,
167 c_right=c_right,
168 sorption=shock1.sorption,
169 )
171 if not merged.satisfies_entropy():
172 msg = (
173 f"Shock merger created non-entropic shock at θ={theta_event:.3f}. "
174 f"This may indicate complex wave interaction requiring special handling."
175 )
176 raise RuntimeError(msg)
178 shock1.deactivate(theta_event)
179 shock2.deactivate(theta_event)
181 return [merged]
184def handle_shock_characteristic_collision(
185 shock: ShockWave,
186 char: CharacteristicWave,
187 theta_event: float,
188 v_event: float,
189) -> list:
190 """Shock catches or is caught by a characteristic.
192 The characteristic concentration modifies one side of the shock:
194 - Shock catches char (shock faster): modifies ``c_right``.
195 - Char catches shock (char faster): modifies ``c_left``.
197 If the resulting shock satisfies entropy it is emitted (compression);
198 otherwise a rarefaction is created (expansion) to preserve mass balance.
199 """
200 s_shock = shock.speed
201 s_char = characteristic_speed(char.concentration, char.sorption)
203 if s_shock > s_char:
204 new_shock = ShockWave(
205 theta_start=theta_event,
206 v_start=v_event,
207 c_left=shock.c_left,
208 c_right=char.concentration,
209 sorption=shock.sorption,
210 )
211 else:
212 new_shock = ShockWave(
213 theta_start=theta_event,
214 v_start=v_event,
215 c_left=char.concentration,
216 c_right=shock.c_right,
217 sorption=shock.sorption,
218 )
220 # Both parents are consumed by the collision on every outcome.
221 shock.deactivate(theta_event)
222 char.deactivate(theta_event)
224 if new_shock.satisfies_entropy():
225 return [new_shock]
227 # Expansion regime: emit a rarefaction whose head is the faster state and
228 # tail the slower state.
229 if s_shock > s_char:
230 c_head = shock.c_left
231 c_tail = char.concentration
232 else:
233 c_head = char.concentration
234 c_tail = shock.c_right
236 s_head = characteristic_speed(c_head, shock.sorption)
237 s_tail = characteristic_speed(c_tail, shock.sorption)
239 if s_head <= s_tail:
240 # Edge case (s_head == s_tail within machine precision): emit nothing.
241 return []
243 return [
244 RarefactionWave(
245 theta_start=theta_event,
246 v_start=v_event,
247 c_head=c_head,
248 c_tail=c_tail,
249 sorption=shock.sorption,
250 )
251 ]
254def handle_shock_rarefaction_collision(
255 shock: ShockWave,
256 raref: RarefactionWave,
257 theta_event: float,
258 v_event: float,
259 boundary_type: str | None,
260) -> list:
261 """Shock interacts with a rarefaction fan (tail or head boundary).
263 Every shock↔rarefaction collision is resolved exactly by a single
264 :class:`~gwtransport.fronttracking.waves.DecayingShockWave` whose trajectory
265 subsumes the fan + shock together, for any
266 :class:`~gwtransport.fronttracking.math.NonlinearSorption`:
268 - **Head collision** (rarefaction head catches the leading shock): the
269 decaying side is the left, ``c_decay_initial = raref.c_head``,
270 ``c_fixed = shock.c_right``, and ``c_fan_tail = raref.c_tail`` (the fan's
271 other boundary, which bounds the decay so partial drying is handled).
272 - **Tail collision** (trailing shock catches the rarefaction tail): the
273 decaying side is the right, ``c_decay_initial = raref.c_tail``,
274 ``c_fixed = shock.c_left``, and ``c_fan_tail = raref.c_head``.
276 The fan is bounded by ``c_fan_tail``: the solver's ``DSW_FAN_EXHAUSTED``
277 event spawns a regular shock once the decaying side reaches it, so partial
278 drying (``raref.c_tail != shock.c_right``) is resolved exactly. If the
279 rarefaction boundary is not faster than the shock (degenerate solver/test
280 input), both waves are deactivated and nothing is emitted.
282 Returns
283 -------
284 list of Wave
285 ``[DecayingShockWave]`` for a physical collision, or ``[]`` for
286 degenerate input.
287 """
288 sorption = raref.sorption
289 # Rarefactions only form for nonlinear isotherms, so the DecayingShockWave's
290 # NonlinearSorption requirement is always met here.
291 assert isinstance(sorption, NonlinearSorption) # noqa: S101
293 if boundary_type == "head":
294 # Rarefaction head catches the shock; decaying side is the left.
295 s_raref_boundary = characteristic_speed(raref.c_head, sorption)
296 if s_raref_boundary <= shock.speed:
297 shock.deactivate(theta_event)
298 raref.deactivate(theta_event)
299 return []
300 c_decay_initial = raref.c_head
301 c_fixed = shock.c_right
302 c_fan_tail = raref.c_tail
303 decay_side = "left"
304 elif boundary_type == "tail":
305 # Trailing shock catches the rarefaction tail; decaying side is the right.
306 s_raref_boundary = characteristic_speed(raref.c_tail, sorption)
307 if shock.speed <= s_raref_boundary:
308 shock.deactivate(theta_event)
309 raref.deactivate(theta_event)
310 return []
311 c_decay_initial = raref.c_tail
312 c_fixed = shock.c_left
313 c_fan_tail = raref.c_head
314 decay_side = "right"
315 else:
316 msg = f"handle_shock_rarefaction_collision: unknown boundary_type {boundary_type!r}"
317 raise RuntimeError(msg)
319 decaying = DecayingShockWave(
320 theta_start=theta_event,
321 v_start=v_event,
322 c_decay_initial=c_decay_initial,
323 c_fixed=c_fixed,
324 c_fan_tail=c_fan_tail,
325 decay_side=decay_side,
326 v_origin=raref.v_start,
327 theta_origin=raref.theta_start,
328 sorption=sorption,
329 )
331 shock.deactivate(theta_event)
332 raref.deactivate(theta_event)
333 return [decaying]
336def handle_rarefaction_characteristic_collision(
337 raref: RarefactionWave,
338 char: CharacteristicWave,
339 theta_event: float,
340 v_event: float,
341 boundary_type: str | None,
342) -> list:
343 """Rarefaction boundary intersects a characteristic.
345 When the characteristic's concentration matches the boundary concentration
346 to within tolerance the characteristic is absorbed; otherwise an
347 informative ``RuntimeError`` is raised, because deactivating it would
348 silently destroy mass.
350 Raises
351 ------
352 RuntimeError
353 If the characteristic's concentration does not match the colliding
354 rarefaction boundary concentration within tolerance, or if
355 ``boundary_type`` is not ``'head'`` or ``'tail'``.
356 """
357 rel_tol = 1e-9
358 abs_tol = 1e-12
359 raref_range = abs(raref.c_head - raref.c_tail)
360 tol = max(rel_tol * raref_range, abs_tol)
362 if boundary_type == "head":
363 boundary_c = raref.c_head
364 elif boundary_type == "tail":
365 boundary_c = raref.c_tail
366 else:
367 msg = f"handle_rarefaction_characteristic_collision: unknown boundary_type {boundary_type!r}"
368 raise RuntimeError(msg)
370 if abs(char.concentration - boundary_c) > tol:
371 msg = (
372 f"Rarefaction-characteristic collision at θ={theta_event:.6f}, V={v_event:.6f} would silently "
373 f"destroy mass: characteristic concentration {char.concentration:.6g} differs from "
374 f"rarefaction {boundary_type} concentration {boundary_c:.6g} by "
375 f"{abs(char.concentration - boundary_c):.3g} (tolerance {tol:.3g}). "
376 f"Proper wave splitting at the rarefaction boundary is required for this case."
377 )
378 raise RuntimeError(msg)
380 char.deactivate(theta_event)
381 return []
384def create_inlet_waves_at_theta(
385 c_prev: float,
386 c_new: float,
387 theta: float,
388 sorption: SorptionModel,
389) -> list:
390 """Emit the wave produced by a step change in inlet concentration.
392 All inlet waves originate at the inlet face ``V = 0``. Wave type is
393 determined by characteristic speed comparison in (V, θ):
395 - ``s_new > s_prev``: compression → shock.
396 - ``s_new < s_prev``: expansion → rarefaction.
397 - equal: contact discontinuity → characteristic.
399 For shocks the entropy condition is verified; if violated, an empty list
400 is returned (mass balance may be affected — a known limitation handled
401 by ``DecayingShockWave``).
402 """
403 if abs(c_new - c_prev) < EPSILON_CONCENTRATION:
404 return []
406 s_prev = characteristic_speed(c_prev, sorption)
407 s_new = characteristic_speed(c_new, sorption)
409 if s_new > s_prev + 1e-15:
410 shock = ShockWave(
411 theta_start=theta,
412 v_start=0.0,
413 c_left=c_new,
414 c_right=c_prev,
415 sorption=sorption,
416 )
417 if not shock.satisfies_entropy():
418 return []
419 return [shock]
421 if s_new < s_prev - 1e-15:
422 return [
423 RarefactionWave(
424 theta_start=theta,
425 v_start=0.0,
426 c_head=c_prev,
427 c_tail=c_new,
428 sorption=sorption,
429 )
430 ]
432 return [
433 CharacteristicWave(
434 theta_start=theta,
435 v_start=0.0,
436 concentration=c_new,
437 sorption=sorption,
438 c_ahead=c_prev,
439 )
440 ]