Coverage for src/gwtransport/fronttracking/interactions.py: 94%
201 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""Wave–wave interaction resolution for multi-front nonlinear-sorption transport.
3The event-driven solver in :mod:`gwtransport.fronttracking.solver` resolves collisions among
4characteristics, shocks and rarefactions with the closed-form helpers in
5:mod:`gwtransport.fronttracking.events`. This module adds the missing interaction classes —
6anything a :class:`~gwtransport.fronttracking.waves.DecayingShockWave` or
7:class:`~gwtransport.fronttracking.waves.DoubleFanShockWave` participates in — via one uniform
8calculus:
10- **Faces and feeders.** Every wave is a set of *faces* (a shock face, a contact, or a
11 rarefaction/fan boundary line). Each face separates a *left* (upstream) from a *right*
12 (downstream) :class:`~gwtransport.fronttracking.waves.Feeder` — a constant state or a
13 bounded self-similar fan.
14- **Universal merge.** When a rear (upstream, faster) face overtakes a front (downstream)
15 face, they merge into a single successor built from ``(rear.left_feeder, front.right_feeder)``.
16 This one rule generates shock↔shock merges, shocks entering a fan (fan-entry), a rarefaction
17 head catching a decaying shock (doubly-fed formation), same-apex decaying-shock annihilation,
18 and every composition thereof.
19- **Lipschitz-safe detection.** First crossings are found by a per-pair speed-bounded march
20 (``|dg/dθ| ≤ Λ_pair``) so no crossing is skipped, with a grazing-minimum bracket for
21 double-crossings inside one step. The bound is per-pair (not a global ``1``) because the
22 percolation conductivity isotherms have ``R < 1`` regions.
24Available functions:
26- :class:`Face` - One separating surface of a wave: the wave it belongs to, a ``role`` tag
27 (``'shock'``, ``'contact'``, ``'head'``, ``'tail'`` or ``'boundary'``), its left (upstream) and right
28 (downstream) :class:`~gwtransport.fronttracking.waves.Feeder`, whether its trajectory is curved, and an
29 upper bound on its speed for the Lipschitz march. ``position(theta)`` evaluates the face, returning ``None``
30 outside the owning wave's active θ-window.
32- :func:`iter_faces` - Enumerate the faces a wave exposes at ``theta``: one contact for a characteristic, one
33 shock face for a shock, head and tail for a rarefaction, and for the fan-fed shocks a curved shock face plus
34 each fan boundary line that is still free. The ``theta`` argument selects the historical state, so a
35 boundary already consumed by a wave entering the fan is not re-exposed.
37- :func:`find_face_crossing` - First θ in ``(θ_start, θ_horizon]`` at which two faces coincide, or ``None``.
38 The gap is marched with a per-pair speed bound so no sign change can be stepped over, and each step also
39 brackets the gap's interior minimum to catch a grazing double crossing.
41- :func:`make_wave_from_feeders` - Build the successor a merge produces from a ``(left, right)`` feeder pair
42 at ``(v, θ)``: two constants give a :class:`~gwtransport.fronttracking.waves.ShockWave` (or ``None`` for a
43 zero jump), one constant and one fan give a :class:`~gwtransport.fronttracking.waves.DecayingShockWave`, and
44 two fans give a :class:`~gwtransport.fronttracking.waves.DoubleFanShockWave`. A feeder whose far boundary is
45 already owned marks the successor's boundary consumed on that side.
47- :func:`resolve_merge` - Resolve a two-face collision into the successor waves. It identifies which face is
48 rear (upstream) and which is front just before the crossing, forms the successor from
49 ``(rear.left, front.right)``, and retires the parents — a wave whose shock, contact or rarefaction face
50 merged is deactivated, while a bare fan boundary line is only marked consumed and its wave lives on.
52This file is part of gwtransport which is released under AGPL-3.0 license.
53See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
54"""
56from dataclasses import dataclass
58import numpy as np
59from scipy.optimize import brentq
61from gwtransport.fronttracking.math import NonlinearSorption, SorptionModel, characteristic_speed
62from gwtransport.fronttracking.waves import (
63 CharacteristicWave,
64 DecayingShockWave,
65 DoubleFanShockWave,
66 Feeder,
67 RarefactionWave,
68 ShockWave,
69 Wave,
70)
72EPSILON_POSITION = 1e-15
73# Largest speed used to size the Lipschitz march when a feeder range touches a saturated
74# state (R = 0, λ = +∞ for van Genuchten-Mualem at S_e = 1). The step floor plus the
75# monotonic-outlet-mass tripwire back this rare case; the sorption isotherms never hit it.
76MAX_FINITE_SPEED = 1e6
77MERGE_MATCH_TOL = 1e-9 # feeder-equality tolerance for degenerate (zero-jump) successors
80@dataclass
81class Face:
82 """One separating surface of a wave, for event detection and the reader sweep.
84 A face carries the wave it belongs to, a role tag, its position as a function of θ, the
85 left (upstream) and right (downstream) feeders, whether its trajectory is curved (a
86 fan-fed shock) and an upper bound on its speed (for the Lipschitz march).
87 """
89 wave: Wave
90 role: str # 'shock' | 'contact' | 'head' | 'tail' | 'boundary'
91 left: Feeder
92 right: Feeder
93 is_curved: bool
94 speed_bound: float
95 line: tuple[float, float, float] | None = None
96 """For a ``boundary`` face: ``(v_apex, theta_apex, speed)`` of the linear characteristic."""
98 def position(self, theta: float) -> float | None:
99 """Face position at ``theta`` (``None`` outside the wave's active θ-window)."""
100 wave = self.wave
101 if not wave.was_active_at(theta):
102 return None
103 if self.role == "boundary" and self.line is not None:
104 v_apex, theta_apex, speed = self.line
105 return v_apex + speed * (theta - theta_apex)
106 if isinstance(wave, RarefactionWave):
107 return wave.head_position_at_theta(theta) if self.role == "head" else wave.tail_position_at_theta(theta)
108 return wave.position_at_theta(theta)
111def _max_characteristic_speed(feeder: Feeder, sorption: SorptionModel) -> float:
112 """Upper bound on ``|λ|`` over a feeder's concentration range (monotone in ``c``)."""
113 if feeder.is_const:
114 s = characteristic_speed(feeder.c_a, sorption)
115 else:
116 s = max(characteristic_speed(feeder.c_a, sorption), characteristic_speed(feeder.c_b, sorption))
117 return min(s, MAX_FINITE_SPEED) if np.isfinite(s) else MAX_FINITE_SPEED
120def _far_bound(feeder: Feeder, sorption: SorptionModel, *, upstream: bool) -> float:
121 """Return the fan boundary concentration on the far (upstream/downstream) side.
123 Upstream (the apex/tail side) is the larger-retardation bound; downstream (the head
124 side) is the smaller-retardation bound. Monotonicity-agnostic (works for the n<1 mirror).
125 """
126 r_a = float(sorption.retardation(feeder.c_a))
127 r_b = float(sorption.retardation(feeder.c_b))
128 if upstream:
129 return feeder.c_a if r_a >= r_b else feeder.c_b
130 return feeder.c_a if r_a <= r_b else feeder.c_b
133def iter_faces(wave: Wave, theta: float) -> list[Face]:
134 """Enumerate the faces of ``wave`` at ``theta`` (shock/contact/boundary).
136 ``theta`` selects the *historical* boundary state: a free fan boundary line is a face
137 only for ``θ`` before it was consumed by a wave entering the fan (mirroring
138 """
139 if isinstance(wave, CharacteristicWave):
140 left = Feeder.constant(wave.concentration)
141 right = Feeder.constant(wave.c_ahead)
142 return [Face(wave, "contact", left, right, is_curved=False, speed_bound=abs(wave.speed()))]
144 if isinstance(wave, ShockWave):
145 left = Feeder.constant(wave.c_left)
146 right = Feeder.constant(wave.c_right)
147 return [Face(wave, "shock", left, right, is_curved=False, speed_bound=abs(wave.speed))]
149 if isinstance(wave, RarefactionWave):
150 fan = Feeder.fan(wave.v_start, wave.theta_start, wave.c_tail, wave.c_head, wave.sorption)
151 head = Face(
152 wave, "head", fan, Feeder.constant(wave.c_head), is_curved=False, speed_bound=abs(wave.head_speed())
153 )
154 tail = Face(
155 wave, "tail", Feeder.constant(wave.c_tail), fan, is_curved=False, speed_bound=abs(wave.tail_speed())
156 )
157 return [head, tail]
159 if isinstance(wave, DecayingShockWave):
160 return _decaying_shock_faces(wave, theta)
162 if isinstance(wave, DoubleFanShockWave):
163 return _double_fan_faces(wave, theta)
165 return []
168def _decaying_shock_faces(wave: DecayingShockWave, theta: float) -> list[Face]:
169 """Shock face (curved) plus the free fan boundary line of a decaying shock at ``theta``."""
170 s = wave.sorption
171 c_lo = min(wave.c_fan_tail, wave.c_decay_initial)
172 c_hi = max(wave.c_fan_tail, wave.c_decay_initial)
173 boundary_free = not wave.fan_boundary_consumed and theta < wave.theta_fan_boundary_consumed
174 fan = Feeder.fan(wave.v_origin, wave.theta_origin, c_lo, c_hi, s, far_boundary_free=boundary_free)
175 fixed = Feeder.constant(wave.c_fixed)
176 speed_bound = max(_max_characteristic_speed(fan, s), _max_characteristic_speed(fixed, s))
177 if wave.decay_side == "left":
178 shock = Face(wave, "shock", fan, fixed, is_curved=True, speed_bound=speed_bound)
179 else:
180 shock = Face(wave, "shock", fixed, fan, is_curved=True, speed_bound=speed_bound)
181 faces = [shock]
183 if boundary_free:
184 tail_c = Feeder.constant(wave.c_fan_tail)
185 line_speed = characteristic_speed(wave.c_fan_tail, s)
186 if np.isfinite(line_speed):
187 # A wave entering through this boundary rides into a fan whose FAR edge is this
188 # shock's own face (owned), so the entrant's successor must not re-expose a
189 # boundary there: the boundary-side fan feeder is far_boundary_free=False.
190 entry_fan = Feeder.fan(wave.v_origin, wave.theta_origin, c_lo, c_hi, s, far_boundary_free=False)
191 if wave.decay_side == "left":
192 # fan is upstream of the shock; its far (upstream) boundary carries c_fan_tail.
193 left, right = tail_c, entry_fan
194 else:
195 left, right = entry_fan, tail_c
196 faces.append(_boundary_line(wave, wave.v_origin, wave.theta_origin, line_speed, left, right))
197 return faces
200def _double_fan_faces(wave: DoubleFanShockWave, theta: float) -> list[Face]:
201 """Shock face (curved) plus the free left/right fan boundary lines of a doubly-fed shock."""
202 s = wave.sorption
203 left_free = not wave.left_boundary_consumed and theta < wave.theta_left_boundary_consumed
204 right_free = not wave.right_boundary_consumed and theta < wave.theta_right_boundary_consumed
205 left_fan = Feeder.fan(
206 wave.left_feeder.v_apex,
207 wave.left_feeder.theta_apex,
208 wave.left_feeder.c_a,
209 wave.left_feeder.c_b,
210 s,
211 far_boundary_free=left_free,
212 )
213 right_fan = Feeder.fan(
214 wave.right_feeder.v_apex,
215 wave.right_feeder.theta_apex,
216 wave.right_feeder.c_a,
217 wave.right_feeder.c_b,
218 s,
219 far_boundary_free=right_free,
220 )
221 speed_bound = max(_max_characteristic_speed(left_fan, s), _max_characteristic_speed(right_fan, s))
222 faces = [Face(wave, "shock", left_fan, right_fan, is_curved=True, speed_bound=speed_bound)]
224 if left_free:
225 c_far = _far_bound(left_fan, s, upstream=True)
226 line_speed = characteristic_speed(c_far, s)
227 if np.isfinite(line_speed):
228 entry_fan = Feeder.fan(
229 left_fan.v_apex, left_fan.theta_apex, left_fan.c_a, left_fan.c_b, s, far_boundary_free=False
230 )
231 faces.append(
232 _boundary_line(
233 wave, left_fan.v_apex, left_fan.theta_apex, line_speed, Feeder.constant(c_far), entry_fan
234 )
235 )
236 if right_free:
237 c_far = _far_bound(right_fan, s, upstream=False)
238 line_speed = characteristic_speed(c_far, s)
239 if np.isfinite(line_speed):
240 entry_fan = Feeder.fan(
241 right_fan.v_apex, right_fan.theta_apex, right_fan.c_a, right_fan.c_b, s, far_boundary_free=False
242 )
243 faces.append(
244 _boundary_line(
245 wave, right_fan.v_apex, right_fan.theta_apex, line_speed, entry_fan, Feeder.constant(c_far)
246 )
247 )
248 return faces
251def _boundary_line(wave: Wave, v_apex: float, theta_apex: float, speed: float, left: Feeder, right: Feeder) -> Face:
252 """Build a linear fan boundary characteristic face (``V = v_apex + speed·(θ − θ_apex)``)."""
253 return Face(
254 wave, "boundary", left, right, is_curved=False, speed_bound=abs(speed), line=(v_apex, theta_apex, speed)
255 )
258def find_face_crossing(
259 face_a: Face, face_b: Face, theta_current: float, theta_horizon: float
260) -> tuple[float, float] | None:
261 """First θ in ``(θ_start, θ_horizon]`` where two faces coincide, else ``None``.
263 ``θ_start = max(θ_current, both faces born)``. The gap ``g(θ) = pos_b − pos_a`` is
264 marched with a per-pair speed bound ``Λ = speed_bound_a + speed_bound_b`` so a step
265 ``Δθ = max(|g|/Λ, floor)`` cannot straddle a sign change unseen; each step also brackets
266 the gap's interior minimum to catch a grazing double-crossing.
268 The loose near-coincidence admission at birth is enabled only for cross-wave pairs (exact
269 coincidence is always reported). Same-wave exhaustion pairs — a wave's shock face against
270 its own fan boundary line — suppress it: a front born a small distance from its own
271 boundary line may be legitimately diverging, and the loose admission would fire a spurious
272 immediate exhaustion.
273 """
274 wave_a, wave_b = face_a.wave, face_b.wave
275 # Loose born-coincidence admission is a cross-wave notion (a newborn wave touching a
276 # different neighbour it must immediately merge with); a wave's own two faces suppress it.
277 allow_born_coincident = wave_a is not wave_b
278 theta_start = max(theta_current, wave_a.theta_start, wave_b.theta_start)
279 if theta_start >= theta_horizon:
280 return None
282 def gap(theta: float) -> float | None:
283 pa = face_a.position(theta)
284 pb = face_b.position(theta)
285 if pa is None or pb is None:
286 return None
287 return pb - pa
289 lam = face_a.speed_bound + face_b.speed_bound
290 lam = lam if np.isfinite(lam) and lam > 0 else MAX_FINITE_SPEED
291 floor = 1e-9 * max(abs(theta_start), 1.0)
293 theta = theta_start
294 g_prev = gap(theta)
295 if g_prev is None:
296 return None
297 # Born-coincident / already-crossed event: when a wave is created by a near-simultaneous
298 # (triple-) collision, it can be born within FP of — or just past — an adjacent wave it
299 # must immediately merge with. The gap is already ≈0 or slightly negative at the search
300 # start, so the forward march never sees a sign change. Detect it directly: if one wave
301 # was just born (θ_start == search start) and the faces coincide within a small position
302 # tolerance, emit the merge now (at the search start).
303 pa0 = face_a.position(theta_start)
304 born_now = abs(wave_a.theta_start - theta_start) <= floor or abs(wave_b.theta_start - theta_start) <= floor
305 coincidence_tol = 1e-3 * max(abs(pa0) if pa0 is not None else 0.0, 1.0)
306 if abs(g_prev) < EPSILON_POSITION or (allow_born_coincident and born_now and abs(g_prev) < coincidence_tol):
307 return (theta, pa0) if pa0 is not None else None
309 while theta < theta_horizon:
310 step = max(abs(g_prev) / lam, floor)
311 theta_next = min(theta + step, theta_horizon)
312 g_next = gap(theta_next)
313 if g_next is None:
314 return None
315 if g_prev * g_next <= 0.0:
316 # The step endpoints straddle zero (or one sits on it), so the root is
317 # bracketed; take the endpoint directly when it is already the root.
318 if abs(g_prev) < EPSILON_POSITION:
319 root = theta
320 elif abs(g_next) < EPSILON_POSITION:
321 root = theta_next
322 else:
323 root = float(brentq(gap, theta, theta_next, xtol=1e-13))
324 v = face_a.position(root)
325 return (root, v) if v is not None else None
326 theta, g_prev = theta_next, g_next
327 return None
330def make_wave_from_feeders(left: Feeder, right: Feeder, v: float, theta: float, sorption: SorptionModel) -> Wave | None:
331 """Build the successor wave a merge produces from ``(left, right)`` feeders at ``(v, θ)``.
333 ``(const, const)`` → :class:`ShockWave` (or ``None`` for a zero jump);
334 ``(const, fan)``/``(fan, const)`` → :class:`DecayingShockWave`;
335 ``(fan, fan)`` → :class:`DoubleFanShockWave`. Fan feeders whose far boundary is not free
336 (``far_boundary_free=False``) mark the successor's boundary consumed on that side.
337 """
338 if left.is_const and right.is_const:
339 if abs(left.c_a - right.c_a) < MERGE_MATCH_TOL:
340 return None
341 shock = ShockWave(theta_start=theta, v_start=v, c_left=left.c_a, c_right=right.c_a, sorption=sorption)
342 return shock if shock.satisfies_entropy() else None
344 if left.is_const or right.is_const:
345 return _make_decaying_shock(left, right, v, theta, sorption)
347 assert isinstance(sorption, NonlinearSorption) # noqa: S101
348 return DoubleFanShockWave(
349 theta_start=theta,
350 v_start=v,
351 left_feeder=left,
352 right_feeder=right,
353 sorption=sorption,
354 left_boundary_consumed=not left.far_boundary_free,
355 right_boundary_consumed=not right.far_boundary_free,
356 )
359def _make_decaying_shock(left: Feeder, right: Feeder, v: float, theta: float, sorption: SorptionModel) -> Wave | None:
360 """Successor for one const + one fan feeder → a decaying shock (or a plain shock).
362 The decaying side evolves from ``c_decay_initial`` (the fan value at the collision)
363 toward ``c_fixed``. The fan only *feeds* that motion while there is a fan edge on the
364 ``c_fixed`` side of ``c_decay_initial``; ``c_fan_tail`` is that edge (the exhaustion
365 target — the decay asymptotes at ``c_fixed`` if it lies before the edge). If no fan edge
366 lies toward ``c_fixed`` (the shock leaves the fan immediately — e.g. a collision exactly
367 at a fan boundary), there is no fan-fed decay and the successor is a plain
368 :class:`ShockWave` ``(c_decay_initial | c_fixed)``.
369 """
370 assert isinstance(sorption, NonlinearSorption) # noqa: S101
371 if left.is_const:
372 fixed_c = left.c_a
373 fan = right
374 decay_side = "right"
375 else:
376 fixed_c = right.c_a
377 fan = left
378 decay_side = "left"
379 c_decay_initial = fan.value(v, theta)
380 if abs(c_decay_initial - fixed_c) < MERGE_MATCH_TOL:
381 return None
383 direction = fixed_c - c_decay_initial
384 edges_toward_fixed = [b for b in (fan.c_a, fan.c_b) if (b - c_decay_initial) * direction > MERGE_MATCH_TOL]
385 if not edges_toward_fixed:
386 # The decay exits the fan at once — no fan-fed side, just a constant-state shock.
387 # decay_side='left' has the decaying (fan) side upstream; 'right' downstream.
388 if decay_side == "left":
389 c_left, c_right = c_decay_initial, fixed_c
390 else:
391 c_left, c_right = fixed_c, c_decay_initial
392 shock = ShockWave(theta_start=theta, v_start=v, c_left=c_left, c_right=c_right, sorption=sorption)
393 return shock if shock.satisfies_entropy() else None
394 # The fan edge farthest along the decay direction is the exhaustion target.
395 c_fan_tail = max(edges_toward_fixed, key=lambda b: (b - c_decay_initial) * direction)
397 return DecayingShockWave(
398 theta_start=theta,
399 v_start=v,
400 c_decay_initial=c_decay_initial,
401 c_fixed=fixed_c,
402 c_fan_tail=c_fan_tail,
403 decay_side=decay_side,
404 v_origin=fan.v_apex,
405 theta_origin=fan.theta_apex,
406 sorption=sorption,
407 fan_boundary_consumed=not fan.far_boundary_free,
408 )
411def resolve_merge(face_a: Face, face_b: Face, theta: float, v: float, sorption: SorptionModel) -> list[Wave]:
412 """Resolve a two-face collision: deactivate/consume parents, return successor waves.
414 Determines the rear (upstream) and front (downstream) face just before the crossing,
415 forms the successor from ``(rear.left, front.right)``, and retires the parents: a whole
416 wave is deactivated when its shock/contact/rarefaction face merges (its fan is absorbed
417 into the successor); a bare fan boundary line is *consumed* while its owning wave lives on.
418 """
419 rear, front = _order_rear_front(face_a, face_b, theta)
420 successors = make_wave_from_feeders(rear.left, front.right, v, theta, sorption)
421 new_waves = [successors] if successors is not None else []
423 _retire_parent(rear, theta)
424 _retire_parent(front, theta)
425 return new_waves
428def _order_rear_front(face_a: Face, face_b: Face, theta: float) -> tuple[Face, Face]:
429 """Return ``(rear, front)`` — the upstream (smaller V just before θ) face first."""
430 eps = 1e-7 * max(abs(theta), 1.0)
431 pa = face_a.position(theta - eps)
432 pb = face_b.position(theta - eps)
433 if pa is None or pb is None:
434 pa = face_a.position(theta)
435 pb = face_b.position(theta)
436 if pa is None or pb is None:
437 return face_a, face_b
438 return (face_a, face_b) if pa <= pb else (face_b, face_a)
441def _retire_parent(face: Face, theta: float) -> None:
442 """Deactivate the wave, or consume just its boundary line if that is the merged face."""
443 wave = face.wave
444 if face.role == "boundary":
445 _consume_boundary(wave, face, theta)
446 return
447 wave.deactivate(theta)
450def _consume_boundary(wave: Wave, face: Face, theta: float) -> None:
451 """Timestamp the crossed fan boundary as consumed at ``theta``; the wave lives on.
453 Retrospective queries at ``θ' < theta`` still see the boundary as free (historical truth),
454 so a later reader query does not retro-erase the boundary before it was actually consumed.
455 """
456 if isinstance(wave, DecayingShockWave):
457 wave.theta_fan_boundary_consumed = theta
458 elif isinstance(wave, DoubleFanShockWave):
459 line = face.line
460 matches_left = (
461 line is not None
462 and abs(line[0] - wave.left_feeder.v_apex) < EPSILON_POSITION
463 and abs(line[1] - wave.left_feeder.theta_apex) < EPSILON_POSITION
464 )
465 if matches_left:
466 wave.theta_left_boundary_consumed = theta
467 else:
468 wave.theta_right_boundary_consumed = theta