Discretization

AbstractLocation

Where a value sits within a cell: Center or Face.

AbstractMaskPolicy

What apply_stencil! does at the edge of the active region: BlankMasked, ShiftWithinRun or ReduceInRun.

A type, like the image and reach conventions in Connectivity: which cells carry a number and which carry masked is a property of the result, so it belongs in the call rather than in a runtime tag.

BlankMasked()

Write masked at a cell that is inactive or whose stencil reads an inactive cell. The default, and the only policy that never invents a value: where the stencil cannot be formed from active data, there is no derivative.

Its cost is a dead band. Every active cell within nodes - 1 of a masked cell is blanked, so a five-point derivative loses two cells either side of every coastline.

Center()

Cell centres — N of them along a direction of N cells.

Face()

Cell boundaries — N+1 of them along a direction of N cells.

Centres do not determine faces: N centres leave the N+1 faces underdetermined by one. The convention here is that faces sit midway between neighbouring centres, with the two outermost placed by linear extrapolation, which is the same rule the curvilinear corner reconstruction uses.

GradientPlan{T}

The geometry of a least-squares gradient, separated from any field: for each cell, the coefficient on each neighbour's difference from it. Built by Connectivity.gradient_plan, applied by gradient!.

apply_stencil! covers the separable case. A CurvilinearGrid has no separable axis to build a scalar stencil along, and a node set's neighbours come from connectivity rather than an index offset, so neither has a gradient without this.

The construction, with tangent-plane displacements Δrₖ from Geometry.project_to_tangent_plane, differences Δfₖ = fₖ - f₀ and weights wₖ = 1/|Δrₖ|²: minimising Σₖ wₖ (∇f·Δrₖ - Δfₖ)² gives A ∇f = b with A = Σₖ wₖ Δrₖ⊗Δrₖ and b = Σₖ wₖ Δrₖ Δfₖ. A depends only on the geometry, so it is inverted once here and the per-neighbour coefficients A⁻¹ wₖ Δrₖ are what the plan stores; each apply is then one dot product per cell and allocates nothing.

What this buys over inverting an index-space Jacobian:

• exact for a linear field on any stencil, however skewed — the least-squares combination cancels the leading truncation term, which a 2×2 Jacobian inverse does not;
• second order on locally symmetric stencils, degrading toward first on strongly skewed cells;
• it reduces to the centred difference where the stencil is separable and orthogonal (A diagonal), so it agrees with apply_stencil! where both apply.

Where A is rank deficient — a tangent direction with no data, at a boundary or beside a mask — that component is zeroed rather than invented, by the pseudo-inverse. Same rule apply_stencil! states for a mask: not determined by the active data, so not produced.

ReduceInRun()

ShiftWithinRun, and where the run cannot hold nodes, use the largest window it can, down to order + 1 samples. masked below that, where no derivative of that order exists.

This trades accuracy order for coverage — a five-point scheme becomes three-point in a strait three cells wide — so it is named rather than reached by fallback. Ask for it when a value everywhere matters more than a uniform order.

Under this policy nodes is a ceiling, not a demand, and that applies to the end of the axis as well as the end of a run: an axis with fewer than nodes samples uses as many as it has instead of raising, and one with fewer than order + 1 is masked throughout. A single-latitude strip, a two-level column and a one-cell-wide channel are ordinary grids, and asking for "second order where the axis allows it" should not require the caller to clamp nodes themselves. The other two policies keep the error, since neither claims to degrade.

ShiftWithinRun()

Shift the stencil to fit inside the run of active samples containing the cell, keeping the full node count — the same thing the stencil already does at the end of a bounded axis, with the end of the active run as the boundary. masked only where the run is shorter than nodes.

The accuracy order is therefore the same everywhere a value is written, which is the property fd_weights exists to preserve. On a run of at least nodes active samples the weights are identical to the unmasked ones, so the interior of an active region is bit-for-bit unchanged.

StencilScratch{T}

The working buffers a degrading apply_stencil! needs to rebuild a window at a mask edge: the Fornberg table and the node list. Build one with stencil_scratch.

One per task, exactly as Connectivity.ball_scratch is — the buffers are written per cell, so chunks cannot share them. A threaded backend therefore allocates its own set per chunk and ignores one passed here.

_stencil_sweep_host!(out, field, indices, weights, mask, masked, Val(dim), nodes, Val(N)) -> out

The host sweep. The index-parallel form exists so one body serves a device launch; on the host it is the wrong shape, and three things it cannot express are worth ~6.6× together:

• iterate the Cartesian range directly, rather than recovering an index per cell from a linear one;
• split the nest at dim, so the stencil row — which depends only on the index along dim — is hoisted out of the contiguous inner loop whenever the differenced direction is not the fastest varying one;
• carry the node count in the type, so the innermost loop has a known trip count, unrolls, and the weights reach registers instead of being re-loaded per node.

The arithmetic and its order are identical to _stencil_cell!, so the two paths agree bit for bit.

_sympinv2(a, b, c, tol) -> (p11, p12, p22)

Pseudo-inverse of the symmetric 2×2 [a b; b c], dropping any eigendirection whose eigenvalue is below tol. Closed form: a 2×2 symmetric eigenproblem has one.

Dropping rather than regularising is the point. A stencil that carries no information along some tangent direction — every neighbour on one line, which happens at a boundary and beside a mask — leaves A singular in that direction, and the gradient there is not determined by the data. Inverting a nudged matrix would answer anyway, with a number governed by the nudge.

_window_start(i, k, lo, hi) -> Int

First index of a k-node window centred on i and shifted to fit inside [lo, hi]. The whole-axis case is lo = 1, hi = n; the masked case is the same expression with the bounds of the active run, which is why both share this.

apply_stencil!(out, field, x, indices, weights, dim; order=1, period=nothing, mask=nothing,
               masked=zero, policy=BlankMasked(), backend=nothing) -> out

Apply a table built by axis_stencils and keep the axis, so any mask policy works.

The table depends on the axis and not on the field, so a caller differencing many fields along one direction should build it once. The bare (indices, weights) form cannot degrade at a mask edge — that needs the axis to rebuild a window from — so it accepts only BlankMasked; this form takes both and serves every policy.

The split is the one the degrade path already makes internally: the precomputed row is used wherever the window is intact, which is every cell away from a mask, and the axis is touched only where a window is actually rebuilt.

Building the table is O(n) against an O(n²) apply, so holding it matters most on small grids — 2.4–13× at n = 48, 10–40% at n = 256, amortized away by n = 1024. The allocation it avoids is there at every size: 49 600 bytes per call at n = 1024.

apply_stencil!(out, field, x, dim; order=1, nodes=order+1, period=nothing,
               mask=nothing, masked=zero) -> out
apply_stencil!(out, field, indices, weights, dim; mask=nothing, masked=zero) -> out

Apply a weight set along direction dim of field, writing out[I] = Σ_q weights[I[dim], q] · field[…, indices[I[dim], q], …].

This is the one field-touching operation here, and it is here because every convention it needs is already settled elsewhere in the package rather than being the caller's to choose: the result sits at the same location as the input, so there is no staggering decision; the stencil shifts inward at a bounded end and wraps on a periodic one, which is fd_weights's stated boundary behaviour and removes any need for a halo. What is not here is anything that does need those choices — a staggered difference, or a multi-direction operator like a divergence or a curl, which additionally needs a result location and a boundary-condition policy.

Pass the axis and an order to have the weights built for you, or precomputed indices/weights from axis_stencils to reuse them across many fields.

With a mask, a cell is written as masked when it is inactive or when its stencil reads an inactive cell — the derivative there is not determined by the active data, so it is not invented. out and field may not alias.

axis_stencils(x, order, nodes; period=nothing) -> (indices, weights)

The order-th derivative's fd_weights at every sample of axis x, as two n × nodes matrices: the axis indices each sample reads, and the weight on each.

One row per sample, so a stretched axis costs nothing extra downstream — the varying weights are already here. Built once and reused by apply_stencil!.

period === nothing shifts the stencil inward at the two ends, exactly as the single-sample fd_weights does. Given a period the stencil stays centred everywhere and wraps, with the wrapped samples' coordinates carried across the seam so the spacing there is the true one.

cell_width(x, i, period=nothing) -> width

The coordinate width of cell i of an axis of cell centres x: the centred width (|h_m| + |h_p|)/2 at an interior cell — and, given a period, at the wrapped boundary too — the one-sided gap to the single neighbour at a genuinely non-periodic boundary, and 1 for a length-1 axis. On a uniform axis every width is the constant step.

Equivalently abs(faces(x)[i+1] - faces(x)[i]), which is what it means: faces places a boundary midway between neighbouring centres, so the width between them is the average of the two adjacent gaps. This form is the one to call per cell, since faces materializes the whole axis.

A width is a physical measure and so is non-negative however the axis is stored — increasing or decreasing, as a dataset holding latitude, depth or pressure levels top-down would. This is the one place that turns a spacing into a length/area/volume contribution, so it is where the abs belongs; local_spacing itself keeps the sign.

A length-1 axis contributes the multiplicative identity, not zero, to a measure that is a product of per-axis widths (Cartesian Δx·Δy), so a degenerate direction reduces an area to a length rather than collapsing the product. The spherical R²cosφ·Δλ·Δφ measure is not a plain product and handles its own singleton case, in the Grids.StructuredGrid constructor.

Grids.cell_width is this on a grid direction, and Grids.cell_widths the whole axis at once.

cell_widths(x, period=nothing) -> AbstractVector

cell_width at every index of an axis at once, for a caller that wants the whole profile rather than one cell — the coordinate widths a separable measure or a flux divergence is weighted by.

A uniform axis gets an Axes.ConstantVector: one number and a length, since every one of its cells has the same width, so nothing is materialized. Anything else is built densely into the same kind of storage as x, by broadcasts over views rather than a scalar loop, so a device-resident axis is widened in place.

Grids.cell_widths is this on a grid direction, taking the period from the grid itself.

centers(f) -> AbstractVector

The N cell centres of an axis of N+1 cell boundaries: the midpoint of each pair.

It inverts faces exactly on a uniform axis. On a stretched one it does not: faces places a boundary midway between two centres, and re-midpointing those boundaries averages neighbouring cell widths rather than recovering the centre. The two are inverse only where the widths are constant.

derivative!(out, field, grid, dim; order=1, nodes=order+1, policy=BlankMasked(),
            masked=zero, active_only=true, backend=nothing) -> out

The derivative with respect to distance along direction dim, rather than with respect to the coordinate: apply_stencil! divided by the metric factor,

∂f/∂sᵈ = (1/hᵈ) · ∂f/∂ξᵈ,    hᵈ = Geometry.scale_factors(geo, p)[d]

On a Cartesian metric every hᵈ is 1 and this is apply_stencil! exactly, at no cost. Anywhere else it is the derivative a physical law is written in, and assembling it from the parts was the one thing every geometry-aware caller had to add — including the coordinate singularity below, which is not theirs to get right.

Where the metric degenerates the derivative does not exist, and masked is written rather than a number invented. Longitude at a pole is the case: h_λ = R cos φ → 0, so 1/h_λ diverges. The test is relative to the geometry's own size and to the precision, |h| ≤ L·√eps(T) — an absolute threshold cannot be right for both, since 1e-12 is below eps(Float32) and in Float32 cos(Float32(π/2)) ≈ -4.4e-8, so a pole row would quietly receive a large finite number instead.

No scale factor in this package depends on longitude, so hᵈ is constant along the first axis whichever direction is differenced. The scaling is applied once per remaining index and swept along that contiguous axis. (It is not generally constant along the differenced direction — on a spheroid h_φ = M(φ) varies with φ — so it is not hoisted that way.)

A divergence or a curl is still the caller's to assemble, needing a result location and a boundary policy this does not choose. Note the flux form when doing so: on a sphere

∇·u = (1/(R cos φ)) [ ∂u_λ/∂λ + ∂(u_φ cos φ)/∂φ ]

so the second term differentiates u_φ cos φ, not u_φ; taking two physical derivatives and adding them is a different, wrong expression.

faces(x) -> AbstractVector

The N+1 cell boundaries of an axis of N cell centres: midpoints of neighbouring centres, with the outermost two extrapolated a half-cell beyond the end centres.

A uniform axis stays uniform — its faces are another Axes.UniformAxis, offset by half a cell — so nothing about the axis's spacing guarantee is lost.

fd_weights!(w, c, nodes, x₀, order) -> w

fd_weights into caller buffers: w holds the length(nodes) weights and c is the length(nodes) × (order+1) recursion table. Both are overwritten.

The allocating form is one of these per call, and a stencil is built once per sample of an axis, so a 4096-sample axis costs ~8000 allocations without this. The degrade path in apply_stencil! needs one per cell near a mask edge, which is the reason it exists.

fd_weights(x, i, order, nodes) -> (indices, weights)

Weights for the order-th derivative at sample i of axis x, using nodes of its samples.

The stencil is centred on i where the axis allows and shifted inward at a boundary, so the accuracy order is the same everywhere — a clipped stencil would silently drop to first order at the two ends. Built on the arbitrary-node form above, so a stretched axis costs nothing extra.

fd_weights(nodes, x₀, order) -> Vector

Finite-difference weights approximating the order-th derivative at x₀ from the values at nodes, by the recursion of Fornberg (1988), Math. Comp. 51, 699–706.

One recursion covers every case: any derivative order, any node count (hence any order of accuracy), any evaluation point — inside the node set or outside it — and arbitrarily spaced nodes. With m nodes the result is exact for polynomials of degree m-1, so accuracy order m - order.

sum(w .* f.(nodes)) is then the derivative estimate. This returns the weights only; applying them to a field is the caller's.

fd_weights([0.0, 1.0, 2.0], 1.0, 1)   # ≈ [-0.5, 0.0, 0.5], the centred first difference

gradient!(g1, g2, field, plan) -> (g1, g2)

Apply a GradientPlan: the two tangent components of ∇field at every cell, written into g1 and g2. One dot product per cell over its neighbours, allocating nothing.

field, g1 and g2 are indexed linearly, so an N-D array of the grid's shape works as is. The components are named by plan.names — (:λ, :φ) on a sphere, (:x, :y) on a plane — and are per unit distance, the tangent plane being metric already.

interpolate(field, grid, p; policy=BlankMasked(), masked=NaN, …) -> value

The value of field at the coordinate p, which is the question observational data asks: a station, a float or a ship track has a coordinate, not a cell index.

interpolation_weights gives this along one axis, and nothing composed them, so a caller had to build the tensor product themselves on a rectilinear grid and had nothing at all on the others.

• StructuredGrid — multilinear, the tensor product of the per-axis weights. A periodic direction interpolates across its seam rather than clamping at the last sample.
• CurvilinearGrid, UnstructuredGrid — a weighted least-squares plane fitted to the k nearest cells in the tangent plane at p, which is exact for a linear field and reproduces a cell's own value at its centre. Falls back to the weighted mean where the fit is rank deficient, that being the part of it the data still determines.

p may be written any way a point is accepted elsewhere.

The mask policies say what an inactive contributor means, as they do for a stencil: BlankMasked — the default — returns masked if any contributor is inactive, and ReduceInRun renormalizes over the active ones. ShiftWithinRun has no meaning here, there being no window to shift, and says so.

interpolation_weights(x, v) -> (i, w)

Linear interpolation on axis x at coordinate v, as the left sample index i and the weight pair w = (w_i, w_{i+1}) with w_i + w_{i+1} == 1.

Weights only: applying them to a field is the caller's loop, and the field's layout is not this module's business. Outside the axis the nearest end is used with weight 1, so the result is a clamp rather than an extrapolation.

jacobian(geometry, point) -> Real

∏ hᵈ from scale_factors: the volume element per unit coordinate volume at point.

lagrange_weights(x, v, nodes) -> (indices, weights)

Lagrange interpolation weights of length(nodes) points on axis x at coordinate v, exact for polynomials up to degree nodes-1 and valid on an arbitrarily spaced axis.

The stencil is centred on v as far as the axis allows and shifted inward at the ends, so the node count is honoured everywhere rather than degrading near a boundary.

local_spacing(x, i, period=nothing) -> (h_m, h_p)

The one-sided coordinate gaps around index i of an axis x: h_m = x[i]-x[i-1] and h_p = x[i+1]-x[i].

This is the primitive a finite-difference operator assembled at the call site is built from — the h_m, h_p of Geometry.nonuniform_first_derivative, and what a staggered difference, a divergence or a curl needs from the grid. Always a scalar subtraction of two already-stored elements, never a heap allocation, so it is safe to call per grid point in a hot loop. On an axis whose spacing is known from its type there is not even a subtraction.

The gaps are signed, so a descending axis reports negative gaps and a derivative stencil keeps the index-versus-coordinate direction: composed with nonuniform_first_derivative this gives the same derivative with respect to the coordinate whichever way the axis is stored. cell_width is where the sign is dropped, being a length rather than a difference.

period, if given (e.g. 2π for a periodic longitude axis), makes the boundary gaps wrap instead of vanishing: at i == 1, h_m is the gap to the unwrapped previous point x[n]-period; at i == n, h_p is the gap to the unwrapped next point x[1]+period. Pass nothing (the default) for a non-periodic axis, where a boundary gap is zero and the caller falls back to a one-sided stencil.

Grids.local_spacing is this on a grid direction, taking the period from the grid itself.

locate(x, v) -> Int

The index of the cell of axis x that contains coordinate v, or 0 when v lies outside.

Cells are the intervals between the axis's faces, so cell i holds f[i] ≤ v < f[i+1] (the last cell includes its far face). O(1) on a uniform axis, from the closed form; O(log n) on a stretched one, by bisection. Both storage orders work.

metric_floor(geometry) -> T

The magnitude below which a scale factor is treated as degenerate: L·√eps(T) for a curved geometry of size L, and 0 for a Cartesian one, whose metric never degenerates.

Relative to both the geometry's size and the element type on purpose — see derivative! for what an absolute constant does in Float32.

nearest_index(x, v) -> Int

The index of the axis sample closest to v, clamped into range. Unlike locate this always returns a valid index, since a nearest sample exists for any v.

Exact ties go to the LOWER index, and the uniform closed form and the general bisection agree on that, so the two paths never disagree. O(1) on a uniform axis, O(log n) on a stretched one.

nodes(x, loc) -> AbstractVector

An axis's sample positions at location loc: x itself at Center, and its faces at Face.

scale_factors(geometry, point) -> NTuple

The metric scale factors hᵈ at point: the physical length of a unit step in each coordinate direction. Cartesian gives 1 in every direction; spherical gives (R cosφ, R) on the surface and (r cosφ, r, 1) with a radius direction.

These are what turns a coordinate derivative into a physical one — ∂/∂sᵈ = (1/hᵈ)·∂/∂ξᵈ — so a divergence or a curl is assembled from these plus fd_weights without this module having to choose a staggering or a boundary condition.

stencil_scratch([T = Float64], order, nodes) -> StencilScratch

Buffers for the degrade path, so a caller taking many derivatives on a masked grid does not allocate them per call. Without one a degrading call allocates a few hundred bytes each time — O(1) in the grid, but per call, so a flux computation taking nine derivatives pays it nine times.

Only the degrading policies need it. An unmasked grid, and any grid under BlankMasked, never rebuilds a window and allocates nothing regardless.