Coverage for src/gwtransport/fronttracking/waves.py: 90%
598 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 21:13 +0000
1"""
2Wave Representation for Front Tracking in (V, θ) coordinates.
4This module implements wave classes for representing characteristics, shocks,
5and rarefaction waves in the front tracking algorithm. Each wave stores its
6formation position in cumulative-flow coordinate ``θ = ∫flow(t') dt'`` and
7knows how to compute its position at any later θ.
9In (V, θ) every wave velocity is a property of the sorption isotherm alone —
10flow does not enter wave dynamics. Time-varying flow is absorbed entirely into
11the θ(t) mapping at the API boundary, so no wave needs recreation when the flow
12rate changes. Attributes named ``theta_*`` are therefore cumulative flow [m³],
13never time; ``v_*`` are volumetric positions [m³] measured from the inlet face
14``V = 0``. Every module in :mod:`gwtransport.fronttracking` follows this convention.
16Available functions:
18- :class:`Wave` - Abstract base of the hierarchy. Holds the formation point ``(v_start, theta_start)``, the
19 ``is_active`` flag and the ``theta_deactivation`` history marker, and requires ``position_at_theta``,
20 ``concentration_left``, ``concentration_right`` and ``concentration_at_point`` from every subclass. Use
21 ``was_active_at(theta)`` for retrospective queries; ``is_active`` describes only the current state.
23- :class:`CharacteristicWave` - A single characteristic line carrying one constant concentration at the
24 flow-free speed ``1/R(c)``. Physically a contact discontinuity in a smooth region: the concentration value
25 travels unchanged, and ``c_ahead`` records the state it separates from downstream.
27- :class:`ShockWave` - A sharp front across which concentration jumps from ``c_left`` to ``c_right``, formed
28 where faster water overtakes slower water. Its ``speed`` is the Rankine-Hugoniot secant
29 ``(c_R − c_L)/(C_T(c_R) − C_T(c_L))``, constant in θ; ``satisfies_entropy`` tests the Lax condition.
31- :class:`RarefactionWave` - An expansion fan spreading between ``c_head`` and ``c_tail``, formed where
32 slower water follows faster water. The interior is self-similar, ``R(c) = (θ − θ_start)/(V − v_start)``, so
33 head and tail move at the characteristic speeds ``1/R(c_head)`` and ``1/R(c_tail)``; construction is
34 rejected when the head is not faster than the tail (that geometry is a compression, hence a shock).
36- :class:`DecayingShockWave` - A shock with a fan on one side, formed when a rarefaction and a shock collide.
37 The fan feeds the ``decay_side``, whose concentration relaxes from ``c_decay_initial`` toward the fan's far
38 bound ``c_fan_tail`` while the other side stays at ``c_fixed``, so the front decelerates along a curved
39 trajectory instead of a straight line. The trajectory uses a closed form where one exists (Freundlich with
40 ``c_fixed = 0`` or ``n = 2``, Langmuir and Brooks-Corey with ``c_fixed = 0``) and a cached quadrature
41 profile otherwise; ``theta_at_fan_exhaustion`` reports the θ at which the fan is spent and the wave hands
42 off to a plain shock.
44- :class:`DoubleFanShockWave` - A shock fed by a self-similar fan on *both* sides, formed when a fan boundary
45 catches a decaying shock whose surviving side is itself a fan. Both side concentrations then relax as the
46 front advances. The trajectory is closed form for Freundlich ``n = 2`` with a shared apex position (the
47 case every inlet-born fan produces) and a cached RK4 spline otherwise; a side ending is detected
48 externally, as the shock face crossing its own fan boundary line.
50- :class:`Feeder` - The boundary state on one side of a front: either a constant concentration or a bounded
51 self-similar fan with apex ``(v_apex, theta_apex)`` spanning ``[c_a, c_b]``, evaluated by inverting
52 ``R = (θ − θ_apex)/(v − v_apex)`` and clamped in ``R``-space so a query past the fan edge reads the plateau
53 value. Feeders are the common currency of the interaction calculus in
54 :mod:`gwtransport.fronttracking.interactions`.
56This file is part of gwtransport which is released under AGPL-3.0 license.
57See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
58"""
60from abc import ABC, abstractmethod
61from dataclasses import dataclass, field
63import numpy as np
64from scipy.interpolate import CubicSpline
65from scipy.optimize import brentq
67from gwtransport.fronttracking.math import (
68 _C_MIN,
69 BrooksCoreyConductivity,
70 FreundlichSorption,
71 LangmuirSorption,
72 NonlinearSorption,
73 SorptionModel,
74 characteristic_speed,
75)
77# Numerical tolerance constants
78EPSILON_POSITION = 1e-15 # Tolerance for checking if two positions are equal
79DECAYING_SHOCK_U_FLOOR = 1e-300 # Lower bracket bound for brentq on Freundlich u-invariant
80DECAYING_SHOCK_BRENTQ_XTOL = (
81 1e-14 # brentq absolute tolerance for monotone θ inversions (exhaustion, outlet, numerical)
82)
83# Cached numerical decay profile (see ``_build_decay_profile``): c-grid resolution, the
84# Gauss-Legendre panel order for the cumulative invariant integral, and the fraction of the
85# c-gap the grid stops short of a secant-speed pole (where ``θ_local → ∞``).
86DECAY_PROFILE_NODES = 6000
87DECAY_PROFILE_GAUSS_ORDER = 10
88DECAY_PROFILE_POLE_FLOOR = 1e-6
89# DoubleFanShockWave numerical-trajectory RK4 substeps per unit of fan age (only used when
90# no closed form applies — distinct fan apex positions or a non-n=2 isotherm). The
91# self-similar fans vary on the scale of their age, so the step is age/this-many.
92DFSW_RK_SUBSTEPS = 512
95@dataclass(frozen=True)
96class Feeder:
97 """One side's boundary state feeding a front: a constant, or a bounded self-similar fan.
99 A ``const`` feeder ignores ``(v, θ)`` and returns its value everywhere. A ``fan``
100 feeder evaluates the self-similar retardation ``R = (θ − θ_apex)/(v − v_apex)`` and
101 inverts it to a concentration, clamped to the fan's physical extent ``[c_a, c_b]``.
102 The clamp is monotonicity-agnostic (it clamps in ``R``-space, so it is correct for
103 both R-decreasing isotherms — Freundlich ``n>1``, Langmuir, Brooks-Corey,
104 van Genuchten-Mualem — and the R-increasing Freundlich ``n<1`` mirror), and it is
105 exactly what makes a fan feeder read the plateau concentration beyond the fan's edge.
107 Feeders are the uniform currency of the interaction calculus: every wave exposes its
108 sides as feeders, event handlers form a successor from ``(rear.left, front.right)``,
109 and the reader evaluates the left feeder of the nearest downstream face.
110 """
112 c_a: float
113 """One physical boundary concentration of the fan (or the constant value)."""
114 c_b: float
115 """The other physical boundary concentration of the fan (unused for a constant)."""
116 is_const: bool
117 """Whether this is a constant state (``True``) or a self-similar fan (``False``)."""
118 v_apex: float = 0.0
119 """Fan apex position [m³] (ignored for a constant)."""
120 theta_apex: float = 0.0
121 """Fan apex cumulative flow [m³] (ignored for a constant)."""
122 sorption: SorptionModel | None = None
123 """Sorption model used to invert ``R`` (required for a fan; ``None`` for a constant)."""
124 far_boundary_free: bool = True
125 """Whether the fan's far edge is a free plateau boundary (a collision line) rather than
126 already terminated by another shock. Propagated through merges so a wave born onto a
127 fan whose far end another wave owns does not re-expose a phantom boundary line."""
129 @classmethod
130 def constant(cls, c: float) -> "Feeder":
131 """Return a constant-concentration feeder."""
132 return cls(c_a=c, c_b=c, is_const=True)
134 @classmethod
135 def fan(
136 cls,
137 v_apex: float,
138 theta_apex: float,
139 c_a: float,
140 c_b: float,
141 sorption: SorptionModel,
142 *,
143 far_boundary_free: bool = True,
144 ) -> "Feeder":
145 """Return a bounded self-similar fan feeder with apex ``(v_apex, theta_apex)`` spanning ``[c_a, c_b]``."""
146 return cls(
147 c_a=c_a,
148 c_b=c_b,
149 is_const=False,
150 v_apex=v_apex,
151 theta_apex=theta_apex,
152 sorption=sorption,
153 far_boundary_free=far_boundary_free,
154 )
156 def value(self, v: float, theta: float) -> float:
157 """Concentration this feeder supplies at ``(v, θ)`` (clamped to the fan extent)."""
158 if self.is_const:
159 return self.c_a
160 sorption = self.sorption
161 assert sorption is not None # noqa: S101 # a fan feeder always carries its sorption model
162 r_a = float(sorption.retardation(self.c_a))
163 r_b = float(sorption.retardation(self.c_b))
164 # c at the smaller-R (faster, "head") and larger-R (slower, "tail"/apex) ends.
165 if r_a <= r_b:
166 r_lo, r_hi, c_at_lo, c_at_hi = r_a, r_b, self.c_a, self.c_b
167 else:
168 r_lo, r_hi, c_at_lo, c_at_hi = r_b, r_a, self.c_b, self.c_a
169 if theta <= self.theta_apex or v <= self.v_apex:
170 return c_at_hi # at/behind the apex: the largest-R (tail) end
171 r = (theta - self.theta_apex) / (v - self.v_apex)
172 if r <= r_lo:
173 return c_at_lo
174 if r >= r_hi:
175 return c_at_hi
176 return float(sorption.concentration_from_retardation(r))
179@dataclass
180class Wave(ABC):
181 """Abstract base class for all wave types in front tracking.
183 All waves share common attributes and must implement methods for
184 computing position and concentration. Waves can be active or inactive
185 (deactivated waves are preserved for history but don't participate in
186 future interactions).
187 """
189 theta_start: float
190 """Cumulative flow at which the wave forms [m³]."""
191 v_start: float
192 """Position at which the wave forms [m³]."""
193 is_active: bool = field(default=True, kw_only=True)
194 """Whether wave is currently active (in the solver's event-loop sense)."""
195 theta_deactivation: float = field(default=float("inf"), kw_only=True)
196 """Cumulative flow at which the wave was deactivated (default ``+∞``).
198 Historical record set by collision handlers when a wave is replaced
199 (e.g., a parent rarefaction superseded by a ``DecayingShockWave``).
200 ``is_active = False`` is the "current state" flag the solver uses for
201 its event loop; ``theta_deactivation`` is the moment in θ-history when
202 the wave stopped contributing. Retrospective queries (any θ in the
203 past) must use ``was_active_at(theta)`` instead of ``is_active`` so
204 that ``compute_domain_mass`` etc. correctly attribute c at v_outlet
205 during the wave's lifetime even after later events have deactivated
206 the wave.
207 """
209 def was_active_at(self, theta: float) -> bool:
210 """Whether the wave was active at cumulative flow ``theta`` (geometric truth).
212 Use for retrospective queries — ``is_active`` reflects only the
213 wave's *current* (post-simulation) state, which is wrong for
214 ``compute_domain_mass`` and similar at θ before a deactivation event.
216 Parameters
217 ----------
218 theta : float
219 Cumulative flow at which to query historical activity [m³].
221 Returns
222 -------
223 bool
224 ``True`` for ``theta_start <= theta < theta_deactivation``.
225 A wave constructed with ``is_active=False`` and no recorded
226 ``theta_deactivation`` (default ``+∞``) is treated as
227 never-active — e.g., synthetic test fixtures that want the
228 wave excluded from dispatch entirely.
229 """
230 if not self.is_active and self.theta_deactivation == float("inf"):
231 return False
232 return self.theta_start <= theta < self.theta_deactivation
234 def deactivate(self, theta: float) -> None:
235 """Mark the wave inactive at cumulative flow ``theta`` (collision handler API).
237 Sets both ``is_active = False`` (solver event-loop flag) and
238 ``theta_deactivation = theta`` (historical record for retrospective
239 ``was_active_at`` queries).
241 Parameters
242 ----------
243 theta : float
244 Cumulative flow at which the wave is deactivated [m³].
245 """
246 self.is_active = False
247 self.theta_deactivation = theta
249 @abstractmethod
250 def position_at_theta(self, theta: float) -> float | None:
251 """Compute wave position at cumulative flow θ.
253 Parameters
254 ----------
255 theta : float
256 Cumulative flow [m³].
258 Returns
259 -------
260 position : float or None
261 Position [m³], or None if θ < θ_start or θ >= theta_deactivation.
262 (Past-θ queries respect the wave's historical lifetime; current-state
263 queries before deactivation behave identically to the ``is_active``
264 check.)
265 """
267 @abstractmethod
268 def concentration_left(self) -> float:
269 """Concentration on the left (upstream) side of the wave."""
271 @abstractmethod
272 def concentration_right(self) -> float:
273 """Concentration on the right (downstream) side of the wave."""
275 @abstractmethod
276 def concentration_at_point(self, v: float, theta: float) -> float | None:
277 """Compute concentration at point (v, θ) if the wave controls it.
279 Returns
280 -------
281 concentration : float or None
282 Concentration [mass/volume] if the wave controls this point, None
283 otherwise.
284 """
287@dataclass
288class CharacteristicWave(Wave):
289 """Characteristic line along which concentration is constant.
291 In smooth regions, concentration travels at speed ``1/R(C)`` in (V, θ)
292 coordinates. Along each characteristic line, the concentration value is
293 constant. This is the fundamental solution element for hyperbolic
294 conservation laws.
296 Examples
297 --------
298 >>> sorption = FreundlichSorption(
299 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
300 ... )
301 >>> char = CharacteristicWave(
302 ... theta_start=0.0, v_start=0.0, concentration=5.0, sorption=sorption
303 ... )
304 >>> speed = char.speed()
305 >>> bool(np.isclose(char.position_at_theta(1000.0), speed * 1000.0))
306 True
307 """
309 concentration: float
310 """Constant concentration carried on the upstream (behind) side [mass/volume]."""
311 sorption: SorptionModel
312 """Sorption model determining the speed."""
313 c_ahead: float = field(default=0.0, kw_only=True)
314 """Concentration on the downstream (ahead) side — the state the contact advances into.
316 A contact separates the carried ``concentration`` (behind, upstream) from ``c_ahead``
317 (ahead, downstream, the pre-existing state). The solver sets it to the previous inlet
318 value; it defaults to ``0`` (the virgin initial condition) for a lone contact. The
319 reader sweep uses it as the contact's downstream feeder."""
320 _speed: float = field(init=False, repr=False, compare=False)
321 """Cached characteristic speed (immutable inputs; set in ``__post_init__``)."""
323 def __post_init__(self) -> None:
324 """Cache the (immutable) characteristic speed once."""
325 self._speed = characteristic_speed(self.concentration, self.sorption)
327 def speed(self) -> float:
328 """Characteristic speed dV/dθ = 1/R(C) (``+∞`` at a saturated state, R = 0)."""
329 return self._speed
331 def position_at_theta(self, theta: float) -> float | None:
332 """Position at cumulative flow θ.
334 ``V(θ) = v_start + speed * (θ - θ_start)``.
335 """
336 if not self.was_active_at(theta):
337 return None
338 return self.v_start + self.speed() * (theta - self.theta_start)
340 def concentration_left(self) -> float:
341 """Concentration on the left (upstream) side; equals the carried value."""
342 return self.concentration
344 def concentration_right(self) -> float:
345 """Concentration on the right (downstream) side; equals the carried value."""
346 return self.concentration
348 def concentration_at_point(self, v: float, theta: float) -> float | None:
349 """Return the carried concentration if the characteristic has reached ``v`` by θ."""
350 v_at_theta = self.position_at_theta(theta)
351 if v_at_theta is None:
352 return None
354 if v_at_theta >= v:
355 return self.concentration
357 return None
360@dataclass
361class ShockWave(Wave):
362 """Shock wave (discontinuity) with jump in concentration.
364 Shocks form when faster water overtakes slower water, creating a sharp
365 front. In (V, θ) the shock speed is given by the Rankine-Hugoniot
366 condition and is independent of flow::
368 dV_s/dθ = (C_R - C_L) / (C_T(C_R) - C_T(C_L))
370 Examples
371 --------
372 >>> sorption = FreundlichSorption(
373 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
374 ... )
375 >>> shock = ShockWave(
376 ... theta_start=0.0,
377 ... v_start=0.0,
378 ... c_left=10.0,
379 ... c_right=2.0,
380 ... sorption=sorption,
381 ... )
382 >>> shock.speed > 0
383 True
384 >>> shock.satisfies_entropy()
385 True
386 """
388 c_left: float
389 """Concentration upstream (behind) shock [mass/volume]."""
390 c_right: float
391 """Concentration downstream (ahead of) shock [mass/volume]."""
392 sorption: SorptionModel
393 """Sorption model."""
394 speed: float = field(init=False)
395 """Shock speed dV/dθ; set in ``__post_init__``."""
397 def __post_init__(self) -> None:
398 """Compute shock speed from Rankine-Hugoniot in (V, θ)."""
399 self.speed = self.sorption.shock_speed(self.c_left, self.c_right)
401 def position_at_theta(self, theta: float) -> float | None:
402 """Position at cumulative flow θ. Shock propagates linearly in θ."""
403 if not self.was_active_at(theta):
404 return None
405 return self.v_start + self.speed * (theta - self.theta_start)
407 def concentration_left(self) -> float:
408 """Upstream concentration of the shock."""
409 return self.c_left
411 def concentration_right(self) -> float:
412 """Downstream concentration of the shock."""
413 return self.c_right
415 def concentration_at_point(self, v: float, theta: float) -> float | None:
416 """Return c_left if upstream of the shock at θ, c_right if downstream.
418 At the exact shock position the average is returned (convention; the
419 shock is infinitesimally thin in practice).
420 """
421 v_shock = self.position_at_theta(theta)
422 if v_shock is None:
423 return None
425 # Position-scaled face width (~1 ULP at all positions), matching
426 # DecayingShockWave.concentration_at_point; a fixed 1e-15 falls below
427 # one ULP for any v_shock > ~1 m³ and degenerates to bit-equality.
428 tol = 1e-15 * max(abs(v_shock), 1.0)
430 if v < v_shock - tol:
431 return self.c_left
432 if v > v_shock + tol:
433 return self.c_right
434 return 0.5 * (self.c_left + self.c_right)
436 def satisfies_entropy(self) -> bool:
437 """Check Lax entropy condition in (V, θ): ``λ_θ(C_L) ≥ s ≥ λ_θ(C_R)``."""
438 return self.sorption.check_entropy_condition(self.c_left, self.c_right, self.speed)
441@dataclass
442class RarefactionWave(Wave):
443 """Rarefaction (expansion fan) with smooth concentration gradient.
445 Rarefactions form when slower water follows faster water, creating an
446 expanding region where concentration varies smoothly. In (V, θ) the
447 solution is self-similar in ``(V - v_start)`` vs ``(θ - θ_start)``::
449 R(C) = (θ - θ_start) / (V - v_start)
451 Head and tail propagate at flow-free speeds ``1/R(C_head)`` and
452 ``1/R(C_tail)``.
454 Raises
455 ------
456 ValueError
457 If head speed <= tail speed (would be a compression, not a rarefaction).
459 Examples
460 --------
461 >>> sorption = FreundlichSorption(
462 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
463 ... )
464 >>> raref = RarefactionWave(
465 ... theta_start=0.0,
466 ... v_start=0.0,
467 ... c_head=10.0,
468 ... c_tail=2.0,
469 ... sorption=sorption,
470 ... )
471 >>> raref.head_speed() > raref.tail_speed()
472 True
473 >>> raref.contains_point(v=150.0, theta=2000.0)
474 True
475 """
477 c_head: float
478 """Concentration at leading edge (faster) [mass/volume]."""
479 c_tail: float
480 """Concentration at trailing edge (slower) [mass/volume]."""
481 sorption: SorptionModel
482 """Sorption model (must be concentration-dependent)."""
483 _head_speed: float = field(init=False, repr=False, compare=False)
484 """Cached head celerity (immutable inputs; set in ``__post_init__``)."""
485 _tail_speed: float = field(init=False, repr=False, compare=False)
486 """Cached tail celerity (immutable inputs; set in ``__post_init__``)."""
488 def __post_init__(self):
489 """Cache head/tail celerities and verify this is a rarefaction (head faster than tail)."""
490 self._head_speed = characteristic_speed(self.c_head, self.sorption)
491 self._tail_speed = characteristic_speed(self.c_tail, self.sorption)
493 if self._head_speed <= self._tail_speed:
494 msg = (
495 f"Not a rarefaction: head_speed={self._head_speed:.6g} <= tail_speed={self._tail_speed:.6g}. "
496 f"This would be a compression (shock) instead."
497 )
498 raise ValueError(msg)
500 def head_speed(self) -> float:
501 """Speed of rarefaction head dV/dθ = 1/R(C_head) (``+∞`` at a saturated state, R = 0)."""
502 return self._head_speed
504 def tail_speed(self) -> float:
505 """Speed of rarefaction tail dV/dθ = 1/R(C_tail) (``+∞`` at a saturated state, R = 0)."""
506 return self._tail_speed
508 def head_position_at_theta(self, theta: float) -> float | None:
509 """Position of rarefaction head at cumulative flow θ."""
510 if not self.was_active_at(theta):
511 return None
512 return self.v_start + self.head_speed() * (theta - self.theta_start)
514 def tail_position_at_theta(self, theta: float) -> float | None:
515 """Position of rarefaction tail at cumulative flow θ."""
516 if not self.was_active_at(theta):
517 return None
518 return self.v_start + self.tail_speed() * (theta - self.theta_start)
520 def position_at_theta(self, theta: float) -> float | None:
521 """Head position (leading edge of rarefaction). Implements abstract Wave method."""
522 return self.head_position_at_theta(theta)
524 def contains_point(self, v: float, theta: float) -> bool:
525 """Return ``True`` if ``(v, θ)`` lies between the fan's tail and head."""
526 if theta <= self.theta_start or theta >= self.theta_deactivation:
527 return False
529 v_head = self.head_position_at_theta(theta)
530 v_tail = self.tail_position_at_theta(theta)
532 if v_head is None or v_tail is None:
533 return False
535 return v_tail <= v <= v_head
537 def concentration_left(self) -> float:
538 """Upstream concentration is the trailing-edge value c_tail."""
539 return self.c_tail
541 def concentration_right(self) -> float:
542 """Downstream concentration is the leading-edge value c_head."""
543 return self.c_head
545 def concentration_at_point(self, v: float, theta: float) -> float | None:
546 """Self-similar concentration inside the fan: ``R(C) = (θ - θ_start)/(v - v_start)``.
548 Outside the fan returns None. For ``ConstantRetardation``, rarefactions
549 don't form (all concentrations travel at the same speed), so this also
550 returns None.
552 Examples
553 --------
554 >>> sorption = FreundlichSorption(
555 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
556 ... )
557 >>> raref = RarefactionWave(0.0, 0.0, 10.0, 2.0, sorption)
558 >>> c = raref.concentration_at_point(v=150.0, theta=2000.0)
559 >>> c is not None
560 True
561 >>> 2.0 <= c <= 10.0
562 True
563 """
564 if abs(v - self.v_start) < EPSILON_POSITION and theta >= self.theta_start:
565 return self.c_tail
567 if not self.contains_point(v, theta):
568 return None
570 r_target = (theta - self.theta_start) / (v - self.v_start)
572 if r_target <= 1.0:
573 return None # Unphysical
575 try:
576 c = self.sorption.concentration_from_retardation(r_target)
577 except NotImplementedError:
578 # ConstantRetardation case — rarefactions don't form
579 return None
581 # contains_point(v, theta) was True, so the point is geometrically inside
582 # the fan. The inverted c may drift by a few ULPs past [c_tail, c_head]
583 # — clamp rather than rejecting so callers at the head/tail boundaries
584 # get the correct boundary concentration.
585 c_lo = min(self.c_tail, self.c_head)
586 c_hi = max(self.c_tail, self.c_head)
587 return min(max(float(c), c_lo), c_hi)
590@dataclass
591class DecayingShockWave(Wave):
592 r"""Merging shock with closed-form (or quadrature) trajectory in θ-space.
594 Formed when a rarefaction fan and a shock collide. The shock then has
595 one side fed by the fan's self-similar profile (the "decay" side) and
596 the other side at the original outer state (the "fixed" side). Valid for
597 any :class:`~gwtransport.fronttracking.math.NonlinearSorption`.
599 Two collision regimes are supported via ``decay_side``:
601 - ``'left'`` (favorable head-collision): the rarefaction's head (faster)
602 catches a leading shock. After collision, the shock's ``c_left`` decays
603 from the rarefaction head value toward ``c_fan_tail`` (the unchanged
604 downstream c_right is ``c_fixed``).
605 - ``'right'`` (unfavorable tail-collision, n<1 mirrored): a trailing shock
606 catches the rarefaction's tail. After collision, the shock's ``c_right``
607 decays from the rarefaction tail value toward ``c_fan_tail`` (the
608 unchanged upstream c_left is ``c_fixed``).
610 The wave is valid only while ``c_decay ∈ (c_fan_tail, c_decay_initial]``;
611 once ``c_decay`` reaches ``c_fan_tail`` the fan is exhausted (see the
612 solver's ``DSW_FAN_EXHAUSTED`` event).
614 **Dispatch.** ``_c_decay_at_theta_local`` is the single dispatch site
615 (position, fan-exhaustion and outlet-crossing all route through it): a
616 closed form is used where one exists, otherwise the per-wave cached numerical
617 profile (:func:`_build_decay_profile`). No combination raises — any
618 :class:`~gwtransport.fronttracking.math.NonlinearSorption` is valid. With
619 ``θ_local := θ − theta_origin`` measured from the rarefaction apex,
620 ``α := ρ_b · k_f / n_por`` for Freundlich, and ``u_d := c_decay^(1/n)``:
622 - Freundlich, ``c_fixed = 0`` (general ``n > 0``, ``n ≠ 1``) — closed form:
623 invariant ``θ_local · u_d^n = K · (n · u_d^(n-1) + α)``,
624 position ``V_s(θ) = v_origin + n · K / u_d(θ)``.
625 - Freundlich, ``c_fixed > 0``, ``n = 2`` (either decay orientation) — closed form:
626 invariant ``(u_d - u_R)² · θ_local = K · (2 u_d + α)`` with ``u_R := c_fixed^(1/2)``,
627 position ``V_s(θ) = v_origin + 2 K · u_d(θ) / (u_d - u_R)²``. The root of the
628 quadratic in ``u_d`` is selected by the decay orientation (``u_d > u_R`` shrinking,
629 ``u_d < u_R`` growing); both are exact (verified to ~1e-14 vs a DOP853 integration).
630 - Langmuir, ``c_fixed = 0`` — closed form:
631 invariant ``θ_local · c_d² = K · ((K_L + c_d)² + a)`` with
632 ``a := ρ_b · s_max · K_L / n_por``,
633 position ``V_s(θ) = v_origin + K · (K_L + c_d)² / c_d²``.
634 - Brooks-Corey, ``c_fixed = 0`` — closed form:
635 invariant ``θ_local ∝ R(c_decay)^{a/(a−1)}`` (``R·S = 1/a`` constant),
636 so ``R(c_d) = R(c0)·(θ_local/θ_local_coll)^{(a−1)/a}``.
637 - Every other ``(isotherm, c_fixed)`` combination (Freundlich ``c_fixed>0,
638 n≠2``, Langmuir/Brooks-Corey ``c_fixed>0``, any van Genuchten) — cached
639 numerical profile (:func:`_build_decay_profile`): the decay-agnostic
640 invariant ``θ_local(c_d) = θ_local_coll · exp(∫ R'/[(1 − R·S)·R] dc)`` with
641 the symmetric secant speed ``S = (c − c_fixed)/(C_T(c) − C_T(c_fixed))``,
642 built once by composite quadrature and inverted for ``c_d(θ)`` by monotone
643 spline interpolation.
645 Every path shares the fan-continuity identity
646 ``V_s = v_origin + θ_local / R(c_decay)``, which ``position_at_theta`` and
647 ``outlet_crossing_theta`` use uniformly across all isotherms.
649 The invariant constant ``K`` (closed-form Freundlich/Langmuir only) is set
650 in ``__post_init__`` from the collision IC ``(theta_start, c_decay_initial)``.
651 ``theta_start``/``v_start`` are the collision coordinates; a fan-consistent
652 construction has ``v_start = v_origin + (theta_start − theta_origin)/R(c_decay_initial)``.
654 See Also
655 --------
656 ShockWave : Linear-θ shock (no decaying side).
657 RarefactionWave : Self-similar expansion fan.
658 """
660 c_decay_initial: float
661 """Concentration on the decaying side at θ=theta_start [mass/volume]. Non-negative.
663 A fully-drained collision value of ``0`` is floored to the shared dry-soil singularity
664 floor ``_C_MIN`` so the retardation and secant-speed evaluations stay finite."""
665 c_fixed: float
666 """Concentration on the non-decaying side [mass/volume]. Non-negative, constant in θ."""
667 c_fan_tail: float
668 """Concentration at the fan's far boundary [mass/volume]; bounds the decay. Non-negative.
670 The wave is valid only while ``c_decay ∈ (c_fan_tail, c_decay_initial]``; at
671 ``c_fan_tail`` the fan is exhausted."""
672 decay_side: str
673 """``'left'`` (favorable head-collision) or ``'right'`` (n<1 mirrored)."""
674 v_origin: float
675 """Position of the rarefaction apex [m³]."""
676 theta_origin: float
677 """Cumulative flow at the rarefaction apex [m³]; strictly less than ``theta_start``."""
678 sorption: NonlinearSorption
679 """Sorption model (any concentration-dependent isotherm)."""
680 K: float = field(init=False)
681 """Invariant constant set in ``__post_init__`` (closed-form Freundlich ``c_fixed=0``/``n≈2`` and Langmuir
682 ``c_fixed=0`` cases; ``nan`` for every numerical case)."""
683 _freundlich_cf: bool = field(init=False, repr=False, compare=False)
684 """Cached Freundlich-closed-form predicate (immutable inputs; set in ``__post_init__``)."""
685 _langmuir_cf: bool = field(init=False, repr=False, compare=False)
686 """Cached Langmuir-closed-form predicate."""
687 _brooks_corey_cf: bool = field(init=False, repr=False, compare=False)
688 """Cached Brooks-Corey ``c_fixed=0`` closed-form predicate."""
689 _numerical: bool = field(init=False, repr=False, compare=False)
690 """Cached predicate: no closed form applies, so the decay routes to the cached numerical profile."""
691 _decay_profile_cache: tuple | None = field(default=None, init=False, repr=False, compare=False)
692 """Lazily-built monotone ``θ_local(c)`` map for the numerical decay path (see ``_decay_profile``)."""
693 fan_boundary_consumed: bool = field(default=False, kw_only=True)
694 """Whether the fan's far-boundary line is owned from birth (a fan-entry successor rides a
695 fan another wave already terminates downstream, so its boundary is never a free face)."""
696 theta_fan_boundary_consumed: float = field(default=float("inf"), kw_only=True)
697 """Cumulative flow at which a *free* boundary line was later consumed by a wave entering
698 the fan. Retrospective reader/event queries treat the boundary as free only for
699 ``θ < theta_fan_boundary_consumed`` (historical truth, mirroring ``was_active_at``)."""
701 def __post_init__(self) -> None:
702 """Validate inputs and compute the closed-form invariant K when applicable."""
703 if self.decay_side not in {"left", "right"}:
704 msg = f"decay_side must be 'left' or 'right', got {self.decay_side!r}"
705 raise ValueError(msg)
706 if self.c_decay_initial < 0.0:
707 msg = f"c_decay_initial must be non-negative, got {self.c_decay_initial}"
708 raise ValueError(msg)
709 # Floor a fully-drained fan tail (c_decay_initial == 0) to _C_MIN so the
710 # retardation and secant-speed evaluations stay finite (package floor convention).
711 self.c_decay_initial = max(self.c_decay_initial, _C_MIN)
712 if self.c_fixed < 0.0:
713 msg = f"c_fixed must be non-negative, got {self.c_fixed}"
714 raise ValueError(msg)
715 if self.c_fan_tail < 0.0:
716 msg = f"c_fan_tail must be non-negative, got {self.c_fan_tail}"
717 raise ValueError(msg)
718 if self.theta_origin >= self.theta_start:
719 msg = (
720 f"theta_origin ({self.theta_origin}) must be strictly less than "
721 f"theta_start ({self.theta_start}); rarefaction apex precedes collision"
722 )
723 raise ValueError(msg)
725 if not isinstance(self.sorption, NonlinearSorption):
726 msg = f"DecayingShockWave requires a NonlinearSorption, got {type(self.sorption).__name__}"
727 raise TypeError(msg)
729 # Classify the decay path once (immutable inputs). Closed forms exist for:
730 # Freundlich c_fixed=0 (general n) or the n≈2 quadratic for either decay orientation
731 # (shrinking c_decay_initial > c_fixed picks the +√ root, growing c_decay_initial <
732 # c_fixed the −√ root; both invert the same (u_d−u_R)² invariant, verified exact to
733 # ~1e-14 vs a DOP853 integration of the fan-fed shock ODE);
734 # Langmuir c_fixed=0; Brooks-Corey c_fixed=0. Everything else is numerical.
735 s = self.sorption
736 self._freundlich_cf = isinstance(s, FreundlichSorption) and (
737 self.c_fixed == 0.0 or bool(np.isclose(s.n, 2.0, rtol=1e-12) and self.c_decay_initial != self.c_fixed)
738 )
739 self._langmuir_cf = isinstance(s, LangmuirSorption) and self.c_fixed == 0.0
740 self._brooks_corey_cf = isinstance(s, BrooksCoreyConductivity) and self.c_fixed == 0.0
741 self._numerical = not (self._freundlich_cf or self._langmuir_cf or self._brooks_corey_cf)
743 # K is the closed-form invariant constant; set only for the Freundlich/Langmuir
744 # closed forms and left NaN for every numerical (and Brooks-Corey) case. The
745 # isinstance guards narrow ``s`` for the typed helpers (the cached predicate already
746 # implies the type; the ``np.isclose`` cost is not re-incurred).
747 self.K = float("nan")
748 if self._freundlich_cf and isinstance(s, FreundlichSorption):
749 self.K = _compute_k_freundlich(
750 s,
751 self.theta_start - self.theta_origin,
752 self.c_decay_initial,
753 self.c_fixed,
754 )
755 elif self._langmuir_cf and isinstance(s, LangmuirSorption):
756 self.K = _compute_k_langmuir(
757 s,
758 self.theta_start - self.theta_origin,
759 self.c_decay_initial,
760 )
762 def _decay_profile(self) -> tuple:
763 """Lazily build & cache the monotone ``θ_local(c)`` map for the numerical decay path.
765 Returns ``(c_of_i, i_max, c_limit_node)``: a ``CubicSpline`` mapping the
766 cumulative invariant ``I = ln(θ_local/θ_local_coll)`` to ``c_decay``, the
767 largest resolved ``I`` (endpoint of the reachable c-range), and the c at
768 that endpoint. Built once per wave (see :func:`_build_decay_profile`).
769 """
770 if self._decay_profile_cache is None:
771 self._decay_profile_cache = _build_decay_profile(
772 self.sorption,
773 self.c_decay_initial,
774 self.c_fixed,
775 self.c_fan_tail,
776 )
777 return self._decay_profile_cache
779 def c_decay_at_theta(self, theta: float) -> float | None:
780 """Concentration on the decaying side at cumulative flow θ.
782 Returns ``None`` for ``θ < theta_start`` or when the wave is inactive;
783 otherwise delegates to the single per-isotherm dispatch in
784 ``_c_decay_at_theta_local``.
785 """
786 if not self.was_active_at(theta):
787 return None
788 return self._c_decay_at_theta_local(theta - self.theta_origin)
790 def position_at_theta(self, theta: float) -> float | None:
791 """Shock position ``V_s(θ)`` via the fan-continuity identity.
793 ``V_s = v_origin + θ_local / R(c_decay)`` for every isotherm. Returns
794 ``None`` for ``θ < theta_start`` or when inactive.
795 """
796 if not self.was_active_at(theta):
797 return None
799 theta_local = theta - self.theta_origin
800 c_d = self._c_decay_at_theta_local(theta_local)
801 return float(self.v_origin + theta_local / float(self.sorption.retardation(c_d)))
803 def theta_at_fan_exhaustion(self) -> float | None:
804 """Cumulative flow θ at which ``c_decay`` reaches ``c_fan_tail``.
806 ``c_decay(θ)`` is strictly monotone from ``c_decay_initial`` toward
807 ``c_fan_tail``, so the exhaustion θ is well-defined. The crossing test is
808 orientation-agnostic: it holds for both the shrinking decay
809 (``c_decay_initial > c_fan_tail``) and the growing decay
810 (``c_decay_initial < c_fan_tail``). Returns ``None`` when ``c_fan_tail``
811 is not strictly between ``c_fixed`` and ``c_decay_initial`` — e.g. full drying
812 (``c_fan_tail == c_fixed``), where the decay asymptotically merges with
813 the fixed state and no finite exhaustion event occurs.
815 Returns
816 -------
817 float or None
818 Cumulative flow θ at exhaustion, or ``None`` if not reached.
819 """
820 # An interior exhaustion needs c_fan_tail strictly between c_fixed and
821 # c_decay_initial (orientation-agnostic via min/max). Full drying
822 # (c_fan_tail == c_fixed) merges asymptotically with no finite crossing —
823 # return None rather than grow the bracket forever (van Genuchten would hang).
824 c_lo = min(self.c_fixed, self.c_decay_initial)
825 c_hi = max(self.c_fixed, self.c_decay_initial)
826 if not (c_lo < self.c_fan_tail < c_hi):
827 return None
829 theta_local_collision = self.theta_start - self.theta_origin
831 if self._numerical:
832 # The numerical forward map saturates AT c_fan_tail (it never crosses it),
833 # so a forward-map bracket cannot see the crossing for either decay
834 # orientation. Evaluate θ_local(c_fan_tail) from the un-clamped invariant
835 # directly: the gate above guarantees c_fan_tail is the reachable limit, so
836 # it is exactly the cached profile's endpoint ``i_max``.
837 _c_of_i, i_max, _c_limit_node = self._decay_profile()
838 return self.theta_origin + theta_local_collision * float(np.exp(i_max))
840 # Closed forms cross c_fan_tail smoothly (always a shrinking decay); invert the
841 # monotone forward map by bracketing (orientation-agnostic — no early return).
842 def f(theta_local: float) -> float:
843 return self._c_decay_at_theta_local(theta_local) - self.c_fan_tail
845 theta_local_exhaust = _invert_monotone_theta_local(
846 f, theta_hi_seed=theta_local_collision, f_seed=f(theta_local_collision)
847 )
848 if theta_local_exhaust is None:
849 return None
850 return self.theta_origin + theta_local_exhaust
852 def _c_decay_at_theta_local(self, theta_local: float) -> float:
853 """Decaying concentration as a function of ``θ_local`` (apex-relative).
855 The SOLE isotherm dispatch site: closed where an exact form exists,
856 otherwise the cached numerical decay profile. The closed forms
857 (Freundlich ``c_fixed=0`` or ``n≈2``; Langmuir ``c_fixed=0``;
858 Brooks-Corey ``c_fixed=0``) are selected by the ``__post_init__``
859 predicates; every other ``(isotherm, c_fixed)`` combination falls through
860 to the per-wave cached invariant profile. Takes ``θ_local`` directly and
861 skips the activity check. ``c_decay_at_theta``, ``position_at_theta``,
862 ``theta_at_fan_exhaustion`` and ``outlet_crossing_theta`` all route
863 through here rather than repeating the dispatch.
864 """
865 theta_local_collision = self.theta_start - self.theta_origin
866 s = self.sorption
867 if self._freundlich_cf and isinstance(s, FreundlichSorption):
868 return _c_decay_freundlich(
869 s, self.K, self.c_decay_initial, self.c_fixed, theta_local_collision, theta_local
870 )
871 if self._langmuir_cf and isinstance(s, LangmuirSorption):
872 return _c_decay_langmuir(s, self.K, theta_local)
873 if self._brooks_corey_cf and isinstance(s, BrooksCoreyConductivity):
874 return _c_decay_brooks_corey(s, self.c_decay_initial, theta_local_collision, theta_local)
876 # Numerical path: invert the cached monotone θ_local(c) map. c_decay stays at
877 # c_decay_initial up to the collision and clamps at the reachable c-limit past it.
878 if theta_local <= theta_local_collision:
879 return self.c_decay_initial
880 c_of_i, i_max, c_limit_node = self._decay_profile()
881 i_target = np.log(theta_local / theta_local_collision)
882 if i_target >= i_max:
883 return float(c_limit_node)
884 return float(c_of_i(i_target))
886 def outlet_crossing_theta(self, v_outlet: float) -> float | None:
887 """Cumulative flow at which ``V_s = v_outlet``.
889 Returns ``None`` if the outlet is upstream of the wave's birth
890 position or no crossing exists in ``(theta_start, +∞)``. The wave's
891 current activity flag is not consulted — callers asking
892 retrospectively about a historical crossing need the answer regardless
893 of subsequent deactivation.
895 The closed-form Freundlich/Langmuir cases invert the fan-continuity
896 identity ``V_s − v_origin = θ_local / R(c_decay)`` analytically (valid
897 only when ``_c_decay_at_theta_local`` itself uses the closed form, so
898 the same conditions are mirrored here); every other case inverts the
899 monotone ``V_s(θ)`` via ``brentq``.
900 """
901 if v_outlet <= self.v_start:
902 return None
904 # V_s is monotonically increasing in θ (positive shock speed); invert
905 # via the fan-continuity identity V_s - v_origin = θ_local / R(c_decay)
906 # combined with the invariant to eliminate u, then solve for θ.
907 s = self.sorption
908 if self._freundlich_cf and isinstance(s, FreundlichSorption):
909 return _outlet_crossing_freundlich(
910 s,
911 self.K,
912 self.c_decay_initial,
913 self.c_fixed,
914 self.v_origin,
915 self.theta_origin,
916 v_outlet,
917 )
918 if self._langmuir_cf and isinstance(s, LangmuirSorption):
919 return _outlet_crossing_langmuir(
920 s,
921 self.K,
922 self.v_origin,
923 self.theta_origin,
924 v_outlet,
925 )
926 return self._outlet_crossing_numerical(v_outlet)
928 def _outlet_crossing_numerical(self, v_outlet: float) -> float | None:
929 """θ at which ``V_s = v_outlet`` for every non-closed-form case.
931 ``V_s(θ) = v_origin + θ_local / R(c_decay(θ))`` is monotone increasing;
932 invert by ``brentq`` on ``θ_local``.
933 """
934 theta_local_collision = self.theta_start - self.theta_origin
936 def f(theta_local: float) -> float:
937 c = self._c_decay_at_theta_local(theta_local)
938 return self.v_origin + theta_local / float(self.sorption.retardation(c)) - v_outlet
940 f_lo = f(theta_local_collision)
941 if f_lo >= 0.0:
942 # Already at/past the outlet at collision; the linear-shock guards
943 # in the solver handle the duplicate-crossing suppression.
944 return self.theta_start
945 # Seed at the collision (f_lo < 0 established above) and let the helper grow the
946 # bracket upward — no dimensional floor, so crossings within θ_local < 1 of the
947 # apex are found (mirrors theta_at_fan_exhaustion's closed-form bracket).
948 theta_local_cross = _invert_monotone_theta_local(f, theta_hi_seed=theta_local_collision, f_seed=f_lo)
949 if theta_local_cross is None:
950 return None
951 return self.theta_origin + theta_local_cross
953 def concentration_left(self) -> float:
954 """Concentration on the left (upstream) side at θ=theta_start.
956 For ``decay_side='left'`` returns the decaying c at the collision
957 moment; for ``decay_side='right'`` returns the fixed side.
958 """
959 return self.c_decay_initial if self.decay_side == "left" else self.c_fixed
961 def concentration_right(self) -> float:
962 """Concentration on the right (downstream) side at θ=theta_start.
964 For ``decay_side='right'`` returns the decaying c at the collision
965 moment; for ``decay_side='left'`` returns the fixed side.
966 """
967 return self.c_decay_initial if self.decay_side == "right" else self.c_fixed
969 def concentration_at_point(self, v: float, theta: float) -> float | None:
970 """Concentration at ``(v, θ)`` if controlled by this decaying shock.
972 Three regions:
974 1. ``v == V_s(θ)`` (within FP): average of decay-side and fixed-side c.
975 2. ``v > V_s(θ)`` (downstream): fixed-side c if ``decay_side='left'``;
976 decay-side c at θ if ``decay_side='right'``.
977 3. ``v < V_s(θ)`` (upstream, inside the fan): the fan's self-similar
978 concentration ``R(c) = (θ − theta_origin)/(v − v_origin)``. Outside
979 the fan — i.e. the decay-side characteristic from the apex hasn't
980 reached v yet, OR the point lies beyond the ``c_fan_tail`` boundary
981 (the fan's far edge) — returns ``None``.
983 Returns ``None`` for ``θ < theta_start`` or inactive waves.
984 """
985 if not self.was_active_at(theta):
986 return None
988 # Compute the decaying-side concentration once and derive V_s from it
989 # (inlining position_at_theta's body) so the shock-face branch can reuse
990 # it instead of re-running the numerical-isotherm root-find.
991 theta_local = theta - self.theta_origin
992 c_d = self._c_decay_at_theta_local(theta_local)
993 v_s = float(self.v_origin + theta_local / float(self.sorption.retardation(c_d)))
995 tol = 1e-15 * max(abs(v_s), 1.0)
997 if abs(v - v_s) < tol:
998 return 0.5 * (c_d + self.c_fixed)
1000 # Region selection depends on decay_side:
1001 # 'left' (favorable n>1, Langmuir): fan extends upstream of V_s
1002 # (v < V_s), c_fixed downstream (v > V_s).
1003 # 'right' (n<1 mirror): fan extends downstream of V_s (v > V_s),
1004 # c_fixed upstream (v < V_s).
1005 if self.decay_side == "left":
1006 v_fan_side = v < v_s - tol
1007 v_fixed_side = v > v_s + tol
1008 else:
1009 v_fan_side = v > v_s + tol
1010 v_fixed_side = v < v_s - tol
1012 if v_fixed_side:
1013 return self.c_fixed
1015 if not v_fan_side:
1016 return None # within tol of shock face — handled above
1018 # Fan-interior: self-similar profile with apex at (v_origin, theta_origin).
1019 if v == self.v_origin:
1020 return None
1021 r_target = (theta - self.theta_origin) / (v - self.v_origin)
1022 if r_target <= 1.0:
1023 return None
1024 try:
1025 c_fan = self.sorption.concentration_from_retardation(r_target)
1026 except NotImplementedError:
1027 return None
1028 c_fan = float(c_fan)
1030 # The fan the DSW controls spans concentrations between the shock face
1031 # (c_decay, ≤ c_decay_initial) and the fan's far boundary c_fan_tail.
1032 # A point past c_fan_tail belongs to whatever lies beyond the fan, not
1033 # to this wave — reject so the fan is not extended past its extent.
1034 c_lo = min(self.c_fan_tail, self.c_decay_initial)
1035 c_hi = max(self.c_fan_tail, self.c_decay_initial)
1036 if c_fan < c_lo - EPSILON_POSITION or c_fan > c_hi + EPSILON_POSITION:
1037 return None
1038 return c_fan
1041@dataclass
1042class DoubleFanShockWave(Wave):
1043 r"""Shock fed by a self-similar fan on BOTH sides (a doubly-fed front).
1045 Formed when a fan boundary (rarefaction head/tail, or another fan-fed shock) catches a
1046 :class:`DecayingShockWave` whose surviving side is itself a fan — the merged front then
1047 has a fan feeder on each side. The shock trajectory solves
1049 .. math::
1050 \frac{dV}{d\theta} = S\bigl(c_L(V,\theta),\, c_R(V,\theta)\bigr), \qquad
1051 R(c_i) = \frac{\theta - \theta_{{\rm apex},i}}{V - v_{{\rm apex},i}},
1053 with ``S`` the Rankine-Hugoniot secant and each ``c_i`` the self-similar value of its fan.
1055 **Closed form (Freundlich ``n = 2``, shared apex position ``v_L = v_R = v_o``).** With
1056 ``u_i = \sqrt{c_i} = A V'/(2(\tau_i - V'))`` (``V' = V - v_o``, ``\tau_i = \theta -
1057 \theta_{{\rm apex},i}``, ``A = \rho_b k_f/n_{por}``), the product ``K = u_L u_R`` is a
1058 first integral of the ODE, and the trajectory is the physical root of
1060 .. math::
1061 (A^2 - 4K)\,V'^2 + 4K(\tau_L + \tau_R)\,V' - 4K\,\tau_L\,\tau_R = 0,
1062 \qquad 0 < V' < \min(\tau_L, \tau_R).
1064 Every fan born at the inlet shares ``v_o = 0``, so inlet-driven inputs use the closed
1065 form. **General fallback** (distinct apex positions, or a non-``n=2`` isotherm): the ODE
1066 is integrated once by fixed-step RK4 (``\Delta\theta = \text{fan age}/DFSW_RK_SUBSTEPS``,
1067 speed-independent since the fans vary on the scale of their age) into a cached monotone
1068 spline. Both paths answer position, side concentrations and outlet crossing uniformly
1069 (side exhaustion is detected externally, as the shock face crossing its own fan boundary).
1071 See Also
1072 --------
1073 DecayingShockWave : One fan side; the doubly-fed front degrades to this on side exhaustion.
1074 Feeder : The bounded-fan side descriptor this wave carries.
1075 """
1077 left_feeder: Feeder
1078 """Left (upstream) side feeder — a fan."""
1079 right_feeder: Feeder
1080 """Right (downstream) side feeder — a fan."""
1081 sorption: NonlinearSorption
1082 """Sorption model (concentration-dependent)."""
1083 left_boundary_consumed: bool = field(default=False, kw_only=True)
1084 """Whether the left fan's far-boundary line is owned from birth (not a free collision face)."""
1085 right_boundary_consumed: bool = field(default=False, kw_only=True)
1086 """Whether the right fan's far-boundary line is owned from birth (not a free collision face)."""
1087 theta_left_boundary_consumed: float = field(default=float("inf"), kw_only=True)
1088 """Cumulative flow at which a free left boundary was later consumed (historical, see DSW)."""
1089 theta_right_boundary_consumed: float = field(default=float("inf"), kw_only=True)
1090 """Cumulative flow at which a free right boundary was later consumed (historical, see DSW)."""
1091 _closed_form: bool = field(init=False, repr=False, compare=False)
1092 """Whether the n=2 shared-apex closed form applies."""
1093 _k: float = field(init=False, repr=False, compare=False)
1094 """Closed-form first integral ``K = u_L·u_R`` (``nan`` for the numerical path)."""
1095 _traj_cache: tuple | None = field(default=None, init=False, repr=False, compare=False)
1096 """Lazily-built ``(theta_grid, CubicSpline)`` for the numerical trajectory."""
1098 def __post_init__(self) -> None:
1099 """Classify closed-form vs numerical and set the first integral ``K``."""
1100 if self.left_feeder.is_const or self.right_feeder.is_const:
1101 msg = "DoubleFanShockWave requires two fan feeders; use DecayingShockWave for a one-fan front"
1102 raise ValueError(msg)
1103 if self.theta_origin_left >= self.theta_start or self.theta_origin_right >= self.theta_start:
1104 msg = "both fan apexes must precede the collision (theta_apex < theta_start)"
1105 raise ValueError(msg)
1106 s = self.sorption
1107 shared_apex = abs(self.left_feeder.v_apex - self.right_feeder.v_apex) < EPSILON_POSITION
1108 self._closed_form = isinstance(s, FreundlichSorption) and bool(np.isclose(s.n, 2.0, rtol=1e-12)) and shared_apex
1109 self._k = float("nan")
1110 if self._closed_form:
1111 v_prime = self.v_start - self.left_feeder.v_apex
1112 tau_l = self.theta_start - self.theta_origin_left
1113 tau_r = self.theta_start - self.theta_origin_right
1114 a = self._alpha()
1115 u_l = a * v_prime / (2.0 * (tau_l - v_prime))
1116 u_r = a * v_prime / (2.0 * (tau_r - v_prime))
1117 self._k = u_l * u_r
1119 @property
1120 def theta_origin_left(self) -> float:
1121 """Cumulative flow at the left fan's apex [m³]."""
1122 return self.left_feeder.theta_apex
1124 @property
1125 def theta_origin_right(self) -> float:
1126 """Cumulative flow at the right fan's apex [m³]."""
1127 return self.right_feeder.theta_apex
1129 def _alpha(self) -> float:
1130 """Freundlich lumped coefficient ``A = ρ_b·k_f/n_por`` (closed-form path only)."""
1131 s = self.sorption
1132 assert isinstance(s, FreundlichSorption) # noqa: S101
1133 return s.bulk_density * s.k_f / s.porosity
1135 def _v_closed(self, theta: float) -> float:
1136 """Closed-form shock position at ``theta`` (n=2 shared apex)."""
1137 v_o = self.left_feeder.v_apex
1138 tau_l = theta - self.theta_origin_left
1139 tau_r = theta - self.theta_origin_right
1140 a2 = self._alpha() ** 2 - 4.0 * self._k
1141 a1 = 4.0 * self._k * (tau_l + tau_r)
1142 a0 = -4.0 * self._k * tau_l * tau_r
1143 if abs(a2) < EPSILON_POSITION:
1144 # A² = 4K: the quadratic degenerates to linear a1·V' + a0 = 0.
1145 v_prime = -a0 / a1
1146 return v_o + v_prime
1147 disc = a1 * a1 - 4.0 * a2 * a0
1148 sqrt_disc = np.sqrt(max(disc, 0.0))
1149 upper = min(tau_l, tau_r)
1150 for sign in (-1.0, 1.0):
1151 v_prime = (-a1 + sign * sqrt_disc) / (2.0 * a2)
1152 if 0.0 < v_prime < upper:
1153 return v_o + v_prime
1154 # Numerical fallback: pick the root closest to the valid open interval.
1155 candidates = [(-a1 + sign * sqrt_disc) / (2.0 * a2) for sign in (-1.0, 1.0)]
1156 v_prime = min(candidates, key=lambda x: abs(x - 0.5 * upper))
1157 return v_o + float(np.clip(v_prime, 0.0, upper))
1159 def _rhs(self, v: float, theta: float) -> float:
1160 """Shock-speed ODE right side ``S(c_L, c_R)`` at ``(v, θ)`` (numerical path)."""
1161 c_l = self.left_feeder.value(v, theta)
1162 c_r = self.right_feeder.value(v, theta)
1163 return float(self.sorption.shock_speed(c_l, c_r))
1165 def _ensure_numerical(self, theta: float) -> tuple:
1166 """Build/extend the cached RK4 trajectory spline out to at least ``theta``."""
1167 age0 = max(min(self.theta_start - self.theta_origin_left, self.theta_start - self.theta_origin_right), 1.0)
1168 step0 = age0 / DFSW_RK_SUBSTEPS
1169 if self._traj_cache is None:
1170 thetas = [self.theta_start]
1171 vs = [self.v_start]
1172 else:
1173 thetas, vs, _ = self._traj_cache
1174 thetas = list(thetas)
1175 vs = list(vs)
1176 # March past ``theta`` and always keep at least two nodes so CubicSpline is well-posed
1177 # (a query at ``theta_start`` alone would otherwise leave a single node).
1178 target = max(theta, self.theta_start + step0)
1179 while thetas[-1] < target:
1180 t0 = thetas[-1]
1181 v0 = vs[-1]
1182 age = min(t0 - self.theta_origin_left, t0 - self.theta_origin_right, age0 + (t0 - self.theta_start))
1183 h = max(age, age0) / DFSW_RK_SUBSTEPS
1184 k1 = self._rhs(v0, t0)
1185 k2 = self._rhs(v0 + 0.5 * h * k1, t0 + 0.5 * h)
1186 k3 = self._rhs(v0 + 0.5 * h * k2, t0 + 0.5 * h)
1187 k4 = self._rhs(v0 + h * k3, t0 + h)
1188 vs.append(v0 + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4))
1189 thetas.append(t0 + h)
1190 theta_arr = np.asarray(thetas)
1191 v_arr = np.asarray(vs)
1192 spline = CubicSpline(theta_arr, v_arr)
1193 self._traj_cache = (theta_arr, v_arr, spline)
1194 return self._traj_cache
1196 def _v_at(self, theta: float) -> float:
1197 """Shock position at ``theta`` (closed form or cached numerical spline)."""
1198 if self._closed_form:
1199 return self._v_closed(theta)
1200 _theta_arr, _v_arr, spline = self._ensure_numerical(theta)
1201 return float(spline(theta))
1203 def position_at_theta(self, theta: float) -> float | None:
1204 """Shock position ``V_s(θ)``; ``None`` for ``θ < theta_start`` or when inactive."""
1205 if not self.was_active_at(theta):
1206 return None
1207 return self._v_at(theta)
1209 def concentration_left(self) -> float:
1210 """Left-side concentration at the collision moment."""
1211 return self.left_feeder.value(self.v_start, self.theta_start)
1213 def concentration_right(self) -> float:
1214 """Right-side concentration at the collision moment."""
1215 return self.right_feeder.value(self.v_start, self.theta_start)
1217 def outlet_crossing_theta(self, v_outlet: float) -> float | None:
1218 """Cumulative flow at which ``V_s = v_outlet`` (monotone, inverted by brentq).
1220 ``θ_local`` is measured from the older of the two fan apexes, so the shared
1221 bracket-then-brentq inversion sees a positive seed.
1222 """
1223 if v_outlet <= self.v_start:
1224 return None
1225 theta_origin = min(self.theta_origin_left, self.theta_origin_right)
1226 theta_local_cross = _invert_monotone_theta_local(
1227 lambda theta_local: self._v_at(theta_origin + theta_local) - v_outlet,
1228 theta_hi_seed=self.theta_start - theta_origin,
1229 f_seed=self.v_start - v_outlet,
1230 )
1231 if theta_local_cross is None:
1232 return None
1233 return theta_origin + theta_local_cross
1235 def concentration_at_point(self, v: float, theta: float) -> float | None:
1236 """Concentration at ``(v, θ)`` if controlled by this doubly-fed shock.
1238 Left of the shock face the left fan controls; right of it the right fan; at the face
1239 the average. Outside both fans' physical extents returns ``None`` (another wave owns
1240 the point).
1241 """
1242 if not self.was_active_at(theta):
1243 return None
1244 v_s = self._v_at(theta)
1245 tol = 1e-15 * max(abs(v_s), 1.0)
1246 if abs(v - v_s) < tol:
1247 return 0.5 * (self.left_feeder.value(v_s, theta) + self.right_feeder.value(v_s, theta))
1248 feeder = self.left_feeder if v < v_s else self.right_feeder
1249 # Only claim the point if it lies within the fan's live self-similar range.
1250 if theta <= feeder.theta_apex or v <= feeder.v_apex:
1251 return None
1252 r = (theta - feeder.theta_apex) / (v - feeder.v_apex)
1253 r_a = float(self.sorption.retardation(feeder.c_a))
1254 r_b = float(self.sorption.retardation(feeder.c_b))
1255 if r < min(r_a, r_b) - EPSILON_POSITION or r > max(r_a, r_b) + EPSILON_POSITION:
1256 return None
1257 return feeder.value(v, theta)
1260def _invert_monotone_theta_local(f, *, theta_hi_seed: float, f_seed: float | None = None) -> float | None:
1261 """Bracket-then-brentq a monotone ``f(θ_local)`` with a sign change above the seed.
1263 Shared by the closed-form branch of ``theta_at_fan_exhaustion``, by
1264 ``DecayingShockWave._outlet_crossing_numerical`` and by
1265 ``DoubleFanShockWave.outlet_crossing_theta``: each inverts a monotone function of
1266 ``θ_local`` whose sign at the collision is already known to differ from its
1267 sign at large ``θ_local``. Geometrically grows ``θ_hi`` (``×2``, ≤200 iters)
1268 from ``theta_hi_seed`` until ``f`` flips sign, then inverts with ``brentq``.
1269 Returns ``None`` if no sign change is bracketed within the iteration budget.
1271 Parameters
1272 ----------
1273 f : callable
1274 Monotone residual ``f(θ_local)``; ``f(theta_hi_seed)`` and the far-field
1275 value must straddle zero. Caller-specific early sentinels
1276 (already-past-outlet) are handled by the caller.
1277 theta_hi_seed : float
1278 ``θ_local`` lower bracket; the search grows ``θ_hi`` from here.
1279 f_seed : float, optional
1280 Pre-evaluated ``f(theta_hi_seed)``; callers that already computed it pass
1281 it to avoid a redundant evaluation. ``None`` recomputes it here.
1283 Returns
1284 -------
1285 float or None
1286 Root ``θ_local`` of ``f``, or ``None`` if not bracketed.
1287 """
1288 if f_seed is None:
1289 f_seed = f(theta_hi_seed)
1290 theta_hi = theta_hi_seed
1291 for _ in range(200):
1292 theta_hi *= 2.0
1293 if f(theta_hi) * f_seed < 0.0:
1294 return float(brentq(f, theta_hi_seed, theta_hi, xtol=DECAYING_SHOCK_BRENTQ_XTOL))
1295 return None
1298def _build_decay_profile(
1299 sorption: NonlinearSorption,
1300 c_decay_initial: float,
1301 c_fixed: float,
1302 c_fan_tail: float,
1303) -> tuple:
1304 r"""Build the per-wave monotone ``θ_local(c)`` map for collisions with no closed form.
1306 Decay-agnostic: the fan-continuity + Rankine-Hugoniot relations do not
1307 depend on which side decays. The secant speed
1308 ``S(c) = (c − c_fixed)/(C_T(c) − C_T(c_fixed))`` is symmetric in
1309 ``(c_decay, c_fixed)``, so the same invariant
1310 ``θ_local(c) = θ_local_coll · exp(I(c))``, ``I(c) = ∫_{c0}^{c} R'/[(1 − R·S)·R] dc``
1311 (``R'`` by central finite difference) holds for Freundlich ``c_fixed>0, n≠2``,
1312 Langmuir ``c_fixed>0``, Brooks-Corey ``c_fixed>0`` and any van Genuchten case
1313 alike. ``I(c)`` is built ONCE by a single vectorised composite Gauss-Legendre
1314 cumulative quadrature over a c-grid from ``c_decay_initial`` to the reachable
1315 limit, then inverted by monotone-spline interpolation. The reachable limit is the
1316 fan tail ``c_fan_tail`` UNLESS ``c_fixed`` lies strictly between
1317 ``c_decay_initial`` and ``c_fan_tail`` — then
1318 the secant speed has a pole at ``c_fixed`` (``R·S → 1``, ``θ_local → ∞``): the
1319 shock asymptotes to the fixed state, so the grid stops a hair short of it and
1320 ``c_decay`` clamps there.
1322 Parameters
1323 ----------
1324 sorption : NonlinearSorption
1325 Sorption model.
1326 c_decay_initial : float
1327 Decaying-side concentration at the collision (``c0``) [mass/volume].
1328 c_fixed : float
1329 Non-decaying-side concentration [mass/volume].
1330 c_fan_tail : float
1331 Concentration at the fan's far boundary [mass/volume]; bounds the decay.
1333 Returns
1334 -------
1335 tuple
1336 ``(c_of_i, i_max, c_limit_node)``: a ``CubicSpline`` mapping the cumulative
1337 invariant ``I = ln(θ_local/θ_local_coll)`` to ``c_decay``, the endpoint
1338 ``I`` of the reachable c-range, and the ``c`` at that endpoint. ``I`` is
1339 collision-independent (``θ_local_coll`` enters only at query time).
1340 """
1341 ct_fixed = float(sorption.total_concentration(c_fixed))
1343 def integrand(c):
1344 c = np.asarray(c, dtype=float)
1345 h = np.maximum(1e-9, 1e-7 * np.abs(c))
1346 r_prime = (np.asarray(sorption.retardation(c + h)) - np.asarray(sorption.retardation(c - h))) / (2.0 * h)
1347 r = np.asarray(sorption.retardation(c))
1348 ct = np.asarray(sorption.total_concentration(c))
1349 secant = (c - c_fixed) / (ct - ct_fixed)
1350 return r_prime / ((1.0 - r * secant) * r)
1352 pole = (c_decay_initial - c_fixed) * (c_fan_tail - c_fixed) < 0.0
1353 c_limit = c_fixed if pole else c_fan_tail
1355 # c-grid from c0 to the reachable limit. Toward a pole (θ_local → ∞) the grid is
1356 # geometric, stopping a fraction DECAY_PROFILE_POLE_FLOOR of the gap short; the
1357 # non-pole grid reaches c_fan_tail exactly (so i_max is the exhaustion integral).
1358 frac = np.linspace(0.0, 1.0, DECAY_PROFILE_NODES)
1359 gap0 = abs(c_decay_initial - c_limit)
1360 gaps = gap0 * DECAY_PROFILE_POLE_FLOOR**frac if pole else gap0 * (1.0 - frac)
1361 c_nodes = c_limit + np.sign(c_decay_initial - c_limit) * gaps
1363 # Cumulative composite Gauss-Legendre integral of the invariant integrand.
1364 x_gl, w_gl = np.polynomial.legendre.leggauss(DECAY_PROFILE_GAUSS_ORDER)
1365 lo = c_nodes[:-1]
1366 hi = c_nodes[1:]
1367 mid = 0.5 * (lo + hi)
1368 half = 0.5 * (hi - lo)
1369 points = mid[:, None] + half[:, None] * x_gl[None, :]
1370 panel = (integrand(points.ravel()).reshape(points.shape) * w_gl[None, :]).sum(axis=1) * half
1371 i_nodes = np.concatenate([[0.0], np.cumsum(panel)])
1373 # Inverse-interpolation precondition: keep the strictly-increasing prefix (the
1374 # near-pole tail can lose monotonicity as the singular integrand outruns the grid).
1375 non_increasing = np.nonzero(np.diff(i_nodes) <= 0.0)[0]
1376 if non_increasing.size:
1377 cut = non_increasing[0] + 1
1378 i_nodes = i_nodes[:cut]
1379 c_nodes = c_nodes[:cut]
1380 c_of_i = CubicSpline(i_nodes, c_nodes)
1381 return c_of_i, float(i_nodes[-1]), float(c_nodes[-1])
1384def _c_decay_brooks_corey(
1385 sorption: BrooksCoreyConductivity,
1386 c_decay_initial: float,
1387 theta_local_collision: float,
1388 theta_local: float,
1389) -> float:
1390 r"""Brooks-Corey ``c_fixed = 0`` closed form for the decaying-side concentration.
1392 For Brooks-Corey with ``c_fixed = 0`` the product ``R·S = 1/a`` is constant
1393 (``a = sorption.a``), so the universal invariant integrates to
1394 ``θ_local ∝ R(c_decay)^{a/(a−1)}``. Inverting,
1395 ``R(c_d) = R(c0)·(θ_local/θ_local_coll)^{(a−1)/a}`` and
1396 ``c_d = concentration_from_retardation(R)``.
1398 Parameters
1399 ----------
1400 sorption : BrooksCoreyConductivity
1401 Sorption model.
1402 c_decay_initial : float
1403 Decaying-side concentration at the collision [mass/volume].
1404 theta_local_collision : float
1405 ``θ_local`` at the collision [m³].
1406 theta_local : float
1407 ``θ_local`` at which to evaluate the decaying concentration [m³].
1409 Returns
1410 -------
1411 float
1412 Decaying-side concentration ``c`` at ``theta_local``.
1413 """
1414 if theta_local <= theta_local_collision:
1415 return c_decay_initial
1416 a = sorption.a
1417 r0 = float(sorption.retardation(c_decay_initial))
1418 r_target = r0 * (theta_local / theta_local_collision) ** ((a - 1.0) / a)
1419 return float(sorption.concentration_from_retardation(r_target))
1422def _compute_k_freundlich(
1423 sorption: FreundlichSorption,
1424 theta_local: float,
1425 c_decay_initial: float,
1426 c_fixed: float,
1427) -> float:
1428 """Closed-form invariant K for Freundlich DecayingShockWave.
1430 Derivation: see plan §"Closed-form derivations". For c_fixed=0,
1431 K = θ_local · u_c^n / (n · u_c^(n-1) + α); for c_fixed>0 (n=2 only),
1432 K = θ_local · (u_c - u_r)^2 / (2 · u_c + α). Here α = ρ_b · k_f / n_por.
1434 Parameters
1435 ----------
1436 sorption : FreundlichSorption
1437 Sorption model.
1438 theta_local : float
1439 Cumulative flow from rarefaction apex to collision [m³].
1440 c_decay_initial : float
1441 Decaying-side concentration at the collision [mass/volume].
1442 c_fixed : float
1443 Non-decaying-side concentration [mass/volume].
1445 Returns
1446 -------
1447 float
1448 Invariant constant K.
1449 """
1450 n = sorption.n
1451 alpha = sorption.bulk_density * sorption.k_f / sorption.porosity
1452 u_d = c_decay_initial ** (1.0 / n)
1454 if c_fixed == 0.0:
1455 return float(theta_local * u_d**n / (n * u_d ** (n - 1.0) + alpha))
1457 # c_fixed > 0, n=2
1458 u_r = c_fixed**0.5
1459 return float(theta_local * (u_d - u_r) ** 2 / (2.0 * u_d + alpha))
1462def _compute_k_langmuir(
1463 sorption: LangmuirSorption,
1464 theta_local: float,
1465 c_decay_initial: float,
1466) -> float:
1467 """Closed-form invariant K for Langmuir DecayingShockWave (c_fixed=0).
1469 K = θ_local · c_d^2 / ((K_L + c_d)^2 + a) with a = ρ_b · s_max · K_L / n_por.
1471 Parameters
1472 ----------
1473 sorption : LangmuirSorption
1474 Sorption model.
1475 theta_local : float
1476 Cumulative flow from rarefaction apex to collision [m³].
1477 c_decay_initial : float
1478 Decaying-side concentration at the collision [mass/volume].
1480 Returns
1481 -------
1482 float
1483 Invariant constant K.
1484 """
1485 return float(theta_local * c_decay_initial**2 / ((sorption.k_l + c_decay_initial) ** 2 + sorption.a_coeff))
1488def _c_decay_freundlich(
1489 sorption: FreundlichSorption,
1490 k_invariant: float,
1491 c_decay_initial: float,
1492 c_fixed: float,
1493 theta_local_collision: float,
1494 theta_local: float,
1495) -> float:
1496 """Invert the Freundlich invariant to get c on the decaying side at θ_local.
1498 For n=2 c_fixed=0 (quadratic in u): closed form
1499 ``u = (K + sqrt(K^2 + K·θ_local·α)) / θ_local``. For general n with
1500 c_fixed=0 (transcendental): brentq on the monotone bracket
1501 ``(tiny, c_decay_initial^(1/n)]``. For n=2 c_fixed>0 (quadratic in u with u_r):
1502 closed form; the root is selected by the decay orientation — the ``+√`` root
1503 (``u > u_r``) for a shrinking decay (``c_decay_initial > c_fixed``), the ``−√``
1504 root (``u < u_r``) for a growing decay (``c_decay_initial < c_fixed``). Both
1505 approach ``u_r`` as ``θ_local → ∞``; the initial side of ``u_r`` fixes the branch.
1507 Returns
1508 -------
1509 float
1510 Decaying-side concentration c at θ_local.
1511 """
1512 n = sorption.n
1513 alpha = sorption.bulk_density * sorption.k_f / sorption.porosity
1515 if c_fixed == 0.0:
1516 # atol=0.0: the default atol=1e-8 would silently route any n within 1e-8 of 2
1517 # to the n=2 closed form instead of the general-n inversion.
1518 if np.isclose(n, 2.0, rtol=1e-12, atol=0.0):
1519 disc = k_invariant * k_invariant + theta_local * k_invariant * alpha
1520 u = (k_invariant + np.sqrt(disc)) / theta_local
1521 return float(u * u)
1522 u_root = _invert_freundlich_cr_zero(k_invariant, c_decay_initial, n, alpha, theta_local_collision, theta_local)
1523 return float(u_root**n)
1525 # n=2, c_fixed > 0. Growing decay (c_decay_initial < c_fixed) starts below u_r and
1526 # takes the −√ root; shrinking decay starts above u_r and takes the +√ root.
1527 u_r = c_fixed**0.5
1528 disc = k_invariant * (theta_local * (2.0 * u_r + alpha) + k_invariant)
1529 sqrt_disc = np.sqrt(disc)
1530 if c_decay_initial < c_fixed:
1531 u = (u_r * theta_local + k_invariant - sqrt_disc) / theta_local
1532 else:
1533 u = (u_r * theta_local + k_invariant + sqrt_disc) / theta_local
1534 return float(u * u)
1537def _invert_freundlich_cr_zero(
1538 k_invariant: float,
1539 c_decay_initial: float,
1540 n: float,
1541 alpha: float,
1542 theta_local_collision: float,
1543 theta_local: float,
1544) -> float:
1545 """Invert ``θ_local · u^n = K · (n·u^(n-1) + α)`` for u via brentq.
1547 Returns
1548 -------
1549 float
1550 Root u of the invariant at θ_local.
1551 """
1552 u_collision = c_decay_initial ** (1.0 / n)
1554 def f(u: float) -> float:
1555 return theta_local * u**n - k_invariant * (n * u ** (n - 1.0) + alpha)
1557 if theta_local <= theta_local_collision:
1558 # Earlier than (or at) collision: c_decay equals c_decay_initial.
1559 return u_collision
1561 u_root = brentq(f, DECAYING_SHOCK_U_FLOOR, u_collision, xtol=1e-15)
1562 return float(u_root) # type: ignore[arg-type]
1565def _c_decay_langmuir(sorption: LangmuirSorption, k_invariant: float, theta_local: float) -> float:
1566 """Invert the Langmuir invariant ``θ_local · c^2 = K · ((K_L+c)^2 + a)`` for c.
1568 Expanded: ``(θ_local - K)·c^2 - 2·K·K_L·c - K·(K_L^2 + a) = 0`` (quadratic
1569 in c). Positive root chosen for c > 0.
1571 Returns
1572 -------
1573 float
1574 Decaying-side concentration c at θ_local.
1575 """
1576 k_l = sorption.k_l
1577 a_coeff = sorption.a_coeff
1578 denom = theta_local - k_invariant
1579 disc = k_invariant * (k_invariant * k_l * k_l + denom * (k_l * k_l + a_coeff))
1580 return float((k_invariant * k_l + np.sqrt(disc)) / denom)
1583def _outlet_crossing_freundlich(
1584 sorption: FreundlichSorption,
1585 k_invariant: float,
1586 c_decay_initial: float,
1587 c_fixed: float,
1588 v_origin: float,
1589 theta_origin: float,
1590 v_outlet: float,
1591) -> float | None:
1592 """θ at which a Freundlich DecayingShockWave reaches v_outlet.
1594 Returns
1595 -------
1596 float or None
1597 Cumulative flow at crossing, or None if no crossing.
1598 """
1599 n = sorption.n
1600 alpha = sorption.bulk_density * sorption.k_f / sorption.porosity
1601 delta_v = v_outlet - v_origin
1603 if c_fixed == 0.0:
1604 u_target = n * k_invariant / delta_v
1605 if u_target <= 0.0:
1606 return None
1607 theta_local = k_invariant * (n * u_target ** (n - 1.0) + alpha) / u_target**n
1608 return float(theta_origin + theta_local)
1610 # n=2, c_fixed > 0: V_s - v_origin = 2·K·u / (u - u_r)^2 ⇒ quadratic in u.
1611 # The two roots multiply to u_r² and sum to 2u_r + 2K/delta_v > 2u_r, so exactly one
1612 # exceeds u_r (the +√ root, shrinking decay) and one lies below (the −√ root, growing
1613 # decay). Select by the decay orientation to match _c_decay_freundlich's branch.
1614 u_r = c_fixed**0.5
1615 b_coef = -(2.0 * delta_v * u_r + 2.0 * k_invariant)
1616 c_coef = delta_v * u_r * u_r
1617 disc = b_coef * b_coef - 4.0 * delta_v * c_coef
1618 if disc < 0:
1619 return None
1620 sqrt_disc = np.sqrt(disc)
1621 if c_decay_initial < c_fixed:
1622 u_target = (-b_coef - sqrt_disc) / (2.0 * delta_v)
1623 else:
1624 u_target = (-b_coef + sqrt_disc) / (2.0 * delta_v)
1625 if u_target <= 0.0:
1626 return None
1627 theta_local = k_invariant * (2.0 * u_target + alpha) / (u_target - u_r) ** 2
1628 return float(theta_origin + theta_local)
1631def _outlet_crossing_langmuir(
1632 sorption: LangmuirSorption,
1633 k_invariant: float,
1634 v_origin: float,
1635 theta_origin: float,
1636 v_outlet: float,
1637) -> float | None:
1638 """θ at which a Langmuir DecayingShockWave reaches v_outlet.
1640 From V_s - v_origin = K·(K_L + c)^2 / c^2 ⇒ (K_L + c)/c = sqrt(Δv/K) =: ratio,
1641 so c = K_L/(ratio - 1). Substitute into the invariant for θ_local.
1643 Returns
1644 -------
1645 float or None
1646 Cumulative flow at crossing, or None if no crossing exists.
1647 """
1648 delta_v = v_outlet - v_origin
1649 ratio = np.sqrt(delta_v / k_invariant)
1650 if ratio <= 1.0:
1651 return None
1652 c_target = sorption.k_l / (ratio - 1.0)
1653 theta_local = k_invariant * ((sorption.k_l + c_target) ** 2 + sorption.a_coeff) / (c_target * c_target)
1654 return float(theta_origin + theta_local)