Coverage for src/gwtransport/_radial_asr_compose.py: 100%

1r"""Flux-resident step response of a constant-Q radial phase in the flushed-volume clock. 

2 

3The single injection/readout primitive the composition engines are built from 

4(:mod:`gwtransport._radial_asr_reuse`): ``G1(S; V') = L^{-1}[ghat_FR(p; V')/p](S)``, the flux-resident 

5step response in flushed volume (FR mode: flux injection at the well, resident detection at volume 

6``V'``). It carries both directions of the transport: 

7 

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 

11 

12 ``f(V') = sum_j cin'_j [G1(S_inj - sigma_j; V') - G1(S_inj - sigma_{j+1}; V')]``. 

13 

142. **Extraction -> arrival.** Each resident parcel at ``V'`` returns to the well with the duality 

15 arrival kernel whose flushed-extraction-volume Laplace transform is the same ``ghat_FR(p; V')`` 

16 (``|Q| h_bar = ghat_FR``). The flow-weighted average over an output (extraction) volume bin 

17 ``[T_i, T_{i+1}]`` is therefore ``[G1(T_{i+1}; V') - G1(T_i; V')]/(T_{i+1}-T_i)``. 

18 

19The flushed-volume FR transfer function is the autonomous form: ``ghat`` depends on the Laplace 

20variable only through ``beta = s/(alpha_L A_0)``, and the S-clock substitution makes 

21``beta = 2 c_geo R p / alpha_L`` (``c_geo = pi b n``, ``p`` conjugate to flushed volume), independent 

22of the flow magnitude. So the kernel is evaluated as ``transfer_function(s = 2 c_geo p, a0 = 1, R)``. 

23 

24This file is part of gwtransport which is released under AGPL-3.0 license. 

25See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details. 

26""" 

27 

28import numpy as np 

29import numpy.typing as npt 

30 

31from gwtransport._radial_asr_dehoog import dehoog_inverse 

32from gwtransport._radial_asr_kernels import transfer_function 

33 

34# Default de Hoog series length and front-anchored scaling margin for the radial-ASR inversions: the FR 

35# step response (here) and the field propagators (_radial_asr_reuse) both import these so the two never 

36# silently desync on de Hoog resolution. 

37_DEHOOG_TERMS = 44 

38_SCALE_MARGIN = 1.3 

39 

40 

41def _fr_step_response( 

42 v_prime: float, 

43 corner_volumes: npt.NDArray[np.floating], 

44 *, 

45 c_geo: float, 

46 r_w: float, 

47 alpha_l: float, 

48 retardation_factor: float, 

49 flow_scale: float, 

50 molecular_diffusivity: float, 

51) -> npt.NDArray[np.floating]: 

52 r"""Flux-resident step response ``G1(S; V') = L^{-1}[ghat_FR(p; V')/p](S)`` at one ``V'``. 

53 

54 For ``D_m = 0`` the flushed-volume Airy kernel depends on the Laplace variable ``p`` only through 

55 ``beta = 2 c_geo R p / alpha_L``, so it is evaluated directly in the flow-free canonical form 

56 (``s = 2 c_geo p``, ``A_0 = 1``) -- exact for arbitrary within-phase variable flow and bit-independent of 

57 ``flow_scale``. For ``D_m > 0`` the kernel depends on ``A_0 = flow_scale / (2 c_geo)`` separately, so the 

58 Laplace variable enters as ``s = flow_scale * p`` and ``flow_scale`` must be the (constant) phase flow 

59 magnitude. 

60 

61 The de Hoog half-period is anchored to the FR arrival-volume mean at ``V'`` 

62 (``mu = R c_geo[(r'+alpha_L)^2 + alpha_L^2 - r_w^2]`` -- the breakthrough front), bounded below by 

63 the requested corner volumes, so the front is resolved even when the output extends far past it. 

64 Corners ``<= 0`` map to ``0`` (no breakthrough yet). 

65 

66 Returns 

67 ------- 

68 ndarray 

69 ``G1(S; V')`` for each ``S`` in ``corner_volumes`` (same shape). 

70 """ 

71 r_p = np.sqrt(r_w**2 + v_prime / c_geo) 

72 mu = retardation_factor * c_geo * ((r_p + alpha_l) ** 2 + alpha_l**2 - r_w**2) 

73 # D_m = 0: the Airy S-clock kernel depends on p only through beta = 2 c_geo R p / alpha_L, so evaluate it 

74 # in the flow-free canonical form (a0 = 1, s = 2 c_geo p) -- routing flow_scale through s and a0 would 

75 # round-trip it and leave ~1-ulp bit-noise the de Hoog QD stage amplifies. D_m > 0: the kernel depends on 

76 # A_0 = flow_scale/(2 c_geo) separately, so use the (constant) phase flow magnitude. 

77 s_mult = 2.0 * c_geo if molecular_diffusivity == 0.0 else flow_scale 

78 a0 = s_mult / (2.0 * c_geo) # 1.0 exactly for D_m = 0 (flow-free), flow_scale/(2 c_geo) for D_m > 0 

79 

80 def f_hat(p: npt.NDArray[np.complexfloating]) -> npt.NDArray[np.complexfloating]: 

81 return ( 

82 transfer_function( 

83 s=s_mult * p, 

84 r=r_p, 

85 r_w=r_w, 

86 alpha_l=alpha_l, 

87 a0=a0, 

88 d_m=molecular_diffusivity, 

89 retardation_factor=retardation_factor, 

90 detect="resident", 

91 ) 

92 / p 

93 ) 

94 

95 cv = np.asarray(corner_volumes, dtype=float) 

96 out = np.zeros_like(cv) 

97 positive = cv > 0.0 

98 if np.any(positive): 

99 scaling = _SCALE_MARGIN * max(mu, float(cv[positive].max())) 

100 out[positive] = dehoog_inverse(f_hat=f_hat, t=cv[positive], n_terms=_DEHOOG_TERMS, scaling=scaling) 

101 return out