API

HelmholtzDecomposition.jl — Helmholtz–Hodge decomposition of velocity fields.

Decomposes a velocity field into rotational (divergence-free), divergent (curl-free) and harmonic parts, in any number of dimensions, on Cartesian and spherical grids.

Why this package exists

On the sphere, filtering velocity Cartesian components does not commute with differential operators (Aluie 2019, Proposition 2). The correct approach — mathematically equivalent to the generalized convolution that does commute — is to filter the scalar Helmholtz potentials (ψ, χ) separately. This package provides the decomposition step.

Solver Extensions (important for performance!)

The base package includes only the SOR iterative solver, which works on any grid but may be orders of magnitude slower than spectral solvers. Load an appropriate extension:

Quick Start

using HelmholtzDecomposition: HelmholtzDecomposition as HD
using FlowGeometries: FlowGeometries as FG
using FFTW: FFTW  # load spectral extension for Cartesian grids

grid = FG.Grids.StructuredGrid(FG.Geometry.CartesianGeometry{Float64}(), xs, ys)
result = HD.helmholtz_decompose(u, v, grid)
# result.u_rot, result.v_rot  — rotational velocity
# result.u_div, result.v_div  — divergent velocity
# result.ψ, result.χ          — scalar potentials

Nothing is exported: every name is reached as HelmholtzDecomposition.name.

References

• Aluie (2019): doi:10.1007/s13137-019-0123-9 — Convolutions on the sphere
• Glötzl & Richters (2023): doi:10.1016/j.jmaa.2023.127138 — n-dimensional Helmholtz potentials
• Buzzicotti et al. (2023): doi:10.1126/sciadv.adi7420 — Global cascade of kinetic energy
• Storer et al. (2022): doi:10.1038/s41467-022-33031-3 — Global energy spectrum

helmholtz_decompose(u, grid; boundary = Neumann(), solver = AutoSolver()) -> HelmholtzResult

Decompose a cell-centred, component-last velocity field on grid.

Builds a HelmholtzPlan and discards it. A caller decomposing more than one field on the same grid should build the plan once and call helmholtz_decompose! — the plan holds the grid's Laplacian coefficients, which are the expensive part.

helmholtz_decompose!(result, u, plan, ws = allocate_workspace(plan)) -> result

Decompose the cell-centred, component-last velocity u into result, in place.

The solver is not an argument here: it is fixed by plan. Accepting one would either be ignored — a silent no-op — or force the plan's transform state to be rebuilt for a different solver, which is the per-call cost the plan exists to remove. Pass solver to plan_helmholtz.

helmholtz_decompose_batch(plan, fields; backend = AutoBackend())

Decompose many fields sharing one grid. The fields are independent, so the batch is the parallel axis: ThreadedBackend (OhMyThreads ext), DistributedBackend (Distributed ext) and MPIBackend (MPI ext) each spread it across their workers; a serial backend maps sequentially. Results come back in input order.

One plan serves the whole batch — its Laplacian coefficients are identical for every field and are the expensive part — while each task takes its own HelmholtzWorkspace, since those buffers are written through.

backend here is the batch axis. A parallel one pins the decomposition inside each field to serial; a serial one leaves the inner loops on the plan's own backend, so a short batch of large fields still parallelises — just along the other axis.

helmholtz_decompose_batch!(batch, fields, plan; backend = AutoBackend()) -> batch

Decompose into a preallocated HelmholtzBatch.

plan_helmholtz(grid; boundary = Neumann(), solver = AutoSolver(), backend = AutoBackend())

Resolve everything a decomposition on grid fixes once: the Laplacian coefficients on the primal and dual grids, which solver will run, that solver's reusable state, and which execution backend the loops inside one decomposition use.

The solver is chosen here rather than per call deliberately. Choosing it reads the mask, which is O(N), and preparing it builds the transform plans; doing both inside helmholtz_decompose! repeats them for each of the P + 1 potentials and again for every field of a batch, which is exactly the work a batch exists to amortize.

backend is the intra-field one — it drives the operator and solver loops, so a single large field is parallel rather than only a batch of them. A batch overrides it with a serial inner backend, because an outer and an inner loop each claiming every thread is slower than either.

HelmholtzPlan

The shareable part of a decomposition on one (grid, boundary, eltype): the Laplacian coefficients on the primal grid and on each dual grid, and the dual grids themselves. Buffers live in a HelmholtzWorkspace, one per task, precisely so that a whole batch can share a single plan without racing.

The geometry behind those coefficients — face areas need the scale factors at the face, which on a curved grid is trigonometry per cell per direction — is invariant across solver iterations, across the P + 1 potentials, and across a whole batch. Rebuilding it per call is the single largest avoidable cost in a decomposition, so it is built here and only read afterwards.

HelmholtzWorkspace

The mutable buffers one decomposition writes through, held apart from the plan.

The split is what makes a batch correct rather than merely possible. A HelmholtzPlan holds only things that are the same for every field on a grid — the Laplacian coefficients on the primal and dual grids — so a batch shares exactly one, which is the whole point, those being the expensive part. A workspace holds the face and corner buffers, which every field writes through, so each task needs its own; sharing one across threads would race on them.

allocate_workspace(plan) -> HelmholtzWorkspace

Buffers for one task working through plan.

HelmholtzResult

Velocity-like fields use the component-last layout (dims..., N); potentials and scalar diagnostics are (dims...).

The array fields are const and the diagnostics are not: an in-place decomposition writes its convergence information into the result it was handed rather than allocating a fresh struct to carry two numbers.

allocate_result(plan) -> HelmholtzResult

A zeroed result matching plan. Array types follow the plan's buffers, so a device-resident plan gives a device-resident result — nothing here hardcodes Array.

HelmholtzBatch

B decompositions on one grid, stored contiguously: each field of a HelmholtzResult gains a trailing batch axis, and batch[b] is a result of views onto slice b.

One allocation per output for the whole batch rather than one per field.

allocate_batch(plan, nfields) -> HelmholtzBatch

streamfunction(result) -> Array

The 2-D streamfunction ψ, on the corner (dual) grid. Only defined in 2D.

velocity_potential(result) -> Array

The scalar velocity potential χ, at cell centres, in any dimension.

vector_potential(result) -> (A1, A2, A3)

The 3-D vector potential, the Hodge dual of the rotation potential: A1 = R_23, A2 = −R_13, A3 = R_12. Only defined in 3D.

count_holes(grid) -> Int

Number of inactive regions fully enclosed by active cells — an estimate of the first Betti number b₁ of the active region, i.e. the dimension of the harmonic subspace.

This is what makes harmonic_fraction readable: a large harmonic part on a domain with count_holes == 0 is a boundary-circulation effect, while on one with holes it is the circulation around them, which no single-valued potential can represent.

The flood fill is FlowGeometries.Connectivity.connected_components, which walks the grid's own wrapping — so on a periodic direction a region running off one side and back on the other is one region, and "touching the boundary" is asked only of directions that actually have one.

face_velocity(ws), face_divergent(ws), face_rotational(ws) -> NTuple{N,Array}

The staggered fields the last decomposition through ws computed: one face-normal array per direction, nfaces(grid, d) long where the direction ends and n where it wraps.

These are where the projection is exact. face_divergent reproduces the input's divergence to round-off and face_rotational is divergence-free to round-off, both discretely; the collocated u_div and u_rot on the result are these interpolated back to cell centres, and that round trip is the entire harmonic residue a smooth field shows — 9.6e-3 at n = 32, falling at second order. A caller who works on the C-grid takes these and has no interpolation error at all.

Valid until the next decomposition through the same workspace, which overwrites them.

face_divergent(ws) -> NTuple{N,Array}

Gχ on faces — see face_velocity.

face_rotational(ws) -> NTuple{N,Array}

−δR on faces, divergence-free to round-off — see face_velocity.

corner_rotation_potential(ws) -> NTuple{P,Array}

R on the corner (dual) grid, one array per rotation pair — see face_velocity.

AbstractBoundaryCondition

Supertype for conditions on a bounded direction. Concrete: Dirichlet, Neumann.

Homogeneous Dirichlet condition (Φ = 0 on the boundary).

Homogeneous Neumann condition (∂Φ/∂n = 0 on the boundary): no flux through the domain edge.

AbstractPoissonSolver

Supertype for Poisson solvers. Concrete subtypes implement

solve_poisson!(Φ, RHS, grid, solver; boundary, kwargs...) -> SolverResult

AutoSolver()

Sentinel for automatic selection. The only thing permitted to choose a solver on the caller's behalf, and it chooses on real capability — geometry, node layout, mask, axis uniformity, the requested boundary condition, and which extensions are loaded.

CGSolver(; max_iter = 1000, rtol = 1e-10)

Conjugate gradients on −L, which is symmetric positive (semi)definite by construction — see Operators.jl. Jacobi-preconditioned, and on a closed problem the iterate and residual are kept orthogonal to the constants.

Works on any grid, any mask, any boundary condition. Stage 6 adds a multigrid preconditioner in place of Jacobi; the outer iteration is unchanged by that.

solve_poisson!(Φ, RHS, grid, solver; boundary, coefficients, state, backend) -> SolverResult

Solve L Φ = RHS in place, with L the operator Operators.jl builds.

coefficients and state are the reusable parts — the Hodge factors of every face, and whatever the solver prepared for this (grid, boundary) pair. A decomposition solves P + 1 right-hand sides and a batch multiplies that by the field count, so both are built once and passed in rather than derived here.

Extensions add methods for their own solver types; the iterative CGSolver is the one that works on any grid, mask and boundary condition.

prepare_solver(solver, grid, boundary) -> state

Whatever solver can compute once for a (grid, boundary) pair and reuse on every solve — returned opaquely and handed back to solve_poisson! as its state keyword.

This exists because the transform plans are the dominant per-call cost and depend on nothing that changes between calls. A decomposition solves P + 1 right-hand sides on one grid and a batch multiplies that by the number of fields, so a solver planning inside solve_poisson! rebuilds the same plan (P + 1) · B times: in 3-D that is 17 FFTW plans per field, and for the non-uniform transforms a plan and a node upload per conjugate-gradient iteration.

nothing is the default and means "nothing to reuse", which is correct for the iterative solvers — their reusable part is the LaplacianCoefficients, which the plan already holds.

select_solver(solver, grid, boundary) -> AbstractPoissonSolver

The concrete solver for this problem, resolved and validated once. Resolution reads the mask, so a caller with several right-hand sides on one grid — which a decomposition always has, P + 1 of them — calls this once and hands the result to each solve.

SolverResult{T}

Convergence diagnostics from a solve: whether it met its tolerance, how many iterations it took (1 for a direct or spectral solve), and the final residual.

supports_boundary(solver, boundary) -> Bool

Whether solver solves the problem boundary describes. A solver that does not is never selected by AutoSolver and errors when named directly — solving a different problem from the one asked for is indistinguishable from a correct answer at the call site.

requires_full_domain(solver) -> Bool

Whether solver needs every cell active. The spectral solvers transform the whole array, so a masked cell would be transformed as though it held data; they set this and are refused a masked grid rather than returning a field that is wrong wherever the mask bites.

requires_uniform_axes(solver) -> Bool

Whether solver needs constant spacing in every direction. An FFT-based solver does; the question is answered from the axis TYPE by FlowGeometries, at no runtime cost.

requires_periodic_domain(solver) -> Bool

Whether solver needs every direction to wrap. An FFT expands in a periodic basis, so it solves the whole-torus problem and nothing else.

This is a question about the grid's topology, not about a boundary condition — which is why there is no Periodic boundary condition to ask instead. A direction that wraps has no boundary to impose a condition on, so on a grid this solver accepts, boundary is vacuous rather than honoured or refused.

register_spectral_solver!(algorithm, solver_type; priority)

Declare that solver_type implements algorithm. Lower priority is tried first: a native implementation takes a lower number than a generic one that would also work.

gradient!(g, χ, grid, bc) -> g

g[d][F] = (χ_above − χ_below) / gap on every face. A face with area but only one cell is a Dirichlet edge, whose ghost value is zero — which is exactly what the condition says.

divergence!(δ, v, grid, bc) -> δ

δ[I] = (1/V_I) Σ_d (A·v)[face above] − (A·v)[face below] — the negative adjoint of gradient! under the cell-measure inner product, which is what makes L = D G symmetric.

laplacian!(out, χ, grid, bc, scratch) -> out

L χ = D G χ, from exactly the two operators above, so a solver inverts the same L the decomposition differentiates with.

apply_laplacian!(out, Φ, grid, c) -> out

out = L Φ, from the prebuilt coefficients. A face at the outer edge of a bounded direction has no cell beyond it and contributes c·(0 − Φ_I) — the zero ghost a Dirichlet condition places there. Under Neumann that face's coefficient is zero and the term vanishes, so one expression serves both.

curl!(W, v, grid, bc) -> W

W_ab = ∂_a v_b − ∂_b v_a on the (a,b) corner, from face-normal velocity.

The circulation v·g is differenced rather than v itself, so what is differenced is the metric-free d; the metric re-enters once, in the division at the end. A corner is included only when all four faces bounding it are open — a loop running half through a mask edge has no zero circulation, and requiring only two faces is what previously broke curl(grad χ) = 0 on a masked grid.

rotational_velocity!(u_rot, R, grid, bc) -> u_rot

u_rot_a = −Σ_b ∂_b R_ab, moving each component of R from its corner onto the a-faces. The adjoint of curl!, so div(u_rot) = 0 holds discretely.

to_faces!(vf, uc, grid, bc) -> vf

Cell-centred velocity (dims..., N) to face-normal velocity, averaging the two cells a face separates. A face of zero area takes nothing; an open boundary face has only one cell to read.

to_centres!(uc, vf, grid, bc) -> uc

Face-normal velocity back to cell centres, averaging a cell's two faces — the adjoint of to_faces!. A cell against a closed face averages only the faces that carry flux.

LaplacianCoefficients{N,T,A,B}

The per-face coefficients c_f = A_f / g_f of L, and the cell measures, evaluated once for a (grid, boundary) pair.

Every iterative solve applies L repeatedly, and each application would otherwise re-evaluate the geometry: face areas involve the scale factors at the face, which on a curved geometry means trigonometry per cell per direction. That is invariant across iterations — and across the P + 1 potentials and the whole batch — so it is reduced once here.

coef[I, d] is the coefficient of the face below cell I along d; diag[I] is the negative sum of a cell's own face coefficients, divided by its measure.

laplacian_coefficients(grid, bc) -> LaplacianCoefficients

The Hodge factor c_F = A_F / g_F of every face, and the cell measures, reduced once for a (grid, boundary) pair.

Every iterative solve applies L repeatedly, and each application would otherwise re-evaluate the geometry — face areas involve the scale factors at the face, i.e. trigonometry per cell per direction on a curved grid. That is invariant across iterations, across the P + 1 potentials, and across a whole batch, so it is computed once here and only read afterwards.

FaceMetrics{N,T,A}

Every face's flux-carrying area and centre-to-centre gap, evaluated once for a (grid, boundary) pair.

face_area and face_gap each read the grid's coordinates and evaluate the geometry's scale factors at the face — on a curved grid, trigonometry per face per direction. None of it depends on the field, yet the operators call them inside their loops, so a decomposition recomputed the same metric 2N + 4P times per field and again for every field of a batch. Reduced here to two array reads.

face_area(grid, F, d, bc, T) -> T

The area that carries flux across face F — the coefficient the operators use, and where the boundary condition lives, because on a flux-form operator that is what a boundary condition is:

• between two active cells → the geometric area;
• against a masked-out cell → 0. No flux crosses, which is at once the no-flux condition, the mask treatment, and the reason D = −G* survives both;
• at the outer edge of a bounded direction → 0 under Neumann, the geometric area under Dirichlet, whose ghost value beyond the edge is zero.

On a covering spherical grid the poles need no special case: the φ-face area carries the cos φ that vanishes there, so the metric closes the surface itself.

face_gap(grid, F, d, T) -> T

Physical centre-to-centre distance across face F: the mean of the two cells' physical widths, or half a cell where the face has only one. Consistent with face_area — a great-circle distance would not be, being shorter than the coordinate line the area is built on.

nfaces(grid, d) -> Int

Number of faces normal to direction d: n if the direction wraps, n + 1 if it ends.

ncorners(grid, d) -> Int

Corners along a staggered direction: n where it wraps, n - 1 where it ends.

A bounded direction has n + 1 faces but only n - 1 corners that carry an unknown. The two outermost have no closed loop of cells around them, so the circulation there is not defined and R is not solved for — it is zero. Leaving them in the array and masking them off says the same thing, but says it as data, and a masked grid is refused by every direct transform. Dropping them says it as a Dirichlet condition on a smaller domain, which a sine transform inverts exactly. Measured at 512²: 98 iterations and 4.6 s became one transform and 19 ms.

corner_offset(grid, d) -> Int

What to add to a corner index in direction d to get the face index it sits on: 1 where the outer pair was dropped, 0 where the direction wraps and nothing was.

face_dims(grid, d) -> NTuple{N,Int}

Shape of the face array normal to direction d.

corner_dims(grid, a, b) -> NTuple{N,Int}

Shape of the array holding a rotation-potential component staggered in both a and b.

allocate_faces(T, grid) -> NTuple{N,Array{T,N}}

One zeroed face array per direction. This is the internal staggered layout; the public API stays collocated, so nothing outside this package sees it.

allocate_corners(T, grid; backend) -> NTuple{P,AbstractArray{T,N}}

One zeroed corner array per rotation pair, shaped by corner_dims.

rotation_pairs(Val(N)) -> NTuple{P,Tuple{Int,Int}}

The pairs (a, b), a < b, lexicographically — the independent components of an antisymmetric 2-tensor, and equally the directions each component is staggered in.

rotation_terms(Val(N)) -> NTuple{N,NTuple{N-1,Tuple{Int,Int,Int}}}

For each direction c, the (e, p, s) triples saying how u_rot_c is assembled: difference component p of R along direction e with sign s. Inverting rotation_pairs this way turns the accumulation u_rot_a −= ∂_b R_ab, u_rot_b += ∂_a R_ab — a scatter, where several pairs write the same face — into a gather, so each face is written exactly once and the loop can be handed to Execution unchanged. Homogeneous, so indexing it with a runtime direction is stable.

n_rotation_components(N) -> Int

N(N-1)/2: 0 in 1D, 1 in 2D, 3 in 3D.

face_coordinates(grid, d) -> axis

Coordinates of the faces normal to direction d — the midpoints between consecutive cell centres, plus the two outer edges where the direction is bounded.

A periodic direction gets n of them and a bounded one n + 1, matching nfaces.

A uniform axis returns a range. Uniformly spaced cells have uniformly spaced faces, but FlowGeometries proves uniformity from the axis TYPE, so returning a Vector here would make every dual grid report itself as stretched — and requires_uniform_axes would then refuse the direct transform on a grid that qualifies for it. That is not a small loss: it left the rotation-potential solve iterative at every size, and it was the whole cost of a 512² decomposition.

corner_mask(grid, a, b) -> Array{Bool}

Which (a,b) corners have a closed loop of active cells around them — the same four-face test curl! applies, hoisted so the dual grid carries it.

A corner without a closed loop is not in the complex: the circulation around it is not defined, so R there is not solved for and contributes nothing back.

dual_grid(grid, a, b, bc) -> StructuredGrid

The grid the (a,b) rotation-potential component lives on: corner coordinates in directions a and b, cell coordinates elsewhere, the primal's geometry and per-direction topology, and corner_mask as its mask.

A bounded direction contributes its interior corners only — see ncorners. The pair it drops is where R is pinned to zero, which the dual solve states as Dirichlet rather than as a mask, so an unmasked domain stays unmasked here and keeps its direct transform.

execution_backend(backend)

The object FlowGeometries.Execution dispatches on, given a ComputationalBackends one.

The two libraries answer different questions: ComputationalBackends names what kind of execution is wanted, while Execution's device methods dispatch on the KernelAbstractions device object itself. A GPUBackend therefore hands over the device it names; serial and threaded backends pass through, FlowGeometries' own ComputationalBackends extension having methods for them. Anything else reaches run_indices unchanged and raises a MethodError there rather than running serially in silence.

resolve_execution_backend(backend) -> concrete backend

The backend a single decomposition's own loops run on, resolved once so nothing below has to ask.

AutoBackend reads the thread count — the defect that made every documented parallel path dead at defaults was resolving it to serial unconditionally. A GPUBackend is checked here, at plan time, against whether FlowGeometries' KernelAbstractions extension is actually loaded, so the error names the package to load rather than surfacing as a MethodError from inside a loop.

DistributedBackend and MPIBackend are refused: they spread fields across processes, which is helmholtz_decompose_batch's axis, not the index space inside one field.

allocate_zeros(backend, T, dims) -> zeroed array

Every buffer this package owns goes through here, so that a device backend gets device memory rather than a host array a kernel cannot reach. The default is a host Array; the KernelAbstractions extension adds the device method.

to_backend(backend, x) -> x on the backend's memory

Move an already-built array — or a struct of them — to where backend executes.

The plan's coefficients and face metrics are built by walking the grid's geometry, which is host work done once; only their results are read in the inner loops. So they are assembled on the host and moved here, rather than every geometry accessor being made device-callable to build them in place. Identity by default; the KernelAbstractions extension routes it through Adapt.

multigrid(grid, bc, T; ω = 0.8, ν = 2, maxlevels = 12) -> MultigridPreconditioner

Build the hierarchy. ω = 0.8 is the usual damped-Jacobi factor for a 2-D five-point Laplacian — undamped Jacobi does not reduce the highest-frequency error at all, which is the error a smoother exists to remove.

MultigridPreconditioner

A hierarchy built once for a (grid, boundary) pair, applied as z ← M⁻¹ r inside conjugate gradients.

MultigridLevel

One rung: the grid, its Laplacian coefficients, and the buffers a cycle needs there.

coarsen(grid, bc) -> grid or nothing

The next grid down: every other coordinate in each direction that still has enough cells to halve.

A coarse cell is active when any of the fine cells it covers is. The alternative — requiring all of them — erodes the domain by a cell per level, so a narrow channel or a coastline would vanish partway down the hierarchy and the correction there would be identically zero.

galerkin_coefficients(fine, fgrid, cgrid) -> LaplacianCoefficients

The coarse operator as R A P rather than a fresh discretization of the coarse grid.

Re-discretizing is what breaks a masked domain. The coarse mask cannot reproduce the fine one — mark a coarse cell active if any child is and the domain grows outward a cell per level, require all and it erodes — so the coarse operator ends up solving a differently-shaped problem and its correction is inconsistent with the fine one. That showed up as the rotation-potential solves running to their iteration cap while the primal solve converged in 13.

For this transfer pair — restriction that averages a cell's children, prolongation that adds a coarse value to each of them — the Galerkin product has a closed form and needs no matrix: the conductance of a coarse face is the sum of the fine conductances crossing it, and a coarse cell's measure is the sum of its children's. Nothing about the geometry is recomputed, so nothing can disagree with the fine level.

restrict!(coarse, fine, cgrid, fgrid, fc, cc)

Measure-weighted average of each coarse cell's fine children.

The weighting is not a refinement — it is what makes R the adjoint of prolongation, and a non-adjoint pair makes the preconditioner non-symmetric, which invalidates conjugate gradients. With prolongation by injection and the coarse measure defined as the sum of its children's,

⟨P c, f⟩_fine = Σ_I Σ_{J∈children(I)} V_J c_I f_J = ⟨c, R f⟩_coarse
⟹  (R f)_I = Σ_J V_J f_J / Σ_J V_J

An unweighted average satisfies that only when every V_J is equal, which is why it worked on a Cartesian grid and failed outright on a sphere, where the measure varies as cos φ: multigrid ran to the iteration cap at every resolution while plain Jacobi converged.

No mask term is needed — an inactive cell's measure is already zero, so it carries no weight.

prolong_add!(fine, coarse, cgrid, fgrid, fc)

Add each coarse value to the fine cells it covers — injection, the adjoint of the measure-weighted restrict!, which is what keeps the cycle symmetric.

smooth!(x, b, level, ω, n)

n sweeps of damped Jacobi on −L, the positive-definite orientation.

Damped Jacobi rather than Gauss–Seidel: it is a pure Ax plus an axpy, so it is the same expression serial, threaded and on a device, and it is symmetric — which Gauss–Seidel is not unless the sweep order is reversed between the pre- and post-smoothing, and a non-symmetric preconditioner breaks conjugate gradients.

vcycle!(levels, ω, ν, backend)

One V-cycle over levels, coarsest last: smooth, restrict the residual, recurse on the tail, prolong the correction, smooth again. Equal pre- and post-smoothing keeps the cycle symmetric.

The recursion is on the tuple's tail, not on an index. Recursing with Val(l + 1) gives the compiler no proof the depth terminates, so it abandons inference at the recursive call and boxes the arguments — 1 KiB per level per cycle, which at 81 conjugate-gradient iterations was most of a megabyte per solve. A shorter tuple type each step terminates structurally.

helmholtz_project_spectral!(û_rot, û_div, û_harm, velocity_hat, ks::NTuple{N})

In-place Leray projection. velocity_hat and the three outputs are component-last spectral arrays of size (kdims..., N); ks holds the per-axis wavenumber vectors. Writes the rotational (divergence-free) part into û_rot, the divergent (curl-free) part into û_div, and the k = 0 mode into û_harm, so that the three sum to velocity_hat.

The k = 0 mode is separated rather than left in û_rot because a constant field is both curl-free and divergence-free: it is the harmonic part, not the rotational one. GPU-compatible (pure broadcast).

helmholtz_project_spectral(velocity_hat, ks::NTuple) -> SpectralCartesianResult

Allocating Leray projection from a component-last spectral array.

helmholtz_potentials_spectral(velocity_hat, ks::NTuple{N}) -> (χ_hat, R_hat)

Compute the spectral scalar velocity potential χ_hat (size (kdims...)) and the rotation-potential components R_hat (component-last, size (kdims..., P), P = N(N-1)/2) from a component-last spectral velocity array. Uses the spectral Poisson inverses χ̂ = −i (k·û)/k² and R̂_ab = −i (k_a û_b − k_b û_a)/k², with the k = 0 mode set to zero.

helmholtz_decompose_spectral(u, grid; kwargs...)
helmholtz_decompose_spectral(u, v, grid; kwargs...)        # 2D convenience
helmholtz_decompose_spectral(u, v, w, grid; kwargs...)     # 3D convenience

Decompose a physical velocity field on grid using a spectral transform, returning a physical HelmholtzResult (CPU) — or, on the GPU path, a (; u_rot, u_div, u_harm) NamedTuple of CuArrays. Requires the appropriate extension (using FFTW, using FastSphericalHarmonics, …). Pass solver= to select among loaded spectral backends.

For raw spectral coefficients, use the lower-level helmholtz_project_spectral.

velocity_norm(U, grid, scratch) -> T

Measure-weighted ‖U‖ = sqrt(∫|U|² dV) over the active cells.

project_out_constant!(Φ, grid, c)

Remove the measure-weighted mean of Φ over the active cells.

A closed problem — every direction periodic, or a Neumann boundary — leaves the constants in L's null space. Krylov iterations must stay orthogonal to that null space or they drift along it, so this is applied to the right-hand side once and to the iterate as it goes.