Coverage for src/gwtransport/fronttracking/handlers.py: 78%
120 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +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.
16"""
18from gwtransport.fronttracking.math import (
19 FreundlichSorption,
20 NonlinearSorption,
21 SorptionModel,
22 characteristic_speed,
23)
24from gwtransport.fronttracking.waves import CharacteristicWave, DecayingShockWave, RarefactionWave, ShockWave
26# Numerical tolerance constants
27# Shared single source for the negligible-concentration-change tolerance; the
28# solver imports this rather than redefining it.
29EPSILON_CONCENTRATION = 1e-15
32def handle_characteristic_collision(
33 char1: CharacteristicWave,
34 char2: CharacteristicWave,
35 theta_event: float,
36 v_event: float,
37) -> list[ShockWave]:
38 """Two characteristics collide → emit a shock.
40 The faster characteristic catches the slower one from behind. By the
41 entropy condition this compressive interaction is always a shock,
42 independently of the sorption regime (Freundlich n>1, n<1, or constant
43 retardation).
45 Parameters
46 ----------
47 char1, char2 : CharacteristicWave
48 Colliding characteristics.
49 theta_event : float
50 Cumulative flow at which the collision occurs [m³].
51 v_event : float
52 Position at which the collision occurs [m³].
54 Returns
55 -------
56 list of ShockWave
57 Single shock created at the collision point.
59 Raises
60 ------
61 RuntimeError
62 If the resulting shock fails the Lax entropy condition.
63 """
64 s1 = characteristic_speed(char1.concentration, char1.sorption)
65 s2 = characteristic_speed(char2.concentration, char2.sorption)
67 if s1 > s2:
68 c_left = char1.concentration
69 c_right = char2.concentration
70 else:
71 c_left = char2.concentration
72 c_right = char1.concentration
74 shock = ShockWave(
75 theta_start=theta_event,
76 v_start=v_event,
77 c_left=c_left,
78 c_right=c_right,
79 sorption=char1.sorption,
80 )
82 if not shock.satisfies_entropy():
83 msg = (
84 f"Characteristic collision created non-entropic shock at θ={theta_event:.3f}, "
85 f"V={v_event:.3f}. c_left={c_left:.3f}, c_right={c_right:.3f}, "
86 f"shock_speed={shock.speed:.6g}"
87 )
88 raise RuntimeError(msg)
90 char1.deactivate(theta_event)
91 char2.deactivate(theta_event)
92 return [shock]
95def handle_shock_collision(
96 shock1: ShockWave,
97 shock2: ShockWave,
98 theta_event: float,
99 v_event: float,
100) -> list[ShockWave]:
101 """Two shocks collide → merge into a single shock connecting outer states.
103 The merged shock has ``c_left`` from the faster (upstream) shock,
104 ``c_right`` from the slower (downstream) shock; its speed is recomputed
105 via Rankine-Hugoniot.
107 Parameters
108 ----------
109 shock1, shock2 : ShockWave
110 Colliding shocks.
111 theta_event, v_event : float
112 Cumulative flow [m³] and position [m³] of the collision.
114 Returns
115 -------
116 list of ShockWave
117 Single merged shock.
119 Raises
120 ------
121 RuntimeError
122 If the merged shock violates the entropy condition.
123 """
124 if shock1.speed > shock2.speed:
125 c_left = shock1.c_left
126 c_right = shock2.c_right
127 else:
128 c_left = shock2.c_left
129 c_right = shock1.c_right
131 merged = ShockWave(
132 theta_start=theta_event,
133 v_start=v_event,
134 c_left=c_left,
135 c_right=c_right,
136 sorption=shock1.sorption,
137 )
139 if not merged.satisfies_entropy():
140 msg = (
141 f"Shock merger created non-entropic shock at θ={theta_event:.3f}. "
142 f"This may indicate complex wave interaction requiring special handling."
143 )
144 raise RuntimeError(msg)
146 shock1.deactivate(theta_event)
147 shock2.deactivate(theta_event)
149 return [merged]
152def handle_shock_characteristic_collision(
153 shock: ShockWave,
154 char: CharacteristicWave,
155 theta_event: float,
156 v_event: float,
157) -> list:
158 """Shock catches or is caught by a characteristic.
160 The characteristic concentration modifies one side of the shock:
162 - Shock catches char (shock faster): modifies ``c_right``.
163 - Char catches shock (char faster): modifies ``c_left``.
165 If the resulting shock satisfies entropy it is emitted (compression);
166 otherwise a rarefaction is created (expansion) to preserve mass balance.
167 """
168 s_shock = shock.speed
169 s_char = characteristic_speed(char.concentration, char.sorption)
171 if s_shock > s_char:
172 new_shock = ShockWave(
173 theta_start=theta_event,
174 v_start=v_event,
175 c_left=shock.c_left,
176 c_right=char.concentration,
177 sorption=shock.sorption,
178 )
179 else:
180 new_shock = ShockWave(
181 theta_start=theta_event,
182 v_start=v_event,
183 c_left=char.concentration,
184 c_right=shock.c_right,
185 sorption=shock.sorption,
186 )
188 if not new_shock.satisfies_entropy():
189 # Expansion regime: emit a rarefaction whose head is the faster
190 # state and tail the slower state.
191 if s_shock > s_char:
192 c_head = shock.c_left
193 c_tail = char.concentration
194 else:
195 c_head = char.concentration
196 c_tail = shock.c_right
198 s_head = characteristic_speed(c_head, shock.sorption)
199 s_tail = characteristic_speed(c_tail, shock.sorption)
201 if s_head > s_tail:
202 raref = RarefactionWave(
203 theta_start=theta_event,
204 v_start=v_event,
205 c_head=c_head,
206 c_tail=c_tail,
207 sorption=shock.sorption,
208 )
209 shock.deactivate(theta_event)
210 char.deactivate(theta_event)
211 return [raref]
212 # Edge case (s_head == s_tail within machine precision): deactivate
213 # and emit nothing.
214 shock.deactivate(theta_event)
215 char.deactivate(theta_event)
216 return []
218 shock.deactivate(theta_event)
219 char.deactivate(theta_event)
220 return [new_shock]
223def handle_shock_rarefaction_collision(
224 shock: ShockWave,
225 raref: RarefactionWave,
226 theta_event: float,
227 v_event: float,
228 boundary_type: str | None,
229) -> list:
230 """Shock interacts with a rarefaction fan (tail or head boundary).
232 Every shock↔rarefaction collision is resolved exactly by a single
233 :class:`~gwtransport.fronttracking.waves.DecayingShockWave` whose trajectory
234 subsumes the fan + shock together, for any
235 :class:`~gwtransport.fronttracking.math.NonlinearSorption`:
237 - **Head collision** (rarefaction head catches the leading shock): the
238 decaying side is the left, ``c_decay_initial = raref.c_head``,
239 ``c_fixed = shock.c_right``, and ``c_fan_tail = raref.c_tail`` (the fan's
240 other boundary, which bounds the decay so partial drying is handled).
241 - **Tail collision** (trailing shock catches the rarefaction tail): the
242 decaying side is the right, ``c_decay_initial = raref.c_tail``,
243 ``c_fixed = shock.c_left``, and ``c_fan_tail = raref.c_head``.
245 The fan is bounded by ``c_fan_tail``: the solver's ``DSW_FAN_EXHAUSTED``
246 event spawns a regular shock once the decaying side reaches it, so partial
247 drying (``raref.c_tail != shock.c_right``) is resolved exactly. If the
248 rarefaction boundary is not faster than the shock (degenerate solver/test
249 input), both waves are deactivated and nothing is emitted.
251 Returns
252 -------
253 list of Wave
254 ``[DecayingShockWave]`` for a physical collision, or ``[]`` for
255 degenerate input.
256 """
257 sorption = raref.sorption
258 # Rarefactions only form for nonlinear isotherms, so the DecayingShockWave's
259 # NonlinearSorption requirement is always met here.
260 assert isinstance(sorption, NonlinearSorption) # noqa: S101
262 if boundary_type == "head":
263 # Rarefaction head catches the shock; decaying side is the left.
264 s_raref_boundary = characteristic_speed(raref.c_head, sorption)
265 if s_raref_boundary <= shock.speed:
266 shock.deactivate(theta_event)
267 raref.deactivate(theta_event)
268 return []
269 c_decay_initial = raref.c_head
270 c_fixed = shock.c_right
271 c_fan_tail = raref.c_tail
272 decay_side = "left"
273 elif boundary_type == "tail":
274 # Trailing shock catches the rarefaction tail; decaying side is the right.
275 s_raref_boundary = characteristic_speed(raref.c_tail, sorption)
276 if shock.speed <= s_raref_boundary:
277 shock.deactivate(theta_event)
278 raref.deactivate(theta_event)
279 return []
280 c_decay_initial = raref.c_tail
281 c_fixed = shock.c_left
282 c_fan_tail = raref.c_head
283 decay_side = "right"
284 else:
285 msg = f"handle_shock_rarefaction_collision: unknown boundary_type {boundary_type!r}"
286 raise RuntimeError(msg)
288 decaying = DecayingShockWave(
289 theta_start=theta_event,
290 v_start=v_event,
291 c_decay_initial=c_decay_initial,
292 c_fixed=c_fixed,
293 c_fan_tail=c_fan_tail,
294 decay_side=decay_side,
295 v_origin=raref.v_start,
296 theta_origin=raref.theta_start,
297 sorption=sorption,
298 )
300 shock.deactivate(theta_event)
301 raref.deactivate(theta_event)
302 return [decaying]
305def handle_rarefaction_characteristic_collision(
306 raref: RarefactionWave,
307 char: CharacteristicWave,
308 theta_event: float,
309 v_event: float,
310 boundary_type: str | None,
311) -> list:
312 """Rarefaction boundary intersects a characteristic.
314 The safe option (b) from the front-tracking rebuild plan: when a
315 characteristic's concentration matches the boundary concentration to
316 within tolerance the characteristic is absorbed; otherwise an
317 informative ``RuntimeError`` is raised because deactivating it would
318 silently destroy mass.
320 Raises
321 ------
322 RuntimeError
323 If the characteristic's concentration does not match the colliding
324 rarefaction boundary concentration within tolerance, or if
325 ``boundary_type`` is not ``'head'`` or ``'tail'``.
326 """
327 rel_tol = 1e-9
328 abs_tol = 1e-12
329 raref_range = abs(raref.c_head - raref.c_tail)
330 tol = max(rel_tol * raref_range, abs_tol)
332 if boundary_type == "head":
333 boundary_c = raref.c_head
334 elif boundary_type == "tail":
335 boundary_c = raref.c_tail
336 else:
337 msg = f"handle_rarefaction_characteristic_collision: unknown boundary_type {boundary_type!r}"
338 raise RuntimeError(msg)
340 if abs(char.concentration - boundary_c) > tol:
341 msg = (
342 f"Rarefaction-characteristic collision at θ={theta_event:.6f}, V={v_event:.6f} would silently "
343 f"destroy mass: characteristic concentration {char.concentration:.6g} differs from "
344 f"rarefaction {boundary_type} concentration {boundary_c:.6g} by "
345 f"{abs(char.concentration - boundary_c):.3g} (tolerance {tol:.3g}). "
346 f"Proper wave splitting at the rarefaction boundary is required for this case."
347 )
348 raise RuntimeError(msg)
350 char.deactivate(theta_event)
351 return []
354def handle_outlet_crossing(wave, theta_event: float, v_outlet: float) -> dict:
355 """Record a wave crossing the outlet boundary.
357 The wave is NOT deactivated — it remains for concentration queries at
358 points between its origin and the outlet. The returned event record
359 holds the cumulative flow ``theta`` at which the crossing occurs; the
360 solver translates this to the user-facing time when appending to
361 ``state.events``.
362 """
363 return {
364 "theta": theta_event,
365 "type": "outlet_crossing",
366 "wave": wave,
367 "location": v_outlet,
368 "concentration_left": wave.concentration_left(),
369 "concentration_right": wave.concentration_right(),
370 }
373def create_inlet_waves_at_theta(
374 c_prev: float,
375 c_new: float,
376 theta: float,
377 sorption: SorptionModel,
378) -> list:
379 """Emit the wave produced by a step change in inlet concentration.
381 All inlet waves originate at the inlet face ``V = 0``. Wave type is
382 determined by characteristic speed comparison in (V, θ):
384 - ``s_new > s_prev``: compression → shock.
385 - ``s_new < s_prev``: expansion → rarefaction.
386 - equal: contact discontinuity → characteristic.
388 For shocks the entropy condition is verified; if violated, an empty list
389 is returned (mass balance may be affected — a known limitation handled
390 by ``DecayingShockWave``).
391 """
392 if abs(c_new - c_prev) < EPSILON_CONCENTRATION:
393 return []
395 c_min = getattr(sorption, "c_min", 0.0)
396 is_n_lt_1 = isinstance(sorption, FreundlichSorption) and sorption.n < 1.0
398 # n<1, c_prev=0 or c_new=0: emit a single CharacteristicWave; clean water
399 # has a well-defined speed since R(0)=1.
400 if (c_prev <= c_min or c_new <= c_min) and is_n_lt_1 and c_min == 0:
401 return [
402 CharacteristicWave(
403 theta_start=theta,
404 v_start=0.0,
405 concentration=c_new,
406 sorption=sorption,
407 )
408 ]
410 s_prev = characteristic_speed(c_prev, sorption)
411 s_new = characteristic_speed(c_new, sorption)
413 if s_new > s_prev + 1e-15:
414 shock = ShockWave(
415 theta_start=theta,
416 v_start=0.0,
417 c_left=c_new,
418 c_right=c_prev,
419 sorption=sorption,
420 )
421 if not shock.satisfies_entropy():
422 return []
423 return [shock]
425 if s_new < s_prev - 1e-15:
426 return [
427 RarefactionWave(
428 theta_start=theta,
429 v_start=0.0,
430 c_head=c_prev,
431 c_tail=c_new,
432 sorption=sorption,
433 )
434 ]
436 return [
437 CharacteristicWave(
438 theta_start=theta,
439 v_start=0.0,
440 concentration=c_new,
441 sorption=sorption,
442 )
443 ]