Internals

AbstractSpectralHelmholtzResult{T}

Abstract supertype for spectral decomposition results.

SpectralCartesianResult{T,A}

Cartesian spectral decomposition result holding the component-last complex Fourier coefficients of the rotational, divergent and harmonic velocity ((kdims..., N)).

u_harm carries the k = 0 mode. On a periodic domain the constant fields are exactly the harmonic subspace — H¹(Tᴺ) has dimension N, and a constant field is both curl-free and divergence-free — so the mean flow belongs to neither of the other two parts.

SpectralSphericalResult{T,V}

Spherical spectral decomposition result holding the spherical-harmonic coefficients of the streamfunction ψ and velocity potential χ, up to degree lmax.

build_cartesian_result(grid, U, velocity_hat, ks, inverse) -> HelmholtzResult

Assemble a physical HelmholtzResult from a component-last spectral velocity array. inverse(spectral_scalar) maps one spectral scalar field ((kdims...)) back to a physical one ((dims...)). Used by the regular-grid spectral extensions (FFTW).

N + P + 2 inverse transforms are taken, not 2N + 2P + 2:

• u_rot is U − u_div − u_harm in physical space, so it costs nothing;
• u_harm is the k = 0 mode, i.e. a constant field, and is the component mean — no transform;
• the rotation tensor W = ΔR is not materialised at all. It was P inverse transforms for a quantity recomputable from R, and no longer has a field on the result.

On this path rotation_potential holds cell-centred arrays, where the finite-difference path holds corner-staggered ones. The field is type-parameterised for exactly this reason, but the staggering does depend on which path produced the result.

CGWorkspace{A}

The four vectors a conjugate-gradient iteration needs, held so the solve allocates nothing.

CGState{W,P}

Everything CGSolver reuses across solves on one (grid, boundary): its vectors and its preconditioner.

The vectors belong here and not in a default argument. helmholtz_decompose! solves P + 1 right-hand sides and a batch multiplies that by the field count, so defaulting them allocates four grid arrays per solve for the lifetime of the program.

_precondition!(z, r, grid, c, preconditioner)

Apply M⁻¹ to the residual: the Jacobi diagonal when there is no preconditioner, a V-cycle when there is (Multigrid.jl adds that method, where the type lives).

_detect_singular(grid, c) -> Bool

Whether the constants lie in L's null space, decided by applying L to the constant field and seeing whether anything comes back.

Asking the boundary condition instead is not enough, and the dual grid is why: its boundary ring of corners has no closed loop of cells around it, so those corners are masked out and every face of the active region is closed — the operator there is singular even under Dirichlet, which a condition-based test reports as nonsingular. Conjugate gradients on a singular system with no null-space projection then drifts along the null space instead of converging, which showed up as a harmonic fraction of 1.7e+32.

L·1 = 0 is exactly the property that matters and costs one operator application at plan time.

_require_domain(solver, grid)

Throw unless grid is a domain solver can solve on. Reads the mask, so it is O(N) and is called once where a solver is chosen — never from solve_poisson!, which a decomposition enters P + 1 times per field and once more per field of a batch.

_require_sampling(solver, grid)

Extra, solver-specific demand on the node layout — a transform defined on one node set and no other says so here. A hook rather than an override of _require_domain: overriding that would silently drop the mask, uniformity and topology checks it also performs.

_spectral_algorithms(geometry_type) -> Tuple

The algorithms worth trying on a geometry, in preference order: the uniform transform first, then the non-uniform one, which subsumes it at greater cost. A candidate that does not apply refuses itself through its own capability check, so this order only decides between several that do.

_spectral_dispatch(u, grid; solver=AutoSolver(), kwargs...)

Resolve a spectral solver (extensions register them; AutoSolver picks the best available) and dispatch to the extension's _decompose_spectral(solver, geometry, u, grid; …). Dispatching on the solver type lets several spectral backends (FFTW + FINUFFT, FSH + NUFSHT) coexist for the same geometry without method clashes. The CUDA extension overrides this for CuArray inputs to take the CUFFT path directly.

cell_above(grid, F, d) -> CartesianIndex or nothing

The cell on the high side of face F along d; nothing at the outer edge of a bounded direction. On a periodic direction every face has both.

cell_below(grid, F, d) -> CartesianIndex or nothing

The cell on the low side of face F along d; nothing at the outer edge of a bounded direction.

face_below(I, d) / face_above(grid, I, d) -> CartesianIndex

The two faces of cell I normal to d, as indices into face array d. On a periodic direction the face above the last cell wraps to face 1.

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

Area of face F normal to d: the product of the physical cell widths in every other direction, evaluated at the face. Purely geometric — whether flux crosses it is face_area's question.

copy_slice!(dest, i, src, j) -> dest

Copy slice j of one batch into slot i of another, fields and diagnostics alike.

_resolve_batch_backend(backend, fields) -> AbstractExecutionBackend

Resolve Auto over the batch axis against real capability: more than one field, more than one thread, and the threading extension actually loaded. A concrete backend passes through and is then honoured or refused — never quietly replaced.

_extension_loaded(name) -> Bool

Whether one of this package's extensions is loaded. hasmethod cannot answer this — the dispatch stubs exist unconditionally, so it is always true.

_unsupported_backend_message(backend, entry) -> String

Why entry cannot honour backend, and what to load so it can. Reaching this is always an error: a backend named explicitly is executed as named or refused, never downgraded.

_grid_wavenumbers(velocity_hat, grid) -> NTuple{N}

Reconstruct the per-axis angular wavenumber vectors for a component-last spectral array on a Cartesian grid. Axis 1 is treated as an rfft axis when its spectral length equals N₁÷2 + 1, otherwise as a full fft axis.

_stack_components(grid, comps...) -> Array

Stack scalar component arrays into the component-last layout (dims..., M) the transforms take.