FreeSurface: prescribed surface datum loses ~7% net flux in parallel (surface-array <-> field round trip)

[lmoresi]

Summary

uw.systems.FreeSurface's consistent (third) solve prescribes u.n = ũ_n so the surface
stays a material boundary. The datum field _un_target carries a ~6-7% net flux on 2
ranks versus ~0.06% in serial — a ~100x degradation that is independent of how the
surface mean is removed, so it is not the demean.

An incompressible interior over a closed no-slip base can neither gain nor lose volume, so
∮ ũ_n ds = 0 is required. In parallel it is badly violated.

Measurement

Isoviscous annulus (r_in=0.5, r_out=1.0, cellSize=0.1), gap-Ra 3e4, rho_g=5e5,
3 steps at advect_scale=10. Net flux = |∮ ũ_n ds| / ∮|ũ_n| ds, both by exact FE
uw.maths.BdIntegral over the surface (not the manager's own nodal sum, so it cannot pass
by construction):

ranks	node-count demean	arc-length demean
1	8.96e-3	5.63e-4
2	7.06e-2	5.75e-2

What is NOT the cause

The ring machinery and the arc-length quadrature are parallel-exact. On np2 the
gathered ring has 66 entries for 64 distinct nodes (the 2 partition-seam vertices appear on
both ranks) and the trapezoidal rule splits their weight correctly:

• sum(weights) = 6.2831853072 = 2*pi*r_out to machine precision
• weighted mean of 1 = 1.000000
• weighted mean of cos(k*theta), k = 1,2,3 = ~1e-16

Pinned by test_freesurface_ring_quadrature_is_exact_in_parallel (runs on any rank count).

So the loss is in the surface-array <-> field round trip across the partition seam: the
demeaned ring array is flux-free by construction, but reading the values back out of
_un_target is not. Suspects, in order:

1. _field_rows matches local surface coords to local field rows with a KD-tree; behaviour
   at seam vertices that are owned by one rank and ghost on the other is unverified.
2. Writing _un_target.array[rows] = u_tilde sets local values only — a ghost/owner scatter
   may overwrite the seam DOF, or the two ranks may disagree on its value.
3. _surface_node_coords returns seam vertices on both ranks (by design, for the ring);
   the field write does not obviously want that duplication.

Why it matters

This is pre-existing — it predates the arc-length demean and affects the penalty path
too. It was never caught because the datum's net flux had not been measured in parallel;
the weak penalty absorbs the incompatibility quietly, where a strong rotated constraint is
over-determined by it. With consistent_constraint="strong" (the serial-validated default,
6.4x better on the material-boundary error and leak-free at 7e-5 in serial) a 7% incompatible
datum in parallel is exactly the over-determined case, so parallel runs are the risk here.

Repro

tests/test_1070_free_surface_plume.py::test_freesurface_prescribed_rate_is_flux_free and
::test_freesurface_strong_constraint_beats_penalty are currently skipped on mpi.size > 1
with this reason; remove the _PARALLEL_DATUM_DEFECT skip to reproduce.

Underworld development team with AI support from Claude Code