Validation

Each calculation is cross-checked against an independent reference. This page lists those oracles, strongest first, with what each one can and cannot catch.

The ordering matters. An oracle that is exact to round-off will catch a defect that a percent-level physical check never notices, and most of the defects found in this package's history were of exactly that kind — a plausible number that was wrong by a constant factor.

The oracles

Every route with a device kernel — the point kernels, the transform engine on every schedule, the non-uniform FFT route and the tensor kernel — is also run on an A100 against the CPU by gpu/test_cuda_gridded_parity.jl and gpu/runtests.jl, counts exact and sums to round-off; the tolerance policy below says what "round-off" means there.

1. Closed-form Fourier modes

For a single mode $u(x) = A\,\hat{e}\cos(k\cdot x + \varphi)$ averaged over a full period,

\[D_{ab}(r) = A^2\, \hat{e}_a \hat{e}_b \left[1 - \cos(k\cdot r)\right],\]

and every odd-order moment vanishes identically. On a periodic grid whose cell count is a multiple of the mode's period the average over cells is the average over a full period, so this holds to round-off rather than to a sampling tolerance.

Distinct grid harmonics are orthogonal over the cells, so a superposition's cross terms cancel exactly and the single-mode form simply adds. Giving each mode two orthonormal polarisations perpendicular to its wavevector makes the field divergence-free and the polarisation sum the transverse projector, so $\sum_p (\hat{e}_p\cdot\hat{r})^2 = 1 - (\hat{k}\cdot\hat{r})^2$. That makes a prescribed spectrum recoverable mode by mode, and the package reproduces it to a relative 1e-12.

One subtlety is worth stating, because it is easy to get wrong when building such a test field: the two polarisations of a single mode share a wavevector, so they are not orthogonal over the cells. Their cross term is traceless — it leaves the trace exact — but it does project onto $\hat{r}$, so it corrupts the longitudinal component. Offsetting their phases by a quarter turn removes it exactly.

2. Transform against direct sweep

On a uniform grid the structure function of any integer order is computable exactly by transform, because the binomial expansion turns the increment moment into a sum of cross-correlations and the correlation theorem evaluates all lags at once. That gives two independent algorithms for one definition, and they must agree to round-off on the same data.

They are genuinely independent: the sweep visits each lag and reduces over cells, while the transform never forms a lag at all until after the inverse transform. A defect in either shows up immediately as disagreement.

Bounded directions are zero-padded to at least $2n-1$ so that the circular correlation equals the linear one; periodic directions are not padded, because there the circular correlation is exactly the sum wanted. Counts are computed as exact integer arithmetic rather than by rounding an inverse transform.

3. Gridded against the unstructured pair loop

The pair loop is the reference implementation: it enumerates pairs, computes each separation, and bins it. On a bounded grid the separations are plain Euclidean, so the pair loop over the same points must give exactly equal counts and matching sums.

There is one systematic difference worth knowing about, and the gridded path has the better behaviour. All pairs sharing a lag have one exact separation, so the sweep bins them identically, while the pair loop recomputes each pair's coordinate difference and those differ by an ulp. When a bin edge falls on an achievable separation the two disagree about the whole shell. The tests place edges between achievable separations so the comparison is unambiguous.

6. Unit invariance

A physical quantity may not change when the unit of length changes. This is a sharper test than it sounds: the Helmholtz decomposition once multiplied a cumulative integral by $r$, which left the energy identity $D_{rot} + D_{div} = D_{LL} + D_{TT}$ passing to 4e-16 — because the spurious terms cancel in that sum — while producing a 33 % divergent signal on a field with no divergent component, and moving that signal by a factor of 2600 under a change from metres to millimetres.

The lesson generalises: a conservation identity that the defect preserves is not a gate. Assert the invariance the defect actually breaks.

7. Frame invariance

The trace of the second-order tensor is the second-order structure function, in any frame. On a curved manifold that is only true if the tensor's components are transported into a common frame per pair, so trace(D_ab) == S2SF is a direct test of the transport. It is silent on a flat geometry, where a raw coordinate difference and a transported increment coincide — which is why the test asserts it on a sphere as well as on a plane.

10–13. The transforms

Each transform is checked against something with a closed form, because a transform that is wrong by a constant returns a curve of exactly the right shape.

The isotropic transform against an analytic spectrum. A Gaussian correlation $C(r) = σ^2 e^{-r^2/2\ell^2}$ gives $S_2(r) = 2σ^2[1 - e^{-r^2/2\ell^2}]$ and the closed-form density $σ^2 \ell^D e^{-k^2\ell^2/2} / (2π)^{D/2}$. Because it decays, the transform is not truncation-limited and the comparison is pointwise: agreement to 6.6e-05, 2.3e-09 and 9.9e-14 in one, two and three dimensions. A wrong kernel, a wrong solid angle, a wrong (2π)^D or a wrong sign would each break that, so one comparison covers all four.

A discrete spectral line will not do for this. Its correlation does not decay, so the truncated transform is a sinc whose sidelobes fall off like 1/k and are cut off by any finite range — worth 1.9 % in one dimension, and refining the wavenumber grid does not move it. Lines are used only to check that peaks land on the right wavenumbers.

The gridded transform against the field's own spectrum. Over the whole lag space nothing is angularly averaged and nothing is radially binned, so this must agree with transforming the field directly, and it does: 6.4e-17, 3.0e-16, 1.3e-16. That is the standard the isotropic route cannot meet on a grid, where the lattice's separations are biased toward its axes.

The flux relation against a closed-form integral. $\int_0^R J_1(Kr)dr = (1 - J_0(KR))/K$, so a constant advective structure function c must give $Π_K = -(c/2)(1 - J_0(KR))$ exactly. That pins the kernel and the prefactor with nothing left to fit — which matters, because a flux wrong by a constant, or by a sign, still looks like a cascade.

The third-order routes against each other. For an isotropic incompressible flow the advective, S3 and L3 structure functions are tied by SF_A = (1/2r)\,d(r\,S3)/dr and S3 = (1/3r^2)\,d(r^3 L3)/dr, so the three flux relations must agree on any analytic family built through those relations. They do only with the boundary terms their integrations by parts leave at the last separation: the tests show the integrals alone miss by tens of percent and can change sign, and the three routes then agree to 10⁻⁶. The enstrophy routes are gated the same way through SF_{Aω} = -∇² SF_{Au}.

The Helmholtz split by linearity. $D_{rot} + D_{div} = D_{LL} + D_{TT}$ exactly and the transform is linear, so the rotational and divergent spectra must sum to the spectrum of the trace, whatever the field is. This is an identity rather than a physical expectation, which is what makes it a gate.

A note on positive-definiteness

A covariance matrix must be positive semi-definite, and one assembled from a sampled covariance function need not be — interpolating a positive-definite kernel does not preserve positive-definiteness. The error falls as the square of the separation spacing:

So two quite different failures land in the same place, and the check distinguishes them: a covariance function that is invalid is off by the scale of the matrix itself, while one that is valid but under-resolved is off by a discretisation error. A coarse covariance that trips the check is reporting that the representation cannot support a valid matrix.

22–27. The later routes

Weights. The weighted statistic is $\sum w_i w_j v_{ij} / \sum w_i w_j$. The gate is a pair loop written in the test with the same weights: the sweep matches it bit for bit, the transform to 1e-12, and weights of one reproduce the unweighted results exactly, so the weighted code path cannot have drifted from the unweighted one.

The non-uniform FFT route. Points placed exactly on the mode grid make the soft-binned route the periodic gridded transform, which is exact, and the two agree to 1e-9 in one, two and three dimensions with masks, weights and multi-fields. Off the grid the one-dimensional kernel identity is written out in the test — the pair sum with the Dirichlet or Gaussian-tapered kernel — and the route matches it to 1e-9; against the hard-binned pair loop the error falls monotonically as the mode count grows. The route is soft-binned by construction and the tests say so: no finite mode count reproduces the hard bins.

The sorted line. For one-dimensional points the pair loop and the sorted route bin every pair with the same digitize call on the same difference, so the counts are equal bit for bit and the sums to 1e-12; the test also holds the route's cost linear in the points (10⁵ points in well under a second where the pair loop would need 5·10⁹ pairs).

Tensors. The transform's symmetric moment store, expanded to the dense tensor, equals the point tensor on the grid's points at orders 2, 3 and 4 on flat grids and on a lat-lon grid, and the joint tensor over angle marginalises to the tensor and, through its trace, to the joint histogram of S2. Odd ranks read a pair canonically in a fixed frame and take no sign in the sphere's geodesic frame; both are tested.

Fits. The forward models are held to closed forms (the one-dimensional bin integral is elementary; the flux model is exact for a piecewise-constant flux and equals the quadrature of the defining integral to 1e-8), the inversions to round trips of the models' own data, the posterior covariance to $σ²(HᵀH)^{-1}$, and the fitted flux to the flux the J₂ transform of the same S3 recovers between the model's jumps — so the forward model of the fits and the inverse relation of the transforms are the same transform.

Backend parity: the tolerance policy

Two backends computing one statistic on one data set must give equal counts and sums that differ by round-off only. "Round-off" is stated as a relative bound on the largest sum, max |Δsum| / max |sum|, and the bound is what a different summation order can produce:

A test that needs a looser bound than these is testing something other than parity, and says what.

Every parity script fixes its draw, so that a row is reproducible and a change in it is attributable to the code. A fixed draw is never what makes a row pass: each row's verdict is one of the bounds above, which hold for any draw. A boolean that happens to come out true on one draw is not a criterion — the Float32 value-axis row above was such a boolean until 2026-09-16, and it could not distinguish a regression from luck. The one exception recorded: NonuniformFFTs chooses a different kernel evaluation on CUDA than it does on the host, and under it a Float64 transform of ours came out near 1e-5 rather than 1e-15. The package's extension therefore names the kernel and the evaluation explicitly on every backend rather than taking a backend's default, which brings the device transform to 1e-15. Why the default evaluation loses that much is not established here, and the difference belongs to that library's choice of defaults, not to a claim about it.

The printed boundary terms of Pearson et al. (2025)

The third-order flux relations (the J₂ and J₃ routes and the enstrophy route from the velocity's advective structure function) come from integrating the J₁ relation by parts at a finite radius $R$, and the boundary terms are essential: the closed-form gate on the power-law family fails without them at every wavenumber, and the estimate can change sign. Appendix B of the paper prints two boundary terms in a form that is dimensionally inconsistent as read — 3 SF_Luu J₁(Kr) in (B4) where a flux $∼ u^3/L$ requires (3/K) SF_Luu J₁(Kr), and K (dSF_Au/dr) J₁ in (B7) where the enstrophy flux $∼ u^3/L^3$ requires (1/K)(dSF_Au/dr) J₁. The consistent forms are the ones that close on the power-law family and make the three energy routes and the three enstrophy routes agree to 1e-6 on analytic isotropic families; they are what the package ships. The paper has no erratum at the time of writing; the difference may be a rendering of a small /K.

What is checked statistically, and why it is not a gate

The isotropic relation

\[D_{LL}(r) = \int_0^\infty E(k)\, f(kr)\,dk, \qquad f(x) = 4\left[\tfrac{1}{3} - \frac{\sin x - x\cos x}{x^3}\right]\]

holds for a three-dimensional isotropic field. The kernel itself is confirmed numerically here: the transverse-projector direction average equals $\tfrac{1}{2} f(k_0 r)$ at every $r$ and every shell radius, to within the spherical quadrature's own error of about 4e-6.

The relation is nonetheless not used as a test assertion, because neither available construction makes it tight:

• A field built from cubic-grid harmonics in a thin shell is 0.4–5 % anisotropic, and that does not improve with shell radius — a shell of 3722 modes still shows 2 %.
• A field built from spherical-quadrature directions evaluated on scattered points converges only as the sampling error allows: 5–8 % at 16 000 points.

In both cases the residual is a property of the mode set and the sampling, not of this package, so an assertion at that tolerance would mostly be testing the test. The exact multi-mode oracle in section 1 covers the same code path to round-off and is used instead.

The same reasoning applies to the inertial-range laws. A synthetic Gaussian field carries no energy flux, so its third-order moments are consistent with zero and the four-fifths law has nothing to recover. What is tested is the inversion itself: given a moment that obeys a law exactly, the corresponding routine returns the prescribed constant, and the routines are not interchangeable — applying the four-fifths law to the scalar moment is off by exactly 5/3, which is asserted so that the two cannot be confused.

Reproducing

Each oracle lives in a targeted test file that can be run on its own:

julia --project=test test/test_known_truth.jl
julia --project=test test/test_gridded_fft.jl
julia --project=test test/test_spherical_geometry.jl
julia --project=test test/test_gridded_weights.jl
julia --project=test test/test_scattered_modes.jl
julia --project=test test/test_sorted_line.jl
julia --project=test test/test_tensor_khm.jl
julia --project=test test/test_fits.jl

The device parity runs need a GPU allocation: sbatch gpu/run_cuda_gridded_parity.sh writes its table to gpu/benchmark_results/, and julia --project=gpu gpu/runtests.jl runs every CUDA suite.