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.

`AbstractJVPBackend`

AbstractJVPBackend

How a linearized operator evaluates its Jacobian–vector product. Selecting the backend is deliberately explicit rather than ambient: the backends differ in accuracy and in which operations they support, so which one you get must not depend on whether some package happens to be loaded.

FiniteDifferenceJVP — the dependency-free default.
EnzymeJVP — exact, and the preferred choice; needs using Enzyme.

`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.

`ArithmeticMean`

ArithmeticMean()

Arithmetic face averaging κ_f = (κ_a + κ_b)/2. The default for diffusion: second-order for smooth coefficients, and imposes no sign restriction on κ.

`BlockField`

BlockField{L,P}(blocks, grid[, generation])
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) and regrid-transfer policy P (default Interpolated; see Conservative/SlopeLimited and with_transfer). 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.

`Conservative`

Conservative()

Mean-preserving regrid-transfer policy: refined leaves are filled by the cell-conservative linear reconstruction u_child = u_parent + Σ_d ξ_d·σ_d (ξ_d = ∓¼), one shared slope per parent cell per dimension — centered in the block interior, one-sided difference at parent-block edges. The volume-weighted mean of the 2ᴺ children equals the parent exactly for any slope, including at block and physical boundaries, so Σ V·u is preserved to roundoff across regrid! (coarsening already is: the 2⁻ᴺ child mean). Exact on linears and second-order like Interpolated. The policy for conserved state on smooth solutions.

See also: SlopeLimited, with_transfer.

`Derivative`

Derivative(grid, dim, order)

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

`Diffusion`

Diffusion(grid, κ, avg)

Compact flux-form variable-coefficient diffusion ∇·(κ∇u). Construct with diffusion, which is what extends κ into its ghost layers; the inner constructor takes κ as given and is the seam for supplying cross-block coefficient ghosts from outside: _slab_op (src/distributed.jl) builds through it with a κ sliced from the global one including its ghosts, and _leaf_op (src/operators/forest.jl) with a block view whose ghosts fill_coefficient_ghosts! filled at construction.

`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.

`EnzymeJVP`

EnzymeJVP()

Evaluate the JVP by forward-mode automatic differentiation, and the transpose by reverse mode. Exact to machine precision, one operator application per product, and — unlike FiniteDifferenceJVP — it supplies a real adjoint, so transpose-needing Krylov methods work on a Jacobian-free Newton–Krylov operator.

Provided by the Enzyme extension: using Enzyme before calling linearize(F, u₀, EnzymeJVP()).

`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)

`FiniteDifferenceJVP`

FiniteDifferenceJVP()

Evaluate the JVP by a central finite difference, (F(u₀+εv) - F(u₀-εv))/2ε. The default, because it needs no dependencies — but it costs two operator applications per product, is accurate only to about √eps, and has no transpose, so adjoint of the resulting Jacobian throws. Prefer EnzymeJVP.

`Gradient`

Gradient(grid)

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

`HarmonicMean`

HarmonicMean()

Harmonic face averaging κ_f = 2κ_aκ_b/(κ_a + κ_b). Conserves flux across a jump in κ, so it is the right choice for piecewise-constant or discontinuous media. Requires κ > 0 — the expression is singular at κ_a = -κ_b and NaN at (0, 0).

Note for inverse problems: the harmonic mean is dominated by the smaller of its two arguments, so ∂κ_f/∂κ_large → 0 as the contrast grows. That weakens identifiability of the high-κ side of an interface. It is a real property of the discretization, not a bug.

`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.

`Interpolated`

Interpolated()

Regrid-transfer policy trait (the default): refined leaves are filled by the linear-exact per-dimension interpolation of the old parent. Second-order and exact on linears, but not mean-preserving — the two children of an interior parent cell use side-biased slopes, leaving the child mean off by ⅛·δ²u per dimension. The right transfer for smooth, non-conserved fields (geometric indicators, coefficients). Attach a policy at construction (scalar_field/vector_field’s transfer keyword) or with with_transfer; regrid! resolves it per field at transfer time, so nothing on an operator hot path ever consults it.

See also: Conservative, SlopeLimited.

`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,P}(data, levels, grid[, generation])

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. L is the location trait and P the regrid-transfer policy (default Interpolated); both survive pack/unpack, and P is consulted only by regrid!, never on an operator path. 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.

Those copies are the price of the Krylov boundary, not of the operator: a prepared operator also takes fields directly through apply!(du, P, u), the explicit time-stepping path.

The traits islinear, isconstant, isselfadjoint and isdiagonal are those of the operator handed to prepare — binding buffers changes nothing about the map. A PreparedOperator is the solver boundary rather than a node of the lazy algebra, so it is not an AbstractOperator: it does not enter +, *, adjoint, or apply; compose first, then prepare.

`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.

`SlopeLimited`

SlopeLimited()

Mean-preserving and bounds-preserving regrid-transfer policy: the same cell-conservative reconstruction as Conservative with the per-dim slope minmod-limited, and dropped to zero at parent-block edges (the transfer is interior-only, so no second slope exists there to limit against). Children never leave the hull of the parent’s neighborhood — no new extrema across a regrid — at the price of first-order transfer at extrema and block edges. Note the edge cost scales with the block: every parent cell on a block face is injected, so a 4×4 block reconstructs only its inner 2×2 and larger blocks shrink that fraction. The policy for conserved state with steep fronts or discontinuities.

See also: Conservative, with_transfer.

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 in-place apply! instead (explicit stepping), or prepare + mul! at the Krylov boundary.

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. The three-argument form takes the grid from x.

This is the explicit time-stepping path; prepare + mul! is the Krylov boundary and stages through flat vectors. *-composed and adjoint trees allocate their intermediates here — prepare those once and call apply!(du, P, u) instead. Like every apply! this is the homogeneous linear part: with inhomogeneous BCs the RHS is L(u) + b, b = boundary_rhs(L, g).

Examples

L = laplacian(g)
du = similar(u)
apply!(du, L, u)                        # du = Δu
interior(u) .+= dt .* interior(du)      # one forward-Euler step

See also: apply, apply_adjoint!, prepare.

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

Apply L over a block-structured forest: exchange inter-block halos once, then run the stencil sweep over every block.

apply!(y::AbstractField, P::PreparedOperator, x::AbstractField, α=true, β=false) -> y
apply!(y::AbstractField, P::PreparedForest, x::AbstractField, α=true, β=false) -> y

Apply a prepared operator at field level: y = α·P(x) + β·y on the interior of y, reusing the scratch prepare bound (so *-composed and adjoint nodes do not reallocate) and skipping the flat staging copies mul! performs. Ghost handling and the homogeneous-BC caveat are those of apply! on an unprepared operator.

prepare once, then apply!(du, P, u) per explicit stage; keep mul! for Krylov. x and y must be on P.grid with the prototype’s element type, and the forest form additionally rejects a regridded forest.

`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.

`assemble_rhs`

assemble_rhs(P, f) -> distributed vector

Solve-ready right-hand side f - boundary_rhs(P) for a distributed prepared operator P (prepare_distributed), assembled entirely partition-locally — nothing global is materialized on one device.

f is either a function of physical coordinates or one Field per partition on the grids local_grids reports. Implemented by package extensions, like prepare_distributed; there is no single-device method, since a serial right-hand side is the one-liner flatten(f) .- flatten(boundary_rhs(L, g)).

This is the whole right-hand side, not a sibling of boundary_rhs — it is the function that consumes one.

Examples

P = prepare_distributed(laplacian(g), 2)
b = assemble_rhs(P, x -> sin(x[1]) * exp(-x[2]))
u, stats = Krylov.cg(P, b)

See also: local_grids, boundary_rhs.

`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).

Evaluated at the global cell index local_range reports, against the global extent origin — so a partition_grid slab and the grid it was cut from return bitwise identical coordinates for the same cell. Recomputing from a slab-local origin instead would round twice and drift by an ulp, which is enough to make a coordinate-assembled right-hand side depend on the partition count. For every undistributed grid and every BlockForest leaf grid first(local_range[d]) == 1, so this is bit-for-bit the plain formula min + (i - 1/2)Δ.

`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]))

`diffusion`

diffusion(g::AbstractGrid, κ::Field; averaging=ArithmeticMean(), check=true) -> Diffusion

Build a matrix-free variable-coefficient diffusion operator ∇·(κ∇u) on g, discretized in compact flux form: fluxes q_{i+½} = κ_{i+½}(u_{i+1} − u_i)/Δ live on faces, with κ averaged to faces by averaging, and the cell balance is (q_{i+½} − q_{i−½})/Δ.

Prefer this to the algebraic composition divergence(g) * scaling(κ) * gradient(g). That composition chains two centered first differences and carries a wide 2Δ stencil: interior rows sample fluxes at I±e and reach solution values at I±2e, skipping adjacent solution cells. The compact form couples adjacent solution cells and removes this odd–even decoupling in u.

Coefficient identifiability is a separate question. With ArithmeticMean, a checkerboard perturbation δκ[i,j] = ε(-1)^(i+j) cancels at every interior face. On periodic grids with even cell counts in every dimension, or with homogeneous Neumann walls, it leaves the entire operator unchanged for every u. Multiple excitations cannot resolve that ambiguity. Dirichlet wall coefficients use the adjacent cell’s κ and can break it; the successful Dirichlet inversion in examples/inverse_diffusion.jl does not establish unique recovery for other boundary conditions or data.

The compact form is also exactly symmetric for real κ under all built-in boundary conditions (the composed form is not), and it has an operator_diagonal (the composed form cannot), which is what makes it usable with the multigrid smoothers.

With κ ≡ c constant this reduces to c * laplacian(g) up to floating-point operation order.

Boundary treatment is exact under the package’s homogeneous ghost fills, which reflect about the face: a Neumann face carries zero flux, so its face coefficient is irrelevant; a Dirichlet face carries κ_I(0 − u_I)/(Δ/2), which the one-sided face coefficient κ_f = κ_I reproduces exactly. That one-sided value is O(Δ) accurate in κ at a wall while interior faces are O(Δ²), so wall-adjacent κ cells have a different sensitivity structure than interior ones — worth knowing when inverting for κ.

κ must be a scalar-eltype Field matching g. The operator stores its own copy with ghost layers extended by fill_coefficient_ghosts!, so whatever κ holds in its ghosts is ignored — and so mutating κ after construction does not affect the operator. Pass check=false to skip the κ > 0 validation that HarmonicMean requires, which is what an optimization loop wants: the operator is rebuilt every objective evaluation and the check would run inside the differentiated region.

Building the leaf inside that region — which is exactly what an inversion loop does — needs Enzyme.set_runtime_activity(Enzyme.Reverse): the operator stores the Const grid alongside the active coefficient, and static activity analysis cannot clear that store before Julia 1.12. See examples/inverse_diffusion.jl, which sets it once on the backend object.

Examples

g = CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (64, 64))
κ = set!(scalar_field(g), x -> 1 + x[1]^2)
L = diffusion(g, κ)                                  # ∇·(κ∇u), compact 5-point
Lj = diffusion(g, κ; averaging=HarmonicMean())       # flux-conserving across jumps

diffusion(bf::BlockForest, κ::AbstractBlockField; averaging=ArithmeticMean(), check=true) -> Diffusion

Build the same compact flux-form operator over a block-structured forest. Semantics, boundary treatment, and the constant-κ reduction to c * laplacian(bf) are exactly the single-grid method’s; only the coefficient’s ghost fill differs, because a forest’s blocks are separate arrays. κ must be a scalar-eltype block field on the same forest as bf — _leaf_op slices it by this forest’s leaf indices, so a coefficient bound to a different topology would alias the wrong leaves.

The operator stores its own copy with every ghost layer filled by fill_coefficient_ghosts!, so mutating κ afterwards does not affect it, and the exchange is paid once here rather than per application.

On a non-uniform forest the operator is discretely conservative: each application rewrites every coarse-side coarse–fine ghost so the coarse stencil’s κ-weighted face flux equals the area-weighted sum of the abutting fine-face fluxes (_cf_flux_rewrite! — one authoritative flux per interface, issue #58), and Σ V·(Lu) telescopes to the net boundary flux to roundoff exactly as on a single grid. The rewrite divides by the coarse face coefficient, so construction additionally validates that avg(κ₁, κ_ghost) is finite and nonzero on every coarse–fine face — with check=false that validation (like the HarmonicMean one) is skipped, and a sign-changing κ under ArithmeticMean can then produce Inf at apply time. Coarse–fine coupling is not symmetric, so isselfadjoint is false (via _selfadjoint_grid) and the adjoint runs the declared transpose — the per-leaf gather plus the rewrite’s exact transpose. operator_diagonal is unavailable on forest leaves, matching laplacian.

Examples

bf = BlockForest(CartesianGrid(((0.0, 1.0), (0.0, 1.0)), (16, 16)); blocksize=(8, 8), maxlevel=2)
κ = set!(scalar_field(bf), x -> 1 + x[1]^2)
L = diffusion(bf, κ)

`diffusion_stencil`

diffusion_stencil(u, κ, I::CartesianIndex, inv_h2::NTuple, avg) -> (uc, div)

Per-cell compact flux-form stencil: returns the center value uc and

div = Σ_d [ κ_{I+½δ}·(u[I+δ] − uc) − κ_{I−½δ}·(uc − u[I−δ]) ] / Δd²

the finite-volume balance of the face fluxes κ∇u, with κ_{I±½δ} = avg(κ[I], κ[I±δ]). This is the single stencil body — the built-in Diffusion leaf calls it, and custom fused operators (§4a of the design) must reuse it so the numerical definition never forks. Ghost layers of both u and κ must be filled first: u by apply_bc!, κ by fill_coefficient_ghosts!.

Every equation couples κ and u at adjacent cells, which is what the wide divergence ∘ scaling ∘ gradient composition does not do — see diffusion.

Examples

uc, div = diffusion_stencil(u.data, κ.data, CartesianIndex(2, 2),
                            inv.(spacing(g) .^ 2), ArithmeticMean())

`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 can be assembled as divergence(g) * scaling(κ) * gradient(g), but prefer the compact flux-form leaf diffusion — chaining two centered differences doubles the stencil width.

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

See also: diffusion, gradient, scaling.

`fill_coefficient_ghosts!` {#fill_coefficient_ghosts!}

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

Extend a coefficient array into its ghost layers: periodic wrap where the dimension is periodic, even (symmetric) mirror at every physical wall, and Interface faces left untouched because those ghosts hold a neighbour’s coefficient and are an external input.

Deliberately not apply_bc!, which would antisymmetrize the coefficient at a Dirichlet wall — the same hazard _average_to_coarse avoids by not using Restriction. The even mirror is what makes a wall face self-average to the adjacent cell’s value, which is the one-sided face coefficient the finite-volume balance calls for.

fill_coefficient_ghosts!(κ::AbstractBlockField, bf::BlockForest) -> κ

Extend a coefficient field into every block’s ghost layers: same-level face copies, coarse→fine injection, fine→coarse volume averaging, and an even mirror at physical walls. The forest counterpart of the single-grid method above, run once by diffusion at construction so the per-application exchange count is unchanged from the Laplacian baseline. A forest has no global coefficient array to window, so unlike a partition slab (_slab_coeff_field) this is a real exchange — legitimate only because κ is constant through a solve.

Deliberately not halo_update!, whose coarse–fine phases are tuned for a solution: the Martin–Cartwright quadratic interpolant assumes a smoothness a material coefficient need not have, and the flux-matching restriction is a stencil on u, not a volume average. Deliberately not apply_bc! either, which antisymmetrizes at a Dirichlet wall and would negate the coefficient there — the same hazard the single-grid method and _average_to_coarse document.

Written as a plain topology walk rather than through the cached ExchangeSchedule descriptors on purpose. It must stay differentiable with respect to κ, and the GhostFill descriptor sweep runs under the halo Enzyme rules (_exchange_storage!, issue #26); those rules report nothing for every derivative slot, so routing a coefficient through them would silently zero its gradient. Recomputing topology here is free — this runs once per operator, never per application.

`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, backend=FiniteDifferenceJVP())

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.

backend selects how the product is evaluated — see AbstractJVPBackend. The default FiniteDifferenceJVP keeps the core dependency-free; EnzymeJVP is exact and additionally provides the transpose.

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

using Enzyme                                 # exact JVP, and adjoint(J) works
Jad = linearize(F, u0, EnzymeJVP())

`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_grids`

local_grids(P) -> Vector{<:AbstractGrid}

The grid each partition of a distributed prepared operator owns, in partition order.

The escape hatch for source data assemble_rhs cannot build from a coordinate function: allocate a Field on one of these, fill it however you like, and pass the vector of fields to assemble_rhs. Each grid records its span of the global grid in local_range, and cell_center on it agrees bitwise with the uncut grid. Implemented by package extensions.

`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 and the global extent verbatim — its position is carried entirely by local_range, which records the global plane range it owns, so cell_center agrees bitwise with the uncut grid. It gets Interface faces on partition cuts (both cut-dimension faces on every partition when the global cut-dimension BC is Periodic). 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.

mul! is the solver boundary: every call stages the flat vector into a halo-padded field and back out again. Explicit integrators do not need that — step at field level with apply! instead.

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 stepping stays at field level:
uf = set!(scalar_field(g), x -> sin(π * x[1]))
du = similar(uf)
dt = 0.4 * spacing(g)[1]^2              # forward-Euler bound
apply!(du, A, uf)                       # or apply!(du, laplacian(g), uf)
interior(uf) .+= dt .* interior(du)

`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 = assemble_rhs(P, x -> sin(x[1]) * exp(-x[2]))
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 the field’s regrid-transfer policy from its old parent, and a coarsened leaf by the conservative 2⁻ᴺ mean of its old children. The policy is per field, resolved from the field’s type at transfer time: Interpolated (the default) is the linear-exact per-dimension interpolation — second-order, right for indicators and coefficients, not mean-preserving; Conservative and SlopeLimited use the cell-conservative reconstruction, preserving Σ V·u to roundoff across every regrid — balance-induced refinements included. Conservation is deliberately not a default: mark conserved state explicitly (the transfer keyword of scalar_field/vector_field, or with_transfer) so indicators, coefficients, and state do not silently share one policy. 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.

!!! note “For ∇·(κ∇u), reach for diffusion” diffusion(g, κ) is the compact flux form: exactly symmetric, with an operator_diagonal, and ~3× faster. divergence(g) * scaling(κ) * gradient(g) also builds a valid ∇·(κ∇u) and is a good demonstration that the algebra composes, but it chains two centered differences: the resulting stencil reaches u[I±2δ] and never u[I±δ], so the grid decouples into interleaved sublattices and κ at a cell never enters that cell’s own equation — which makes a κ-inversion fit disjoint halves of the data.

Examples

g = CartesianGrid(((0.0, 1.0),), (64,))
κ = set!(scalar_field(g), x -> 1 + x[1]^2)
H = laplacian(g) - scaling(κ)                    # Helmholtz-type: Δu − κ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]))

`shares_exchange`

shares_exchange(L::AbstractOperator) -> Bool

Whether L’s action on a BlockForest is one stencil sweep over an already-exchanged input: it reads x’s interiors and ghosts, writes only y’s interior, and never writes into x or an intermediate needing its own exchange. An Added whose operands all claim this runs a single halo_update!/apply_bc! for the set. Defaults to false; stencil leaves opt in, Added/Scaled propagate, Composed/AdjointOp never claim it. A forgotten declaration costs a redundant exchange, never a wrong result.

`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}

`with_transfer`

with_transfer(f, policy) -> field

The same field — shared storage, same location trait and generation — carrying policy (Interpolated, Conservative, or SlopeLimited) as its regrid-transfer policy. The way to mark an existing field (e.g. conserved state built with the plain constructors) before handing it to regrid!; no data is copied.

Types

AbstractBC

AbstractField

AbstractGrid

AbstractJVPBackend

AbstractOperator

Added

AdjointOp

Advection

ArithmeticMean

BlockField

BlockForest

Keyword Arguments

Examples

CartesianGrid

Arguments

Keyword Arguments

Examples

Center

Chebyshev

Composed

Conservative

Derivative

Diffusion

Dirichlet

Divergence

EnzymeJVP

Field

Examples

FiniteDifferenceJVP

Gradient

HarmonicMean

IdentityOp

Interpolated

Jacobi

Laplacian

LinearizedOp

MultigridPreconditioner

Examples

MultigridSolver

Neumann

PackedBlockField

Periodic

PreparedForest

PreparedOperator

Prolongation

Restriction

Scaled

ScalingOp

SelfAdvection

SlopeLimited

Functions

advection

Examples

apply

Examples

apply! {#apply!}

Examples

apply_adjoint! {#apply_adjoint!}

apply_bc! {#apply_bc!}

assemble_rhs

Examples

balance! {#balance!}

boundary_conditions

boundary_rhs

Examples

cell_center

coarsen

Examples

coarsen! {#coarsen!}

component

derivative

Examples

derivative_stencil

Examples

diffusion

Examples

Examples

diffusion_stencil

Examples

`AbstractBC`

`AbstractField`

`AbstractGrid`

`AbstractJVPBackend`

`AbstractOperator`

`Added`

`AdjointOp`

`Advection`

`ArithmeticMean`

`BlockField`

`BlockForest`

`CartesianGrid`

`Center`

`Chebyshev`

`Composed`

`Conservative`

`Derivative`

`Diffusion`

`Dirichlet`

`Divergence`

`EnzymeJVP`

`Field`

`FiniteDifferenceJVP`

`Gradient`

`HarmonicMean`

`IdentityOp`

`Interpolated`

`Jacobi`

`Laplacian`

`LinearizedOp`

`MultigridPreconditioner`

`MultigridSolver`

`Neumann`

`PackedBlockField`

`Periodic`

`PreparedForest`

`PreparedOperator`

`Prolongation`

`Restriction`

`Scaled`

`ScalingOp`

`SelfAdvection`

`SlopeLimited`

`advection`

`apply`

`apply!` {#apply!}

`apply_adjoint!` {#apply_adjoint!}

`apply_bc!` {#apply_bc!}

`assemble_rhs`

`balance!` {#balance!}

`boundary_conditions`

`boundary_rhs`

`cell_center`

`coarsen`

`coarsen!` {#coarsen!}

`component`

`derivative`

`derivative_stencil`

`diffusion`

`diffusion_stencil`

`dimension`

`divergence`

`fill_coefficient_ghosts!` {#fill_coefficient_ghosts!}

`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`