feat(operators): anisotropic diffusion ∇·(D∇u) with a full tensor D

Start with src/operators/diffusion.jl, src/operators/scaling.jl, src/Fields.jl, src/boundaries.jl, and the existing operator_diagonal implementation in src/operators/diagonal.jl. Run the CartesianGrid diffusion and adjoint tests before tracing the axis-aligned stencil path. Done means the scoped Diffusion leaf supports tensor fields with arithmetic face averaging on CartesianGrid, has an exact materialized-grid adjoint and updated diagonal, and includes the hostile anisotropy test and benchmark guidance described here.

From Dennis Ogiermann, 2026-08-14: "With AI we can add anisotropic 3D finite difference stencils (from the work by Nick Trefethen)."

Which anisotropy — needs Dennis

That sentence reads two ways, and they are different pieces of work:

• Physical anisotropy — a full conductivity/diffusivity tensor D, so cardiac fiber orientation enters the operator. This is what the issue below assumes.
• Numerical anisotropy — stencils whose truncation error is isotropic, so a wave on a Cartesian grid does not propagate faster along the axes than along the diagonals. That is closer to Trefethen's known line (group velocity and numerical anisotropy in difference schemes), and it means a wider isotropic-error Laplacian: 9-point in 2D, 19- or 27-point in 3D.

Both are worth having and neither implies the other. Filing under the first reading because that is what monodomain needs to express its own physics; if Dennis meant the second, this issue splits.

What's missing

κ is scalar-only, and it is enforced at both entry points:

• src/operators/diffusion.jl:545-549 — eltype(κ.data) <: Number || throw(ArgumentError("diffusion coefficient must be a scalar-eltype Field, …")), with the same guard on the BlockForest method at :608-612
• src/operators/scaling.jl:47-54 — same guard on scaling(κ::Field), and on scaling(κ::AbstractBlockField) at :55-62

So the fused Diffusion leaf is closed to a tensor, and so is the composed escape hatch divergence(g) * scaling(K) * gradient(g). There is no tensor_field constructor, and _scalar_eltype handles Number and SVector only (src/Fields.jl:247-248). This contradicts DESIGN.md:145-149, which already promises a tensor field over Array{SMatrix{…}} with "the backing array's eltype is the only thing that changes."

The face-averaging policies are also scalar binary ops — ArithmeticMean and HarmonicMean at src/operators/diffusion.jl:9-25 are (a+b)/2 and 2ab/(a+b).

Why

Cardiac conduction is fiber-aligned: D = σ_l ff^T + σ_t (I − ff^T), with σ_l/σ_t commonly 3–10. Every monodomain result this package produces — including examples/niederer_benchmark.jl, whose whole point is a standard cardiac verification problem — is currently restricted to isotropic conduction because the operator cannot hold the tensor. That is a physics gap, not an ergonomics one.

Sketch

g = CartesianGrid(((0.0, 1.0), (0.0, 1.0), (0.0, 1.0)), (n, n, n); bc = …)
D = tensor_field(g)                       # Field over Array{SMatrix{3,3,T,9}}
set!(D, x -> σ_t * I + (σ_l - σ_t) * (f(x) * f(x)'))
L = diffusion(g, D)                       # ∇·(D∇u), cross terms included

Staged, because the third stage is where the actual work is:

1. tensor_field + _scalar_eltype(::SMatrix) + lift the two eltype guards. Mostly plumbing — with one exception that must land before any guard is lifted: adjoint_operator(::ScalingOp{<:Field}) returns Field(conj.(κ.data), κ.grid) (src/operators/scaling.jl:78-82, and the AbstractBlockField method at :83-91 does the same per block). Elementwise conj is the adjoint of a scalar coefficient, not of a matrix one — for an SMatrix κ it must be adjoint.(κ) (the pointwise transpose for real κ). Lift the guard without that and scaling(K) with a non-symmetric tensor silently gets a wrong adjoint, which is exactly what "adjoints are declared, never assumed" forbids. Open PR #107's generic scaling (GHScalingOp in its title, scaling(f, args...) in the diff) is not a shortcut here: it requires f to return a real scalar, trusts the user on that, and declares itself self-adjoint unconditionally (isselfadjoint(S::ScalingOp) = true, adjoint_operator(S::ScalingOp) = S), so it must not become the tensor-κ path. (That PR also renames today's single-coefficient op to CoeffScaling, an alias of its new ScalingOp; if it merges first, the conj → adjoint fix moves to adjoint_operator(::CoeffScaling{<:Field}) unchanged.)
2. Tensor-aware face averaging. Arithmetic generalizes componentwise. Harmonic does not — the scalar 2ab/(a+b) has several inequivalent matrix generalizations, and which one preserves flux continuity across a fiber-orientation jump is a real question, not a naming choice.
3. The cross-derivative terms ∂/∂xᵢ(D_ij ∂u/∂xⱼ) for i ≠ j. Everything below is about this.

Constraints

• No leaf reads a diagonal neighbour today, and the infrastructure says so out loud. src/transfer.jl:514-518: "Per-leaf stencil adjoints must leave corner ghosts exactly zero (true for axis-aligned stencils, whose transposes never scatter there) … A future cross-derivative leaf needs a corner-aware exchange." Same statement for slabs at src/distributed.jl:650-655, and for the packed forest at DESIGN.md:425 ("last dim wins; no axis-aligned stencil reads them"). So the forest and distributed corner ghosts are currently garbage by design, and a cross-derivative leaf on those paths reads them.
• A plain CartesianGrid is fine already. Corner ghosts are filled consistently as ghost-of-ghost values at src/boundaries.jl:107-108. That makes CartesianGrid-only the obvious first slice, with BlockForest / slab support and its corner-aware exchange and corner-aware adjoint fold as an explicit follow-up.
• Adjoint. The cross terms have to be transposed exactly, including how they touch corner cells near a boundary. ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ on a materialized small grid is the gate, per the package's standing rule that adjoints are declared and never assumed.
• operator_diagonal. src/operators/diagonal.jl:125-135 handles Diffusion on CartesianGrid only; the tensor form changes the diagonal and multigrid depends on it.
• Monotonicity. A standard 7-point-plus-cross-terms discretization of a strongly anisotropic tensor loses the M-matrix property once the off-diagonal entries dominate, which breaks the discrete maximum principle and degrades Jacobi/Gauss–Seidel smoothing exactly where multigrid needs it most. This deserves a named test with a deliberately hostile σ_l/σ_t and a fiber direction at 45° to the grid, not a footnote.
• The ntuple trap. CLAUDE.md records that a stencil reading two arrays must not unroll dimensions with ntuple(Val(N)) — the closure stops inlining, vectorization dies, and a 256² sweep went from 33 µs to 260 µs with identical numerics. A tensor stencil reads more arrays than that. Unroll by recursion over Val(D), following _diff_axes in src/operators/diffusion.jl:377-383, and benchmark before believing anything.

Scope note

CartesianGrid only, Diffusion leaf only, arithmetic face averaging only, for the first slice. Forest/distributed corner exchange, harmonic tensor averaging, and the wider isotropic-error stencils from the second reading above are all separate follow-ups.

—
🤖 Claude added a transpose warning before anyone learns it the hard way.
_{claude-opus-5-5[1m] (high) · claude-code 2.1.285 · unreviewed}