feat(grids): finite differences on non-trivial geometry — pick the seam

[kylebeggs]

From Dennis Ogiermann, 2026-08-14: "There are papers from Pras Pathmanathan showing how to adopt FD to non-trivial geometries."

Where this sits relative to #61

#61 is the first-order answer; this is the second-order one. A binary active-cell mask makes a cell either in or out, which puts the domain boundary on a staircase — and a staircase boundary caps the geometric error at first order no matter how formally accurate the interior stencil is. For the hemisphere-in-a-box case that motivates #61 that may be perfectly adequate. For a curved tissue or ventricle boundary where the boundary flux is the physics, it is not, and the difference shows up as a convergence-rate ceiling rather than as a visible bug.

This issue is to decide, with numbers, which seam MatrixFreeOperators.jl adopts for non-box domains — not to build one speculatively.

What exists today

Nothing. There are exactly two AbstractGrid subtypes and both are axis-aligned boxes:

• CartesianGrid (src/Grids.jl:42) — uniform derived spacing (src/Grids.jl:67), no per-cell metric, no cell-volume or face-area accessor.
• BlockForest (src/BlockForest.jl:33) — a forest over the same box, whose constructor requires an exact root tiling (src/BlockForest.jl:49-54). AMR changes resolution, never participation.

Boundary conditions are per-dimension (low, high) face pairs attached to the grid (src/Grids.jl:56, src/boundaries.jl:12-71), implemented as ghost mirroring through _bc_sign (src/boundaries.jl:69-71) — six faces of a box, not per-cell and not per-face-fragment. The only mask in src/ is bounds-masking inside adjoint gathers (_maskedget, src/operators/abstract.jl:212), which is unrelated.

DESIGN.md is silent on this. It is not in §1 locked decisions, not in §9 forward-looking seams, not in §10 open decisions. The only relevant sentence is DESIGN.md:9-11 — "scoped initially to structured grids but built so any grid amenable to matrix-free application fits later" — plus the convention-based grid interface at DESIGN.md:113-114. #54 explicitly lists "unstructured, curvilinear, cut-cell, or embedded-boundary mesh framework" as out of scope, which is the right call for #54 and is exactly the boundary this issue is asking to revisit deliberately.

Candidate seams to evaluate

1. Staircase active-cell mask (feat(grids): active-cell mask so stencil operators can skip inactive cells #61). Cheapest, already scoped, keeps every stencil array-level. First-order geometric error.
2. Shortley–Weller / cut-cell one-sided stencils. Recovers second order for Dirichlet data on a curved boundary by shortening the arm that crosses it. Needs a per-cell distance-to-boundary and a per-cell stencil.
3. Generalized finite differences from a local neighbour set — weights solved per node against a Taylor system, typically via a pseudo-inverse. This is the family Dennis is pointing at. Most general, most expensive per cell, and the one that most clearly stops being "a stencil".
4. Level-set / embedded boundary with a signed-distance field carried as a Field and volume/aperture fractions derived from it. Composes best with AMR, and shares machinery with the geometry accessors feat(fv): add conservative finite-volume operators without sacrificing matrix-free AMR #54 already proposes.

The tensions each option has to survive

These are the package invariants any answer must not quietly break, and they are the reason this is an investigation rather than a PR:

• Array-level leaf bodies. CLAUDE.md requires leaf bodies stay broadcast/slicing so they are device-agnostic and AD-friendly with no per-backend code. A per-cell variable stencil is a direct hit on that, and the KernelAbstractions @kernel escape hatch is supposed to be per-operator, not the default.
• Halo width. Every current stencil is compact — one ghost layer. Options 2–4 can widen the neighbour set, and halo width is baked into grid construction and into the exchange schedule.
• The adjoint is the gate. Whatever apply! does, apply_adjoint! must remain its exact discrete transpose. A shortened boundary arm has to be shortened identically on the way back, or ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ fails. That test, not a plot, decides whether an implementation is correct.
• operator_diagonal becomes per-cell. Multigrid's smoothers and the coarsest-level LU (src/multigrid.jl:288) consume it, and the current implementations assume a diagonal that varies only by BC face.
• Conditioning. Small cut cells produce a small diagonal and a badly conditioned operator; that is the classic cut-cell failure mode and it lands on the Krylov path, not on the operator.
• Don't abstract ahead of evidence. DESIGN.md:862-870 — an abstraction is "deferred until a second topology backend … actually exists … Do not introduce it speculatively." This issue is allowed to produce a recommendation and one concrete implementation, not an AbstractGeometry layer.

Deliverable

• A short written comparison of the four options against the constraints above, naming the recommended one and why.
• A convergence study that makes the argument concrete rather than rhetorical: Poisson on a disk (2D) and a sphere (3D) with a manufactured analytic solution, measuring observed order for the staircase mask versus one chosen alternative on the same grids.
• The adjoint dot-product identity holding for whatever is implemented, on a materialized small problem.
• A DESIGN.md §10 entry recording the decision, so the next person does not re-open it from scratch.

Out of scope

• Building an unstructured or curvilinear mesh framework. The question is how far a structured grid can be pushed, not how to stop being one.
• Masked BlockForest or masked distributed slabs. Uniform CartesianGrid first, same boundary feat(grids): active-cell mask so stencil operators can skip inactive cells #61 and feat(fv): add conservative finite-volume operators without sacrificing matrix-free AMR #54 drew.
• Moving or deforming geometry.

Open question — needs Dennis

The exact Pathmanathan references. A search turned up a generalized-finite-difference bidomain method on arbitrary irregular meshes that matches the description, but I could not confirm the author list, and there is more than one plausible line of work here. Worth getting the citations before spending a day on the wrong family.

---

_{🤖 Claude drew the staircase, measured it, and concluded stairs are not a curve.}

[kylebeggs]

A fifth candidate seam, and evidence on two of the existing four

From Kaboudian, Cherry & Fenton, Sci. Adv. 5:eaav6019 (2019) and the same authors' Computing in Cardiology 50 (2023). Both solve exactly this problem — monodomain on a real heart in a Cartesian box — and between them they use two different seams.

5. Phase-field — the option this issue is missing

Paper 1 handles irregular heart geometry with the phase-field method (Fenton, Cherry, Karma & Rappel, Chaos 15:013502, 2005). A smooth indicator φ ∈ [0,1] is carried as a field, and the equation is rewritten

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

which enforces no-flux on the φ level set with no cut-cell geometry anywhere.

In our algebra that is, today, with no new leaf:

scaling(inv.(φ)) * diffusion(g, φ .* κ)

Against this issue's own constraint list, it is the only candidate that clears every row:

constraint	phase-field
array-level leaf bodies	✅ unchanged — φ is a coefficient field, the stencil is the existing Diffusion one
halo width	✅ still one ghost layer, still compact
adjoint is the gate	✅ composition of two already-declared adjoints; nothing new to derive
operator_diagonal per-cell	✅ falls out of the Diffusion diagonal times 1/φ
cut-cell conditioning	✅ no small cut cells exist — this is the failure mode it avoids
don't abstract ahead of evidence	✅ zero new abstraction; it is a usage, not a type

More usefully, it collapses the framing: #61's binary mask is the φ ∈ {0,1} special case. #61 and this issue are one field at two smoothness settings, not two features. That is worth deciding deliberately, because it changes what #61 should build — an active::Field{Bool} and a φ::Field{T} want to be the same seam.

Two honest caveats, so this doesn't read as a free lunch:

• φ does not replace the mask, it needs one. 1/φ blows up outside tissue, so φ must be floored or the domain masked anyway. Phase-field is feat(grids): active-cell mask so stencil operators can skip inactive cells #61 plus a smoothing, not an alternative to it.
• Accuracy is governed by the interface width ξ, which must span a few cells. It trades staircase geometric error for a resolution requirement, and the standard critique is that resolving ξ eats the grid savings. Whether that is a net win is exactly the convergence study this issue already asks for — Poisson on a disk/sphere with a manufactured solution, observed order, staircase vs. phase-field on the same grids. That comparison is now the cheapest one to run, because both sides are expressible with existing leaves.

Evidence on option 1: the group Dennis pointed at ships the staircase

Paper 2 is from the FDA Division of Biomedical Physics and cites Pathmanathan's own credibility-assessment work (Galappaththige, Gray, Costa, Niederer & Pathmanathan, PLoS Comput Biol 18(10):e1010541, 2022). For patient-specific cardiac EP on real anatomy at 0.026 cm they use a voxelized staircase with mirrored ghost cells — not generalized finite differences.

Read that as option 1 is shippable for this application class, not option 1 is right. Activation-time metrics are far less sensitive to boundary geometric error than a Poisson flux is, and cardiac EP is a wave-arrival problem where the boundary mostly says "don't leak." For the microvascular O₂ case in #61, where the boundary flux is the quantity, that transfer is not automatic. It narrows the question rather than closing it.

It also narrows the "needs Dennis" item at the bottom of this issue. Pathmanathan appears in this literature as a verification & credibility author on voxelized-Cartesian cardiac work, not as the author of a generalized-FD scheme. So the sharper question for Dennis is: are you thinking of the FDA voxelization line, or of a distinct generalized-FD-on-irregular-meshes paper? — rather than an open-ended request for citations.

A concrete citation for option 4

Paper 2's zero-flux reference is Chai, Luo, Wang, Zheng & Fan, "Imposing mixed Dirichlet–Neumann–Robin boundary conditions on irregular domains in a level set / ghost fluid based finite difference framework," Computers & Fluids 214:104772 (2021) (DOI).

That is option 4 with all three BC types worked out, in a finite-difference framework, and it is simultaneously the missing prior art for #62 (Robin). If option 4 is ever evaluated seriously, that is the paper to start from.

Suggested revision to the deliverable

The convergence study this issue specifies is now cheaper and more informative as a three-way: staircase mask vs. phase-field vs. one sharp-interface option (Shortley–Weller), on the same disk/sphere grids. The first two need no new leaves, so the cost is dominated by the third — which is also the one most likely to be dropped once the first two are measured.

---

_{🤖 Claude went looking for a curve and came back with a smooth indicator function.}