Discretization

AbstractLocation

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

AbstractStencilPlan{T}

A finite-difference weight set for one axis. UniformStencilPlan or TabulatedStencilPlan, chosen by the axis's Axes.spacing_trait.

Build one with stencil_plan and apply it with Operators.apply_stencil!.

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.

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.

TabulatedStencilPlan{T,K,MI,MW}

The weight set of a stretched axis: one row of node indices and weights per sample, as the n × K pair of matrices axis_stencils returns.

UniformStencilPlan{T,K,MI,MW}

The weight set of a uniformly spaced axis: one centred row of K weights, plus the rows at the two ends where a bounded window shifts inward.

Every interior sample reads field[j - left + q - 1] weighted by weights[q], the same K numbers throughout, so they reach the sweep as a tuple in registers with no per-cell gather. Storage is O(K²), independent of the axis length. A wrapping axis has no shifted rows at all, every sample being centred.

left and right count the nodes before and after the sample and are equal only for odd K: a 4-node window reads one sample back and two forward, so two rows shift at the top end and one at the bottom.

The weights are built from the spacing. Mathematically they are the row axis_stencils holds, and the table carries per-row round-off from reconstructing each window's coordinates, at ~1e-14 relative on a 5-node first derivative; one exact row applied everywhere carries none.

_uniform_row(T, h, s, Val(K), order) -> NTuple{K,T}

The weights of a row whose first node sits s steps of size h from the sample it is evaluated at.

s = -left is the centred row. Only the offsets enter, so this is the exact weight set for every sample sharing that offset — on a uniform axis, all of them away from a shifted end.

_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 passes the bounds of the active run.

axis_length(plan) -> Int

The number of samples on the axis plan describes. A field differenced with it has to have that many along the differenced direction.

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

axis_stencils into two caller-owned n × nodes matrices.

The table is O(n·k), so a caller rebuilding it — at each step of a stretching mesh, say — owns the matrices and passes them in. Pass a stencil_scratch as well and the call allocates nothing: the remaining buffers are the Fornberg table and the node list, O(k) each, allocated per call without one.

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]): faces places a boundary midway between neighbouring centres, so the width between them is the average of the two adjacent gaps. Call this form per cell; faces materializes the whole axis.

A width is a physical measure, so it is non-negative however the axis is stored — increasing or decreasing, as a dataset holding latitude, depth or pressure levels top-down is. This is where a spacing becomes a length, area or volume contribution, and where the abs is applied; local_spacing keeps the sign.

A length-1 axis contributes the multiplicative identity to a measure that is a product of per-axis widths (Cartesian Δx·Δy), so a degenerate direction reduces an area to a length. 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: 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, every one of its cells having the same width. Anything else is built densely into the same kind of storage as x, by broadcasts over views, so a device-resident axis is widened on its device.

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

centers!(out, f) -> out

centers into a caller-owned vector of length N, for an axis of N+1 faces.

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. The two are inverse only where the widths are constant.

derivative_order(plan) -> Int

Which derivative plan's weights compute.

faces!(out, x) -> out

faces into a caller-owned vector of length N+1.

The form to use when staggering repeatedly — a moving mesh, a column solver stepping in time — where the N+1 faces are otherwise allocated per call.

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 allocates two arrays per call, and a stencil is built once per sample of an axis. The degrade path in apply_stencil! rebuilds one per cell near a mask edge, and holds these buffers across all of them.

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, holding the accuracy order at the two ends as well as the interior. 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

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. Outside the axis the nearest end is used with weight 1, so a coordinate beyond the ends clamps.

jacobian(geometry, point) -> Real

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

lagrange_weights!(w, x, v, nodes) -> (indices, w)

lagrange_weights into a caller-owned vector of length nodes.

The indices come back as a range, a stack value, so with w supplied the call allocates nothing. This is the form to use when interpolating many points along one axis.

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 holds at a boundary as well as in the interior.

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 drops the sign, a width being a length.

period, if given (e.g. 2π for a periodic longitude axis), makes the boundary gaps wrap: 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.

It scales with both the geometry's size and the element type. At Float32 on Earth's radius, cos(Float32(π/2)) ≈ -4.4e-8 puts h_λ at the pole around 0.28 m, above any fixed small constant.

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.

nnodes(plan) -> Int

How many samples each row of plan reads.

nodes(x, loc) -> AbstractVector

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

plan_row(plan, j) -> (nodes, weights)

Row j of plan, as two K-tuples: the axis indices that sample reads and the weight on each.

The generic accessor, for a path that wants a row whatever form the plan took. The uniform sweep does not go through it — the whole point there is that the interior row is not fetched per cell.

scale_factors(geometry, point) -> NTuple
scale_factors(geometry, point, Val(N)) -> NTuple{N}

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.

They turn a coordinate derivative into a physical one, ∂/∂sᵈ = (1/hᵈ)·∂/∂ξᵈ, so a divergence or a curl is assembled from these plus fd_weights at the call site, where the staggering and the boundary condition are chosen.

Val(N) names how many components point has, giving a concrete NTuple{N} from a point whose length is a runtime value.

stencil_plan(x, order, nodes; period = nothing) -> AbstractStencilPlan

The order-th derivative's weights along axis x, held for reuse.

A uniform axis gives a UniformStencilPlan and a stretched one a TabulatedStencilPlan; which it is comes from the axis's type through Axes.spacing_trait, never from inspecting its values.

A wrapping uniform axis takes the uniform form only when the period given is the one the axis's own spacing implies, n·h. This checks the argument against the axis type, and reads no coordinate: any other period puts a seam gap other than h at the wrap, which the tabulated form describes.

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.