Coverage for src/gwtransport/_radial_asr_dehoog.py: 95%
59 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1r"""Vectorized double-precision de Hoog numerical Laplace inversion.
3This private module provides :func:`dehoog_inverse`, an implementation of the de Hoog, Knight &
4Stokes (1982) accelerated Fourier-series method for inverting a Laplace transform
5:math:`\bar f(s) \to f(t)`. It is the foundational numerical primitive of the exact radial
6advection-dispersion module: the per-phase transfer functions (:mod:`gwtransport._radial_asr_kernels`)
7are known in closed form only in the Laplace domain, and the bin-level observable needs the
8real-time (here: real-flushed-volume) kernel and its antiderivatives, obtained by inverting
9:math:`\hat g/s` and :math:`\hat g/s^2`.
11Why de Hoog and not Talbot
12--------------------------
13The radial ASR observables include step-like resident profiles and sharp breakthroughs. The Talbot
14(deformed-contour) family is built for smooth, monotone-decaying responses and produces large
15spurious oscillations -- even negative concentrations -- on these step-like inputs. The de Hoog
16quotient-difference / Pade acceleration of the Bromwich-Fourier series handles them reliably. This
17is verified in the test suite against ``mpmath.invertlaplace`` on both smooth and step-like
18transforms; Talbot is intentionally not offered.
20Algorithm
21---------
22For a real evaluation time ``t`` the inverse is the Bromwich integral, discretized as the
23Fourier series (Crump, 1976)
25.. math::
27 f(t) \approx \frac{e^{\gamma t}}{T}\Big[\tfrac12 \bar f(\gamma)
28 + \sum_{k=1}^{\infty}\operatorname{Re}\big\{\bar f(\gamma + i k\pi/T)\,e^{i k\pi t/T}\big\}\Big],
30with abscissa :math:`\gamma = \alpha - \ln(\text{tol})/(2T)` placed to the right of every
31singularity of :math:`\bar f` (``alpha`` = the largest real part of any singularity; ``0`` for the
32decaying radial kernels, whose poles lie on the non-positive real axis). de Hoog et al. accelerate
33the truncated power series :math:`\sum_k a_k z^k`, :math:`z=e^{i\pi t/T}`,
34:math:`a_k=\bar f(\gamma+ik\pi/T)` (with :math:`a_0` halved), by its continued-fraction (Pade)
35representation whose coefficients come from the quotient-difference (QD) recurrence, plus a tail
36remainder estimate for the last convergent. The continued fraction is evaluated by the standard
37three-term recurrence.
39This file is part of gwtransport which is released under AGPL-3.0 license.
40See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
41"""
43from collections.abc import Callable
45import numpy as np
46import numpy.typing as npt
49def dehoog_inverse(
50 *,
51 f_hat: Callable[[npt.NDArray[np.complexfloating]], npt.NDArray[np.complexfloating]],
52 t: npt.ArrayLike,
53 n_terms: int = 24,
54 scaling: float | None = None,
55 alpha: float = 0.0,
56 tol: float = 1e-9,
57) -> npt.NDArray[np.floating]:
58 r"""Invert a Laplace transform -- optionally a whole batch of them -- at an array of times.
60 The transform ``f_hat`` is sampled once on ``2 * n_terms + 1`` complex Bromwich nodes; the
61 continued-fraction acceleration and its evaluation are vectorized over ``t`` and over any trailing
62 batch axes ``f_hat`` carries, so one call inverts every entry of a propagator matrix in a single
63 QD/continued-fraction pass. Only ``t > 0`` is meaningful; ``t <= 0`` returns ``0`` (the inverse of a
64 one-sided transform vanishes for negative time, and the series is undefined at ``t = 0``).
66 Parameters
67 ----------
68 f_hat : callable
69 Vectorized Laplace transform :math:`\bar f(s)`. Receives a complex ``ndarray`` of shape
70 ``(2 * n_terms + 1,)`` and returns a complex ``ndarray`` whose leading axis is ``2 * n_terms + 1``,
71 optionally with trailing batch axes (e.g. all ``n_quad**2`` entries of a propagator matrix).
72 t : array-like
73 Evaluation times (same clock as the inverse, e.g. flushed volume). 1-D or scalar.
74 n_terms : int, optional
75 Number of acceleration terms ``M``; the transform is evaluated at ``2*M + 1`` nodes.
76 Default 24.
77 scaling : float, optional
78 The half-period ``T`` of the Bromwich-Fourier approximation. The result is accurate for
79 ``0 < t < 2T``, best for ``t <~ T``. Defaults to ``2 * max(t)`` (so every requested time
80 sits in the well-resolved first half).
81 alpha : float, optional
82 Largest real part of any singularity of ``f_hat`` (the Bromwich abscissa is placed to its
83 right). Default 0 (correct for the decaying radial kernels and for ``f_hat`` with poles only
84 at ``s <= 0``, including the ``1/s`` and ``1/s^2`` antiderivative factors).
85 tol : float, optional
86 Target relative accuracy controlling the abscissa offset. Default ``1e-9``.
88 Returns
89 -------
90 ndarray
91 Real inverse ``f(t)``, shape ``t.shape + batch`` where ``batch`` is the trailing shape of
92 ``f_hat``'s output (a scalar ``t`` with no batch yields a 0-d array).
94 Raises
95 ------
96 ValueError
97 If ``f_hat``'s returned array does not have leading axis ``2 * n_terms + 1``.
99 Notes
100 -----
101 The QD recurrence is run on 1-D arrays of decreasing length (the standard rhombus rules of de
102 Hoog, Knight & Stokes 1982); the continued fraction is then evaluated by the three-term
103 recurrence ``A_n = A_{n-1} + d_n z A_{n-2}``, ``B_n = B_{n-1} + d_n z B_{n-2}`` with a final
104 remainder term that sums the truncated tail.
106 References
107 ----------
108 de Hoog, F. R., Knight, J. H., & Stokes, A. N. (1982). An improved method for numerical
109 inversion of Laplace transforms. SIAM Journal on Scientific and Statistical Computing, 3(3),
110 357-366.
111 """
112 t_arr = np.asarray(t, dtype=float)
113 scalar_input = t_arr.ndim == 0
114 t_flat = np.atleast_1d(t_arr)
115 n_t = t_flat.size
117 m = int(n_terms)
118 n_nodes = 2 * m + 1
119 t_max = float(np.max(t_flat)) if t_flat.size else 1.0
120 big_t = float(scaling) if scaling is not None else 2.0 * t_max
121 if big_t <= 0.0:
122 big_t = 1.0 # all t <= 0: the result is masked to 0 below, this only keeps the nodes finite
123 gamma = alpha - np.log(tol) / (2.0 * big_t)
125 # Sample the transform on the Bromwich nodes s_k = gamma + i k pi / T, k = 0 .. 2M. ``f_hat`` may
126 # carry a trailing batch shape; the QD/continued-fraction pass below broadcasts over it, so one call
127 # inverts a whole batch (e.g. every entry of a propagator matrix) at once.
128 k = np.arange(n_nodes)
129 s = gamma + 1j * k * np.pi / big_t
130 a = np.asarray(f_hat(s), dtype=complex)
131 if a.shape[0] != n_nodes:
132 msg = f"f_hat must return leading axis {n_nodes}, got {a.shape}"
133 raise ValueError(msg)
134 batch = a.shape[1:]
135 nb = len(batch)
136 a = a.copy()
137 # A column is *underflow-degenerate* when its transform stays FINITE but decays to the double-precision
138 # floor -- identically zero (a decoupled azimuthal mode with no source) or underflowing toward zero at
139 # the high-frequency nodes (a heavily damped far-field propagator entry: rest / Airy / Riccati). Both
140 # drive the quotient-difference recurrence into 0/0 or x/0 and poison the column with NaN, yet the
141 # physical inverse is ~0. Such columns are detected here (finite transform) and their non-finite output
142 # is zeroed below. A column whose transform OVERFLOWED (inf/nan already in ``a`` -- an ill-scaled kernel)
143 # is a genuine breakdown and is deliberately NOT masked: it must still surface as NaN so a real failure
144 # can never masquerade as a physical zero.
145 finite_transform = np.all(np.isfinite(a), axis=0)
146 a[0] *= 0.5
148 # Quotient-difference algorithm -> continued-fraction coefficients d[0 .. 2M] (per batch entry; the
149 # rhombus rules slice the leading node axis and broadcast over the trailing batch). The degenerate
150 # columns above form 0/0 / x/0 here and in the continued fraction, so every division is evaluated under
151 # errstate; their non-finite output is masked to zero after the assembly.
152 d = np.empty((n_nodes, *batch), dtype=complex)
153 d[0] = a[0]
154 with np.errstate(all="ignore"): # degenerate columns form 0/0, x/0 and over/underflow; masked below
155 q = a[1:] / a[:-1] # q_1^{(i)}, length 2M
156 d[1] = -q[0]
157 e = q[1:] - q[:-1] # e_1^{(i)} = e_0^{(i+1)} + q_1^{(i+1)} - q_1^{(i)}, length 2M-1
158 d[2] = -e[0]
159 for r in range(2, m + 1):
160 q = q[1:-1] * e[1:] / e[:-1] # q_r^{(i)}, length 2M-2r+2
161 d[2 * r - 1] = -q[0]
162 e = e[1:-1] + q[1:] - q[:-1] # e_r^{(i)}, length 2M-2r+1
163 d[2 * r] = -e[0]
165 # Evaluate the continued fraction at z = exp(i pi t / T) by the three-term recurrence
166 # A_n = A_{n-1} + d_n z A_{n-2}, B_n = B_{n-1} + d_n z B_{n-2}, broadcasting time (leading axis) against
167 # the batch (trailing axes). The loop stops one coefficient early (n up to 2M-1) so the de Hoog tail
168 # remainder can replace the bare last coefficient d_{2M} in the final convergent -- the "improved"
169 # estimate of the truncated tail. Keeping the two trailing convergents avoids any division forming it.
170 time_shape = (n_t, *([1] * nb)) # broadcast the leading time axis against the trailing batch
171 z = np.exp(1j * np.pi * t_flat / big_t).reshape(time_shape) # (n_t, 1..1)
172 a_pp = np.zeros((n_t, *batch), dtype=complex) # A_{-1}
173 b_pp = np.ones((n_t, *batch), dtype=complex) # B_{-1}
174 a_p = np.broadcast_to(d[0], (n_t, *batch)) # A_0 (read-only view; only ever read below, never mutated)
175 b_p = np.ones((n_t, *batch), dtype=complex) # B_0
176 with np.errstate(all="ignore"): # degenerate columns propagate Inf/NaN through the CF; masked below
177 for n in range(1, n_nodes - 1):
178 dz = d[n] * z # d[n] (batch) broadcasts against z (time, 1..1) -> (n_t, *batch)
179 a_pp, a_p = a_p, a_p + dz * a_pp
180 b_pp, b_p = b_p, b_p + dz * b_pp
181 # a_p, b_p hold the (2M-1) convergent; a_pp, b_pp the (2M-2) convergent.
182 rem = 0.5 * (1.0 + z * (d[n_nodes - 2] - d[n_nodes - 1]))
183 rem = -rem * (1.0 - np.sqrt(1.0 + z * d[n_nodes - 1] / (rem * rem)))
184 a_final = a_p + rem * a_pp
185 b_final = b_p + rem * b_pp
186 result = (np.exp(gamma * t_flat).reshape(time_shape) / big_t) * np.real(a_final / b_final)
187 result = np.where(t_flat.reshape(time_shape) > 0.0, result, 0.0)
188 # Zero the underflow-degenerate columns (finite transform, non-finite QD output). Overflowed columns
189 # (non-finite transform) are left as NaN so a genuine breakdown surfaces rather than reads as zero.
190 result = np.where(~np.isfinite(result) & finite_transform, 0.0, result)
191 if scalar_input:
192 return result[0]
193 return result.reshape(t_arr.shape + batch)