API

Types

`AbstractBC`

AbstractBC

Supertype for boundary conditions attached to grid faces.

apply_bc! enforces only the homogeneous part of a boundary condition, so every linear operator satisfies L(0) = 0; inhomogeneous boundary data enters a solve separately through boundary_rhs.

`AbstractField`

AbstractField

Supertype for fields the operator algebra and solver boundary act on. Field is the single-grid case; BlockField is the block-structured (forest) case. The shared contract is similar/zero_ghosts! and the flat-vector boundary (flatten/flat_to_interior!/interior_to_flat!/flat_length); operators are otherwise written against single-grid Fields and reused per block.

`AbstractGrid`

AbstractGrid{N}

Supertype for all N-dimensional grids. Operators are authored once against the grid interface (spacing, interior, halo_update!, …); single-device, distributed, and adaptive grids differ only in what the grid object is and what halo_update! does.

`AbstractOperator`

AbstractOperator

Supertype of all matrix-free operators. A concrete operator provides:

apply!(y, L, x, grid, α, β) — the action y = α·L(x) + β·y, writing the interior of y (ghost layers of x are scratch: filled by halo_update! and apply_bc! before the stencil reads neighbors)
traits islinear, isconstant, isselfadjoint, isdiagonal (default false — leaves opt in)
for linear operators: an adjoint via adjoint_operator and apply_adjoint!

Operators compose lazily with +, -, *, and scalar scaling, and are applied with apply, L(x), or L * x. For Krylov solvers, wrap with prepare to get an allocation-free mul!.

`Added`

Added(a, b)

Lazy operator sum a + b. Produced by a + b and a - b.

`AdjointOp`

AdjointOp(L)

Lazy adjoint wrapper: applying it calls apply_adjoint! of the wrapped operator. Produced by adjoint(L) for leaves without a cheaper adjoint expression.

`Advection`

Advection(grid, velocity)

Matrix-free advection leaf v·∇. Linear when the velocity is a prescribed Field (passive transport — a field-valued coefficient does not make an operator nonlinear); nonlinear when the velocity is SelfAdvection (u·∇u). Construct with advection.

`BlockField`

BlockField{L}(blocks, grid)
BlockField(blocks, grid)

Field over a BlockForest: one halo-padded array per leaf block, indexed in the forest’s Morton (storage) order, with location trait L (default Center). The vector-of-blocks layout is the simplest correct storage; a packed contiguous buffer can replace it later without touching operators. Block storage is tied to the leaf set at allocation time — the forest’s regrid generation is stamped into the field, and any use after a refine!/coarsen! that changed the leaf set throws; allocate a fresh field after a regrid.

`BlockForest`

BlockForest(base::CartesianGrid; blocksize, maxlevel) -> BlockForest

Block-structured adaptive grid: a Forest of fixed-size leaf-blocks laid over the physical domain of base. Each leaf is an ordinary CartesianGrid of blocksize cells with the same halo as base; refining a block replaces it with 2ᴺ children at half the spacing. Leaf grids carry the Interface boundary on every face; inter-block ghosts are filled by halo_update! and physical domain faces by the forest-level apply_bc! face pass from base’s boundary conditions (kept on the forest).

base ncells must be divisible by blocksize (the quotient is the root tiling). Operators built on a BlockForest run the existing per-CartesianGrid stencil code unchanged on every leaf.

Keyword Arguments

blocksize::NTuple{N,Int}: cells per block per dimension
maxlevel::Int: maximum refinement level (root blocks are level 0)

Examples

base   = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16);
                       bc=((Dirichlet(), Dirichlet()), (Dirichlet(), Dirichlet())))
forest = BlockForest(base; blocksize=(8, 8), maxlevel=4)
refine!(forest, x -> sum(abs2, x .- 0.5) < 0.05)   # refine a disc near the center

See also: refine!, coarsen!, leaves.

`CartesianGrid`

CartesianGrid(extent, ncells; bc, halo, device)

Uniform cell-centered Cartesian grid.

Arguments

extent::NTuple{N,Tuple{T,T}}: physical (min, max) per dimension
ncells::NTuple{N,Int}: interior cell counts per dimension

Keyword Arguments

bc: per-dimension (low, high) boundary-condition pairs (default: homogeneous Dirichlet on every face)
halo: ghost-layer width per dimension (default: 1 per dimension)
device: KernelAbstractions backend used for field allocation (default: CPU())

Cell spacing is derived from extent and ncells. Cell centers sit at min + (i - 1/2)Δ for interior cell i.

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))                       # 1-D, Dirichlet
g = CartesianGrid(((0.0, 1.0), (0.0, 2.0)), (32, 64);
                  bc = ((Periodic(), Periodic()),
                        (Neumann(), Neumann())))              # 2-D, mixed BCs

`Center`

Center

Location trait for cell-centered (collocated) fields — the only location in v1. Staggered layouts add face locations as new trait types.

`Chebyshev`

Chebyshev(degree=3)

Chebyshev polynomial smoother configuration for MultigridPreconditioner: a fixed degree-degree polynomial in D⁻¹·A targeting the upper part of the spectrum [λmax/4, λmax], with λmax estimated by power iteration at construction. Needs only repeated apply! plus operator_diagonal; fixed coefficients keep the V-cycle symmetric.

`Composed`

Composed(a, b)

Lazy operator composition (a ∘ b)(x) = a(b(x)). Produced by a * b.

`Derivative`

Derivative(grid, dim, order)

Matrix-free partial-derivative leaf ∂^order/∂x_dim^order. Construct with derivative.

`Dirichlet`

Dirichlet(value=0)

Dirichlet boundary condition u = value at the face. The homogeneous part fills ghost cells by antisymmetric mirroring (ghost = -interior); value contributes only through boundary_rhs.

`Divergence`

Divergence(grid)

Matrix-free divergence leaf — a rank-changer mapping an SVector{N}-valued field to a scalar field. Construct with divergence.

`Field`

Field(data, grid)
Field{L}(data, grid)

Halo-padded field on grid with location trait L (default Center). data must have size padded_size(grid). The element type carries the tensor rank: a scalar field stores numbers, a vector field stores SVectors — see scalar_field and vector_field.

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))
u = Field(zeros(padded_size(g)), g)

`Gradient`

Gradient(grid)

Matrix-free gradient leaf — a rank-changer mapping a scalar field to an SVector{N}-valued field. Construct with gradient.

`IdentityOp`

IdentityOp()

Grid-free identity operator. Resolves its size and element type from a composition sibling or the field it is applied to. Construct with identity_op.

`Jacobi`

Jacobi(ω=2//3, sweeps=2)

Weighted-Jacobi smoother configuration for MultigridPreconditioner: x ← x + ω·D⁻¹·(b - A·x) repeated sweeps times, with D extracted once per level via operator_diagonal. Applied identically before and after each coarse-grid correction, which keeps the V-cycle symmetric.

`Laplacian`

Laplacian(grid)

Matrix-free Laplacian (∇²) leaf bound to grid. Construct with laplacian.

`LinearizedOp`

LinearizedOp

Matrix-free Jacobian J = ∂F/∂u of an operator F, frozen at a state u₀. Its action is the Jacobian–vector product J·v = ∂/∂ε F(u₀ + εv)|₀, evaluated by a central finite-difference JVP. Linear by construction — the operator handed to Krylov in Jacobian-free Newton–Krylov. Build with linearize; refresh the frozen state in place with linearize!.

`MultigridPreconditioner`

MultigridPreconditioner(L; smoother=Jacobi(), cycle=:V, levels=:auto)

Geometric-multigrid V-cycle preconditioner for a linear operator on a uniform CartesianGrid. Each mul!(z, M, r) runs one V-cycle from a zero initial guess, so M acts as a fixed linear approximation of L⁻¹ — pass it as the preconditioner M to a Krylov solver.

Levels are built by coarsen-ing the grid (levels=:auto stops at ≤64 DOFs, odd cell counts, or 10 levels; pass an Integer to force a depth); each level rediscretizes L on its grid, transfers residuals/corrections with restriction/prolongation, smooths with smoother (Jacobi or Chebyshev) symmetrically before and after the coarse correction, and solves the coarsest level with a dense LU factored at construction. For a symmetric L the cycle is symmetric, so Krylov.cg applies.

Stateful and single-threaded, like prepare — construct one per concurrent solve.

Examples

g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
L = -laplacian(g)
A = prepare(L)
M = MultigridPreconditioner(L)
b = flatten(set!(scalar_field(g), x -> sinpi(x[1]) * sinpi(x[2])))
u, stats = Krylov.cg(A, b; M)

`MultigridSolver`

MultigridSolver(L; smoother=Jacobi(), cycle=:V, levels=:auto, maxiter=200)

Standalone multigrid solver: the stationary iteration u ← u + M·(b - A·u) with M a MultigridPreconditioner of L and A = prepare(L). Run it with solve.

`Neumann`

Neumann(flux=0)

Neumann boundary condition ∂u/∂n = flux at the face (outward normal). The homogeneous part fills ghost cells by symmetric mirroring (ghost = interior); flux contributes only through boundary_rhs.

`PackedBlockField`

PackedBlockField{L}(data, levels, grid)

Packed twin of BlockField: every leaf’s halo-padded block stored in one contiguous (blocksize .+ 2halo ..., nleaves) array (data, leaves in Morton order along the trailing dimension), with the per-leaf refinement levels in a device-resident vector (levels) — the geometry SoA forest-native kernels read; spacing derives from the forest’s root spacing and the level. Behind the shared AbstractBlockField interface every operator runs on packed storage via the per-leaf fallback sweep (each leaf is a trailing-dim view), and operators with a forest-native kernel sweep all leaves in a single launch. Convert with pack / unpack; BlockField remains the reference layout for AD, Reactant, and regrid! (re-pack after a regrid).

`Periodic`

Periodic()

Periodic boundary condition. Must be paired on both faces of a dimension.

`PreparedForest`

PreparedForest

Block-forest counterpart of PreparedOperator: the flat mul!/size/ eltype boundary for an operator over a BlockForest. Beyond the padded scratch fields it carries a scratch BlockField for accumulating adjoint sweeps (adjscratch), and prepare warms the grid’s per-generation halo-exchange schedule (see halo_update!), so the inter-block ghost copies and the physical-BC face passes run as flat descriptor loops with no topology queries. Every leaf grid is one concrete all-Interface type, so the leaf sweep is type-stable and rebuilding a leaf grid allocates nothing; the residual cost is the per-leaf stencil apply! itself (allocation-free only when inlined into a single mul!); see the allocation-free-kernel follow-up. The prepared operator is tied to the forest’s regrid generation; a refine!/coarsen!/balance! after prepare invalidates it and mul! throws — re-run prepare on the new forest.

`PreparedOperator`

PreparedOperator

A linear operator bound to pre-allocated scratch fields, exposing the flat mul!/size/eltype interface Krylov solvers need. Built with prepare; after warm-up, mul! runs without steady-state allocations on single-grid Fields. The BlockField path uses the sibling PreparedForest.

The flat vectors span interior DOFs only — ghost cells are determined by boundary conditions and halo_update!, never solver unknowns. Each mul! copies the flat vector into a halo-padded scratch field, applies the operator, and copies the interior back out fused with the α/β axpby, so the solver’s vectors are never mutated by halo or BC fills.

mul! writes only the interior of xpad, deliberately leaving ghost cells as it found them. The distributed path (prepare_distributed) depends on that same discipline throughout: it drives the operator tree itself so it can fill Interface ghost slabs from a neighbor exchange between applies, and anything that zeroed ghosts on the way in would silently discard exchanged data.

`Prolongation`

Prolongation(coarse, fine)

Coarse-to-fine transfer operator on a 2:1 cell-centered grid pair: per-dim linear interpolation blending each fine cell’s coarse parent (weight 3/4) with the neighbor on the child’s side (weight 1/4). Boundary children read the coarse homogeneous ghost fill, so the interpolant is consistent with the grid’s boundary conditions (Dirichlet wall children get u₁/2, Neumann get u₁). Adjoint of Restriction up to scaling: Pᵀ = 2^N · R. Construct with prolongation.

`Restriction`

Restriction(fine, coarse)

Fine-to-coarse transfer operator on a 2:1 cell-centered grid pair, defined as the scaled transpose of Prolongation: R = 2⁻ᴺ·Pᵀ — full weighting, per-dim interior weights (1/8, 3/8, 3/8, 1/8). Sharing P’s kernels makes the adjoint identity exact by construction. Input ghost layers are scratch (zeroed before the gather), so only fine interior values contribute and the affine boundary lift is identically zero. Construct with restriction.

`Scaled`

Scaled(op, α)

Lazy scalar multiple α·op. Produced by α * L, L * α, -L, and L / α.

`ScalingOp`

ScalingOp(coeff)

Pointwise multiplication by a coefficient — a Number or a scalar-eltype Field κ(x). The diagonal, parameter-carrying leaf: its coefficient field is a differentiable operator parameter (material coefficients for inverse problems). Construct with scaling.

`SelfAdvection`

SelfAdvection()

Velocity marker for advection selecting self-advection u·∇u: the velocity is the advected state itself, making the operator nonlinear.

Functions

`advection`

advection(g::AbstractGrid, velocity::Field) -> Advection
advection(g::BlockForest, velocity::AbstractBlockField) -> Advection
advection(g::AbstractGrid, ::SelfAdvection) -> Advection

Build a matrix-free advection operator v·∇ bound to g, with second-order central differences. A prescribed velocity must be an SVector-valued field with one component per grid dimension — on a BlockForest a block field on that same forest (matching the applied field’s layout lets the forest-native kernel engage); the resulting operator is linear in the advected input and acts componentwise on scalar or vector inputs. With SelfAdvection the operator computes u·∇u of a vector field and is nonlinear — see linearize for its Jacobian.

Examples

g = CartesianGrid(((0.0, 2π), (0.0, 2π)), (32, 32);
                  bc=((Periodic(), Periodic()), (Periodic(), Periodic())))
v = set!(vector_field(g), x -> SVector(1.0, 0.0))
A = advection(g, v)                  # linear: passive transport by v
B = advection(g, SelfAdvection())    # nonlinear: u·∇u

`apply`

apply(L::AbstractOperator, x::AbstractField) -> AbstractField

Allocating application L(x). The pure path used by autodiff; hot loops should use prepare + mul! or in-place apply! instead.

Examples

g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
u = set!(scalar_field(g), x -> sin(x[1]))
Δu = apply(laplacian(g), u)        # equivalently laplacian(g)(u) or laplacian(g) * u

`apply!` {#apply!}

apply!(y::AbstractField, L::AbstractOperator, x::AbstractField, g::AbstractGrid, α=true, β=false) -> y

Apply the operator in place: y = α·L(x) + β·y on the interior of y. Ghost layers of x are treated as scratch (overwritten with halo/BC fills); ghost layers of y are left untouched by the accumulating form.

`apply_adjoint!` {#apply_adjoint!}

apply_adjoint!(x̄::AbstractField, L::AbstractOperator, ȳ::AbstractField, g::AbstractGrid, α=true, β=false) -> x̄

Apply the adjoint of a linear operator in place: x̄ = α·Lᵀ(ȳ) + β·x̄. Ghost layers of ȳ are treated as scratch (zeroed — only interior values are adjoint inputs, matching the flat Krylov boundary).

apply_adjoint!(x̄::AbstractBlockField, L, ȳ::AbstractBlockField, g::BlockForest, α=true, β=false) -> x̄

Adjoint action over the forest. For a self-adjoint operator (e.g. the Laplacian on a uniform forest, whose same-level halo copy couples both neighbors symmetrically) this is the forward action. Otherwise it is the exact transpose of stencil ∘ BC fill ∘ halo exchange: per-leaf stencil-transpose gathers (leaving all ghost cotangents in place), the forest-level fold_bc! folding physical ghosts, then halo_update_adjoint! folding interface ghosts into the neighbor interiors.

`apply_bc!` {#apply_bc!}

apply_bc!(data, g::AbstractGrid) -> data

Fill all ghost layers of the halo-padded array data with the homogeneous part of the grid’s boundary conditions (periodic wrap, Dirichlet antisymmetric mirror, Neumann symmetric mirror). Dimensions are filled in order 1:N, so corner ghosts are consistent ghost-of-ghost values.

`balance!` {#balance!}

balance!(forest::Forest) -> forest

Enforce the 2:1 balance invariant across faces: no leaf has a face-adjacent leaf more than one level finer. Iterates to a fixpoint, refining the coarser leaf of any violating pair. This is what reduces every coarse–fine interface to the three tractable cases (same level / one-coarser / one-finer).

balance!(bf::BlockForest) -> bf

Re-establish the 2:1 balance invariant across leaf faces.

`boundary_conditions`

boundary_conditions(g::AbstractGrid) -> NTuple{N,Tuple}

Per-dimension (low, high) boundary-condition pairs.

`boundary_rhs`

boundary_rhs(L::AbstractOperator, g::AbstractGrid) -> AbstractField
boundary_rhs(L::AbstractOperator, x_proto::AbstractField) -> AbstractField

Boundary lift of the affine split L_full(x) = L(x) + b: the contribution of inhomogeneous boundary data (Dirichlet values, Neumann fluxes) that apply! deliberately omits so that L stays linear (L(0) = 0). Assemble once per solve and fold into the right-hand side: the discrete problem L_full(u) = f becomes L·u = f - b. The grid form assumes a scalar input field; pass a prototype field for vector inputs. On a BlockForest the lift is assembled by the forest-level inhomogeneous face pass — only physical domain faces contribute; Interface faces stay homogeneous.

Examples

g = CartesianGrid(((0.0, 1.0),), (64,); bc=((Dirichlet(1.0), Dirichlet(2.0)),))
L = laplacian(g)
b = boundary_rhs(L, g)
rhs = flatten(f) .- flatten(b)        # solve  prepare(L) \ rhs  with Krylov

`cell_center`

cell_center(g::AbstractGrid, I::CartesianIndex) -> SVector

Physical coordinates of the center of cell I (in halo-padded index space).

`coarsen`

coarsen(g::CartesianGrid) -> CartesianGrid

The next-coarser grid in a 2:1 multigrid hierarchy: half the cells per dimension over the same extent, so the spacing doubles. Boundary conditions, halo width, and device carry over verbatim. Requires an even cell count in every dimension.

Examples

g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
gc = coarsen(g)          # 32×32, spacing doubled

`coarsen!` {#coarsen!}

coarsen!(forest::Forest, should_coarsen) -> forest

Replace each complete family of 2ᴺ sibling leaves — all present and all satisfying should_coarsen(key) — with their parent, then re-balance. Families that are incomplete or only partially flagged are left untouched.

coarsen!(bf::BlockForest, predicate) -> bf

Coarsen each complete family of 2ᴺ sibling leaves whose centers all satisfy predicate(center::SVector), then re-establish 2:1 balance.

`component`

component(f::Field, d::Integer) -> Field

Extract component d of a vector field as a new (allocated) scalar field.

`derivative`

derivative(g::AbstractGrid, dim::Integer; order=1) -> Derivative

Build a matrix-free partial derivative ∂/∂x_dim (or ∂²/∂x_dim² with order=2) bound to g, discretized with second-order central differences. Acts componentwise on any element type.

Examples

g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
u = set!(scalar_field(g), x -> sin(x[1]))
∂u = derivative(g, 1) * u        # ≈ cos

`derivative_stencil`

derivative_stencil(u, I::CartesianIndex, dim, order, inv_h) -> (uc, du)

Per-cell second-order central derivative stencil along dimension dim: returns the center value uc and, for order == 1, du = (u[I+δ] - u[I-δ]) / 2Δ or, for order == 2, du = (u[I-δ] - 2uc + u[I+δ]) / Δ², where inv_h = 1/Δ. The single stencil body shared by the built-in Derivative leaf and custom fused operators. Ghost layers of u must be filled before calling.

Examples

uc, du = derivative_stencil(u.data, CartesianIndex(3, 3), 1, 1, inv(spacing(g)[1]))

`dimension`

dimension(g::AbstractGrid) -> Int

Spatial dimension of the grid.

`divergence`

divergence(g::AbstractGrid) -> Divergence

Build a matrix-free divergence operator (∇⋅) bound to g, discretized with second-order central differences. Maps an SVector{N}-valued vector field to a scalar field — one of the two rank-changing leaves (with gradient). Variable-coefficient diffusion composes as divergence(g) * scaling(κ) * gradient(g).

Examples

g = CartesianGrid(((0.0, 2π), (0.0, 2π)), (32, 32);
                  bc=((Periodic(), Periodic()), (Periodic(), Periodic())))
v = set!(vector_field(g), x -> SVector(sin(x[1]), cos(x[2])))
divv = divergence(g) * v

`flat_to_interior!` {#flat_to_interior!}

flat_to_interior!(f::Field, v::AbstractVector, α=true, β=false) -> f

Fused axpby copy-in of the flat interior vector v (as produced by flatten) into the interior of f: interior(f) = α * v + β * interior(f) in one broadcast. Ghost cells are untouched — the mirror of interior_to_flat!, and the reason the distributed path can stage exchanged ghosts before a local apply.

`flatten`

flatten(f::Field) -> AbstractVector

Copy the interior of f into a flat vector of scalars — the Krylov-facing representation. Spans interior DOFs only (ghost cells are never solver unknowns); SVector elements are flattened component-fastest.

`fold_bc!` {#fold_bc!}

fold_bc!(data, g::AbstractGrid) -> data

Exact discrete adjoint of apply_bc!: fold each ghost value back into its mirror/wrap source cell with the same sign, then zero all ghost layers. Dimensions are folded in reverse order N:1, transposing the fill order exactly.

fold_bc!(x̄::AbstractBlockField, g::BlockForest) -> x̄

Exact discrete adjoint of the forest-level apply_bc!: fold each physical ghost back into its mirror source with the same sign, then zero it — dimensions in reverse order N:1, transposing the fill. Runs after the per-leaf adjoint gathers (which leave ghost cotangents in place) and before halo_update_adjoint!.

`gradient`

gradient(g::AbstractGrid) -> Gradient

Build a matrix-free gradient operator (∇) bound to g, discretized with second-order central differences. Maps a scalar field to an SVector{N}-valued vector field — one of the two rank-changing leaves (with divergence).

!!! note gradient is also exported by Enzyme; qualify as MatrixFreeOperators.gradient when both are loaded.

Examples

g = CartesianGrid(((0.0, 2π), (0.0, 2π)), (32, 32);
                  bc=((Periodic(), Periodic()), (Periodic(), Periodic())))
u = set!(scalar_field(g), x -> sin(x[1]) * sin(x[2]))
∇u = gradient(g) * u             # SVector{2}-valued Field

`halo_update!` {#halo_update!}

halo_update!(x, g::AbstractGrid) -> x

Fill ghost layers with neighbor data. The single distributed seam: a no-op on single-device grids; distributed grids overload it to exchange halos. Operators call this before any stencil that reads neighbor cells.

halo_update!(x::AbstractBlockField, g::BlockForest) -> x

Fill each leaf block’s interface ghosts from its neighbors, in three phases over the precomputed per-generation ExchangeSchedule: same-level slab copies, then coarse→fine quadratic interpolation, then fine→coarse flux-matching restriction (which reads the interpolation-filled fine ghosts). Domain-boundary faces are left to the forest-level apply_bc! face pass. Must run once over the whole forest before any stencil sweep.

On a PackedBlockField with a GPU backend the phases run as batched descriptor kernels over the flattened device schedule — bit-identical to the reference loops except corner ghost cells of the copy phase (the batched per-dim launches let the last dim win where the host order interleaves), which no axis-aligned stencil reads and adjoints require to be zero.

`halo_width`

halo_width(g::AbstractGrid) -> NTuple{N,Int}

Ghost-layer width per dimension.

`identity_op`

identity_op() -> IdentityOp

Build the identity operator, e.g. for shifted systems A - λ*identity_op().

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))
H = laplacian(g) - 4 * identity_op()             # Helmholtz-style shift

`interior`

interior(g::AbstractGrid) -> CartesianIndices

Indices of the interior (owned, non-halo) cells in halo-padded index space.

interior(f::Field) -> SubArray

View of the interior (owned, non-halo) cells of the field.

`interior_to_flat!` {#interior_to_flat!}

interior_to_flat!(v::AbstractVector, f::Field, α=true, β=false) -> v

Fused axpby copy-out of the field interior into the flat vector: v = α * interior(f) + β * v in one broadcast.

`isconstant`

isconstant(L::AbstractOperator) -> Bool

Whether the operator’s parameters and coefficients are time-invariant. Defaults to false; leaves opt in.

`isdiagonal`

isdiagonal(L::AbstractOperator) -> Bool

Whether L acts pointwise (a diagonal operator). Enables cheap Jacobi-type smoothers. Defaults to false.

`islinear`

islinear(L::AbstractOperator) -> Bool

Whether L is a linear map of its input field. Defaults to false; linear leaves opt in. Gates adjoint and use as a Krylov linear map — see linearize for the nonlinear path.

`isselfadjoint`

isselfadjoint(L::AbstractOperator) -> Bool

Whether ⟨L*x, y⟩ = ⟨x, L*y⟩ holds, including boundary contributions. Defaults to false — boundary conditions break self-adjointness easily, so leaves opt in only when their BC handling makes it exact.

`laplacian`

laplacian(g::AbstractGrid) -> Laplacian

Build a matrix-free Laplacian operator (∇²) bound to g, discretized with second-order central differences. Acts componentwise on any element type, so the same operator serves scalar and SVector-valued fields. Exactly self-adjoint under the homogeneous ghost fills of all built-in boundary conditions.

Examples

g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
u = set!(scalar_field(g), x -> sin(x[1]))
Δu = laplacian(g) * u            # ≈ -u

`laplacian_7pt_noflux`

laplacian_7pt_noflux(u, c, i, j, k, nx, ny, nz, sy, sz, ihx, ihy, ihz) -> (uc, lap)

Flat, ghost-free 3D 7-point no-flux (homogeneous Neumann) Laplacian on a column-major vector u at linear index c with subscripts (i, j, k); boundary cells drop the missing face flux. sy/sz are the y/z strides (nx, nx*ny) and ihx/ihy/ihz the inverse squared spacings. Returns the center value uc alongside lap, since fused callers need both.

The sanctioned escape-hatch stencil for fused custom operators (§4a): unlike laplacian_stencil it needs no halo, so it composes a differential stencil with pointwise terms in a single @kernel over the flat solver state. Numerically identical to the ghost-based Laplacian leaf under homogeneous Neumann — (u[c+1] - uc)·ih equals the leaf’s (u[I+δ] - 2uc + ghost)·ih with the mirror ghost = uc — a parity locked by tests.

Examples

uc, lap = laplacian_7pt_noflux(u, c, i, j, k, nx, ny, nz, nx, nx * ny, ihx, ihy, ihz)

`laplacian_stencil`

laplacian_stencil(u, I::CartesianIndex, inv_h2::NTuple) -> (uc, lap)

Per-cell second-order central Laplacian stencil: returns the center value uc and the Laplacian lap = Σ_d (u[I-δd] - 2uc + u[I+δd]) / Δd². This is the single stencil body — the built-in Laplacian leaf calls it, and custom fused operators (§4a of the design) must reuse it so the numerical definition never forks. Ghost layers of u must be filled before calling.

Examples

uc, lap = laplacian_stencil(u.data, CartesianIndex(2, 2), inv.(spacing(g) .^ 2))

`leaves`

leaves(bf::BlockForest)

Iterator over (key, leaf_grid) pairs for every leaf, in Morton (storage) order. Block storage in a BlockField over bf is indexed in the same order. Leaf grids are all-Interface; physical BCs live on bf.bc.

`linearize`

linearize(F::AbstractOperator, u0::Field) -> LinearizedOp

Linearize the (possibly nonlinear) operator F at the state u0, returning a linear matrix-free Jacobian operator whose apply!/mul! is the JVP ∂/∂ε F(u0 + εv)|₀. Owns a copy of u0; see linearize! for in-place refresh inside Newton–Krylov loops.

Examples

g = CartesianGrid(((0.0, 2π),), (64,); bc=((Periodic(), Periodic()),))
F = advection(g, SelfAdvection())            # nonlinear u·∇u
u0 = set!(vector_field(g), x -> SVector(sin(x[1])))
J = linearize(F, u0)                         # linear: v ↦ (∂F/∂u)|_{u0} · v
P = prepare(J, u0)                           # Krylov-ready JFNK Jacobian

`linearize!` {#linearize!}

linearize!(J::LinearizedOp, u::Field) -> J

Refresh the frozen linearization state of J to u in place, reusing the operator and its scratch fields. Safe inside Newton–Krylov loops because Krylov solves are never differentiated through — sensitivities come from implicit-function-theorem adjoints on the solution, not the iteration.

`local_size`

local_size(g::AbstractGrid) -> NTuple{N,Int}

Interior (owned, non-halo) cell counts per dimension.

`ncomponents`

ncomponents(f::Field) -> Int

Number of components of the field’s element type: 1 for scalar fields, N for SVector{N}-valued fields.

`operator_diagonal`

operator_diagonal(L::AbstractOperator) -> Number | Field

Exact diagonal of L as a linear map over interior DOFs, including boundary contributions — a Number when the diagonal is uniform, a scalar Field (ghost entries zero) otherwise. Powers the multigrid smoothers. There is no generic fallback: leaves declare their diagonal explicitly, so a missing declaration errors instead of degrading to a wrong diagonal.

`pack`

pack(f::BlockField) -> PackedBlockField

Copy f into packed contiguous storage: one (blocksize .+ 2halo ..., nleaves) array on the same device, leaves in the same Morton order. The packed field is what the forest-native kernel sweeps consume; prepare on a packed prototype yields packed scratch, so the prepared mul! runs the single-launch path. Coefficient fields (scaling, advection) follow the same layout rule: their kernels engage when the coefficient is packed too — prepare packs BlockField coefficients under a packed prototype automatically, and any remaining layout mismatch degrades to the per-leaf reference sweep.

Examples

u  = set!(scalar_field(bf), x -> sin(π * x[1]))
P  = prepare(laplacian(bf), pack(u))

`padded_size`

padded_size(g::AbstractGrid) -> NTuple{N,Int}

Array size per dimension including ghost layers on both faces.

`partition_grid`

partition_grid(g::CartesianGrid{N}, nparts::Integer) -> Vector{CartesianGrid{N}}

Split g into nparts slab partitions along dimension N — the memory-contiguous dimension, so each partition owns a contiguous range of global flat interior DOFs. Cell counts split evenly with the remainder going to the first partitions.

Each local grid keeps the global spacing verbatim, gets Interface faces on partition cuts (both cut-dimension faces on every partition when the global cut-dimension BC is Periodic), and records its global plane range in local_range. Interface ghost slabs are filled by a distributed exchange (e.g. the MDLA extension), never by apply_bc!.

Every slab must have at least halo planes along the cut dimension (twice that for a periodic cut into exactly two partitions, where both of a partition’s ghost stacks come from the same neighbor), so each ghost slab has a single owner and ghost requests are duplicate-free. partition_grid(g, 1) returns [g] unchanged.

Examples

g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (32, 32))
parts = partition_grid(g, 2)    # two 32×16 slabs cut along dimension 2

`prepare`

prepare(L::AbstractOperator, x::Field) -> PreparedOperator
prepare(L::AbstractOperator) -> PreparedOperator

Walk the operator tree once, allocating the scratch buffers every node needs, and return a PreparedOperator — or a PreparedForest for a BlockField prototype. mul! is then allocation-free in steady state on a single grid; on a block forest the residual cost is the per-leaf stencil apply. x is a prototype of the input field (contents are ignored); the one-argument form assumes a scalar field on the operator’s grid.

The prepared operator is stateful and single-threaded — prepare once per concurrent solve. Requires islinear(L); linearize nonlinear operators first with linearize.

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))
A = prepare(laplacian(g))
b = flatten(set!(scalar_field(g), x -> sin(π * x[1])))
u, stats = Krylov.minres(A, b)

# explicit time stepping (OrdinaryDiffEq-style RHS closure):
f!(du, u, p, t) = mul!(du, A, u)

`prepare_distributed`

prepare_distributed(L::AbstractOperator, nparts::Integer; devices=nothing)

Prepare L for a distributed multi-device solve: partition its grid into nparts slabs (partition_grid), build one prepared operator per partition on its own device, and return a distributed prepared operator exposing mul!/size/eltype over device-partitioned vectors, with ghost slabs exchanged between partitions around each local apply.

Implemented by package extensions; the function has no methods until one is loaded. The MDLA extension (load MultiDeviceLinearAlgebra.jl, CUDA.jl, and Krylov.jl) maps slabs onto one CUDA device each — devices optionally picks which (0-indexed, unique) — and returns an operator over MDLA MultiDeviceVectors.

Examples

using MultiDeviceLinearAlgebra, CUDA, Krylov
g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (256, 256))
P = prepare_distributed(laplacian(g), 2)
b = MultiDeviceVector(flatten(f), P.spec)
u, stats = Krylov.cg(P, b)

`prolongation`

prolongation(coarse::CartesianGrid, fine::CartesianGrid) -> Prolongation

Build the coarse-to-fine linear interpolation operator for a 2:1 grid pair. The grids must share extent, boundary-condition kinds, and device.

Examples

gf = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
gc = coarsen(gf)
P = prolongation(gc, gf)      # coarse -> fine
R = restriction(gf, gc)       # fine -> coarse, R = 2⁻ᴺ·Pᵀ

`refine!` {#refine!}

refine!(forest::Forest, should_refine) -> forest

Replace every leaf key (below maxlevel) for which should_refine(key) is true with its 2ᴺ children, then re-establish 2:1 balance via balance!.

refine!(bf::BlockForest, predicate) -> bf

Refine every leaf whose center satisfies predicate(center::SVector) (and is below maxlevel), then re-establish 2:1 balance. Fields must be (re)allocated after a regrid — block storage is tied to the leaf set at allocation time.

`regrid!` {#regrid!}

regrid!(u::BlockField, more::BlockField...; refine, coarsen=Returns(false))
    -> u′ | (u′, more′...)

Adapt the forest to the current solution and carry the field(s) across the regrid. Evaluates the criteria on every leaf of u — each receives the leaf block as an ordinary Field (data plus leaf grid), so a criterion is typically a reduction like b -> maximum(abs, interior(b)) > τ — then edits the topology in a single pass (refine-marked leaves split, complete fully-marked sibling families coarsen, refine wins conflicts), re-establishes 2:1 balance, and returns freshly allocated field(s) on the new leaf set. Refine marks at maxlevel and coarsen marks on level-0 or incomplete families are silently ignored, matching refine!/coarsen!.

Solution transfer is interior-only and keyed by leaf identity: an unchanged leaf is copied, a refined leaf is filled by per-dimension linear interpolation from its old parent (one-sided at the parent-block edges), and a coarsened leaf by the conservative 2⁻ᴺ mean of its old children. Ghost cells of the returned fields are zero — they are scratch, filled by halo_update!/apply_bc! on the next operator application — so boundary data re-enters through the next solve.

Additional fields passed as more ride the same transfer (results are returned in input order). Because a regrid invalidates every field allocated on the old leaf set, any field needed afterwards — coefficients, or a precomputed indicator — must be passed here; compute indicator fields before calling.

When the marks change nothing, the forest generation is untouched and the input field(s) are returned unchanged — existing prepared operators remain valid. After a real regrid, stale fields and PreparedForests throw on use; re-run prepare on the returned field.

Examples

base = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (32, 32);
                     bc=((Dirichlet(), Dirichlet()), (Dirichlet(), Dirichlet())))
bf   = BlockForest(base; blocksize=(8, 8), maxlevel=3)
u    = scalar_field(bf)

# canonical adaptive solve loop
for cycle in 1:ncycles
    L = laplacian(bf)
    P = prepare(L, u)
    b = flatten(f) .- flatten(boundary_rhs(L, u))
    sol, stats = Krylov.cg(P, b)
    flat_to_interior!(u, sol)

    η = magnitude_of(gradient(bf) * u)          # indicator, computed pre-regrid
    η, u = regrid!(η, u; refine  = b -> maximum(interior(b)) > τ,
                         coarsen = b -> maximum(interior(b)) < τ / 10)
end

See also: refine!, coarsen!, prepare.

`restriction`

restriction(fine::CartesianGrid, coarse::CartesianGrid=coarsen(fine)) -> Restriction

Build the fine-to-coarse full-weighting operator for a 2:1 grid pair. The grids must share extent, boundary-condition kinds, and device.

Examples

gf = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
R = restriction(gf)                  # to coarsen(gf)
P = prolongation(coarsen(gf), gf)    # its adjoint partner: Rᵀ = 2⁻ᴺ·P

`scalar_field`

scalar_field(g::AbstractGrid, T=eltype(spacing(g))) -> Field

Allocate a zeroed cell-centered scalar field on the grid’s device.

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))
u = scalar_field(g)
v = scalar_field(g, Float32)

`scaling`

scaling(κ) -> ScalingOp

Build a pointwise scaling operator x ↦ κ ⊙ x from a Number or a scalar-eltype coefficient field (Field on a single grid, BlockField/ PackedBlockField on a forest — the coefficient’s layout should match the field the operator is applied to for the forest-native kernel to engage). Diagonal and (for real coefficients) self-adjoint — the natural Jacobi-smoother target. Variable-coefficient diffusion composes as divergence(g) * scaling(κ) * gradient(g).

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))
κ = set!(scalar_field(g), x -> 1 + x[1]^2)
K = divergence(g) * scaling(κ) * gradient(g)     # ∇·(κ∇u)

`set!` {#set!}

set!(f::Field, fun) -> f

Set the interior of f to fun(x) evaluated at cell centers, where x is the SVector of physical coordinates.

Examples

g = CartesianGrid(((0.0, 2π),), (64,))
u = set!(scalar_field(g), x -> sin(x[1]))

`solve`

solve(s::MultigridSolver, b::AbstractVector; rtol=1e-8, atol=0) -> u

Solve A·u = b by multigrid V-cycle iteration on flat interior-DOF vectors (the flatten layout), from a zero initial guess, until ‖b - A·u‖ ≤ max(rtol·‖b‖, atol) or maxiter cycles. Warns and returns the current iterate if the tolerance is not reached.

Note: other packages (CommonSolve/SciML) also export a solve; qualify as MatrixFreeOperators.solve when both are loaded.

Examples

u = solve(MultigridSolver(-laplacian(g)), flatten(f); rtol=1e-10)

`spacing`

spacing(g::AbstractGrid) -> NTuple{N}

Cell spacing per dimension.

`unpack`

unpack(f::PackedBlockField) -> BlockField

Copy packed storage back into the vector-of-blocks reference layout — the BlockField that AD, Reactant, and regrid! consume.

`vector_field`

vector_field(g::AbstractGrid{N}, T=eltype(spacing(g))) -> Field

Allocate a zeroed cell-centered vector field with element type SVector{N,T}. Operators written generically over the element type act componentwise on it.

Examples

g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (32, 32))
v = vector_field(g)        # eltype SVector{2,Float64}

Types

AbstractBC

AbstractField

AbstractGrid

AbstractOperator

Added

AdjointOp

Advection

BlockField

BlockForest

Keyword Arguments

Examples

CartesianGrid

Arguments

Keyword Arguments

Examples

Center

Chebyshev

Composed

Derivative

Dirichlet

Divergence

Field

Examples

Gradient

IdentityOp

Jacobi

Laplacian

LinearizedOp

MultigridPreconditioner

Examples

MultigridSolver

Neumann

PackedBlockField

Periodic

PreparedForest

PreparedOperator

Prolongation

Restriction

Scaled

ScalingOp

SelfAdvection

Functions

advection

Examples

apply

Examples

apply! {#apply!}

apply_adjoint! {#apply_adjoint!}

apply_bc! {#apply_bc!}

balance! {#balance!}

boundary_conditions

boundary_rhs

Examples

cell_center

coarsen

Examples

coarsen! {#coarsen!}

component

derivative

Examples

derivative_stencil

Examples

dimension

divergence

Examples

flat_to_interior! {#flat_to_interior!}

flatten

fold_bc! {#fold_bc!}

gradient

Examples

halo_update! {#halo_update!}

halo_width

identity_op

Examples

interior

interior_to_flat! {#interior_to_flat!}

isconstant

isdiagonal

islinear

`AbstractBC`

`AbstractField`

`AbstractGrid`

`AbstractOperator`

`Added`

`AdjointOp`

`Advection`

`BlockField`

`BlockForest`

`CartesianGrid`

`Center`

`Chebyshev`

`Composed`

`Derivative`

`Dirichlet`

`Divergence`

`Field`

`Gradient`

`IdentityOp`

`Jacobi`

`Laplacian`

`LinearizedOp`

`MultigridPreconditioner`

`MultigridSolver`

`Neumann`

`PackedBlockField`

`Periodic`

`PreparedForest`

`PreparedOperator`

`Prolongation`

`Restriction`

`Scaled`

`ScalingOp`

`SelfAdvection`

`advection`

`apply`

`apply!` {#apply!}

`apply_adjoint!` {#apply_adjoint!}

`apply_bc!` {#apply_bc!}

`balance!` {#balance!}

`boundary_conditions`

`boundary_rhs`

`cell_center`

`coarsen`

`coarsen!` {#coarsen!}

`component`

`derivative`

`derivative_stencil`

`dimension`

`divergence`

`flat_to_interior!` {#flat_to_interior!}

`flatten`

`fold_bc!` {#fold_bc!}

`gradient`

`halo_update!` {#halo_update!}

`halo_width`

`identity_op`

`interior`

`interior_to_flat!` {#interior_to_flat!}

`isconstant`

`isdiagonal`

`islinear`

`isselfadjoint`

`laplacian`

`laplacian_7pt_noflux`

`laplacian_stencil`

`leaves`

`linearize`

`linearize!` {#linearize!}

`local_size`

`ncomponents`

`operator_diagonal`

`pack`

`padded_size`

`partition_grid`

`prepare`

`prepare_distributed`

`prolongation`

`refine!` {#refine!}

`regrid!` {#regrid!}