Coverage for src/gwtransport/fronttracking/waves.py: 93%
411 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +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 θ.
9The change from (V, t) to (V, θ) makes every wave velocity a property of the
10sorption isotherm alone — flow no longer enters into wave dynamics. Time-
11varying flow is absorbed entirely into the θ(t) mapping at the API boundary;
12no wave needs recreation when the flow rate changes.
14This file is part of gwtransport which is released under AGPL-3.0 license.
15See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
16"""
18from abc import ABC, abstractmethod
19from dataclasses import dataclass, field
21import numpy as np
22from scipy.interpolate import CubicSpline
23from scipy.optimize import brentq
25from gwtransport.fronttracking.math import (
26 _C_MIN,
27 BrooksCoreyConductivity,
28 FreundlichSorption,
29 LangmuirSorption,
30 NonlinearSorption,
31 SorptionModel,
32 characteristic_speed,
33)
35# Numerical tolerance constants
36EPSILON_POSITION = 1e-15 # Tolerance for checking if two positions are equal
37DECAYING_SHOCK_U_FLOOR = 1e-300 # Lower bracket bound for brentq on Freundlich u-invariant
38DECAYING_SHOCK_BRENTQ_XTOL = (
39 1e-14 # brentq absolute tolerance for monotone θ inversions (exhaustion, outlet, numerical)
40)
41# Cached numerical decay profile (see ``_build_decay_profile``): c-grid resolution, the
42# Gauss-Legendre panel order for the cumulative invariant integral, and the fraction of the
43# c-gap the grid stops short of a secant-speed pole (where ``θ_local → ∞``).
44DECAY_PROFILE_NODES = 6000
45DECAY_PROFILE_GAUSS_ORDER = 10
46DECAY_PROFILE_POLE_FLOOR = 1e-6
49@dataclass
50class Wave(ABC):
51 """Abstract base class for all wave types in front tracking.
53 All waves share common attributes and must implement methods for
54 computing position and concentration. Waves can be active or inactive
55 (deactivated waves are preserved for history but don't participate in
56 future interactions).
58 Parameters
59 ----------
60 theta_start : float
61 Cumulative flow at which the wave forms [m³].
62 v_start : float
63 Position at which the wave forms [m³].
64 is_active : bool, optional
65 Whether wave is currently active. Default True.
66 """
68 theta_start: float
69 """Cumulative flow at which the wave forms [m³]."""
70 v_start: float
71 """Position at which the wave forms [m³]."""
72 is_active: bool = field(default=True, kw_only=True)
73 """Whether wave is currently active (in the solver's event-loop sense)."""
74 theta_deactivation: float = field(default=float("inf"), kw_only=True)
75 """Cumulative flow at which the wave was deactivated (default ``+∞``).
77 Historical record set by collision handlers when a wave is replaced
78 (e.g., a parent rarefaction superseded by a ``DecayingShockWave``).
79 ``is_active = False`` is the "current state" flag the solver uses for
80 its event loop; ``theta_deactivation`` is the moment in θ-history when
81 the wave stopped contributing. Retrospective queries (any θ in the
82 past) must use ``was_active_at(theta)`` instead of ``is_active`` so
83 that ``compute_domain_mass`` etc. correctly attribute c at v_outlet
84 during the wave's lifetime even after later events have deactivated
85 the wave.
86 """
88 def was_active_at(self, theta: float) -> bool:
89 """Whether the wave was active at cumulative flow ``theta`` (geometric truth).
91 Use for retrospective queries — ``is_active`` reflects only the
92 wave's *current* (post-simulation) state, which is wrong for
93 ``compute_domain_mass`` and similar at θ before a deactivation event.
95 Parameters
96 ----------
97 theta : float
98 Cumulative flow at which to query historical activity [m³].
100 Returns
101 -------
102 bool
103 ``True`` for ``theta_start <= theta < theta_deactivation``.
104 A wave constructed with ``is_active=False`` and no recorded
105 ``theta_deactivation`` (default ``+∞``) is treated as
106 never-active — e.g., synthetic test fixtures that want the
107 wave excluded from dispatch entirely.
108 """
109 if not self.is_active and self.theta_deactivation == float("inf"):
110 return False
111 return self.theta_start <= theta < self.theta_deactivation
113 def deactivate(self, theta: float) -> None:
114 """Mark the wave inactive at cumulative flow ``theta`` (collision handler API).
116 Sets both ``is_active = False`` (solver event-loop flag) and
117 ``theta_deactivation = theta`` (historical record for retrospective
118 ``was_active_at`` queries).
120 Parameters
121 ----------
122 theta : float
123 Cumulative flow at which the wave is deactivated [m³].
124 """
125 self.is_active = False
126 self.theta_deactivation = theta
128 @abstractmethod
129 def position_at_theta(self, theta: float) -> float | None:
130 """Compute wave position at cumulative flow θ.
132 Parameters
133 ----------
134 theta : float
135 Cumulative flow [m³].
137 Returns
138 -------
139 position : float or None
140 Position [m³], or None if θ < θ_start or θ >= theta_deactivation.
141 (Past-θ queries respect the wave's historical lifetime; current-state
142 queries before deactivation behave identically to the ``is_active``
143 check.)
144 """
146 @abstractmethod
147 def concentration_left(self) -> float:
148 """Concentration on the left (upstream) side of the wave."""
150 @abstractmethod
151 def concentration_right(self) -> float:
152 """Concentration on the right (downstream) side of the wave."""
154 @abstractmethod
155 def concentration_at_point(self, v: float, theta: float) -> float | None:
156 """Compute concentration at point (v, θ) if the wave controls it.
158 Returns
159 -------
160 concentration : float or None
161 Concentration [mass/volume] if the wave controls this point, None
162 otherwise.
163 """
166@dataclass
167class CharacteristicWave(Wave):
168 """Characteristic line along which concentration is constant.
170 In smooth regions, concentration travels at speed ``1/R(C)`` in (V, θ)
171 coordinates. Along each characteristic line, the concentration value is
172 constant. This is the fundamental solution element for hyperbolic
173 conservation laws.
175 Parameters
176 ----------
177 theta_start : float
178 Formation cumulative flow [m³].
179 v_start : float
180 Starting position [m³].
181 concentration : float
182 Constant concentration carried [mass/volume].
183 sorption : SorptionModel
184 Sorption model determining the speed.
185 is_active : bool, optional
186 Activity status. Default True.
188 Examples
189 --------
190 >>> sorption = FreundlichSorption(
191 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
192 ... )
193 >>> char = CharacteristicWave(
194 ... theta_start=0.0, v_start=0.0, concentration=5.0, sorption=sorption
195 ... )
196 >>> speed = char.speed()
197 >>> bool(np.isclose(char.position_at_theta(1000.0), speed * 1000.0))
198 True
199 """
201 concentration: float
202 """Constant concentration carried [mass/volume]."""
203 sorption: SorptionModel
204 """Sorption model determining the speed."""
205 _speed: float = field(init=False, repr=False, compare=False)
206 """Cached characteristic speed (immutable inputs; set in ``__post_init__``)."""
208 def __post_init__(self) -> None:
209 """Cache the (immutable) characteristic speed once."""
210 self._speed = characteristic_speed(self.concentration, self.sorption)
212 def speed(self) -> float:
213 """Characteristic speed dV/dθ = 1/R(C) (``+∞`` at a saturated state, R = 0)."""
214 return self._speed
216 def position_at_theta(self, theta: float) -> float | None:
217 """Position at cumulative flow θ.
219 ``V(θ) = v_start + speed * (θ - θ_start)``.
220 """
221 if not self.was_active_at(theta):
222 return None
223 return self.v_start + self.speed() * (theta - self.theta_start)
225 def concentration_left(self) -> float:
226 """Concentration on the left (upstream) side; equals the carried value."""
227 return self.concentration
229 def concentration_right(self) -> float:
230 """Concentration on the right (downstream) side; equals the carried value."""
231 return self.concentration
233 def concentration_at_point(self, v: float, theta: float) -> float | None:
234 """Return the carried concentration if the characteristic has reached ``v`` by θ."""
235 v_at_theta = self.position_at_theta(theta)
236 if v_at_theta is None:
237 return None
239 if v_at_theta >= v:
240 return self.concentration
242 return None
245@dataclass
246class ShockWave(Wave):
247 """Shock wave (discontinuity) with jump in concentration.
249 Shocks form when faster water overtakes slower water, creating a sharp
250 front. In (V, θ) the shock speed is given by the Rankine-Hugoniot
251 condition and is independent of flow::
253 dV_s/dθ = (C_R - C_L) / (C_T(C_R) - C_T(C_L))
255 Parameters
256 ----------
257 theta_start : float
258 Formation cumulative flow [m³].
259 v_start : float
260 Formation position [m³].
261 c_left : float
262 Concentration upstream (behind) shock [mass/volume].
263 c_right : float
264 Concentration downstream (ahead of) shock [mass/volume].
265 sorption : SorptionModel
266 Sorption model.
267 is_active : bool, optional
268 Activity status. Default True.
269 speed : float, optional
270 Shock speed dV/dθ. Computed from Rankine-Hugoniot in ``__post_init__``.
272 Examples
273 --------
274 >>> sorption = FreundlichSorption(
275 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
276 ... )
277 >>> shock = ShockWave(
278 ... theta_start=0.0,
279 ... v_start=0.0,
280 ... c_left=10.0,
281 ... c_right=2.0,
282 ... sorption=sorption,
283 ... )
284 >>> shock.speed > 0
285 True
286 >>> shock.satisfies_entropy()
287 True
288 """
290 c_left: float
291 """Concentration upstream (behind) shock [mass/volume]."""
292 c_right: float
293 """Concentration downstream (ahead of) shock [mass/volume]."""
294 sorption: SorptionModel
295 """Sorption model."""
296 speed: float = field(init=False)
297 """Shock speed dV/dθ; set in ``__post_init__``."""
299 def __post_init__(self) -> None:
300 """Compute shock speed from Rankine-Hugoniot in (V, θ)."""
301 self.speed = self.sorption.shock_speed(self.c_left, self.c_right)
303 def position_at_theta(self, theta: float) -> float | None:
304 """Position at cumulative flow θ. Shock propagates linearly in θ."""
305 if not self.was_active_at(theta):
306 return None
307 return self.v_start + self.speed * (theta - self.theta_start)
309 def concentration_left(self) -> float:
310 """Upstream concentration of the shock."""
311 return self.c_left
313 def concentration_right(self) -> float:
314 """Downstream concentration of the shock."""
315 return self.c_right
317 def concentration_at_point(self, v: float, theta: float) -> float | None:
318 """Return c_left if upstream of the shock at θ, c_right if downstream.
320 At the exact shock position the average is returned (convention; the
321 shock is infinitesimally thin in practice).
322 """
323 v_shock = self.position_at_theta(theta)
324 if v_shock is None:
325 return None
327 # Position-scaled face width (~1 ULP at all positions), matching
328 # DecayingShockWave.concentration_at_point; a fixed 1e-15 falls below
329 # one ULP for any v_shock > ~1 m³ and degenerates to bit-equality.
330 tol = 1e-15 * max(abs(v_shock), 1.0)
332 if v < v_shock - tol:
333 return self.c_left
334 if v > v_shock + tol:
335 return self.c_right
336 return 0.5 * (self.c_left + self.c_right)
338 def satisfies_entropy(self) -> bool:
339 """Check Lax entropy condition in (V, θ): ``λ_θ(C_L) ≥ s ≥ λ_θ(C_R)``."""
340 return self.sorption.check_entropy_condition(self.c_left, self.c_right, self.speed)
343@dataclass
344class RarefactionWave(Wave):
345 """Rarefaction (expansion fan) with smooth concentration gradient.
347 Rarefactions form when slower water follows faster water, creating an
348 expanding region where concentration varies smoothly. In (V, θ) the
349 solution is self-similar in ``(V - v_start)`` vs ``(θ - θ_start)``::
351 R(C) = (θ - θ_start) / (V - v_start)
353 Head and tail propagate at flow-free speeds ``1/R(C_head)`` and
354 ``1/R(C_tail)``.
356 Parameters
357 ----------
358 theta_start : float
359 Formation cumulative flow [m³].
360 v_start : float
361 Formation position [m³].
362 c_head : float
363 Concentration at leading edge (faster) [mass/volume].
364 c_tail : float
365 Concentration at trailing edge (slower) [mass/volume].
366 sorption : SorptionModel
367 Sorption model (must be concentration-dependent).
368 is_active : bool, optional
369 Activity status. Default True.
371 Raises
372 ------
373 ValueError
374 If head speed <= tail speed (would be a compression, not a rarefaction).
376 Examples
377 --------
378 >>> sorption = FreundlichSorption(
379 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
380 ... )
381 >>> raref = RarefactionWave(
382 ... theta_start=0.0,
383 ... v_start=0.0,
384 ... c_head=10.0,
385 ... c_tail=2.0,
386 ... sorption=sorption,
387 ... )
388 >>> raref.head_speed() > raref.tail_speed()
389 True
390 >>> raref.contains_point(v=150.0, theta=2000.0)
391 True
392 """
394 c_head: float
395 """Concentration at leading edge (faster) [mass/volume]."""
396 c_tail: float
397 """Concentration at trailing edge (slower) [mass/volume]."""
398 sorption: SorptionModel
399 """Sorption model (must be concentration-dependent)."""
400 _head_speed: float = field(init=False, repr=False, compare=False)
401 """Cached head celerity (immutable inputs; set in ``__post_init__``)."""
402 _tail_speed: float = field(init=False, repr=False, compare=False)
403 """Cached tail celerity (immutable inputs; set in ``__post_init__``)."""
405 def __post_init__(self):
406 """Cache head/tail celerities and verify this is a rarefaction (head faster than tail)."""
407 self._head_speed = characteristic_speed(self.c_head, self.sorption)
408 self._tail_speed = characteristic_speed(self.c_tail, self.sorption)
410 if self._head_speed <= self._tail_speed:
411 msg = (
412 f"Not a rarefaction: head_speed={self._head_speed:.6g} <= tail_speed={self._tail_speed:.6g}. "
413 f"This would be a compression (shock) instead."
414 )
415 raise ValueError(msg)
417 def head_speed(self) -> float:
418 """Speed of rarefaction head dV/dθ = 1/R(C_head) (``+∞`` at a saturated state, R = 0)."""
419 return self._head_speed
421 def tail_speed(self) -> float:
422 """Speed of rarefaction tail dV/dθ = 1/R(C_tail) (``+∞`` at a saturated state, R = 0)."""
423 return self._tail_speed
425 def head_position_at_theta(self, theta: float) -> float | None:
426 """Position of rarefaction head at cumulative flow θ."""
427 if not self.was_active_at(theta):
428 return None
429 return self.v_start + self.head_speed() * (theta - self.theta_start)
431 def tail_position_at_theta(self, theta: float) -> float | None:
432 """Position of rarefaction tail at cumulative flow θ."""
433 if not self.was_active_at(theta):
434 return None
435 return self.v_start + self.tail_speed() * (theta - self.theta_start)
437 def position_at_theta(self, theta: float) -> float | None:
438 """Head position (leading edge of rarefaction). Implements abstract Wave method."""
439 return self.head_position_at_theta(theta)
441 def contains_point(self, v: float, theta: float) -> bool:
442 """Return ``True`` if ``(v, θ)`` lies between the fan's tail and head."""
443 if theta <= self.theta_start or theta >= self.theta_deactivation:
444 return False
446 v_head = self.head_position_at_theta(theta)
447 v_tail = self.tail_position_at_theta(theta)
449 if v_head is None or v_tail is None:
450 return False
452 return v_tail <= v <= v_head
454 def concentration_left(self) -> float:
455 """Upstream concentration is the trailing-edge value c_tail."""
456 return self.c_tail
458 def concentration_right(self) -> float:
459 """Downstream concentration is the leading-edge value c_head."""
460 return self.c_head
462 def concentration_at_point(self, v: float, theta: float) -> float | None:
463 """Self-similar concentration inside the fan: ``R(C) = (θ - θ_start)/(v - v_start)``.
465 Outside the fan returns None. For ``ConstantRetardation``, rarefactions
466 don't form (all concentrations travel at the same speed), so this also
467 returns None.
469 Examples
470 --------
471 >>> sorption = FreundlichSorption(
472 ... k_f=0.01, n=2.0, bulk_density=1500.0, porosity=0.3
473 ... )
474 >>> raref = RarefactionWave(0.0, 0.0, 10.0, 2.0, sorption)
475 >>> c = raref.concentration_at_point(v=150.0, theta=2000.0)
476 >>> c is not None
477 True
478 >>> 2.0 <= c <= 10.0
479 True
480 """
481 if abs(v - self.v_start) < EPSILON_POSITION and theta >= self.theta_start:
482 return self.c_tail
484 if not self.contains_point(v, theta):
485 return None
487 r_target = (theta - self.theta_start) / (v - self.v_start)
489 if r_target <= 1.0:
490 return None # Unphysical
492 try:
493 c = self.sorption.concentration_from_retardation(r_target)
494 except NotImplementedError:
495 # ConstantRetardation case — rarefactions don't form
496 return None
498 # contains_point(v, theta) was True, so the point is geometrically inside
499 # the fan. The inverted c may drift by a few ULPs past [c_tail, c_head]
500 # — clamp rather than rejecting so callers at the head/tail boundaries
501 # get the correct boundary concentration.
502 c_lo = min(self.c_tail, self.c_head)
503 c_hi = max(self.c_tail, self.c_head)
504 return min(max(float(c), c_lo), c_hi)
507@dataclass
508class DecayingShockWave(Wave):
509 r"""Merging shock with closed-form (or quadrature) trajectory in θ-space.
511 Formed when a rarefaction fan and a shock collide. The shock then has
512 one side fed by the fan's self-similar profile (the "decay" side) and
513 the other side at the original outer state (the "fixed" side). Valid for
514 any :class:`~gwtransport.fronttracking.math.NonlinearSorption`.
516 Two collision regimes are supported via ``decay_side``:
518 - ``'left'`` (favorable head-collision): the rarefaction's head (faster)
519 catches a leading shock. After collision, the shock's ``c_left`` decays
520 from the rarefaction head value toward ``c_fan_tail`` (the unchanged
521 downstream c_right is ``c_fixed``).
522 - ``'right'`` (unfavorable tail-collision, n<1 mirrored): a trailing shock
523 catches the rarefaction's tail. After collision, the shock's ``c_right``
524 decays from the rarefaction tail value toward ``c_fan_tail`` (the
525 unchanged upstream c_left is ``c_fixed``).
527 The wave is valid only while ``c_decay ∈ (c_fan_tail, c_decay_initial]``;
528 once ``c_decay`` reaches ``c_fan_tail`` the fan is exhausted (see the
529 solver's ``DSW_FAN_EXHAUSTED`` event).
531 **Dispatch.** ``_c_decay_at_theta_local`` is the single dispatch site
532 (position, fan-exhaustion and outlet-crossing all route through it): a
533 closed form is used where one exists, otherwise the per-wave cached numerical
534 profile (:func:`_build_decay_profile`). No combination raises — any
535 :class:`~gwtransport.fronttracking.math.NonlinearSorption` is valid. With
536 ``θ_local := θ − theta_origin`` measured from the rarefaction apex,
537 ``α := ρ_b · k_f / n_por`` for Freundlich, and ``u_d := c_decay^(1/n)``:
539 - Freundlich, ``c_fixed = 0`` (general ``n > 0``, ``n ≠ 1``) — closed form:
540 invariant ``θ_local · u_d^n = K · (n · u_d^(n-1) + α)``,
541 position ``V_s(θ) = v_origin + n · K / u_d(θ)``.
542 - Freundlich, ``c_fixed > 0``, ``n = 2`` and ``c_decay_initial > c_fixed``
543 — closed form: invariant ``(u_d - u_R)² · θ_local = K · (2 u_d + α)``
544 with ``u_R := c_fixed^(1/2)``,
545 position ``V_s(θ) = v_origin + 2 K · u_d(θ) / (u_d - u_R)²``.
546 (The ``c_decay_initial < c_fixed`` mirror falls through to numerical.)
547 - Langmuir, ``c_fixed = 0`` — closed form:
548 invariant ``θ_local · c_d² = K · ((K_L + c_d)² + a)`` with
549 ``a := ρ_b · s_max · K_L / n_por``,
550 position ``V_s(θ) = v_origin + K · (K_L + c_d)² / c_d²``.
551 - Brooks-Corey, ``c_fixed = 0`` — closed form:
552 invariant ``θ_local ∝ R(c_decay)^{a/(a−1)}`` (``R·S = 1/a`` constant),
553 so ``R(c_d) = R(c0)·(θ_local/θ_local_coll)^{(a−1)/a}``.
554 - Every other ``(isotherm, c_fixed)`` combination (Freundlich ``c_fixed>0,
555 n≠2``, Langmuir/Brooks-Corey ``c_fixed>0``, any van Genuchten) — cached
556 numerical profile (:func:`_build_decay_profile`): the decay-agnostic
557 invariant ``θ_local(c_d) = θ_local_coll · exp(∫ R'/[(1 − R·S)·R] dc)`` with
558 the symmetric secant speed ``S = (c − c_fixed)/(C_T(c) − C_T(c_fixed))``,
559 built once by composite quadrature and inverted for ``c_d(θ)`` by monotone
560 spline interpolation.
562 Every path shares the fan-continuity identity
563 ``V_s = v_origin + θ_local / R(c_decay)``, which ``position_at_theta`` and
564 ``outlet_crossing_theta`` use uniformly across all isotherms.
566 The invariant constant ``K`` (closed-form Freundlich/Langmuir only) is set
567 in ``__post_init__`` from the collision IC ``(theta_start, c_decay_initial)``.
569 Parameters
570 ----------
571 theta_start : float
572 Cumulative flow at which the merged wave forms (collision θ) [m³].
573 v_start : float
574 Position at which the merged wave forms [m³]. Should equal
575 ``v_origin + (V_s) at θ=theta_start`` for a fan-consistent
576 construction.
577 c_decay_initial : float
578 Concentration on the decaying side at θ=theta_start [mass/volume].
579 Must be non-negative; a fully-drained collision value of ``0`` is
580 floored to the shared dry-soil singularity floor ``_C_MIN`` so the
581 retardation and secant-speed evaluations stay finite (issue #222).
582 c_fixed : float
583 Concentration on the non-decaying side [mass/volume]. Constant in θ.
584 Non-negative.
585 c_fan_tail : float
586 Concentration at the fan's far boundary [mass/volume]. The wave is
587 valid only while ``c_decay ∈ (c_fan_tail, c_decay_initial]``; at
588 ``c_fan_tail`` the fan is exhausted. Non-negative.
589 decay_side : str
590 ``'left'`` or ``'right'``. See class docstring.
591 v_origin : float
592 Position of the rarefaction apex [m³].
593 theta_origin : float
594 Cumulative flow at the rarefaction apex [m³]. Must satisfy
595 ``theta_origin < theta_start``.
596 sorption : NonlinearSorption
597 Sorption model (any concentration-dependent isotherm).
598 is_active : bool, optional
599 Activity flag. Default True.
601 See Also
602 --------
603 ShockWave : Linear-θ shock (no decaying side).
604 RarefactionWave : Self-similar expansion fan.
605 """
607 c_decay_initial: float
608 """Concentration on the decaying side at θ=theta_start [mass/volume]."""
609 c_fixed: float
610 """Concentration on the non-decaying side [mass/volume]."""
611 c_fan_tail: float
612 """Concentration at the fan's far boundary [mass/volume]; bounds the decay."""
613 decay_side: str
614 """``'left'`` (favorable head-collision) or ``'right'`` (n<1 mirrored)."""
615 v_origin: float
616 """Position of the rarefaction apex [m³]."""
617 theta_origin: float
618 """Cumulative flow at the rarefaction apex [m³]."""
619 sorption: NonlinearSorption
620 """Sorption model (any concentration-dependent isotherm)."""
621 K: float = field(init=False)
622 """Invariant constant set in ``__post_init__`` (closed-form Freundlich ``c_fixed=0``/``n≈2`` and Langmuir
623 ``c_fixed=0`` cases; ``nan`` for every numerical case)."""
624 _freundlich_cf: bool = field(init=False, repr=False, compare=False)
625 """Cached Freundlich-closed-form predicate (immutable inputs; set in ``__post_init__``)."""
626 _langmuir_cf: bool = field(init=False, repr=False, compare=False)
627 """Cached Langmuir-closed-form predicate."""
628 _brooks_corey_cf: bool = field(init=False, repr=False, compare=False)
629 """Cached Brooks-Corey ``c_fixed=0`` closed-form predicate."""
630 _numerical: bool = field(init=False, repr=False, compare=False)
631 """Cached predicate: no closed form applies, so the decay routes to the cached numerical profile."""
632 _decay_profile_cache: tuple | None = field(default=None, init=False, repr=False, compare=False)
633 """Lazily-built monotone ``θ_local(c)`` map for the numerical decay path (see ``_decay_profile``)."""
635 def __post_init__(self) -> None:
636 """Validate inputs and compute the closed-form invariant K when applicable."""
637 if self.decay_side not in {"left", "right"}:
638 msg = f"decay_side must be 'left' or 'right', got {self.decay_side!r}"
639 raise ValueError(msg)
640 if self.c_decay_initial < 0.0:
641 msg = f"c_decay_initial must be non-negative, got {self.c_decay_initial}"
642 raise ValueError(msg)
643 # Floor a fully-drained fan tail (c_decay_initial == 0, #222) to _C_MIN so the
644 # retardation and secant-speed evaluations stay finite (package floor convention).
645 self.c_decay_initial = max(self.c_decay_initial, _C_MIN)
646 if self.c_fixed < 0.0:
647 msg = f"c_fixed must be non-negative, got {self.c_fixed}"
648 raise ValueError(msg)
649 if self.c_fan_tail < 0.0:
650 msg = f"c_fan_tail must be non-negative, got {self.c_fan_tail}"
651 raise ValueError(msg)
652 if self.theta_origin >= self.theta_start:
653 msg = (
654 f"theta_origin ({self.theta_origin}) must be strictly less than "
655 f"theta_start ({self.theta_start}); rarefaction apex precedes collision"
656 )
657 raise ValueError(msg)
659 if not isinstance(self.sorption, NonlinearSorption):
660 msg = f"DecayingShockWave requires a NonlinearSorption, got {type(self.sorption).__name__}"
661 raise TypeError(msg)
663 # Classify the decay path once (immutable inputs). Closed forms exist for:
664 # Freundlich c_fixed=0 (general n) or the n≈2 quadratic when the decaying side
665 # starts above the fixed side (the +√ / (u_d−u_R)² branch was derived only for
666 # c_decay_initial > c_fixed — the mirror routes to the numerical profile);
667 # Langmuir c_fixed=0; Brooks-Corey c_fixed=0. Everything else is numerical.
668 s = self.sorption
669 self._freundlich_cf = isinstance(s, FreundlichSorption) and (
670 self.c_fixed == 0.0 or bool(np.isclose(s.n, 2.0, rtol=1e-12) and self.c_decay_initial > self.c_fixed)
671 )
672 self._langmuir_cf = isinstance(s, LangmuirSorption) and self.c_fixed == 0.0
673 self._brooks_corey_cf = isinstance(s, BrooksCoreyConductivity) and self.c_fixed == 0.0
674 self._numerical = not (self._freundlich_cf or self._langmuir_cf or self._brooks_corey_cf)
676 # K is the closed-form invariant constant; set only for the Freundlich/Langmuir
677 # closed forms and left NaN for every numerical (and Brooks-Corey) case. The
678 # isinstance guards narrow ``s`` for the typed helpers (the cached predicate already
679 # implies the type; the ``np.isclose`` cost is not re-incurred).
680 self.K = float("nan")
681 if self._freundlich_cf and isinstance(s, FreundlichSorption):
682 self.K = _compute_k_freundlich(
683 s,
684 self.theta_start - self.theta_origin,
685 self.c_decay_initial,
686 self.c_fixed,
687 )
688 elif self._langmuir_cf and isinstance(s, LangmuirSorption):
689 self.K = _compute_k_langmuir(
690 s,
691 self.theta_start - self.theta_origin,
692 self.c_decay_initial,
693 )
695 def _decay_profile(self) -> tuple:
696 """Lazily build & cache the monotone ``θ_local(c)`` map for the numerical decay path.
698 Returns ``(c_of_i, i_max, c_limit_node)``: a ``CubicSpline`` mapping the
699 cumulative invariant ``I = ln(θ_local/θ_local_coll)`` to ``c_decay``, the
700 largest resolved ``I`` (endpoint of the reachable c-range), and the c at
701 that endpoint. Built once per wave (see :func:`_build_decay_profile`).
702 """
703 if self._decay_profile_cache is None:
704 self._decay_profile_cache = _build_decay_profile(
705 self.sorption,
706 self.c_decay_initial,
707 self.c_fixed,
708 self.c_fan_tail,
709 )
710 return self._decay_profile_cache
712 def c_decay_at_theta(self, theta: float) -> float | None:
713 """Concentration on the decaying side at cumulative flow θ.
715 Returns ``None`` for ``θ < theta_start`` or when the wave is inactive;
716 otherwise delegates to the single per-isotherm dispatch in
717 ``_c_decay_at_theta_local``.
718 """
719 if not self.was_active_at(theta):
720 return None
721 return self._c_decay_at_theta_local(theta - self.theta_origin)
723 def position_at_theta(self, theta: float) -> float | None:
724 """Shock position ``V_s(θ)`` via the fan-continuity identity.
726 ``V_s = v_origin + θ_local / R(c_decay)`` for every isotherm. Returns
727 ``None`` for ``θ < theta_start`` or when inactive.
728 """
729 if not self.was_active_at(theta):
730 return None
732 theta_local = theta - self.theta_origin
733 c_d = self._c_decay_at_theta_local(theta_local)
734 return float(self.v_origin + theta_local / float(self.sorption.retardation(c_d)))
736 def theta_at_fan_exhaustion(self) -> float | None:
737 """Cumulative flow θ at which ``c_decay`` reaches ``c_fan_tail``.
739 ``c_decay(θ)`` is strictly monotone from ``c_decay_initial`` toward
740 ``c_fan_tail``, so the exhaustion θ is well-defined. The crossing test is
741 orientation-agnostic: it holds for both the shrinking decay
742 (``c_decay_initial > c_fan_tail``) and the growing decay
743 (``c_decay_initial < c_fan_tail``). Returns ``None`` when ``c_fan_tail``
744 is not strictly between ``c_fixed`` and ``c_decay_initial`` — e.g. full drying
745 (``c_fan_tail == c_fixed``), where the decay asymptotically merges with
746 the fixed state and no finite exhaustion event occurs.
748 Returns
749 -------
750 float or None
751 Cumulative flow θ at exhaustion, or ``None`` if not reached.
752 """
753 # An interior exhaustion needs c_fan_tail strictly between c_fixed and
754 # c_decay_initial (orientation-agnostic via min/max). Full drying
755 # (c_fan_tail == c_fixed) merges asymptotically with no finite crossing —
756 # return None rather than grow the bracket forever (van Genuchten would hang).
757 c_lo = min(self.c_fixed, self.c_decay_initial)
758 c_hi = max(self.c_fixed, self.c_decay_initial)
759 if not (c_lo < self.c_fan_tail < c_hi):
760 return None
762 theta_local_collision = self.theta_start - self.theta_origin
764 if self._numerical:
765 # The numerical forward map saturates AT c_fan_tail (it never crosses it),
766 # so a forward-map bracket cannot see the crossing for either decay
767 # orientation. Evaluate θ_local(c_fan_tail) from the un-clamped invariant
768 # directly: the gate above guarantees c_fan_tail is the reachable limit, so
769 # it is exactly the cached profile's endpoint ``i_max``.
770 _c_of_i, i_max, _c_limit_node = self._decay_profile()
771 return self.theta_origin + theta_local_collision * float(np.exp(i_max))
773 # Closed forms cross c_fan_tail smoothly (always a shrinking decay); invert the
774 # monotone forward map by bracketing (orientation-agnostic — no early return).
775 def f(theta_local: float) -> float:
776 return self._c_decay_at_theta_local(theta_local) - self.c_fan_tail
778 theta_local_exhaust = _invert_monotone_theta_local(
779 f, theta_hi_seed=theta_local_collision, f_seed=f(theta_local_collision)
780 )
781 if theta_local_exhaust is None:
782 return None
783 return self.theta_origin + theta_local_exhaust
785 def _c_decay_at_theta_local(self, theta_local: float) -> float:
786 """Decaying concentration as a function of ``θ_local`` (apex-relative).
788 The SOLE isotherm dispatch site: closed where an exact form exists,
789 otherwise the cached numerical decay profile. The closed forms
790 (Freundlich ``c_fixed=0`` or ``n≈2``; Langmuir ``c_fixed=0``;
791 Brooks-Corey ``c_fixed=0``) are selected by the ``__post_init__``
792 predicates; every other ``(isotherm, c_fixed)`` combination falls through
793 to the per-wave cached invariant profile. Takes ``θ_local`` directly and
794 skips the activity check. ``c_decay_at_theta``, ``position_at_theta``,
795 ``theta_at_fan_exhaustion`` and ``outlet_crossing_theta`` all route
796 through here rather than repeating the dispatch.
797 """
798 theta_local_collision = self.theta_start - self.theta_origin
799 s = self.sorption
800 if self._freundlich_cf and isinstance(s, FreundlichSorption):
801 return _c_decay_freundlich(
802 s, self.K, self.c_decay_initial, self.c_fixed, theta_local_collision, theta_local
803 )
804 if self._langmuir_cf and isinstance(s, LangmuirSorption):
805 return _c_decay_langmuir(s, self.K, theta_local)
806 if self._brooks_corey_cf and isinstance(s, BrooksCoreyConductivity):
807 return _c_decay_brooks_corey(s, self.c_decay_initial, theta_local_collision, theta_local)
809 # Numerical path: invert the cached monotone θ_local(c) map. c_decay stays at
810 # c_decay_initial up to the collision and clamps at the reachable c-limit past it.
811 if theta_local <= theta_local_collision:
812 return self.c_decay_initial
813 c_of_i, i_max, c_limit_node = self._decay_profile()
814 i_target = np.log(theta_local / theta_local_collision)
815 if i_target >= i_max:
816 return float(c_limit_node)
817 return float(c_of_i(i_target))
819 def outlet_crossing_theta(self, v_outlet: float) -> float | None:
820 """Cumulative flow at which ``V_s = v_outlet``.
822 Returns ``None`` if the outlet is upstream of the wave's birth
823 position or no crossing exists in ``(theta_start, +∞)``. The wave's
824 current activity flag is not consulted — callers asking
825 retrospectively about a historical crossing need the answer regardless
826 of subsequent deactivation.
828 The closed-form Freundlich/Langmuir cases invert the fan-continuity
829 identity ``V_s − v_origin = θ_local / R(c_decay)`` analytically (valid
830 only when ``_c_decay_at_theta_local`` itself uses the closed form, so
831 the same conditions are mirrored here); every other case inverts the
832 monotone ``V_s(θ)`` via ``brentq``.
833 """
834 if v_outlet <= self.v_start:
835 return None
837 # V_s is monotonically increasing in θ (positive shock speed); invert
838 # via the fan-continuity identity V_s - v_origin = θ_local / R(c_decay)
839 # combined with the invariant to eliminate u, then solve for θ.
840 s = self.sorption
841 if self._freundlich_cf and isinstance(s, FreundlichSorption):
842 return _outlet_crossing_freundlich(
843 s,
844 self.K,
845 self.c_fixed,
846 self.v_origin,
847 self.theta_origin,
848 v_outlet,
849 )
850 if self._langmuir_cf and isinstance(s, LangmuirSorption):
851 return _outlet_crossing_langmuir(
852 s,
853 self.K,
854 self.v_origin,
855 self.theta_origin,
856 v_outlet,
857 )
858 return self._outlet_crossing_numerical(v_outlet)
860 def _outlet_crossing_numerical(self, v_outlet: float) -> float | None:
861 """θ at which ``V_s = v_outlet`` for every non-closed-form case.
863 ``V_s(θ) = v_origin + θ_local / R(c_decay(θ))`` is monotone increasing;
864 invert by ``brentq`` on ``θ_local``.
865 """
866 theta_local_collision = self.theta_start - self.theta_origin
868 def f(theta_local: float) -> float:
869 c = self._c_decay_at_theta_local(theta_local)
870 return self.v_origin + theta_local / float(self.sorption.retardation(c)) - v_outlet
872 f_lo = f(theta_local_collision)
873 if f_lo >= 0.0:
874 # Already at/past the outlet at collision; the linear-shock guards
875 # in the solver handle the duplicate-crossing suppression.
876 return self.theta_start
877 # Seed at the collision (f_lo < 0 established above) and let the helper grow the
878 # bracket upward — no dimensional floor, so crossings within θ_local < 1 of the
879 # apex are found (mirrors theta_at_fan_exhaustion's closed-form bracket).
880 theta_local_cross = _invert_monotone_theta_local(f, theta_hi_seed=theta_local_collision, f_seed=f_lo)
881 if theta_local_cross is None:
882 return None
883 return self.theta_origin + theta_local_cross
885 def concentration_left(self) -> float:
886 """Concentration on the left (upstream) side at θ=theta_start.
888 For ``decay_side='left'`` returns the decaying c at the collision
889 moment; for ``decay_side='right'`` returns the fixed side.
890 """
891 return self.c_decay_initial if self.decay_side == "left" else self.c_fixed
893 def concentration_right(self) -> float:
894 """Concentration on the right (downstream) side at θ=theta_start.
896 For ``decay_side='right'`` returns the decaying c at the collision
897 moment; for ``decay_side='left'`` returns the fixed side.
898 """
899 return self.c_decay_initial if self.decay_side == "right" else self.c_fixed
901 def concentration_at_point(self, v: float, theta: float) -> float | None:
902 """Concentration at ``(v, θ)`` if controlled by this decaying shock.
904 Three regions:
906 1. ``v == V_s(θ)`` (within FP): average of decay-side and fixed-side c.
907 2. ``v > V_s(θ)`` (downstream): fixed-side c if ``decay_side='left'``;
908 decay-side c at θ if ``decay_side='right'``.
909 3. ``v < V_s(θ)`` (upstream, inside the fan): the fan's self-similar
910 concentration ``R(c) = (θ − theta_origin)/(v − v_origin)``. Outside
911 the fan — i.e. the decay-side characteristic from the apex hasn't
912 reached v yet, OR the point lies beyond the ``c_fan_tail`` boundary
913 (the fan's far edge) — returns ``None``.
915 Returns ``None`` for ``θ < theta_start`` or inactive waves.
916 """
917 if not self.was_active_at(theta):
918 return None
920 # Compute the decaying-side concentration once and derive V_s from it
921 # (inlining position_at_theta's body) so the shock-face branch can reuse
922 # it instead of re-running the numerical-isotherm root-find.
923 theta_local = theta - self.theta_origin
924 c_d = self._c_decay_at_theta_local(theta_local)
925 v_s = float(self.v_origin + theta_local / float(self.sorption.retardation(c_d)))
927 tol = 1e-15 * max(abs(v_s), 1.0)
929 if abs(v - v_s) < tol:
930 return 0.5 * (c_d + self.c_fixed)
932 # Region selection depends on decay_side:
933 # 'left' (favorable n>1, Langmuir): fan extends upstream of V_s
934 # (v < V_s), c_fixed downstream (v > V_s).
935 # 'right' (n<1 mirror): fan extends downstream of V_s (v > V_s),
936 # c_fixed upstream (v < V_s).
937 if self.decay_side == "left":
938 v_fan_side = v < v_s - tol
939 v_fixed_side = v > v_s + tol
940 else:
941 v_fan_side = v > v_s + tol
942 v_fixed_side = v < v_s - tol
944 if v_fixed_side:
945 return self.c_fixed
947 if not v_fan_side:
948 return None # within tol of shock face — handled above
950 # Fan-interior: self-similar profile with apex at (v_origin, theta_origin).
951 if v == self.v_origin:
952 return None
953 r_target = (theta - self.theta_origin) / (v - self.v_origin)
954 if r_target <= 1.0:
955 return None
956 try:
957 c_fan = self.sorption.concentration_from_retardation(r_target)
958 except NotImplementedError:
959 return None
960 c_fan = float(c_fan)
962 # The fan the DSW controls spans concentrations between the shock face
963 # (c_decay, ≤ c_decay_initial) and the fan's far boundary c_fan_tail.
964 # A point past c_fan_tail belongs to whatever lies beyond the fan, not
965 # to this wave — reject so the fan is not extended past its extent.
966 c_lo = min(self.c_fan_tail, self.c_decay_initial)
967 c_hi = max(self.c_fan_tail, self.c_decay_initial)
968 if c_fan < c_lo - EPSILON_POSITION or c_fan > c_hi + EPSILON_POSITION:
969 return None
970 return c_fan
973def _invert_monotone_theta_local(f, *, theta_hi_seed: float, f_seed: float | None = None) -> float | None:
974 """Bracket-then-brentq a monotone ``f(θ_local)`` with a sign change above the seed.
976 Shared by the closed-form branch of ``theta_at_fan_exhaustion`` and by
977 ``_outlet_crossing_numerical``: both invert a monotone function of
978 ``θ_local`` whose sign at the collision is already known to differ from its
979 sign at large ``θ_local``. Geometrically grows ``θ_hi`` (``×2``, ≤200 iters)
980 from ``theta_hi_seed`` until ``f`` flips sign, then inverts with ``brentq``.
981 Returns ``None`` if no sign change is bracketed within the iteration budget.
983 Parameters
984 ----------
985 f : callable
986 Monotone residual ``f(θ_local)``; ``f(theta_hi_seed)`` and the far-field
987 value must straddle zero. Caller-specific early sentinels
988 (already-past-outlet) are handled by the caller.
989 theta_hi_seed : float
990 ``θ_local`` lower bracket; the search grows ``θ_hi`` from here.
991 f_seed : float, optional
992 Pre-evaluated ``f(theta_hi_seed)``; callers that already computed it pass
993 it to avoid a redundant evaluation. ``None`` recomputes it here.
995 Returns
996 -------
997 float or None
998 Root ``θ_local`` of ``f``, or ``None`` if not bracketed.
999 """
1000 if f_seed is None:
1001 f_seed = f(theta_hi_seed)
1002 theta_hi = theta_hi_seed
1003 for _ in range(200):
1004 theta_hi *= 2.0
1005 if f(theta_hi) * f_seed < 0.0:
1006 return float(brentq(f, theta_hi_seed, theta_hi, xtol=DECAYING_SHOCK_BRENTQ_XTOL))
1007 return None
1010def _build_decay_profile(
1011 sorption: NonlinearSorption,
1012 c_decay_initial: float,
1013 c_fixed: float,
1014 c_fan_tail: float,
1015) -> tuple:
1016 r"""Build the per-wave monotone ``θ_local(c)`` map for collisions with no closed form.
1018 Decay-agnostic: the fan-continuity + Rankine-Hugoniot relations do not
1019 depend on which side decays. The secant speed
1020 ``S(c) = (c − c_fixed)/(C_T(c) − C_T(c_fixed))`` is symmetric in
1021 ``(c_decay, c_fixed)``, so the same invariant
1022 ``θ_local(c) = θ_local_coll · exp(I(c))``, ``I(c) = ∫_{c0}^{c} R'/[(1 − R·S)·R] dc``
1023 (``R'`` by central finite difference) holds for Freundlich ``c_fixed>0, n≠2``,
1024 Langmuir ``c_fixed>0``, Brooks-Corey ``c_fixed>0`` and any van Genuchten case
1025 alike. ``I(c)`` is built ONCE by a single vectorised composite Gauss-Legendre
1026 cumulative quadrature over a c-grid from ``c_decay_initial`` to the reachable
1027 limit, then inverted by monotone-spline interpolation — replacing the former
1028 quad-inside-brentq-inside-brentq scalar solve (~1000× fewer integrand evals
1029 across a record). The reachable limit is the fan tail ``c_fan_tail`` UNLESS
1030 ``c_fixed`` lies strictly between ``c_decay_initial`` and ``c_fan_tail`` — then
1031 the secant speed has a pole at ``c_fixed`` (``R·S → 1``, ``θ_local → ∞``): the
1032 shock asymptotes to the fixed state, so the grid stops a hair short of it and
1033 ``c_decay`` clamps there.
1035 Parameters
1036 ----------
1037 sorption : NonlinearSorption
1038 Sorption model.
1039 c_decay_initial : float
1040 Decaying-side concentration at the collision (``c0``) [mass/volume].
1041 c_fixed : float
1042 Non-decaying-side concentration [mass/volume].
1043 c_fan_tail : float
1044 Concentration at the fan's far boundary [mass/volume]; bounds the decay.
1046 Returns
1047 -------
1048 tuple
1049 ``(c_of_i, i_max, c_limit_node)``: a ``CubicSpline`` mapping the cumulative
1050 invariant ``I = ln(θ_local/θ_local_coll)`` to ``c_decay``, the endpoint
1051 ``I`` of the reachable c-range, and the ``c`` at that endpoint. ``I`` is
1052 collision-independent (``θ_local_coll`` enters only at query time).
1053 """
1054 ct_fixed = float(sorption.total_concentration(c_fixed))
1056 def integrand(c):
1057 c = np.asarray(c, dtype=float)
1058 h = np.maximum(1e-9, 1e-7 * np.abs(c))
1059 r_prime = (np.asarray(sorption.retardation(c + h)) - np.asarray(sorption.retardation(c - h))) / (2.0 * h)
1060 r = np.asarray(sorption.retardation(c))
1061 ct = np.asarray(sorption.total_concentration(c))
1062 secant = (c - c_fixed) / (ct - ct_fixed)
1063 return r_prime / ((1.0 - r * secant) * r)
1065 pole = (c_decay_initial - c_fixed) * (c_fan_tail - c_fixed) < 0.0
1066 c_limit = c_fixed if pole else c_fan_tail
1068 # c-grid from c0 to the reachable limit. Toward a pole (θ_local → ∞) the grid is
1069 # geometric, stopping a fraction DECAY_PROFILE_POLE_FLOOR of the gap short; the
1070 # non-pole grid reaches c_fan_tail exactly (so i_max is the exhaustion integral).
1071 frac = np.linspace(0.0, 1.0, DECAY_PROFILE_NODES)
1072 gap0 = abs(c_decay_initial - c_limit)
1073 gaps = gap0 * DECAY_PROFILE_POLE_FLOOR**frac if pole else gap0 * (1.0 - frac)
1074 c_nodes = c_limit + np.sign(c_decay_initial - c_limit) * gaps
1076 # Cumulative composite Gauss-Legendre integral of the invariant integrand.
1077 x_gl, w_gl = np.polynomial.legendre.leggauss(DECAY_PROFILE_GAUSS_ORDER)
1078 lo = c_nodes[:-1]
1079 hi = c_nodes[1:]
1080 mid = 0.5 * (lo + hi)
1081 half = 0.5 * (hi - lo)
1082 points = mid[:, None] + half[:, None] * x_gl[None, :]
1083 panel = (integrand(points.ravel()).reshape(points.shape) * w_gl[None, :]).sum(axis=1) * half
1084 i_nodes = np.concatenate([[0.0], np.cumsum(panel)])
1086 # Inverse-interpolation precondition: keep the strictly-increasing prefix (the
1087 # near-pole tail can lose monotonicity as the singular integrand outruns the grid).
1088 non_increasing = np.nonzero(np.diff(i_nodes) <= 0.0)[0]
1089 if non_increasing.size:
1090 cut = non_increasing[0] + 1
1091 i_nodes = i_nodes[:cut]
1092 c_nodes = c_nodes[:cut]
1093 c_of_i = CubicSpline(i_nodes, c_nodes)
1094 return c_of_i, float(i_nodes[-1]), float(c_nodes[-1])
1097def _c_decay_brooks_corey(
1098 sorption: BrooksCoreyConductivity,
1099 c_decay_initial: float,
1100 theta_local_collision: float,
1101 theta_local: float,
1102) -> float:
1103 r"""Brooks-Corey ``c_fixed = 0`` closed form for the decaying-side concentration.
1105 For Brooks-Corey with ``c_fixed = 0`` the product ``R·S = 1/a`` is constant
1106 (``a = sorption.a``), so the universal invariant integrates to
1107 ``θ_local ∝ R(c_decay)^{a/(a−1)}``. Inverting,
1108 ``R(c_d) = R(c0)·(θ_local/θ_local_coll)^{(a−1)/a}`` and
1109 ``c_d = concentration_from_retardation(R)``.
1111 Parameters
1112 ----------
1113 sorption : BrooksCoreyConductivity
1114 Sorption model.
1115 c_decay_initial : float
1116 Decaying-side concentration at the collision [mass/volume].
1117 theta_local_collision : float
1118 ``θ_local`` at the collision [m³].
1119 theta_local : float
1120 ``θ_local`` at which to evaluate the decaying concentration [m³].
1122 Returns
1123 -------
1124 float
1125 Decaying-side concentration ``c`` at ``theta_local``.
1126 """
1127 if theta_local <= theta_local_collision:
1128 return c_decay_initial
1129 a = sorption.a
1130 r0 = float(sorption.retardation(c_decay_initial))
1131 r_target = r0 * (theta_local / theta_local_collision) ** ((a - 1.0) / a)
1132 return float(sorption.concentration_from_retardation(r_target))
1135def _compute_k_freundlich(
1136 sorption: FreundlichSorption,
1137 theta_local: float,
1138 c_decay_initial: float,
1139 c_fixed: float,
1140) -> float:
1141 """Closed-form invariant K for Freundlich DecayingShockWave.
1143 Derivation: see plan §"Closed-form derivations". For c_fixed=0,
1144 K = θ_local · u_c^n / (n · u_c^(n-1) + α); for c_fixed>0 (n=2 only),
1145 K = θ_local · (u_c - u_r)^2 / (2 · u_c + α). Here α = ρ_b · k_f / n_por.
1147 Parameters
1148 ----------
1149 sorption : FreundlichSorption
1150 Sorption model.
1151 theta_local : float
1152 Cumulative flow from rarefaction apex to collision [m³].
1153 c_decay_initial : float
1154 Decaying-side concentration at the collision [mass/volume].
1155 c_fixed : float
1156 Non-decaying-side concentration [mass/volume].
1158 Returns
1159 -------
1160 float
1161 Invariant constant K.
1162 """
1163 n = sorption.n
1164 alpha = sorption.bulk_density * sorption.k_f / sorption.porosity
1165 u_d = c_decay_initial ** (1.0 / n)
1167 if c_fixed == 0.0:
1168 return float(theta_local * u_d**n / (n * u_d ** (n - 1.0) + alpha))
1170 # c_fixed > 0, n=2
1171 u_r = c_fixed**0.5
1172 return float(theta_local * (u_d - u_r) ** 2 / (2.0 * u_d + alpha))
1175def _compute_k_langmuir(
1176 sorption: LangmuirSorption,
1177 theta_local: float,
1178 c_decay_initial: float,
1179) -> float:
1180 """Closed-form invariant K for Langmuir DecayingShockWave (c_fixed=0).
1182 K = θ_local · c_d^2 / ((K_L + c_d)^2 + a) with a = ρ_b · s_max · K_L / n_por.
1184 Parameters
1185 ----------
1186 sorption : LangmuirSorption
1187 Sorption model.
1188 theta_local : float
1189 Cumulative flow from rarefaction apex to collision [m³].
1190 c_decay_initial : float
1191 Decaying-side concentration at the collision [mass/volume].
1193 Returns
1194 -------
1195 float
1196 Invariant constant K.
1197 """
1198 return float(theta_local * c_decay_initial**2 / ((sorption.k_l + c_decay_initial) ** 2 + sorption.a_coeff))
1201def _c_decay_freundlich(
1202 sorption: FreundlichSorption,
1203 k_invariant: float,
1204 c_decay_initial: float,
1205 c_fixed: float,
1206 theta_local_collision: float,
1207 theta_local: float,
1208) -> float:
1209 """Invert the Freundlich invariant to get c on the decaying side at θ_local.
1211 For n=2 c_fixed=0 (quadratic in u): closed form
1212 ``u = (K + sqrt(K^2 + K·θ_local·α)) / θ_local``. For general n with
1213 c_fixed=0 (transcendental): brentq on the monotone bracket
1214 ``(tiny, c_decay_initial^(1/n)]``. For n=2 c_fixed>0 (quadratic in u with u_r):
1215 closed form with positive sign chosen to give u > u_r as θ_local → ∞.
1217 Returns
1218 -------
1219 float
1220 Decaying-side concentration c at θ_local.
1221 """
1222 n = sorption.n
1223 alpha = sorption.bulk_density * sorption.k_f / sorption.porosity
1225 if c_fixed == 0.0:
1226 if np.isclose(n, 2.0, rtol=1e-12):
1227 disc = k_invariant * k_invariant + theta_local * k_invariant * alpha
1228 u = (k_invariant + np.sqrt(disc)) / theta_local
1229 return float(u * u)
1230 u_root = _invert_freundlich_cr_zero(k_invariant, c_decay_initial, n, alpha, theta_local_collision, theta_local)
1231 return float(u_root**n)
1233 # n=2, c_fixed > 0
1234 u_r = c_fixed**0.5
1235 disc = k_invariant * (theta_local * (2.0 * u_r + alpha) + k_invariant)
1236 u = (u_r * theta_local + k_invariant + np.sqrt(disc)) / theta_local
1237 return float(u * u)
1240def _invert_freundlich_cr_zero(
1241 k_invariant: float,
1242 c_decay_initial: float,
1243 n: float,
1244 alpha: float,
1245 theta_local_collision: float,
1246 theta_local: float,
1247) -> float:
1248 """Invert ``θ_local · u^n = K · (n·u^(n-1) + α)`` for u via brentq.
1250 Returns
1251 -------
1252 float
1253 Root u of the invariant at θ_local.
1254 """
1255 u_collision = c_decay_initial ** (1.0 / n)
1257 def f(u: float) -> float:
1258 return theta_local * u**n - k_invariant * (n * u ** (n - 1.0) + alpha)
1260 if theta_local <= theta_local_collision:
1261 # Earlier than (or at) collision: c_decay equals c_decay_initial.
1262 return u_collision
1264 u_root = brentq(f, DECAYING_SHOCK_U_FLOOR, u_collision, xtol=1e-15)
1265 return float(u_root) # type: ignore[arg-type]
1268def _c_decay_langmuir(sorption: LangmuirSorption, k_invariant: float, theta_local: float) -> float:
1269 """Invert the Langmuir invariant ``θ_local · c^2 = K · ((K_L+c)^2 + a)`` for c.
1271 Expanded: ``(θ_local - K)·c^2 - 2·K·K_L·c - K·(K_L^2 + a) = 0`` (quadratic
1272 in c). Positive root chosen for c > 0.
1274 Returns
1275 -------
1276 float
1277 Decaying-side concentration c at θ_local.
1278 """
1279 k_l = sorption.k_l
1280 a_coeff = sorption.a_coeff
1281 denom = theta_local - k_invariant
1282 disc = k_invariant * (k_invariant * k_l * k_l + denom * (k_l * k_l + a_coeff))
1283 return float((k_invariant * k_l + np.sqrt(disc)) / denom)
1286def _outlet_crossing_freundlich(
1287 sorption: FreundlichSorption,
1288 k_invariant: float,
1289 c_fixed: float,
1290 v_origin: float,
1291 theta_origin: float,
1292 v_outlet: float,
1293) -> float | None:
1294 """θ at which a Freundlich DecayingShockWave reaches v_outlet.
1296 Returns
1297 -------
1298 float or None
1299 Cumulative flow at crossing, or None if no crossing.
1300 """
1301 n = sorption.n
1302 alpha = sorption.bulk_density * sorption.k_f / sorption.porosity
1303 delta_v = v_outlet - v_origin
1305 if c_fixed == 0.0:
1306 u_target = n * k_invariant / delta_v
1307 if u_target <= 0.0:
1308 return None
1309 theta_local = k_invariant * (n * u_target ** (n - 1.0) + alpha) / u_target**n
1310 return float(theta_origin + theta_local)
1312 # n=2, c_fixed > 0: V_s - v_origin = 2·K·u / (u - u_r)^2 ⇒ quadratic in u.
1313 # The plus-sqrt root always satisfies u > u_r for K, alpha, delta_v > 0:
1314 # the quadratic's roots multiply to u_r² and sum to 2u_r + 2K/delta_v > 2u_r,
1315 # so exactly one root exceeds u_r and it is the plus-sqrt branch. The
1316 # minus-sqrt root is the unphysical companion.
1317 u_r = c_fixed**0.5
1318 b_coef = -(2.0 * delta_v * u_r + 2.0 * k_invariant)
1319 c_coef = delta_v * u_r * u_r
1320 disc = b_coef * b_coef - 4.0 * delta_v * c_coef
1321 if disc < 0:
1322 return None
1323 u_target = (-b_coef + np.sqrt(disc)) / (2.0 * delta_v)
1324 theta_local = k_invariant * (2.0 * u_target + alpha) / (u_target - u_r) ** 2
1325 return float(theta_origin + theta_local)
1328def _outlet_crossing_langmuir(
1329 sorption: LangmuirSorption,
1330 k_invariant: float,
1331 v_origin: float,
1332 theta_origin: float,
1333 v_outlet: float,
1334) -> float | None:
1335 """θ at which a Langmuir DecayingShockWave reaches v_outlet.
1337 From V_s - v_origin = K·(K_L + c)^2 / c^2 ⇒ (K_L + c)/c = sqrt(Δv/K) =: ratio,
1338 so c = K_L/(ratio - 1). Substitute into the invariant for θ_local.
1340 Returns
1341 -------
1342 float or None
1343 Cumulative flow at crossing, or None if no crossing exists.
1344 """
1345 delta_v = v_outlet - v_origin
1346 ratio = np.sqrt(delta_v / k_invariant)
1347 if ratio <= 1.0:
1348 return None
1349 c_target = sorption.k_l / (ratio - 1.0)
1350 theta_local = k_invariant * ((sorption.k_l + c_target) ** 2 + sorption.a_coeff) / (c_target * c_target)
1351 return float(theta_origin + theta_local)