feat(fv): FluxDivergence primitive — geometry accessors, adjoint inner product, conservative advection #54

Status (2026-09-30): the diffusion execution stack (#56) has landed; what remains is geometry accessors, the adaptive-cell adjoint inner product, the FluxDivergence decision, and conservative advection — see §5.

Motivation

MatrixFreeOperators.jl currently provides second-order cell-centered finite-difference leaves. For constant-coefficient diffusion on a uniform Cartesian grid, the compact Laplacian is already algebraically identical to a cell-centered finite-volume balance of centered face fluxes. The distinction becomes consequential for variable coefficients, discontinuities, conservative advection, and AMR: the package needs face-consistent numerical fluxes and conservative state transfer, not merely centered derivatives evaluated at cell centers.

Finite volume does not conflict with the package's matrix-free architecture. A finite-volume operator can apply

(L u)ᵢ = -(1 / Vᵢ) Σ_f A_f F̂_f(u⁻, u⁺)

directly into a cell residual without assembling a matrix or necessarily materializing a global face-flux field. The existing operator algebra, prepare/mul!, device abstraction, BlockForest topology, 2:1 balance, packed storage, and Krylov integration should remain intact.

This issue establishes the conservative finite-volume foundation and carries one scalar diffusion operator through the full execution stack. #48 is the first discretization slice: compact variable-coefficient diffusion that removes the current collocated checkerboard decoupling.

Proposal

1. Define finite-volume geometry and state semantics

• Add cell-volume and face-area accessors for CartesianGrid and BlockForest leaves. Keep them compile-time/simple scalar geometry on the current orthogonal Cartesian grids.
• Decide whether cell averages require a distinct CellAverage location trait or whether Center explicitly permits both point samples and cell averages. Do not introduce public face-field storage until a concrete operator requires it.
• Represent the semi-discrete operator as a flux divergence composed with an explicit diagonal inverse-volume/mass operator where that separation improves adjoint and solver semantics.
• Document whether adjoint remains the Euclidean transpose of the flat-vector map, how the volume-weighted physical inner product is expressed, and which form is appropriate for Krylov methods on non-uniform forests.

2. Implement a fused conservative diffusion leaf

• Implement the compact face-flux form from Compact flux-form diffusion: fix the collocated checkerboard decoupling #48, forming face coefficients and fluxes on demand in one cell sweep.
• Support an explicit coefficient-averaging policy. Arithmetic averaging is appropriate for smooth coefficients; harmonic averaging is needed for flux continuity across jumps.
• Reuse the existing compact Laplacian stencil for constant coefficients so the finite-difference and finite-volume forms cannot drift numerically.
• Preserve the linear/affine boundary split. Neumann data enters as a prescribed boundary flux; Dirichlet data defines the exterior face state or equivalent ghost reconstruction.
• Supply the declared adjoint, operator_diagonal, parameter-gradient path, prepare support, and packed GPU kernel required of a first-class leaf.

3. Make AMR conservation explicit

• Compute one authoritative physical flux at every coarse-fine interface. The coarse residual must receive the area-weighted sum/mean of the abutting fine-face fluxes rather than an independently evaluated coarse flux.
• Reuse the existing exchange schedule and its flux-matching machinery where it remains valid for the compact diffusion flux; factor a face-flux descriptor seam where general numerical fluxes need different data.
• Keep one forest halo exchange per operator application. A finite-volume leaf must not add a second exchange merely because its derivation uses faces.
• For the current single-rate whole-forest application, inject the fine flux directly into the adjacent coarse residual. Defer persistent flux registers/refluxing until level subcycling is implemented.
• Preserve exact transpose behavior through coarse-fine coupling, including the packed-field path.

4. Make regridding conservative for conserved fields

• Retain the existing conservative child mean when coarsening.
• Replace or supplement refinement interpolation with mean-preserving prolongation: the volume-weighted mean of new children must equal the old parent value exactly, including at block and physical boundaries.
• Add a slope-limited prolongation policy before using this path for discontinuous solutions; the current linear-exact interpolation may remain available for non-conserved smooth fields.
• Make conservation semantics explicit per field so geometric indicators, coefficients, and conserved state do not all inherit the same transfer policy accidentally.

5. Stage the rollout

The diffusion slice was tracked in #56 (closed) and has landed:

• Land Compact flux-form diffusion: fix the collocated checkerboard decoupling #48 as the uniform-grid conservative diffusion leaf — feat(operators): compact flux-form diffusion leaf #60, including multigrid rediscretization (src/multigrid.jl) and operator_diagonal (src/operators/diagonal.jl).
• Distributed parity on partition slabs — feat(distributed): compact diffusion leaf on partition slabs #66.
• Extend the leaf through BlockForest and conservative coarse-fine coupling — feat(amr): build the compact diffusion leaf on a BlockForest #74, feat(amr): conservative coarse–fine flux for the forest diffusion leaf #75.
• Add PackedBlockField GPU parity, the forest-native kernel, and conservative regrid! (Conservative/SlopeLimited transfer policies set per field via with_transfer) — perf(amr): packed diffusion kernels and conservative regrid transfer #76.

What remains is the part that generalizes beyond diffusion:

• Geometry/mass semantics: cell-volume and face-area accessors on CartesianGrid and BlockForest leaves, the explicit inverse-volume/mass operator, and the CellAverage-vs-Center decision (§1).
• Decide which inner product adjoint documents on adaptive cells (§1).
• Decide whether a reusable FluxDivergence parameterized by a numerical flux becomes the public primitive, with diffusion lowering to it. Compact flux-form diffusion: stage the Diffusion leaf from CartesianGrid through forest and GPU #56 shipped the concrete leaf first because a user-supplied flux callable cannot satisfy the declared-adjoint invariant until there is a second flux to generalize from; the anisotropic tensor leaf in feat(operators): anisotropic diffusion ∇·(D∇u) with a full tensor D #70 is the natural second flux.
• Use the resulting face-flux contract to add conservative advection (Advection is still the non-conservative central v·∇ in src/operators/advection.jl) and scope higher-resolution reconstructions.

Success criteria

• Constant-coefficient diffusion is bitwise-identical, or differs only by documented floating-point operation order, from the existing compact laplacian on uniform grids.
• Variable-coefficient diffusion uses a compact nearest-neighbor stencil and does not exhibit the even/odd checkerboard decoupling described in Compact flux-form diffusion: fix the collocated checkerboard decoupling #48.
• For periodic or zero-flux boundaries, Σᵢ Vᵢ(Lu)ᵢ is zero to roundoff on uniform and non-uniform forests. With prescribed boundary fluxes, it equals the net boundary flux to roundoff.
• Refinement and coarsening preserve Σᵢ Vᵢuᵢ to roundoff for fields using conservative transfer, including balance-induced refinements.
• The discrete adjoint satisfies the package's declared inner-product contract on Field, BlockField, and PackedBlockField; materialized small problems agree with the corresponding transpose.
• CPU/GPU, packed/reference, and single-/multi-device paths agree within the package's existing numerical tolerances.
• Prepared applications remain allocation-free on the established CPU hot paths and retain a leaf-count-independent launch structure on packed GPUs.
• The constant-coefficient finite-volume path is within 5% of the existing compact Laplacian in representative uniform and adaptive benchmarks. Any larger regression is explained by measured additional work rather than accepted as an inherent "finite-volume tax."
• The adapted Niederer-style benchmark performs no more halo exchanges per diffusion application than it does before this work; results include a phase breakdown comparable to meta(amr): forest performance — the coarse–fine exchange, not regridding, is the AMR tax #49.

Out of scope

• Replacing or removing the existing finite-difference leaves.
• Building an unstructured, curvilinear, cut-cell, or embedded-boundary mesh framework.
• Level-subcycled AMR time integration, persistent flux registers, or refluxing across timesteps.
• A complete hyperbolic solver suite, Riemann-solver catalog, WENO/TVD reconstruction framework, or positivity-preserving limiter stack. Those should follow after the scalar flux contract is proven.
• Materializing public face fields solely to make the implementation look conventionally finite-volume; fused ephemeral face fluxes remain the default unless profiling justifies stored fluxes.
• Forest geometric multigrid, which is not current package functionality.

Open questions

• Should cell averages use a distinct CellAverage trait, or is that a semantic mode of the existing Center storage?
• Should the public primitive be a problem-specific conservative_diffusion leaf, a reusable FluxDivergence parameterized by a numerical-flux callable, or both with the former lowering to the latter?
• Is the existing coarse-fine ghost construction the right execution mechanism for general fluxes, or should the exchange schedule gain explicit face-flux descriptors while retaining the same topology cache?
• Face-trace (DG) vs corner (feat(operators): anisotropic diffusion ∇·(D∇u) with a full tensor D #70) vs face-flux (this issue): one exchange mechanism or three? Carried over from the DG investigation in Investigate whether the matrix-free operator algebra extends to discontinuous Galerkin #71, which found that DG face terms need a face-trace exchange, a cross-derivative tensor leaf needs a corner-aware exchange, and a general FluxDivergence may need face-flux descriptors.
• Which inner product should user-facing adjoint document on adaptive cells: the existing Euclidean flat-vector transpose, a volume-weighted physical adjoint, or the Euclidean transpose with an explicit mass-operator conversion?
• How should conservative transfer policy attach to a field without adding runtime branching to operator hot paths?

Related: #48 (compact flux-form diffusion), #56 (the diffusion slice, closed), #70 (anisotropic tensor diffusion, the second flux), #71 (DG investigation, closed), #49 (measured AMR exchange cost), #52 (traffic-reduction and exchange-overlap work).

—
🤖 Claude trimmed the checklist down to what's actually left and adopted #71's best question.
_{claude-opus-5-5[1m] (high) · claude-code 2.1.285 · unreviewed}