Fix MPI ghost-cell bug in FormFunctionLocal

[awickert]

Summary

When WTM is run with more than one MPI process, FormFunctionLocal computes
transmissivity from global DMDA vectors (topo_vec, fdepth_vec, ksat_vec)
and accesses neighbour indices my_T[j][i±1] at subdomain edges.
DMDAVecGetArray on a global vector covers only the locally owned range,
so those neighbour accesses fall outside valid memory at every processor
boundary. In practice PETSc zero-initialises these vectors, so the out-of-
range reads return zero values for fdepth and ksat. With fdepth = 0,
depthIntegratedTransmissivity returns zero, the harmonic-mean face
transmissivity at the boundary collapses to zero, and an implicit no-flow wall
is imposed across every MPI subdomain edge.

Each subdomain then converges correctly — to the solution of a physically wrong
problem in which no groundwater crosses processor boundaries. The solver does
not report any error, and results look plausible at a glance because the
imposed walls are invisible in the output and groundwater diffusion is a local
process. However, large-scale flow paths that cross a processor boundary are
suppressed: the water table is anomalously high on the upstream side and
anomalously low on the downstream side of each internal boundary, with the
error decaying with distance from the wall. Single-process runs are unaffected.

Changes

src/CreateSNES.hpp / src/DMDA_array_pack.hpp / src/DMDA_array_pack.cpp
Replace T_vec (global) with T_local (local ghost vector). Add
topo_local, fdepth_local, ksat_local as local ghost vectors.
A new scatter_static_fields() function populates the global topo/fdepth/ksat
vectors from ArrayPack and scatters them to the local ghost vectors before
DMDA_Array_Pack acquires its GetArray locks — which would otherwise block
the scatter. Each DMGlobalToLocalBegin/End pair is completed sequentially
(the DMDA shares one PetscSF across all scatter operations, so interleaved
calls corrupt internal state).

src/transient_groundwater.cpp
FormFunctionLocal now uses the local ghost vectors and computes 1/T over
the full ghost range so all neighbour lookups are valid under MPI.
After each solve, arp.wtd is assembled across all ranks via
MPI_Allreduce(SUM) so that every rank holds the complete field before output.

src/WTM.cpp
Output is written by rank 0 only, avoiding a file race condition when multiple
ranks called saveGDAL on the same path.

src/CreateSNES.cpp
SNESSetType(SNESANDERSON) is set as the compiled-in default. The binary
previously required -snes_type anderson to be passed explicitly at runtime
(as in the benchmark scripts); without it PETSc defaulted to Newton, which
diverges immediately without a hand-coded Jacobian.

tests/ghost_cell/
Integration test: runs a 12×3 equilibrium domain with a factor-of-10 ksat
discontinuity placed at the MPI processor boundary, then compares 1-process
and 2-process outputs. With the fix they agree to machine precision (~10⁻¹⁴ m).
Without the fix the two halves are hydrologically decoupled and diverge by
O(10 m). Run with cd tests/ghost_cell && ./run_test.sh.

Caveats

The MPI_Allreduce output-gathering approach in transient_groundwater.cpp
works correctly when starting from wtd = 0 (each rank's non-owned cells are
zero and the sum assembles cleanly). For multi-cycle runs the same logic holds
because non-owned cells are zeroed before each reduction. This should be
reviewed if support for non-zero supplied initial water tables (supplied_wt)
under MPI is needed in future.

🤖 Generated with Claude Code