feat(distributed): ImplicitGlobalGrid halo backend for multi-node CartesianGrid

[kylebeggs]

Context

DESIGN.md §10.4 is now resolved: MDLA first (done), ImplicitGlobalGrid.jl for multi-node next.

MDLA is implemented in ext/MatrixFreeOperatorsMDLAExt.jl behind the _dist_scatter! / _dist_reduce! seam and covers single-node multi-GPU with the adjoint carried through. It cannot go multi-node, and it is CUDA-only. IGG is production-proven GPU-aware MPI halo exchange over a Cartesian MPI topology — exactly our uniform CartesianGrid case, multi-node, CPU+GPU.

This issue is to investigate first, then implement — the investigation gate below is real and may change the shape of the implementation.

Scope

• In: uniform CartesianGrid, multi-node, CPU + GPU, behind the existing distributed seam, as ext/MatrixFreeOperatorsIGGExt.jl. No new core dependency.
• Out: multi-node AMR (BlockForest). IGG assumes one fixed uniform Cartesian decomposition; the forest's halo_update! is intra-forest and a distributed forest needs a separate answer. Track that separately if/when it comes up.
• Not adopted: ParallelStencil.jl. We take IGG's halo only. @init_parallel_stencil(backend, precision, ndims) is module-level global state fixed at load, which is incompatible with keeping GPU support in a weak-dep extension; PS also brings no operator algebra, no linear-solve path, and no AMR. Our hot-stencil escape hatch remains KernelAbstractions @kernel (DESIGN §1.A).

Investigation gate (do this first)

IGG owns its own grid construction (init_global_grid) and array allocation (@zeros). The load-bearing question:

Can update_halo! be driven against MFO-owned Field storage, or is only the topology / MPI-Cartesian-communicator part reusable, with the exchange itself written directly on MPI.jl?

Prototype both far enough to answer it before committing to the extension's shape. If only the topology is reusable, the honest outcome may be "depend on MPI.jl directly and borrow IGG's decomposition logic" — which is a fine result for this issue, just a different one.

Secondary questions for the same prototype:

• How init_global_grid's global state interacts with partition_grid / prepare_distributed, which already compute local-vs-global ranges. One of the two has to own the decomposition; they must not disagree.
• Whether the adjoint contract survives: _dist_reduce! must remain the exact transpose of _dist_scatter! across ranks, including the ∂/∂κ scatter-add over owned cells. This is the part with no equivalent anywhere in the IGG/PS ecosystem, and it is the easiest thing to get silently wrong.

Calibrate the expected payoff — do not over-invest in overlap

For the driving workload (reaction–diffusion where a large local ODE dominates each node — cardiac monodomain, ~65 states and ~255 exp per node per step) only the scalar field crosses a boundary while the full state churns over the volume. For a 256³ subdomain that is roughly a 2800:1 byte ratio, so even a naive blocking exchange should cost low single-digit percent of a step.

IGG is being adopted for correctness and multi-node CPU+GPU portability, not because the halo is a bottleneck. Get a correct blocking exchange working and measured before reaching for @hide_communication-style overlap. If a benchmark later shows the exchange is actually material, that is the trigger to revisit — not an assumption to build on.

Done when

• Investigation gate answered and written down (IGG arrays vs. MFO storage)
• ext/MatrixFreeOperatorsIGGExt.jl implements _dist_scatter! / _dist_reduce!
• Adjoint identity ⟨Lx,y⟩ = ⟨x,Lᵀy⟩ verified across ranks, not just within one
• Parameter-gradient (∂/∂κ) parity against the single-rank result
• Multi-node test gated like test/multigpu (see test(ci): officially instantiate test/multigpu (>=3-GPU distributed tests) #30), skipped by default
• Halo cost measured on a realistic domain and recorded, so the overlap question is settled by data

---

_{🤖 Filed by Claude — Kyle supplied the judgement, I supplied the typing.}

[kylebeggs]

Blocked-by note: the IGG extension should implement its backend through the explicit selection API proposed in #47 (prepare_distributed(L, IGGBackend(...))), not by adding a second prepare_distributed method alongside MDLA's.

As things stand, prepare_distributed is exported from the core but its only method lives inside ext/MatrixFreeOperatorsMDLAExt.jl, so a second extension defining the same (operator, integer) signature would be ambiguous. #47 moves the entry point into the core behind an AbstractDistributedBackend tag, mirroring the existing AbstractJVPBackend pattern. Worth landing that first, or at least deciding its shape before this extension is written.

---

_{🤖 Claude again — connecting two issues so future-you doesn't have to.}

[kylebeggs]

Note from the PR #66 review (slice 2c, the Diffusion leaf on slabs): the slab path costs no coefficient communication only because the process running _slab_op holds the global κ — MDLA adapts the whole tree to the host before slicing. No IGG rank will hold it, so this backend owes one setup-time exchange of κ's Interface (cut-plane) ghosts.

DESIGN.md §7 now carries an Owed by the next backend paragraph with the shape that fits the existing seams: build each rank's leaf through the inner constructor from a rank-local κ (fill_coefficient_ghosts! mirrors the physical ghosts and leaves Interface ones untouched), then run one _dist_scatter!-shaped exchange over those planes at prepare time. Reusing _slab_op(::Diffusion) per rank with a rank-local κ would not do — it either throws on an out-of-range view or silently leaves the cut-plane ghosts at zero or the mirror (κ_I/2 face coefficients under ArithmeticMean, 0 under HarmonicMean), and no current guard catches that. Worth a checkbox under Done when.

---

_{🤖 Beep boop — Claude left this note; Kyle is not holding the global κ either.}