API Reference

DivergentSecondOrderStructureFunctionType()

Helmholtz-derived 2D divergent second-order component. This is not a pairwise operator; compute it from binned L2SF/T2SF with helmholtz_decompose_2d.

FullVectorStructureFunctionType{NF}

Generic full-vector norm-power operator, $||δu||^NF$. This is not used for S3SF, whose conventional definition is ThirdOrderStructureFunctionType.

HelmholtzDecomposition2DType()

Derived quantity describing the 2D isotropic Helmholtz decomposition into rotational and divergent second-order components.

LongitudinalTransverseComponentThirdOrderStructureFunctionType()

Per-component variant of L1T2SF, $δu_L * ||δu_t||² / (D - 1)$.

ProjectedStructureFunctionType{NL,NT}

Parametric type representing structure function operators with longitudinal (NL) and transverse (NT) contributions.

(sf::ProjectedStructureFunctionType{NL,NT})(δu, r̂)

Compute the structure function kernel for longitudinal/transverse components.

• NL : power of longitudinal component δu_l
• NT : power of transverse component ||δu_t||

RotationalSecondOrderStructureFunctionType()

Helmholtz-derived 2D rotational second-order component. This is not a pairwise operator; compute it from binned L2SF/T2SF with helmholtz_decompose_2d.

SecondOrderStructureFunctionType()

Full-vector second-order structure function, $S2SF = ||δu||²$.

ThirdOrderStructureFunctionType()

Third-order scalar flux structure function, $S3SF = δu_L * ||δu||² = L3SF + L1T2SF$. It is intentionally not $||δu||³$.

TransverseComponentSecondOrderStructureFunctionType()

Per-component transverse second-order structure function, $||δu_t||² / (D - 1)$. This is distinct from T2SF, which stores the total transverse energy.

_sf_raw(sf, δu, dx, r2)

Evaluate sf from raw pair geometry: separation dx, its squared length r2 = dx⋅dx, and δu.

Identical to sf(δu, dx/√r2) for every operator. The second-order specializations below are polynomials in p = δu⋅dx and ‖δu‖² over a power of r², so they need no sqrt; odd orders and odd transverse orders need r̂ and fall through to the generic method.

get_structure_function_type(order::Int, mode::Symbol)

Map an integer order and a mode symbol (e.g., :longitudinal, :transverse, :scalar/total) to a specific operator instance.

norm2(x)

Compute the sum of squares of elements of x. Faster than norm(x)^2 for small vectors, used for transverse components.

order(sf::AbstractStructureFunctionType)

Returns the order of the structure function.

GPUSFWorkspace(backend, distance_bins; kind=:sf1d)

Reusable GPU device histogram workspace. Requires the GPUExt extension (using KernelAbstractions). Pass to gpu_calculate_structure_function(!) or slice drivers to avoid per-call device allocation.

See also reset_histogram!, release!.

Serial pair-loop accumulation into native $(6, n_bins)$ buffers (no allocation).

Serial pair-loop accumulation into native $(6, n_bins, n_val)$ buffers (no allocation).

_assert_counts_representable(CT, n_points)

Throw unless the worst-case pair count n_points*(n_points-1)÷2 fits in CT.

Every pair can land in one bin, so that product is the only safe bound. UInt32 saturates at N = 92682, past which the counter wraps silently.

_balanced_index_chunks(N, k) -> Vector of k index-lists

Split 1:N into k balanced outer-index lists for the triangular pair loop (work ∝ N-i). Round-robin assignment (i ≡ w (mod k)) gives each chunk a mix of cheap/expensive indices.

_bin_average!(out, sums, counts)
_bin_average(sums, counts)

Per-bin mean sums ./ counts with the empty-bin guard count == 0 → NaN. The cast uses eltype(out) (so Float32 stays Float32, Float64 stays Float64). Elementwise: works for 1D vectors, 2D matrices, and batched (n_bins, batch...) arrays whose sums/counts share a shape. The allocating form returns a fresh array of eltype(sums). This is the single canonical averaging used by _finalize.

Bytes in one full-width (sums, counts) accumulator of shape dims.

_bl_accum_pool(ws, make_accum, widths) -> Vector

One accumulator per task, of the given batch widths, drawn from the workspace when there is one and allocated fresh otherwise. Always built outside the parallel region, so tasks only read their slot.

Accumulator axes after the leading batch axis, per workspace kind.

_bl_add_permuted!(dest, src_bl, perm)

Add the batch-leading accumulator src_bl into dest through the permutation perm.

! entry points across the package accumulate into the caller's buffers; zeroing belongs to the non-mutating wrappers. PermutedDimsArray is a lazy view, so this fuses the permute with the add and allocates nothing.

_bl_batch_chunk_count(accum_bytes, B, n_tasks) -> Int

How many chunks to split the batch axis into, given the full-width accumulator size in bytes.

Every task holds an accumulator, so the live footprint is n_tasks * accum_bytes — that is what turns into GC and page-fault time, and it is what a batch chunk divides. Splitting b costs pair geometry, recomputed once per batch chunk, so the right answer is the SMALLEST split that brings the footprint under budget. Kernels whose footprint already fits (the 1D batch accumulator is tens of KiB) keep Bc == 1 and never pay geometry twice; the single-pass 2D accumulator at large B and many threads reaches hundreds of MiB and is split until it fits.

Contiguous batch-axis chunks; _bl_batch_chunk_count picks n.

_bl_executor(backend) -> executor

The batch-leading executor a local backend runs with, as a value. Backends that compose over a local inner backend (MPI) look it up here; extensions add methods for the backends they provide.

_bl_n_tasks(backend) -> Int

How many tasks a batch call on backend splits into, and therefore how many accumulators a CPUSFWorkspace must hold. A threaded backend falls back to the serial executor when OhMyThreads is not loaded, so the count follows the executor that will actually run.

_bl_partition(B, n_tasks, accum_bytes) -> (bchunks, n_ichunks)

The batch-axis chunks and the number of outer-index chunks each is split into. Shared by the threaded executor and by CPUSFWorkspace so a workspace cannot be sized for a partition the executor will not use.

The full-width reduction accumulator, from the workspace when there is one.

_pair_dims(metric, D) -> (Val{W}, Val{D})

Coordinate and velocity widths as type parameters, for the kernels that take x/u as Arrays (whose axis-1 lengths are values, not type parameters). D is already restricted to 2 or 3 by _validate_spatial_dimension, so dispatching on those two makes both results concrete; W comes from the metric's geometry, which is what defines it.

_partial_2d_sums_counts(inner, sf_type, x_vecs, u_vecs, distance_bins, value_bins, ilist; kwargs...)

Partial 2D-joint sums/counts over an explicit outer-index list ilist, the 2D analogue of _partial_sums_counts. Euclidean D ∈ {2,3} takes the SIMD compute/scatter kernel; other metrics fall back to the scalar per-i kernel. Returns (sums, counts).

_partial_single_pass_1d(x, u, distance_bins, ilist; distance_metric, count_eltype)

Six-invariant partial sums/counts over an explicit outer-index list, for one distributed worker or MPI rank. Euclidean D ∈ {2,3} takes the SIMD kernel; other metrics use the scalar loop. The accumulator element type comes from the inputs.

_partial_single_pass_2d(x, u, distance_bins, value_bins, ilist; distance_metric, count_eltype)

Six-invariant 2D joint partial sums/counts over an explicit outer-index list, for one distributed worker or MPI rank. The accumulator element type comes from the inputs.

_partial_sums_counts(inner, sf_type, x_vecs, u_vecs, distance_bins, ilist; kwargs...)

Partial 1D sums/counts over an explicit outer-index list ilist (each i contributes pairs (i, j>i)). Used by the distributed driver to give each worker a balanced share; inner selects how the worker computes its share locally. This generic method runs SERIALLY for any backend; the OhMyThreads extension adds a ::CB.AbstractThreadedBackend method that threads over ilist (enabling hybrid distributed+threaded). Returns a StructureFunctionSumsAndCounts.

_pf_2d_simd_pairs!(sums2d, counts2d, sf, xc, uc, dist_be, val_be, ::Val{D}, distbuf, valbuf, irange)

2D-joint point-field SIMD compute/scatter kernel over outer indices irange: @simd over j>i computes distance + SF value into buffers (no scatter ⇒ vectorizes), then a scalar loop digitizes both axes and scatters into the (dist,value) cell. Loop lives in this one kernel (called once per range) so the @simd vectorizes; shared by serial + threaded.

_pf_2d_simd_partial!(sums2d, counts2d, sf, x_vecs, u_vecs, dist_be, val_be, ::Val{D}, ilist)

Run _pf_2d_simd_pairs! over an explicit outer-index list, materializing the contiguous component vectors and scratch buffers this worker needs. Shared by the serial, distributed and MPI drivers, whose inputs arrive as strided views.

_pf_simd_pairs!(output, counts, sf, xc, uc, plan, ::Val{D}, r2buf, valbuf, idxbuf, irange)

Accumulate pairs (i, j>i) for i in irange into output/counts.

The @simd half writes r², the SF value, and the approximate bin index to buffers; the scalar half corrects the index and scatters straight into output/counts, skipping out-of-range bins.

The i-loop and the inner @simd must stay in this function body; factoring the inner loop into a per-i helper stops it vectorizing.

_pf_simd_partial!(sums, counts, sf, x_vecs, u_vecs, dist_be, ::Val{D}, ilist)

Run _pf_simd_pairs! over an explicit outer-index list, materializing the contiguous component vectors and scratch buffers this worker needs. Shared by the distributed, MPI and hybrid drivers, whose inputs arrive as strided views.

_pf_simd_run!(output, counts, sf, x_vecs, u_vecs, dist_be, ::Val{D}) -> mutates buffers

Point-field 1D (Euclidean) via the SIMD compute/scatter split. Materializes contiguous per-component vectors (so consecutive j are unit-stride → packed loads), then for each i: @simd over j>i computes distance + SF value into buffers (no scatter ⇒ vectorizes), and a short scalar loop digitizes + scatters into the histogram. Val{D} keeps it type-stable.

_pf_sp_simd_pairs!(sums, counts, xc, uc, dist_be, ::Val{D}, distbuf, duLbuf, dn2buf, irange)

Single-pass (6 invariants) point-field SIMD compute/scatter kernel over outer indices irange. For each i: @simd over j>i computes distance, du_L = du·r̂, and |du|² into buffers (contiguous components ⇒ packed loads, no scatter ⇒ vectorizes), then a scalar loop digitizes and scatters the 6 invariants. Like _pf_simd_pairs!, the loop must live in this one kernel (not a per-i helper) for the @simd to vectorize. Shared by serial + threaded.

_single_pass_collection_1d(sums, counts, distance_bins, ::Type{OT})

Wrap the stacked 1D single-pass (sums, counts) into a NamedTuple keyed by invariant (S2, L2, T2, S3, L3, L1T2), each value a single-operator result of representation OT (default the averaged StructureFunction; pass StructureFunctionSumsAndCounts for raw). Entries are zero-copy views into the stacked accumulator and share the (identical-by-construction) counts row. For point-field input (stacked with the Helmholtz rows) a :helmholtz => HelmholtzDecomposition2D entry is appended.

_single_pass_collection_2d(sums, counts, distance_bins, value_bins, ::Type{OT})

Wrap the stacked 2D single-pass (sums, counts) (shape (6, n_dist, n_val, aux...)) into a NamedTuple keyed by invariant, each value a StructureFunction2DSumsAndCounts view into the stacked accumulator. Unlike 1D, the per-cell counts genuinely differ per invariant (each invariant's value lands in a different value-bin), so counts are taken per-invariant. The 2D joint histogram has no averaged representation, so OT must be StructureFunction2DSumsAndCounts.

_sp1d_derive_rows!(sums, counts)

Fill the two derived invariant rows and replicate the shared count of a (SINGLE_PASS_N, n_bins) accumulator.

T2 = S2 - L2 and L1T2 = S3 - L3 hold for every pair, and a bin is a sum, so the pair loops store only rows 1, 2, 4, 5 — four stores per pair instead of six — and these two are differenced once per call. The count is identical across all six rows, so it is accumulated into row 1 and broadcast here.

Uses =, never +=, so it is idempotent: correct whether a kernel runs once or many times over the same buffer, and correct under threaded/partial reduction because the derivation is linear (Σ(S2-L2) = ΣS2 - ΣL2). That is why it belongs at the end of each pair kernel rather than at the callers' assembly points, of which there are more than twenty.

_sp1d_pairs!(sums, counts, x, u, dist_be, ::Val{W}, ::Val{D}, metric, n_bins, ilist)

Six-invariant scalar pair loop over the outer indices ilist, for geometries with no SIMD fast path. W is the coordinate width and D the velocity width; both are type parameters so the SVectors stay concrete inside the loop.

Accumulate single-pass 2D pairs for outer indices ilist into the caller's sums/counts.

Fill the interleaved accumulator from pairs whose outer index is in ilist. Euclidean D ∈ {2,3} takes the SIMD compute/scatter split; other metrics or dimensions take the scalar loop. The Val{D} branch is a function barrier: D must be a type parameter inside the loop, or SVector{D} builds its type per point.

_sp2d_histogram(OT, n_bins, n_val) -> Array{OT,4}

The single-pass 2D accumulator, laid out (sum|count, invariant, value_bin, distance_bin).

Each pair writes all six invariants at ONE distance bin but six different value bins, so putting the value axis inside the distance axis keeps a pair's six updates inside one distance slab, and interleaving sum with count puts each invariant's two updates on one cache line: 6 lines touched per pair instead of 12, and the cost stops scaling with histogram size.

_sp2d_pairs!(h, x, u, dist_be, plan, value_bins, ::Val{D}, metric, n_bins, n_val, ilist)

Single-pass 2D scalar pair loop over the outer indices ilist, for non-Euclidean metrics. Specialized on the spatial dimension D so the SVectors are concrete.

Scatter the six invariants of one pair into their cells of the interleaved accumulator.

_sp2d_simd_pairs!(h, xc, uc, plan, value_bins, ::Val{D}, keybuf, duLbuf, dn2buf, idxbuf, n_val, irange)

Single-pass 2D point-field SIMD compute/scatter kernel over outer indices irange, the 2D analogue of _pf_sp_simd_pairs!. The @simd half computes distance, du_L and |du|² into buffers; the scalar half derives the six invariants from those two scalars and scatters each into its own (distance, value) cell. Shared by serial + threaded.

_sp2d_simd_partial!(h, x, u, dist_be, value_bins, ::Val{D}, n_val, ilist)

Run _sp2d_simd_pairs! over an explicit outer-index list, with this worker's buffers.

Add the interleaved accumulator into the caller's (6, n_bins, n_val) sums/counts.

Run _pf_sp_simd_pairs! over an explicit outer-index list, with this worker's buffers.

_tensor_bin_average(sums, counts, ::Val{P})

Tensor analogue of _bin_average: counts (indexed by (bin, aux...)) broadcasts over the P leading component axes of sums (shape (D×P..., n_bins, aux...)). Same empty-bin guard (count == 0 → NaN) and eltype preservation. Used by _finalize to average a tensor result.

_validate_array_shape(x, u, distance_metric) -> AbstractFieldShape

Axis-1 of u is the velocity dimension D; axis-1 of x is however many coordinates the metric's geometry needs to locate a point, which is not always D. On a sphere a point takes two coordinates whether or not the velocity carries a third, radial, component — the shell radius belongs to the geometry, not to each point — so x is (2, N) while u may be (3, N).

Throw unless the staged coordinate width matches what geom needs.

The batch entry points take raw (D, N, B…) arrays rather than going through _validate_array_shape, so this is where a mismatched x is caught — before the kernels load it under @inbounds.

Throw unless ws holds accumulators of this kind and shape. Checked before they are used.

Throw unless ws was built for this input shape. Checked before the transpose buffers are used.

Accumulator axes after the batch axis, as the workspace was built for.

append_helmholtz_rotational_divergent_rows(sums, counts, distance_bins)

Append Helmholtz-derived rotational/divergent rows to native six-row single-pass sums and counts. Rows 1 through 6 are copied unchanged; rows 7 and 8 are rotational and divergent second-order components.

auxiliary_joint2d!(sums, counts, sf_type, x, u, distance_bins, value_bins)

Single-type joint 2D batch (serial); output (n_dist, n_val, batch…).

auxiliary_shared_positions!(sums, counts, x_mat, u_batch, sf_type, distance_bins; strip_width=32)

Fixed geometry batch: x is (N_dims, N), u has trailing batch dims.

auxiliary_structure_function!(sums, counts, sf_type, x, u, distance_bins; strip_width=32)

Dispatch CPU batch calculation based on x rank.

auxiliary_varying_positions!(sums, counts, x_batch, u_batch, sf_type, distance_bins)

Varying geometry batch: x and u have matching trailing batch dims.

calculate_structure_function_2d_batch!(sums, counts, sf_type, x, u, distance_bins, value_bins; backend=..., ...)

Batch 2D joint histograms over (N_dims, N_points, T); outputs (n_dist, n_val, T).

calculate_structure_function_batch!(sums, counts, sf_type, x, u, distance_bins; backend=..., workspace=nothing, ...)

Batch structure functions over the third dimension of matrix inputs (N_dims, N_points, T). Host sums, counts must have shape (NB, T) where NB = length(distance_bins) - 1.

GPU: fully implemented via GPUBackend when KernelAbstractions is loaded. CPU backends: not yet implemented (use a loop over t with calculate_structure_function!).

calculate_structure_functions_single_pass!(sums, counts, x, u, distance_bins; backend=CB.SerialBackend(), kwargs...)

Accumulate into pre-allocated $(6, n_bins)$ buffers using the requested execution backend.

calculate_structure_functions_single_pass(x, u, distance_bins; backend=CB.AutoBackend(),
                                          output_type=StructureFunction, kwargs...)

Compute the six native invariant structure functions (S2, L2, T2, S3, L3, L1T2) in one pair pass, returned as a NamedTuple keyed by invariant. Each entry is a single-operator result of the requested output_type (default the averaged StructureFunction; pass StructureFunctionSumsAndCounts for the raw sums+counts). For point-field input a :helmholtz entry (a HelmholtzDecomposition2D) is included.

calculate_structure_functions_single_pass_2d(x, u, distance_bins, value_bins; backend=CB.AutoBackend(),
                                             output_type=StructureFunction2DSumsAndCounts, kwargs...)

Compute the six invariant 2D joint structure-function histograms in one pass, returned as a NamedTuple keyed by invariant (S2, L2, T2, S3, L3, L1T2). Each entry is a StructureFunction2DSumsAndCounts view into the stacked accumulator (the 2D joint histogram has no averaged form, so output_type must be StructureFunction2DSumsAndCounts).

calculate_structure_functions_single_pass_2d_batch!(sums, counts, x, u, distance_bins, value_bins; backend=..., ...)

Batch six invariant distance × value joint histograms over (N_dims, N_points, T); outputs (6, NB, n_val, T). Pass shared bin types or NTuple{6,...}; use Tuple(v...) if you have a length-6 vector of bin objects.

calculate_structure_functions_single_pass_batch!(sums, counts, x, u, distance_bins; backend=..., ...)

Batch six invariant 1D distance histograms over (N_dims, N_points, T); outputs (6, NB, T).

Loop-over-slice gold reference for batch parity.

gpu_calculate_structure_function(...)

GPU-accelerated structure function calculation. Requires loading KernelAbstractions.jl to activate the GPUExt extension. The backend can be KernelAbstractions.CPU() (for testing parity) or any GPU backend like CUDABackend() from CUDA.jl.

Device histogram buffers are UInt32; count_eltype (default UInt32) selects the host count type after download. Pass workspace=GPUSFWorkspace(...) to reuse device histogram buffers across repeated calls (see GPUSFWorkspace).

This stub exists so the extension can legally extend this function.

gpu_calculate_structure_function!(output_sums, output_counts, ...)

In-place GPU structure function reduction. Requires the GPUExt extension. Accumulates into caller-owned output_sums and output_counts (same contract as serial_calculate_structure_function! / threaded_calculate_structure_function!). Reuses host buffers only; pass workspace=GPUSFWorkspace(...) to reuse device histogram buffers across calls.

gpu_calculate_structure_function_2d(sf_type, backend, x_mat, u_mat, distance_bins, value_bins; kwargs...)

GPU 2D joint structure function (distance × SF value histogram) for one sf_type. Requires loading KernelAbstractions.jl to activate the GPUExt extension. Device counts are UInt32; count_eltype selects the host matrix type after download.

GPU 2D joint slice batch; outputs (n_dist, n_val, T). Requires GPUExt.

gpu_calculate_structure_function_batch!(sums, counts, sf_type, backend, x, u, distance_bins; workspace=nothing, ...)

GPU slice batch over (N_dims, N_points, T); host outputs (NB, T). Requires GPUExt.

gpu_calculate_structure_functions_single_pass_2d(backend, x, u, distance_bins, value_bins; kwargs...)

Six invariant native distance × value joint histograms on a KernelAbstractions backend. Requires loading KernelAbstractions.jl to activate the GPUExt extension. Device counts are UInt32; count_eltype selects the host array type after download.

gpu_calculate_structure_functions_single_pass_2d!(sums, counts, backend, x, u, distance_bins, value_bins; kwargs...)

In-place GPU 2D single-pass accumulation. Requires the GPUExt extension.

GPU single-pass 2D slice batch; outputs (6, NB, n_val, T). Requires GPUExt.

GPU single-pass slice batch; outputs (6, NB, T). Requires GPUExt.

gpu_device_caps(backend) -> GPUDeviceCaps

Capabilities of backend. The default is deliberately the universal floor: a backend with no override behaves exactly as the package did before device querying existed, so an unknown or future backend degrades to "correct and portable" rather than to "assumes an A100". StructureFunctionsCUDAExt overrides this with the real device attributes.

gpu_dynamic_smem_budget(caps; target_blocks_per_sm = 2) -> Int

Bytes a dynamic-shared kernel should use per block. Expressed in device-relative terms — the opt-in ceiling, and the per-SM pool divided by an occupancy target — so the same rule sizes correctly on any part rather than encoding one device's byte count. target_blocks_per_sm is the only free parameter and is dimensionless.

Try the CUDA fast 1D launch (N-body broadcast + privatized shared histogram). Returns true if handled, false to fall back to the portable KA tiled kernel. Overridden by StructureFunctionsCUDAExt.

Try the CUDA fast 2D launch (N-body broadcast + privatized histogram, dynamic shared for large single-pass). Returns true if handled, false to fall back to the portable KA tiled kernel. Overridden by StructureFunctionsCUDAExt.

gpu_static_smem_budget(caps) -> Int

Bytes a portable (static @localmem) kernel may use on this device: the static cap, further limited if the device offers less than it.

helmholtz_decompose_2d(distance_bins, sums, counts)

Run the 2D isotropic Helmholtz decomposition from native single-pass sums/counts. Rows 2 and 3 must contain L2SF and T2SF.

helmholtz_decompose_2d(distance_bins, L2_sums, L2_counts, T2_sums, T2_counts)

Run the 2D isotropic Helmholtz decomposition using the trapezoidal rule over binned longitudinal/transverse second-order structure functions. This implements the cumulative integral equations described by Lindborg (JFM 2015) and Bühler, Callies, and Ferrari (JFM 2014).

256-aligned joint compile width (capped at max). Implemented in GPU extension.

Exact joint histogram cell count n_dist × n_val. Implemented in GPU extension.

Compile-time joint 2D shared-histogram width SF_GPU_MAX_2D_HIST (4096). Implemented in GPU extension.

marginalize_sp2d_then_append_helmholtz_rows(sums_6, counts_6, distance_bins)

Marginalize six native invariant 2D joint histograms to 1D, then append rotational/divergent rows via append_helmholtz_rotational_divergent_rows. Returns $(sums_8, counts_8)$.

minmax_i(i, x_vecs, distance_metric)

Calculate the min and max distances from point i to all other points j != i.

Release device buffers held by a GPUSFWorkspace (optional explicit free).

Zero device histogram buffers in a GPUSFWorkspace before the next launch.

Zero every accumulator held by a CPUSFWorkspace; the drivers do this per call.

resolve_auto_backend(shape, threaded_available, distributed_available) -> AbstractExecutionBackend

What AutoBackend resolves to. threaded_available/distributed_available are zero-argument predicates because each entry family probes a different dispatch.

serial_calculate_structure_functions_single_pass!(sums, counts, x, u, distance_bins)

Six-invariant-type single-pass 1D batch (serial).

serial_calculate_structure_functions_single_pass(x, u, distance_bins, sums, counts)

Zero $sums$/$counts$ then accumulate six invariant native 1D structure functions on one thread. For allocation-free reuse, prefer calculate_structure_functions_single_pass!.

serial_calculate_structure_functions_single_pass_2d!(sums, counts, x, u, distance_bins, value_bins)

Six-invariant-type SP2D batch (serial); output (6, n_dist, n_val, batch…).

Velocity components on axis 1 of u.

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 rather than being restated as a shape parameter.

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 (measured <1% of one pass). Zero-cost type tag.

CPUSFWorkspace{kind}(x, u, distance_bins[, value_bins]; backend = SerialBackend(), count_eltype = UInt32)

Reusable CPU scratch for the batch drivers: the batch-leading working copies of x/u, one accumulator per task, and the full-width reduction accumulator. Pass to a batch entry point as workspace = ... to make repeated calls allocation-free.

Without it a batch call allocates one accumulator per task — tens of MiB at large B — on every call. That does not slow the fastest call, but the resulting GC pauses land on some calls, so the typical call is far slower than the best one and the spread is wide. Reusing the buffers removes the pauses, not the peak.

kind matches GPUSFWorkspace: :sf1d, :joint2d, :single_pass or :single_pass_2d. It is a type parameter, so the accumulator rank is fixed at construction and every field is concretely typed.

Built from the same arguments as the call it serves — including backend, which fixes how many accumulators it holds — and checked against them on use: a workspace built for a different shape, element type or backend is a hard error, never a silent reallocation. Composes with BatchLeading — when the input is already batch-leading the transpose buffers have length zero and are never touched.

Not thread-safe: one workspace serves one call at a time, exactly like GPUSFWorkspace.

GPUDeviceCaps

What a GPU backend can actually offer, queried rather than assumed, so shared-memory strategy is chosen per device instead of per hardcoded constant. Shared memory per block differs by an order of magnitude across parts a user may run on (V100 96 KiB, L40S 100 KiB, A100 163 KiB, later parts more), and the right accumulation strategy differs with it.

smem_per_block is the opt-in maximum, reachable only by dynamic shared memory that a kernel explicitly requests; smem_per_sm bounds how many blocks stay resident and is what makes "use every byte" the wrong default.

Static shared memory is capped far lower and independently — see [GPU_SMEM_STATIC_MAX]. KernelAbstractions' @localmem lowers to a static allocation (CuStaticSharedArray on CUDA) and its launch path passes no shmem, so a kernel written once for every backend is bound by the static cap. Reaching smem_per_block therefore takes a backend-specialized kernel that declares dynamic shared memory and opts in at launch — which is what the vendor fast paths do, with the portable kernel remaining as the correctness fallback for backends that have none.

Position/velocity input a batch driver accepts: a plain (D, N, B…) array, or a BatchLeading wrapper around a (B, D, N) one.

Largest static shared allocation a block may declare. Measured, not assumed: on an A100 a static @localmem of 48 KiB compiles and 64 KiB fails ptxas, while dynamic shared reaches the full 163 KiB opt-in. This is an architectural limit rather than a per-device one, so exceeding a device's opt-in maximum is a separate check.

Shared memory every CUDA-class device provides without opting in.

Value-axis specification for calculate_structure_functions_single_pass_2d.

SINGLE_PASS_OPERATORS

The six native single-pass invariants in stacked-row order, keyed by short name. Used to label each row of the single-pass (sums, counts) accumulator when building the result collection.

AbstractStructureFunction

Abstract base type for all structure function result containers.

HelmholtzDecomposition2D(distance_bins, rotational_sums, rotational_counts,
                         divergent_sums, divergent_counts,
                         longitudinal_values, transverse_values)

Result of the 2D isotropic Helmholtz decomposition computed from binned longitudinal/transverse second-order structure functions. Rotational and divergent components are derived quantities, not pairwise operators.

StructureFunction{FT, OT, BT, VT}

Metadata-rich result object for structure function calculations.

• operator: The specific operator used (e.g., L2SF).
• distance: The coordinate container (can be bin midpoints, edges, or point distances).
• values: The computed structure function values.

StructureFunction2DSumsAndCounts{FT, OT, BT, VT, MT}

2D Joint-Probability intermediate result container containing raw sums and counts matrices. Useful for analyzing the PDF of structure function values across separation distance bins.

• operator: The specific operator used (e.g., L2SF).
• distance_bins: 1D container of distance bin edges.
• value_bins: 1D container of structure function value bin edges.
• sums: 2D matrix of accumulated exact values of shape (Ndistancebins, Nvaluebins).
• counts: 2D matrix of accumulated contribution counts of shape (Ndistancebins, Nvaluebins).

StructureFunctionSumsAndCounts{FT, OT, BT, VT}

Intermediate result object containing raw sums and per-bin pair counts (integer). Useful for aggregating measurements before final averaging.

StructureFunctionTensor(order, distance_bins, values)

Averaged binned tensor structure function — the mean tensor $D_{i…}(r)$ per distance bin (the sums ./ counts reduction of a StructureFunctionTensorSumsAndCounts, with the empty-bin guard count == 0 → NaN). For tensor order P and spatial dimension D, values has leading axes (D, D, ..., D, n_bins, auxiliary...) with P repeated component axes. This is the default result of calculate_structure_function_tensor; pass output_type = StructureFunctionTensorSumsAndCounts for the raw accumulator.

StructureFunctionTensorSumsAndCounts(order, distance_bins, sums, counts)

Raw binned tensor structure-function result. For tensor order P and spatial dimension D, sums has leading axes (D, D, ..., D, n_bins, auxiliary...) with P repeated component axes. counts has shape (n_bins, auxiliary...).

marginalize(sf2d::StructureFunction2DSumsAndCounts)

Sum $sums$ and $counts$ over the value-bin axis to produce a 1D $StructureFunctionSumsAndCounts$ (mass-conserving reduction of a 2D joint histogram).

AbstractBinEdges{T} <: AbstractVector{T}

Supertype for all custom, high-performance bin edge collections in StructureFunctions.jl.

Why AbstractBinEdges Exists

In structure function calculations over large datasets, the spatial separation distance $r$ for each of the $O(N^2)$ point pairs must be mapped to its corresponding distance bin (index). Using a standard sorted vector of bin edges requires a binary search (searchsortedfirst), which has $O(\log B)$ complexity where $B$ is the number of bins.

For large $N$, this binary search becomes the dominant CPU bottleneck, causing high branch mispredictions and cache misses. Subtypes of AbstractBinEdges bypass the standard binary search by implementing custom Base.searchsortedfirst overrides that execute in $O(1)$ time:

• LinearBinEdges utilizes Fused Multiply-Add (FMA) arithmetic for uniformly-spaced bins.
• LogBinEdges maps physical queries via log(q) then the same FMA path on the log-space grid.

Wrapping standard arrays in these subtypes allows digitize to execute 5x to 15x faster, resolving the primary computational bottleneck in the package.

Bin edges and backends

• Plain AbstractVector bin edges are valid everywhere but use generic searchsortedfirst / binary search.
• Pass LinearBinEdges, LogBinEdges, or InfPaddedBinEdges for O(1) CPU digitize and matching GPU tiled kernels (see StructureFunctionsKernelAbstractionsExt).
• InfPaddedBinEdges adds implicit catch-all under/overflow bins; do not vcat(-Inf, …, Inf) manually.
• Classic serial/threaded pair-loop paths normalize once via BinEdges; single-pass APIs expect callers to choose edge types explicitly (no auto-wrap on GPU).

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.

BinEdges(edges::AbstractVector{T})

Generic wrapper for arbitrary sorted vectors of bin edges. Bypasses range-specific optimizations but conforms to the AbstractBinEdges interface.

Behavior

• If constructed with an AbstractRange (e.g. StepRange or StepRangeLen), it automatically promotes and returns a LinearBinEdges wrapper to enable O(1) FMA indexing.
• Otherwise, it wraps the vector and delegates to standard $O(\log N)$ binary search methods.

BinEdges(edges)

Normalize flat edge input to AbstractBinEdges for hot-loop digitize:

• existing AbstractBinEdges → unchanged
• AbstractRange → LinearBinEdges
• other AbstractVector → wrapped via the default BinEdges struct constructor (generic binary search)

InfPaddedBinEdges(edges::AbstractVector{T})

Wrapper that implicitly prepends $-\infty$ (or typemin(T)) and appends $+\infty$ (or typemax(T)) to an existing bin edge collection.

Why InfPaddedBinEdges Exists

Structure function distance bins are defined as half-open intervals $(r_i, r_{i+1}]$. When mapping a distance $r$ to a bin, any query value $r < \text{first}(edges)$ or $r > \text{last}(edges)$ is out-of-bounds. Instead of checking for these out-of-bound cases manually using branches in inner loops, InfPaddedBinEdges embeds the infinite endpoints implicitly:

• The first element is treated as typemin(T) ($-\infty$).
• The last element is treated as typemax(T) ($+\infty$).

This guarantees that every valid positive separation distance maps to a valid index without allocating actual padding elements in memory or copying the array.

Prevention of Double Padding

The constructor checks if the input array already has infinite endpoints. If they exist, it trims them before wrapping to prevent nested padding (e.g. $[-\infty, -\infty, ...]$).

LinearBinEdges(edges::AbstractRange{T})

High-performance wrapper for uniformly-spaced ranges (linear spacing).

Mathematical Theory

A standard binary search takes $O(\log B)$ steps. With a uniformly spaced range of bin edges $v_i = v_1 + (i-1)\delta$, searchsortedfirst(v, x) is the smallest index $i$ with $v_i \ge x$: $ i^*(x) = \min{ i : v1 + (i-1)\delta \ge x } = \left\lfloor \frac{x - v1}{\delta} \right\rfloor + 1 = \lceil \frac{x - v_1}{\delta} + 1 \rceil . $ The discrete operator is ceiling / floor+1, not round (which answers a different question).

Precompute inv_step = 1/δ and at query time one FMA gives $t = (x-v_1)/\delta$: t = muladd(x, inv_step, -v_1 * inv_step) then idx = floor(Int, t) + 1.

Float boundary correction

t and reconstructed edges are inexact. After clamping idx, compare the reconstructed left edge u_idx = v_1 + (idx-1)δ; if u_idx < x, return idx + 1. One comparison fixes the at-most-one-bin FP error without a second downward correction.

Performance

Bypasses the Twice-Precision arithmetic of Julia's standard StepRangeLen search. Reduces lookup time from ~46 ns to ~3 ns (a 15x speedup), completely eliminating the linear binning bottleneck.

LogBinEdges(edges::AbstractVector{T})
LogBinEdges_from_log_edges(log_edges)

Log-spaced (geometric) bin edges: uniform grid in log-space.

Digitize semantics

log_edges is the authoritative grid. A physical query q > 0 maps to searchsortedfirst(LinearBinEdges(log_edges), log(q)). See docs/UNIFORM_BIN_DIGITIZE.md and benchmark/LOG_BIN_EDGES_BENCHMARK.md.

LogBinEdges from a physical vector builds the same log grid from range(log(first), log(last); length). LogBinEdges_from_log_edges accepts the log grid directly.

Performance

One log(q) (~6 ns) plus O(1) FMA digitize (~2.5 ns) on the log grid. See benchmark log.

Arbitrary sorted edges: binary search on the precomputed squared edges (skips the sqrt).

Implicit ±Inf catch-all bins around an inner plan; every pair lands somewhere.

Uniform-in-r edges: squares of a uniform grid are not uniform, so this keeps one sqrt and the existing O(1) FMA search.

Log-uniform edges: index from _fast_log2(r²) and one FMA, then corrected.

LogBinEdges_from_log_edges(log_edges)

Build log-uniform bins from log-space edges (u with physical edge exp(u)).

_fast_log2(x)

log2(x) for finite x > 0, max error ~4e-8 (Float64). Exponent extract plus an odd series on a mantissa recentred to [1/√2, √2); vectorizes, unlike the scalar libm log. Approximate — the bin is decided by squared_digitize's correction, not by this.

Map 1-based upper-triangle pair index to (i, j) with i < j for N points.

Trailing auxiliary-axis shape of u after leading (D, N).

View (D, N) field slice at linear auxiliary index b.

Product of trailing auxiliary-axis dimensions.

digitize_key(plan, r2)

The quantity the plan's scalar half compares against, computed in the vectorized half from r².

r² for the log and general plans; √r² for the linear plan, whose grid is uniform in r.

has_vector_index(plan) -> Bool

Whether the plan's index is branch-free, and so worth computing in the vectorized half of a pair kernel. False for the linear and general plans, whose searchsortedfirst branches would de-vectorize the whole @simd body.

n_histogram_bins(edges::AbstractVector) -> Int

Number of histogram bins for flat edges (length == N + 1 → N bins).

Bins covered by the plan (digitize results in 1:n_bins are in range).

Public alias for tests and callers enumerating upper-triangle pairs.

physical_edges_vector(bins::LogBinEdges) -> Vector

Materialize physical bin edges exp(log_edges[i]) for display or generic consumers. Not used on the digitize hot path.

squared_approx_index(plan, r2) -> Int32

Approximate searchsortedfirst index for r², branch-free, for the vectorized half of a pair kernel. Meaningful only when has_vector_index; other plans return 0.

squared_bin(plan, key, i) -> Int

The bin, in the scalar half, from digitize_key. When has_vector_index, i is the precomputed approximate index and this only corrects it; otherwise i is ignored.

squared_correct(plan, r2, i) -> Int

Walk i to the exact searchsortedfirst(sqedges, r²). 0 or 1 step for random separations; a loop rather than a fixed step because within a few ulps of an edge more can be needed.

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 digitize.

squared_digitize_plan(edges) -> AbstractSquaredDigitizePlan

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