Coverage for src/gwtransport/_radial_asr_compose.py: 100%
39 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +0000
1r"""Composition of constant-Q radial phases into the well observable (single cycle).
3For a single inject-then-extract cycle at one well (no intermediate flow reversal) the extracted
4flux concentration is built grid-free from the per-phase kernels of :mod:`gwtransport._radial_asr_kernels`
5(the KB Sec. 10a pipeline), with everything carried in the flushed-volume clock so that arbitrary
6within-phase variable flow is exact for ``D_m = 0`` (the S-clock convolution theorem):
81. **Injection -> resident profile.** A piecewise-constant injected deviation ``cin'`` (concentration
9 minus background) over injection volume bins ``[sigma_j, sigma_{j+1}]`` leaves, after flushing the
10 total injected volume ``S_inj``, the resident profile
12 ``f(V') = sum_j cin'_j [G1(S_inj - sigma_j; V') - G1(S_inj - sigma_{j+1}; V')]``,
14 where ``G1(S; V') = L^{-1}[ghat_FR(p; V')/p](S)`` is the flux-resident step response in flushed
15 volume (FR mode: flux injection at the well, resident detection at volume ``V'``).
172. **Extraction -> echo.** Each resident parcel at ``V'`` returns to the well with the duality arrival
18 kernel whose flushed-extraction-volume Laplace transform is the same ``ghat_FR(p; V')`` (KB Sec. 7,
19 ``|Q| h_bar = ghat_FR``). The flow-weighted average over an output (extraction) volume bin
20 ``[T_i, T_{i+1}]`` is therefore ``[G1(T_{i+1}; V') - G1(T_i; V')]/(T_{i+1}-T_i)``, and superposing
21 over the resident profile gives the bin-averaged echo. The whole map is the weight matrix
23 ``W_{ij} = int_0^inf [G1(S_inj-sigma_j;V') - G1(S_inj-sigma_{j+1};V')]
24 * [G1(T_{i+1};V') - G1(T_i;V')]/(T_{i+1}-T_i) dV'``, ``cout'_i = sum_j W_{ij} cin'_j``.
26The flushed-volume FR transfer function is the autonomous form: ``ghat`` depends on the Laplace
27variable only through ``beta = s/(alpha_L A_0)``, and the S-clock substitution makes
28``beta = 2 c_geo R p / alpha_L`` (``c_geo = pi b n``, ``p`` conjugate to flushed volume), independent
29of the flow magnitude. So the kernel is evaluated as ``transfer_function(s = 2 c_geo p, a0 = 1, R)``.
31This file is part of gwtransport which is released under AGPL-3.0 license.
32See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
33"""
35import numpy as np
36import numpy.typing as npt
38from gwtransport._radial_asr_dehoog import dehoog_inverse
39from gwtransport._radial_asr_kernels import transfer_function
41# Default de Hoog series length and front-anchored scaling margin for the radial-ASR inversions: the FR
42# step response (here) and the field propagators (_radial_asr_reuse) both import these so the two never
43# silently desync on de Hoog resolution.
44_DEHOOG_TERMS = 44
45_SCALE_MARGIN = 1.3
48def _fr_step_response(
49 v_prime: float,
50 corner_volumes: npt.NDArray[np.floating],
51 *,
52 c_geo: float,
53 r_w: float,
54 alpha_l: float,
55 retardation_factor: float,
56 flow_scale: float,
57 molecular_diffusivity: float,
58) -> npt.NDArray[np.floating]:
59 r"""Flux-resident step response ``G1(S; V') = L^{-1}[ghat_FR(p; V')/p](S)`` at one ``V'``.
61 The flushed-volume Laplace variable ``p`` enters the transfer function as ``s = flow_scale * p``
62 with ``A_0 = flow_scale / (2 c_geo)``. For ``D_m = 0`` the flow magnitude cancels (``ghat``
63 depends only on ``beta = 2 c_geo p / alpha_L``), so the response is the flow-independent S-clock
64 kernel -- exact for arbitrary within-phase variable flow. For ``D_m > 0`` the kernel depends on
65 ``A_0`` separately, so ``flow_scale`` must be the (constant) phase flow magnitude.
67 The de Hoog half-period is anchored to the FR arrival-volume mean at ``V'``
68 (``mu = R c_geo[(r'+alpha_L)^2 + alpha_L^2 - r_w^2]``, KB Sec. 7 -- the breakthrough front),
69 bounded below by the requested corner volumes, so the front is resolved even when the output
70 extends far past it. Corners ``<= 0`` map to ``0`` (no breakthrough yet).
72 Returns
73 -------
74 ndarray
75 ``G1(S; V')`` for each ``S`` in ``corner_volumes`` (same shape).
76 """
77 r_p = np.sqrt(r_w**2 + v_prime / c_geo)
78 mu = retardation_factor * c_geo * ((r_p + alpha_l) ** 2 + alpha_l**2 - r_w**2)
80 def f_hat(p: npt.NDArray[np.complexfloating]) -> npt.NDArray[np.complexfloating]:
81 return (
82 transfer_function(
83 s=flow_scale * p,
84 r=r_p,
85 r_w=r_w,
86 alpha_l=alpha_l,
87 a0=flow_scale / (2.0 * c_geo),
88 d_m=molecular_diffusivity,
89 retardation_factor=retardation_factor,
90 inject="flux",
91 detect="resident",
92 )
93 / p
94 )
96 cv = np.asarray(corner_volumes, dtype=float)
97 out = np.zeros_like(cv)
98 positive = cv > 0.0
99 if np.any(positive):
100 scaling = _SCALE_MARGIN * max(mu, float(cv[positive].max()))
101 out[positive] = dehoog_inverse(f_hat=f_hat, t=cv[positive], n_terms=_DEHOOG_TERMS, scaling=scaling)
102 return out
105def single_cycle_echo_matrix(
106 *,
107 inj_volume_edges: npt.NDArray[np.floating],
108 ext_volume_edges: npt.NDArray[np.floating],
109 c_geo: float,
110 r_w: float,
111 alpha_l: float,
112 inj_flow_scale: float,
113 ext_flow_scale: float,
114 retardation_factor: float = 1.0,
115 molecular_diffusivity: float = 0.0,
116 n_quad: int = 240,
117) -> npt.NDArray[np.floating]:
118 r"""Echo weight matrix ``W`` (``cout' = W @ cin'``) for one inject-then-extract cycle, one disk.
120 The injection builds the resident profile with the FR kernel at the injection flow magnitude and
121 the extraction reads it out at the extraction flow magnitude. For ``D_m = 0`` the flow magnitudes
122 cancel (the S-clock kernel is flow-independent), so arbitrary within-phase variable flow is exact;
123 for ``D_m > 0`` the flow scales must be the (constant) per-phase magnitudes.
125 Parameters
126 ----------
127 inj_volume_edges : ndarray
128 Cumulative flushed-volume edges of the injection bins (length ``n_inj + 1``), increasing.
129 ext_volume_edges : ndarray
130 Cumulative flushed-volume edges of the extraction (output) bins (length ``n_ext + 1``),
131 increasing, measured from the start of extraction.
132 c_geo : float
133 Geometry constant ``pi b n`` so that ``V(r) = c_geo (r^2 - r_w^2)`` (m^3 per m^2).
134 r_w : float
135 Well radius (m).
136 alpha_l : float
137 Longitudinal dispersivity (m).
138 inj_flow_scale, ext_flow_scale : float
139 Flow magnitudes ``|Q|`` (m^3/day) of the injection and extraction phases (only relevant when
140 ``molecular_diffusivity > 0``; ignored to within rounding for ``D_m = 0``).
141 retardation_factor : float, optional
142 Linear retardation ``R >= 1``. Default 1.
143 molecular_diffusivity : float, optional
144 Molecular diffusivity (m^2/day). Default 0 (Airy S-clock; exact for variable flow).
145 n_quad : int, optional
146 Number of Gauss-Legendre nodes for the ``V'`` superposition integral. Default 240.
148 Returns
149 -------
150 ndarray, shape (n_ext, n_inj)
151 ``W`` with ``cout'_i = sum_j W_{ij} cin'_j`` (deviation from background).
152 """
153 sigma = np.asarray(inj_volume_edges, dtype=float) - inj_volume_edges[0] # 0 .. S_inj
154 s_inj = sigma[-1]
155 t_edges = np.asarray(ext_volume_edges, dtype=float) - ext_volume_edges[0] # 0 .. T_end
156 dt = np.diff(t_edges)
158 # V' quadrature window: the resident profile f(V') = G1(S_inj; V') has its front at the retarded solute
159 # radius r_front (where the FR arrival mean equals S_inj) and a dispersive tail of breakthrough width
160 # ~sqrt(alpha_L r_front). A flat 4 S_inj/R margin is Peclet-blind and truncates that tail at low Peclet
161 # (r_front/alpha_L << 1), losing mass (int f dV' < S_inj/R). Mirror the reuse engine's 12-sigma advective
162 # reach plus a 6-sigma molecular reach over the cycle duration -- the same grid that conserves the mass
163 # to ~1e-5. D_m dispatch is inside transfer_function (Airy for D_m = 0, Riccati log-derivative for D_m > 0).
164 r_front = np.sqrt(r_w**2 + s_inj / (retardation_factor * c_geo))
165 total_time = s_inj / inj_flow_scale + t_edges[-1] / ext_flow_scale
166 r_max = r_front + 12.0 * np.sqrt(alpha_l * r_front + alpha_l**2) + 6.0 * np.sqrt(molecular_diffusivity * total_time)
167 v_max = c_geo * (r_max**2 - r_w**2)
168 nodes, weights = np.polynomial.legendre.leggauss(n_quad)
169 v_nodes = 0.5 * v_max * (nodes + 1.0)
170 v_weights = 0.5 * v_max * weights
172 inj_corners = s_inj - sigma # G1 evaluated at S_inj - sigma_j (descending, last is 0)
173 w = np.zeros((len(dt), len(sigma) - 1))
174 for v, vw in zip(v_nodes, v_weights, strict=True):
175 g_inj = _fr_step_response(
176 v,
177 inj_corners,
178 c_geo=c_geo,
179 r_w=r_w,
180 alpha_l=alpha_l,
181 retardation_factor=retardation_factor,
182 flow_scale=inj_flow_scale,
183 molecular_diffusivity=molecular_diffusivity,
184 )
185 g_ext = _fr_step_response(
186 v,
187 t_edges,
188 c_geo=c_geo,
189 r_w=r_w,
190 alpha_l=alpha_l,
191 retardation_factor=retardation_factor,
192 flow_scale=ext_flow_scale,
193 molecular_diffusivity=molecular_diffusivity,
194 )
195 f_contrib = g_inj[:-1] - g_inj[1:] # length n_inj: resident-profile contribution of each cin bin
196 ext_avg = (g_ext[1:] - g_ext[:-1]) / dt # length n_ext: bin-averaged arrival per output bin
197 w += vw * np.outer(ext_avg, f_contrib)
198 # Retardation amplitude: the mobile profile integrates to S_inj/R (sorbed mass is immobile) and the
199 # arrival kernel's CDF plateaus at 1, so the bare readout under-recovers by 1/R. Each extracted
200 # mobile parcel mobilizes its sorbed companion (total solute = R x mobile), so scale the readout by R.
201 return retardation_factor * w