Stokes_Constrained: parallel-correctness residuals (partition-dependent [p,λ] gauge + constraint-assembly drift + lossy knockout)

[lmoresi]

Summary

While root-causing the #248 Zhong-benchmark velocity catastrophe (now fixed by #242 — it was a broken SphericalShellInternalBoundary Lower/Internal label bug), three smaller but real parallel-correctness issues in Stokes_Constrained (in-saddle multiplier free-slip) were isolated with reliable integral diagnostics. Tracking them here for a focused fix.

Method note: point evaluation (uw.function.evaluate) is unreliable for serial-vs-parallel field comparison here (gave 134%-of-|v| pointwise noise while volume L2 differed 0.4%). Use partition-independent integrals only.

Ruled out: velocity rigid-rotation nullspace (angular momentum ∫r×v dV reproducible to ~1e-5 serial vs np=8); iterative under-convergence (serial & np=8 each bit-identical at tol 1e-7 vs 1e-11 — both fully converge, to different answers).

1. [p,λ] / topography gauge is partition-inconsistent (the nullspace bug)

The combined constant-(pressure, multiplier) gauge (_build_block_gauge_nullspace_vector) is a true velocity null mode, but its level is partition-dependent: mean pressure = 3.333 (serial) vs 2.983 (np=8), 10.5%. Because topography is read from the multiplier, this corrupts topography across ranks unless mean-stripped. set_pressure_gauge("Upper", 0.0) (#250, global BdIntegral callback) pins it consistently → 0.4%, bit-identical on velocity.
Fix: make the gauge partition-reproducible by construction (default-on pin, or enforce topography(reference="mean")); add a parallel regression asserting meanP/topography match serial to ~1e-6. Highest priority — topography correctness, fix already exists.

2. Residual ~0.4% velocity partition-dependence in constraint/augmentation assembly (separate, NOT a nullspace)

With the gauge pinned, velocity L2/|Ut|/|Ub| still differ 0.36/0.37/0.41% serial vs np=8 (np=8 bit-identical with/without the pin). Fully converged ⇒ the assembled operator/RHS differs ~0.4% by partition — locus is the add_constraint_bc / augmented-Lagrangian (augmentation_base=1e4) boundary assembly, not a nullspace. 3D-specific (2D test_1063 is bit-identical to 1e-9; boundary areas match exactly, so it's the v·n / augmentation term at inter-rank seams on the 2-D free-slip surface).
Fix: find and fix the partition-dependent constraint/augmentation assembly; goal = converged velocity partition-independent to round-off.

3. Interior-multiplier knockout is lossy (docstring says "lossless")

_constrain_interior_multipliers_in_section (gated by _reduce_interior_multiplier, default True) pins interior multiplier DOFs. Disabling it changes the serial answer by 2% (L2/Ut) to 5% (Ub) — so the pinned DOFs are not inert. DOF count is partition-identical (5170 at np=1/2/4/8), so not gross loss. (Off makes parallel spread worse — 1.7% at CMB — so it's a distinct issue from #2, consistent with block-error trade-off.)
Fix: determine whether the boundary-trace closure drops DOFs that matter (or "inert interior multiplier" is the wrong model); fix or re-document. Also: monolithic MUMPS LU on the constrained system gives wrong physics (|Ut|=4e-3 vs validated 1e-2) and segfaults at np=8 — capture as known-bad or fix.

Reproduction

~/+Simulations/repro_248_internal_load.py (lifts the test_1064 Zhong setup; reports volume L2, |Ut|/|Ub|, angular momentum, meanP — all partition-independent). Args: <method> <write|check> <solver> <tol> <knockout on|off> <gauge none|pin>. Build/run in the amr-dev pixi env; parallel needs --with-mpi.

Related: #248 (velocity half resolved by #242; topography half open), #244 (constrained slow in parallel).

cc @gthyagi

Underworld development team with AI support from Claude Code

[lmoresi]

#265 merged to development (fgmres + unpreconditioned norm, rotation-gauge projection, guarded auto pressure gauge, knockout audit). Serial correctness verified. Remaining to close: the np=8 partition-independence numbers (tangential velocity 1.2% → 2.5e-14 were measured on a remote box) — pending re-run on the many-core machine. Keep open until that lands.

Underworld development team with AI support from Claude Code

[lmoresi]

PR #265 — Parallel verification report

Environment

• Branch: development @ <COMMIT> (PR Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 merged 2026-06-24 04:19 UTC)
• PETSc: 3.25.2 (conda-forge, unpatched — dev pixi env)
• petsc4py: 3.25.2
• MPI: MPICH 5.0.1
• Machine: 48-core Linux (this /remote-control session)

Headline

PR #265 achieves partition-independence for the external-boundary
constrained free-slip case at the tolerances claimed. test_1063 passes once
its GOLDEN reference is regenerated under PETSc 3.25 (the existing GOLDEN was
baked under 3.24 and shifts ~2-7e-5).

The 3D SphericalShellInternalBoundary (Zhong) case still shows a
partition-dependent residual at the ~3e-6 level in velocity — exactly the
"known PETSc limitation" PR #265 documents but does not fix. The fix lives in
petsc-custom/patches/plexfem-internal-boundary-ownership-fix.patch, which is
applied to the Underworld-bundled PETSc 3.24/3.25 but not present in
conda-forge PETSc.

1. test_1063_constrained_freeslip_parallel.py — annulus, external boundary

Re-derived serial reference on PETSc 3.25.2 (test currently asserts vs PETSc
3.24-era GOLDEN, hence the test failures below are GOLDEN drift, not solver
regressions):

diag	np=1 (new ref)	np=2	np=4	np=8	max rel diff vs np=1
iso L2	6.194403143796e-01	6.194403144137e-01	6.194403143988e-01	6.194403144573e-01	1.3e-10
ti L2	3.926249050981e-01	3.926249050014e-01	3.926249050691e-01	3.926249048306e-01	6.8e-10
gauge L2	6.194402844557e-01	6.194402844556e-01	6.194402844539e-01	6.194402844573e-01	2.9e-12
gauge meanP	-5.7e-11	-3.8e-10	-1.9e-10	-4.4e-11	(~1e-10 = 0)
iso topo	3.785960214024e+01	(via test)	3.785959770353e+01	3.785960050885e+01	1.2e-7
ti topo	3.707615257402e+01	(via test)	3.707615390743e+01	3.707616948635e+01	4.6e-7
gauge topo	3.781685363609e+01	(via test)	3.781685353448e+01	3.781685364877e+01	3.4e-9

• velocity L2 partition-independent to < 1e-9 across np=1,2,4,8 (target met);
• raw mean pressure pinned to ~1e-10 ≈ 0 (auto pressure gauge fires, partition-reproducible);
• mean-stripped topography partition-independent to < 1e-6 (target met).

Recommendation (action item): regenerate GOLDEN/GOLDEN_GAUGE in
tests/parallel/test_1063_constrained_freeslip_parallel.py from PETSc 3.25.2
serial (values above). Without this, the test reports a failure that is
actually a PETSc-version GOLDEN drift, not a solver regression.

2. Zhong 3D SphericalShellInternalBoundary — internal-boundary natural BC

Setup matches tests/test_1064_constrained_spherical_shell_response.py: cellSize
1/8, qdegree=2, degree=1, Stokes_Constrained with constraint_bc on Upper +
Lower, natural BC Y_l=2 ê_r on the internal boundary, tol 1e-7. Default
solver (PR #265's fgmres + unpreconditioned norm).

Wall-clock and partition scaling (single MPI job at a time, MPICH 5.0.1):

np	wall-clock	speedup	efficiency	converged?
1	4435.2 s	1.00	1.00	yes (SNES=5)
8	798.7 s	5.55	0.69	yes (SNES=5)
16	444.6 s	9.98	0.62	yes (SNES=5)
32	250.8 s	17.7	0.55	yes (SNES=5)
48	—	—	—	SEGV on rank 33 (see below)

diag	np=1 (serial)	np=8	np=16	np=32
velocity L2 √∫v·v	1.806922688e-2	1.806905456e-2	1.806799697e-2	1.806873740e-2
surface |Ut|	1.005469276e-2	1.005464134e-2	1.005421623e-2	1.005452825e-2
CMB |Ub|	1.180547582e-2	1.180520050e-2	1.180427036e-2	1.180459942e-2
angular momentum ∫r×v	9.917e-05	9.921e-05	9.926e-05	9.921e-05
raw mean pressure meanP	3.64094e-3	3.64121e-3	3.64039e-3	3.64173e-3

Relative differences vs serial:

np	rel(vL2)	rel(Ut)	rel(Ub)	abs(meanP)
8	9.5e-6	5.1e-6	2.3e-5	2.7e-7
16	1.3e-5	4.7e-5	1.0e-4	5.5e-7
32	2.7e-5	1.6e-5	7.4e-5	7.8e-6

• Serial values reproduce the Zhong reference (Ut=1.006e-2, Ub=1.186e-2) to
  better than 0.1% — solver correctness is good.
• All non-serial partition residuals (vL2, Ut, Ub) stay in 1e-6…1e-4 —
  exactly the documented "known PETSc limitation" PR Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 leaves behind. They
  do NOT grow systematically with np, consistent with a fixed-magnitude
  internal-facet over-counting that the elliptic solve smooths. Pre-Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265
  numbers (issue Stokes_Constrained: parallel-correctness residuals (partition-dependent [p,λ] gauge + constraint-assembly drift + lossy knockout) #254) showed ~0.4% on velocity / |Ut| and ~10.5% on meanP
  — improvement of two to three orders of magnitude.
• Auto pressure gauge: surface-mean is pinned on the first constraint boundary,
  so volume-mean meanP we report here is gauge-orthogonal. Absolute drift
  stays in 1e-7…1e-5 across np=8/16/32, vs the pre-Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 ~10.5%.
• The remaining ~1e-5 in velocity will only fall to round-off once the PETSc
  patch is in the build (see "PETSc patch status" below).

np=48 SEGV

[33]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation,
  probably memory access out of range
Abort(59) on node 33 (rank 33 in comm 0)

Rank 33 crashes during solve setup (after mesh built / vars created / solver configured, before solve done). Most likely an empty / near-empty partition
on this fine decomposition of a small mesh (cellSize 1/8 → roughly 2-3k cells,
spread over 48 ranks ≈ ~50 cells/rank, easily yielding ranks with no boundary
facets). Not introduced by #265 (no nearby changes to partitioning), but worth
filing as a separate issue if anyone routinely runs Constrained at np ≳ 32 on
small meshes.

PETSc patch status (the user's question)

petsc-custom/patches/plexfem-internal-boundary-ownership-fix.patch was
generated against PETSc 3.24.5. Tried patch --dry-run -p1 against PETSc
3.25.2 source pulled from https://gitlab.com/petsc/petsc -b v3.25.2:

• Hunk 1 (ghost-facet filter in DMPlexComputeBdIntegral) — FAILED
  because PETSc 3.25.2 already has this fix upstream (see plexfem.c:2876).
• Hunks 2-4 (support[0] → support[key.part] in
  DMPlexComputeBdResidualSingleByKey) — apply with line-number offsets;
  the buggy code is still present at plexfem.c:5002-5008.
• Hunk 5 (ghost-facet filter in DMPlexComputeBdResidual_Internal) — applies
  with fuzz; PETSc 3.25.2 already has this fix upstream (plexfem.c:5139).
• Hunk 6 (ghost-facet filter in DMPlexComputeBdJacobian_Internal) —
  FAILED; this fix is still missing in 3.25.2.

So the patch is still needed for PETSc 3.25 but is smaller now:

• the support[0] → support[key.part] fix in DMPlexComputeBdResidualSingleByKey
  (plexfem.c:5002-5045 in 3.25.2),
• the same fix in DMPlexComputeBdJacobianSingleByKey (plexfem.c:6026-6080 in
  3.25.2),
• the ghost-facet filter at the top of DMPlexComputeBdJacobian_Internal.

Hunks 1 and 5 of the old patch are redundant against 3.25 and can be dropped.
A refreshed plexfem-internal-boundary-ownership-fix-v3.25.patch should be a
straightforward port; the upstream MR draft at
petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md is
unchanged in intent.

Action items

1. Update tests/parallel/test_1063_constrained_freeslip_parallel.py GOLDENs
   to PETSc 3.25 (separate small PR). Values from this session above.
2. Refresh the PETSc patch for 3.25 and either:
   • (preferred) push it upstream (the MR description already exists at
     petsc-custom/patches/PETSC_MR_DESCRIPTION_INTERNAL_BOUNDARY_OWNERSHIP.md);
   • (interim) carry the patched custom-PETSc build through CI for the
     constrained tests, or ship a patched conda-forge PETSc.
3. Add test_1063 to the CI test script — currently scripts/test.sh --p 2
   only runs tests/parallel/test_075*py. The constrained-parallel test that
   Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265 added is not in CI.
4. File a follow-up for the np=48 SEGV on Stokes_Constrained (default
   field-split) — likely a small-mesh fine-decomposition issue, not a Stokes_Constrained parallel correctness: gauge, convergence, knockout, rotation-gauge #265
   regression.

[lmoresi]

Triage note (2026-07-23): scope narrowed to the single remaining item

Of the three residuals this issue tracked, #265 (merged 2026-06-24) fixed the partition-dependent [p,λ] gauge and the lossy knockout, and the 2026-07-04 verification report above confirms partition-independence for the external-boundary constrained case at the claimed tolerances (tangential velocity 1.2% → ~1e-14).

What remains open here: the ~3e-6 partition-dependent velocity residual on the 3D SphericalShellInternalBoundary (Zhong) configuration — the documented upstream PETSc internal-boundary limitation, related to the internal-facet orientation work in #327 (UW3-side fix in PR #361; upstream DMPlexOrientLabel petsc4py support in MR !9418). Keeping this open as the tracking issue for that last item only.

Underworld development team with AI support from Claude Code

[lmoresi]

Closed into #651, which consolidates the Stokes_Constrained backlog into one place with the status of each thread re-measured rather than copied forward. Nothing from this issue is dropped — the open items, the reproduction command, the ruled-out list and the method note about evaluate all carry across.

Underworld development team with AI support from Claude Code