Internals

AbstractSeparableSchedule

A pair enumeration with at least one uniform direction. The cells factor into slabs, one per index on the remaining directions, and within a pair of slabs every pair of cells is named by an integer lag along the uniform directions. Both sweeps run over (slab, slab, lag): the direct one reduces the field over each lag, the transform produces every lag of a slab pair at once.

A schedule provides uniform_axes, n_slabs, enumerated_pairs, separable_layout, lag_limits, lag_transport and lag_frame.

UniformLagSchedule(dims, spacing, periodic)

Lag enumeration for a uniform rectilinear grid of dims cells, with constant spacing along each direction, each direction wrapping or not per periodic.

A lag is an integer offset vector naming every pair (cell, cell + lag), so the separation, its direction and its distance bin are properties of the lag alone. Directions are positional and match the trailing axes of the field the schedule is swept over.

spacing is signed, as a descending axis reports it; only its magnitude enters a separation.

RectilinearLagSchedule(uniform, enumerated, axis_order)

Lag enumeration for a rectilinear grid with some uniform directions and some given as coordinate lists. uniform describes the uniform directions, enumerated holds one coordinate vector per remaining direction, and axis_order[k] is the field direction at position k of (uniform..., enumerated...).

Every pair of enumerated indices is a pair of slabs, within which the pairs are named by lags along the uniform directions exactly as on a uniform grid; an enumerated direction never wraps.

ZonalLagSchedule(lats, n_lon, dlon, radius, lon_periodic)

Pair enumeration for a lat-lon grid: one latitude per entry of lats (radians), n_lon cells per row spaced dlon apart in longitude (radians), on a sphere of the given radius.

The latitude axis keeps whatever vector type it arrives as — a range stays a range — so a grid whose axis is a formula is not materialised to carry it.

The field it sweeps is stored (component, longitude, latitude) with components in the local (east, north[, radial]) basis, which is the basis the transport matrices are written in, so no per-pair conversion to ambient coordinates happens. Each latitude row is one slab and longitude is the uniform direction, with the pair frame applied per lag (see FrameTransport).

ScatteredPairs(points, metric)

Pair enumeration for a grid with no structure to exploit: the points themselves, and the metric that measures between them.

A pixelized sphere, a curvilinear mesh, a node set, a grid with no uniform direction — none of these share a separation between many pairs, so there is nothing to hoist and the pairs are enumerated. That is what the unstructured path does, with culling, and this schedule routes to it.

IdentityTransport()
FrameTransport()

How the two ends of a pair are compared. Under IdentityTransport the increment is the plain difference u(J, i+h) - u(I, i) and the lag's displacement is its direction. Under FrameTransport each end is first expressed in the pair's own frame, B·u(J, i+h) - A·u(I, i) with W×W matrices fixed by the lag, and the longitudinal direction is ê₁.

The transport of a flat schedule: the increment is the plain difference and the lag's displacement its direction; see AbstractLagTransport.

The transport of a curved schedule: each end is expressed in the pair's own frame before differencing; see AbstractLagTransport.

uniform_axes(schedule) -> UniformLagSchedule

The uniform directions, in the order they run within a slab.

n_slabs(schedule) -> Int

How many slabs the cells factor into; one for a schedule whose every direction is uniform.

separable_layout(schedule, data, valid, weights) -> (data, valid, weights)

The packed field, its validity and its pair weights with slab I occupying columns (I-1)·N_u + 1 : I·N_u, N_u the cells of a slab, and the uniform directions running column-major within it.

lag_limits(schedule, I, J, r_max) -> NTuple{Du, Int}
lag_limits(schedule, r_max) -> NTuple{Du, Int}

The largest |h_d| along each uniform direction at which a pair of the two slabs can still lie within r_max; the two-argument form bounds it over every slab pair.

uniform_lag_box(schedule) -> Bool

Whether lag_limits(schedule, I, J, r_max) is the same box for every slab pair, so a device launch indexes its work items by division on one box. A schedule that answers false carries one box per pair, and the launch is indexed through the prefix sum of their volumes.

lag_transport(schedule) -> AbstractLagTransport

See IdentityTransport.

NoWeights()

Every pair counts once, as an indexable value: the unweighted sweep is the same kernel with the weight folded away, as AllValid folds the mask away.

AllValid()

Every cell holds a usable datum, as an indexable value.

Distinct from a grid's own mask, which says which cells exist: a cell can exist and still hold nothing, and it is the field that decides. field_validity combines the two. Spelling the complete-field case as a value that answers true keeps one kernel — the test folds away, so a complete field pays nothing, and no separate masked kernel can drift from the unmasked one.

field_validity(u[, cell_mask]) -> BitVector or AllValid

Which cells of u hold a usable datum: every component finite, and — where a cell_mask is given, as a grid carries one — the cell marked as existing. u is (components, cells...) or a Fields.

Returns AllValid when nothing is excluded, so a complete field costs nothing downstream.

BatchLeading(u)

Wrap a velocity/position array that is already stored batch-leading, shape (B, D, N) (batch axis innermost/contiguous — the CPU-optimal SoA layout). CPU batch kernels then run zero-copy. A plain (D, N, B...) array (the default contract, GPU-optimal) is transposed once internally. Zero-cost type tag.

isotropic_kernel(::Val{D}, x)

Angular average of cos(k·r) over the directions of r in D dimensions, at x = kr.

cos on a line, J₀ on a plane and sin(x)/x in a volume. This is the kernel of the isotropic Fourier transform, so a single mode of amplitude A and wavenumber k₀ has S₂(r) = A²[1 - isotropic_kernel(Val(D), k₀ r)] after angular averaging.

AbstractMissingLagPolicy

What a lag-space spectrum does with a lag no pair of held cells names. RefuseMissingLags() throws, since the structure function is undefined there. ZeroDeviationAtMissingLags() sets the autocovariance to zero at that lag, which is what a taper does at the lags beyond its reach.

AbstractFieldShape{D}

Which array-rank pattern a validated (x, u) pair forms, with D the number of velocity components on axis 1 of u. Backends dispatch on this, so D reaches the kernels as a static method parameter.

• PointField: u is (D, N).
• SharedPositionField: u is (D, N, auxiliary...) — one set of positions, many fields.
• VaryingPositionField: x and u both carry matching auxiliary axes.

The coordinate count on axis 1 of x is deliberately not here: it is fixed by the metric's geometry (see HelperFunctions.coordinate_width), so it is read off that geometry once D is static.

CullingPolicy

Whether a calculation may skip pairs that cannot land in a finite bin. Culling never changes the result, so this selects only whether sorting the points is worth paying for.

• AutoCulling — cull when it removes pairs (the default)
• AlwaysCulling — cull whenever pairs can be skipped, and error when they cannot
• NoCulling — sweep every pair

Cull only when the cutoff is small enough to remove pairs; see CullingPolicy.

Cull whenever pairs can be skipped, and error when they cannot; see CullingPolicy.

Never cull; see CullingPolicy.

CellGrid{D, FT, IV, OV}

Uniform cell decomposition of a point set, with points sorted by cell id, so every cell is one contiguous run: cell_run(grid, c) in the permuted order.

Only occupied cells are stored, as the sorted ids cell_ids and their run boundaries run_starts. The cell-id space is unbounded — dims grows as (extent/cutoff)^D, reaching 2.3 GB of ids for a flat 3-D set at r_max/L = 0.003 — while occupancy never exceeds the point count, so storing the occupied cells keeps both the memory and the block enumeration O(N) at any resolution.

The index and offset arrays are type parameters, not Vectors, so one type serves a host build and a device-resident copy.

AbstractSecondAxisSource

What the second axis of a joint histogram bins.

The first axis is always the separation. The second is a choice: the operator's own value, giving the distribution of the structure function at each separation; or the separation's direction, giving the structure function resolved by angle. Both are digitized the same way into the same histogram, so a source is where the binned quantity comes from, not a new kind of output.

InvariantValueAxis()

Bin the operator's own value: the joint distribution of separation and structure-function value.

SeparationAngleAxis(reference_axis)

Bin the angle between the pair's separation and reference_axis, giving S(r, θ) — the structure function resolved by direction, which is what an anisotropic flow needs.

The angle is folded so that a pair and its reverse give the same value, which they must: swapping the two ends flips the separation, and no structure function distinguishes the two. In two dimensions the signed azimuth is taken modulo π, giving [0, π); in three or more the polar angle to the axis is taken from |cos|, giving [0, π/2].

reference_axis need not be normalized. A pair separated along the reference axis has angle exactly 0, and bins are (lo, hi], so bins over this axis start just below zero — axis_bounds does — or an axis-aligned pair falls in no bin.

axis_bounds(source) -> (lo, hi)

The range the source's quantity spans, for callers building bins over it. The operator's value has no bound the source knows.

is_polynomial_operator(sf) -> Bool

Whether sf(δu, r̂) is a homogeneous polynomial in the packed increment δu with coefficients that depend only on r̂, so that Σ_pairs sf(δu, r̂) is moment_contract of the increment moment tensor.

moment_contract(sf, M::SymmetricMoments, r̂, ::Val{V}, ::Val{K})

Σ_pairs sf(δu, r̂) from the increment moment tensor M[i₁, …, i_p] = Σ_pairs Π_k δu[i_k] of the pairs sharing the unit separation r̂. δu is packed as V vector fields of width D = length(r̂) followed by K scalar fields, so M has W = V * D + K components. Operators that are not polynomials in δu have no method.

SymmetricMoments{W, P}(data::SVector)

The rank-P symmetric increment moment tensor over W packed components, stored as its binomial(W + P - 1, P) independent entries M[j₁ ≤ … ≤ j_P] in the order of symmetric_indices.

symmetric_indices(Val(W), Val(P)) -> NTuple{N, NTuple{P, Int}}

The sorted multi-indices j₁ ≤ … ≤ j_P over 1:W in storage order: entry n has n - 1 = Σ_k binomial(j_k + k - 2, k).

order(sf::AbstractStructureFunctionType)

Returns the order of the structure function.

scalar_order(sf) -> Int

The total power of scalar-field increments in sf.

Reading a pair from its other end leaves every vector quantity unchanged — δu, r̂ and n̂ all flip together — and negates each scalar increment, so sf changes sign under that reading exactly when this is odd. Such an operator is evaluated with the pair read from the lower to the upper end along the first coordinate that separates its two points (on a sphere: from south to north, and along one parallel from west to east); where neither end comes first, the two readings are averaged and the moment is zero.

AbstractSquaredDigitizePlan

Maps a squared separation r² to the bin index digitize(r, edges) would give.

Edges are strictly positive, so e_{i-1} < r ≤ e_i iff e_{i-1}² < r² ≤ e_i². The index is approximate; squared_digitize corrects it against the true squared edges.

squared_digitize_plan(edges) -> AbstractSquaredDigitizePlan

Build the r² digitize plan for edges, once per call (never in the pair loop).

squared_digitize(plan, r2) -> Int

Exact digitize(r, edges) computed from r² alone. Out-of-range gives 0 (below) or n_bins + 1 (above), matching HelperFunctions.digitize.