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

[kylebeggs]

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:249-253 — eltype(κ.data) <: Number || throw(ArgumentError("diffusion coefficient must be a scalar-eltype Field, …"))
• src/operators/scaling.jl:47-54 — same guard on scaling(κ::Field)

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:134-135). 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.
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:343-347: "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:694, and for the packed forest at DESIGN.md:421 ("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:124-133 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, 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.

---

_{🤖 Filed by Claude, who read every stencil in the package and found them all suspiciously axis-aligned.}

[kylebeggs]

Corroboration on the "which anisotropy" question, and a data-ingestion note

From Kaboudian, Cherry & Fenton, Sci. Adv. 5:eaav6019 (2019) and the same authors' Computing in Cardiology 50 (2023).

The first reading is the right one to build

Paper 1 writes the cable equation as ∂V/∂t = ∇·(D∇V) − ΣI_ion and states plainly that "D is the diffusion tensor that contains the tissue's structure and rotational anisotropy." Paper 2 solves the same equation on voxelized patient anatomy. Neither treats the tensor as an extension — it is the baseline formulation for cardiac EP, and the scalar-κ case is the simplification.

That doesn't resolve what Dennis meant, but it does settle the priority between the two readings: physical anisotropy is the one that gates real monodomain results, which is the call this issue already made. Numerical-anisotropy (wide isotropic-error stencils) stays a genuine but separate want.

The ingestion path is better than the sketch suggests

The set!(D, x -> σ_t*I + (σ_l - σ_t)*(f(x)*f(x)')) sketch in the body assumes an analytic fiber field. In practice the papers do something that suits a Cartesian grid much better: fiber orientation is imported directly from diffusion-tensor MRI as image data onto the voxel grid, with no meshing or remeshing — the same route by which the segmentation itself arrives. Paper 1 is explicit that this is a reason to prefer the Cartesian representation.

So tensor_field(g) on a CartesianGrid is populate-able from real per-patient data by array assignment, not just from closures. That is worth stating in the issue because it changes what the acceptance test should look like: the interesting case is a measured, spatially varying, non-axis-aligned tensor field, which is exactly the regime where stage 3's cross-derivative terms and the M-matrix concern bite. The Strocchi 24-patient cohort is openly published and Carp2Cartesian converts the meshes to voxels, so that fixture is a download.

Sequencing note against #69

#69 now has a fifth candidate seam on the table — the phase-field method, where irregular geometry is imposed by carrying a smooth indicator φ and solving

∂u/∂t = (1/φ) ∇·(φ D ∇u) − reaction

Note the φD product. If phase-field wins there, the coefficient that Diffusion carries becomes φD — a scalar times a tensor — and the two features meet inside the same field. Not a conflict, but a reason to know #69's answer before finalizing stage 2's tensor face-averaging policy: harmonic averaging of φD across a geometric interface is a different question from harmonic averaging of D across a fiber-orientation jump, and the issue currently only anticipates the second.

Concretely: stage 1 (tensor_field + lifting the two eltype guards) is unaffected and can proceed independently. Stages 2 and 3 want #69 resolved first.

---

_{🤖 Claude checked whether the tensor was a luxury; the cardiac literature says it is the floor.}