feat(fv): add conservative finite-volume operators without sacrificing matrix-free AMR

[kylebeggs]

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

• Land Compact flux-form diffusion: fix the collocated checkerboard decoupling #48 as the uniform-grid conservative diffusion leaf.
• Add geometry/mass semantics and conservation identities on CartesianGrid.
• Extend the leaf through BlockForest, coarse-fine coupling, and conservative regrid!.
• Add PackedBlockField GPU parity and the forest-native kernel.
• Add uniform-grid multigrid rediscretization/diagonal support and distributed parity.
• Use the resulting face-flux contract to scope follow-on conservative advection and 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?
• 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), #49 (measured AMR exchange cost), #52 (traffic-reduction and exchange-overlap work).

---

_{🤖 Claude left the flux ledger balanced; Kyle remains innocent of conversing with his own issue tracker.}

[kylebeggs]

Cross-reference: the diffusion slice of this programme now has its own tracker, #56, with the execution stack broken into stages — #48 (uniform CartesianGrid), #57 (distributed slabs), #58 (BlockForest + coarse–fine flux), #59 (packed-GPU kernels + conservative regrid). That covers step 1 of the rollout here plus the execution work steps 3–5 describe, for one scalar operator.

What stays with this issue is the part that generalizes beyond diffusion, and which #56 deliberately does not decide:

• cell-volume and face-area accessors, and the volume-weighted vs Euclidean inner-product question for adjoint on adaptive cells;
• whether cell averages need a distinct CellAverage location trait or are a semantic mode of Center;
• whether the public primitive should become a reusable FluxDivergence parameterized by a numerical flux, 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, on the grounds that a user-supplied flux callable cannot satisfy the declared-adjoint invariant until there is a second flux to generalize from;
• conservative advection and higher-resolution reconstructions.

Two results from #48 that are inputs here. The compact form is exactly symmetric for real κ under every built-in boundary condition, so the "which inner product should adjoint document" question can be answered on a uniform grid without ambiguity — it only becomes interesting once cell volumes differ. And constant-κ parity with the existing compact laplacian holds to ~1 ulp but not bitwise, because the flux form multiplies by κ per face before summing; the success criterion in this issue should read as a tolerance, not an equality.

---

_{🤖 Claude threading the issue graph — Kyle is not holding a séance with his own tracker.}

[kylebeggs]

Tracker housekeeping, 2026-09-09: #56 (the diffusion slice) is closed — all four stages landed via #60, #66, #74/#75 and #76. Its two review follow-ups were re-homed by theme rather than kept here: #77 (slab Diffusion adjoint fast path) sits under the distributed tracker #79, and #78 (batch the coarse–fine flux rewrite) under the AMR-performance tracker #49.

What remains open in this issue is unchanged: geometry accessors, the volume-weighted vs Euclidean inner-product question for adjoint on adaptive cells, whether FluxDivergence becomes the public primitive, and conservative advection. #71 (DG investigation) depends on the FluxDivergence answer and cross-references here.

---

_{🤖 Beep boop — Claude tidied the flux ledger's footnotes; Kyle is not auditing himself.}