Compact flux-form diffusion: fix the collocated checkerboard decoupling

[kylebeggs]

Problem

The advertised variable-coefficient diffusion divergence(g) * scaling(κ) * gradient(g) chains two collocated centered first differences, producing the wide 2h-stencil Laplacian. A centered first difference annihilates the grid-Nyquist mode, so on a collocated grid this composition has the classic checkerboard pathology: the equation at cell I samples the flux κ∇u only at neighbors I±e, never at I, and the misfit Hessian FᵀF of a κ-inversion exactly decouples the even and odd (i+j)-parity sublattices. Nothing in the data ever couples two adjacent κ pixels — only the regularizer does.

Measured while iterating on examples/inverse_diffusion.jl (143fe79). Exact regularized least-squares minimizers (forward map is linear in κ, n=48, 1% noise):

scenario	rel-err	jaggedness (truth = 0.0046)
1 excitation, λ=5e-4	27.7%	0.415
1 excitation, λ=1e-2 (20× prior)	20.1%	0.307
3 excitations, λ=5e-4	3.5%	0.051

A 20× stronger prior barely helps; the pathology is structural. (The example now uses three excitations, which shrinks each sublattice's variance until the interleaved halves agree — a right-sized fix for the example, not a cure.)

Scope: forward applies are unaffected — the wide stencil is consistent and second-order on smooth fields. The decoupling bites when energy is injected at the grid scale and nothing removes it: noisy parameter inversion w.r.t. κ (the flagship AD use case), and null-mode-sensitive solves. Same disease as collocated incompressible flow (cured there by MAC staggering or Rhie–Chow interpolation).

The true fix: a compact flux-form stencil

Store fluxes at faces: q_{i+½} = κ_{i+½}(u_{i+1} − u_i)/h, then (∇·q)i = (q{i+½} − q_{i−½})/h, with κ averaged to faces (arithmetic for smooth κ, harmonic to conserve flux across jumps). Expanded, this is the compact 5-point variable-coefficient Laplacian: every equation couples κ and u at adjacent cells, so the checkerboard mode is simply not in the null space.

Two routes:

Option A — staggered rank-changers (the composable, durable fix)

gradient: Center → face-located flux components; scaling(κ) averages κ to faces; divergence: faces → Center. The algebra D * S(κ) * G survives. This is exactly what the Field location trait L was reserved for (DESIGN.md: L ∈ {Center, XFace, YFace, …}, v1 Center-only; open question §10.2 recommends collocated v1 with staggered purely additive). Real infrastructure: face fields, BCs and ghost conventions per location, declared adjoints for the new leaves, and the AD invariants re-verified. Also what incompressible flow will eventually want.

Option B — fused diffusion(g) leaf on collocated storage (the contained fix)

A single leaf applying the flux form directly from cell-centered κ and u (face κ formed on the fly). Small to write, immediately unblocks accurate κ-inversion. Costs: gives up the advertised composition for exactly this operator; needs its own declared adjoint including boundary contributions and the four standard checks; the κ-gradient path must stay rule-free per the invariant that custom rules never fire on paths carrying operator parameters. Face-averaging choice (arithmetic vs harmonic) should be an explicit argument, since it changes the adjoint.

Recommendation

Option B near-term (it is a leaf, matching the existing leaf/combinator architecture, and inverse problems in κ are the flagship use case), with Option A folded into the §10.2 staggered-grids decision when that gets resolved — B does not preclude A, and the fused leaf remains useful afterwards as the fast path.

---

_{🤖 Filed by Claude while Kyle was busy not talking to himself in his own issue tracker.}

[kylebeggs]

Restructured: this is now stage 1 of #56, which tracks carrying the compact flux-form leaf through the rest of the execution stack — distributed slabs (#57), BlockForest (#58), and packed-GPU kernels plus conservative regrid (#59).

The scope of this issue is unchanged and stays as filed: the leaf on a uniform CartesianGrid. Option A (staggered rank-changers) remains deferred to the DESIGN.md §10.2 decision and is not part of any stage here.

One correction to this issue's own recommendation, worth recording since it shaped the implementation: the plan to give the fused leaf "its own declared adjoint including boundary contributions" understated the result. The compact form is exactly symmetric for real κ — asymmetry measures as a hard zero, not a tolerance — across Dirichlet/Neumann/periodic and mixed combinations, because each shared face contributes the same κ_f to both A[I,J] and A[J,I] while a wall face contributes only to A[I,I]. That needs just two things: the averaging function is symmetric in its arguments, and the periodic neighbour map is an involution. So the adjoint short-circuits to the forward action on an all-physical grid, and the gather is only needed once Interface faces appear in #58.

Also worth flagging for anyone who later attempts Option A: the operator factors as L = −Bᵀ diag(κ_f) B, but at a Dirichlet wall that row of B is √2/Δ, not 2/Δ — the geometric mean of the one-sided gradient weight and the divergence weight. A staggered G/D pair that assumes D = −Gᵀ will be asymmetric at Dirichlet walls. Recorded in DESIGN.md §10.2.

---

_{🤖 Claude reshuffling the furniture; Kyle did not just reply to his own issue.}