Insert essential boundary values in the rotated solve's field copy-back

[lmoresi]

Fixes #497.

Essential-BC DOFs are absent from the global vector, so the rotated path's per-field globalToLocal scatter left them at zero in the output fields wherever the Dirichlet datum was non-zero. The solve itself was always correct — the residual path inserts the values internally — but every field-based diagnostic (projection, uw.maths.Integral, uw.function.evaluate, renders) read a garbage boundary strip. Measured on the split-node fault study: far-field sigma_xy of 0.79 for a true 1.00, and a projected stress row off by 3x one fault length out. Every existing rotated test used homogeneous Dirichlet on the other walls, so zero was accidentally the right value — which is how this survived.

Fix: petsc_dm_insert_boundary_values in cython/petsc_discretisation.pyx — wrapping DMPlexInsertBoundaryValues, which petsc4py does not expose; the same insertion the consistent-boundary-flux paths already use internally (cf. #407/#411) — called per field after the scatter in _finalize_rotated_solution.

Not affected (both live in the global vector): the rotated boundary's own prescribed datum (add_rotated_freeslip_bc(conds != 0, ...), the free-surface path), and the sigma_nn / dynamic-topography reaction recovery.

Regression test: test_rotated_solve_fields_carry_inhomogeneous_dirichlet_walls in tests/test_1018_rotated_freeslip.py — an inhomogeneous lid with rotated free-slip sides, asserting the FIELD carries the wall values; it fails without the fix. The change is a verbatim cherry-pick of 08f8603 from feature/fault-split-node, where the full rotated suite (19 tests) and the fault-contact suites run green with it.

The tell that found it, for the record: the divergence theorem — Integral(2 e_xy) dA against its own boundary term.

Underworld development team with AI support from Claude Code

[lmoresi]

Adversarial self-review, scope and hazards:

• Scope: three files; the only behavioural change is one insertion call per field in _finalize_rotated_solution, on the rotated path only. The standard SNES path is untouched (it never had the gap — PETSc's own copy-back inserts).
• Could the insertion be wrong for homogeneous BCs? No: inserting zeros over zeros is idempotent, and the pre-existing rotated suite (18 tests, all homogeneous-wall) passes bit-identically with the fix in place.
• Could it double-apply a rotated datum? No: the rotated datum rows are velocity DOFs present in the global vector (constrained algebraically, not eliminated), so DMPlexInsertBoundaryValues does not touch them; the nonlinear-datum test asserts those field values and passes.
• Parallel: the insertion is per-field on the sub-DM local vector, rank-local by construction; verified inside the fault-contact work at np=2 where the affected fields feed collective diagnostics.
• What we did not do: audit every historical rotated result for inhomogeneous-Dirichlet contamination. The scan is cheap where suspected — the divergence-theorem check flags it immediately — and no production study we know of combined the two.

Underworld development team with AI support from Claude Code

[Copilot]

Pull request overview

This PR fixes a rotated free-slip solve copy-back defect where inhomogeneous essential (Dirichlet) boundary values were missing from the per-field local vectors after scattering from the composite global solution, leading to incorrect field-based diagnostics near Dirichlet boundaries.

Changes:

• Adds a Cython shim petsc_dm_insert_boundary_values() wrapping DMPlexInsertBoundaryValues() (not exposed by petsc4py) to insert essential BC values into local vectors.
• Updates the rotated solve finalization to call this insertion after the global-to-local scatter for each field.
• Adds a regression test ensuring inhomogeneous Dirichlet wall values are preserved in the output field after a rotated solve.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 1 comment.

File	Description
tests/test_1018_rotated_freeslip.py	Adds a regression test that fails if rotated copy-back drops inhomogeneous Dirichlet values.
src/underworld3/utilities/rotated_bc.py	Inserts essential BC values into each field local vector after scatter in _finalize_rotated_solution.
src/underworld3/cython/petsc_discretisation.pyx	Introduces petsc_dm_insert_boundary_values() as a Python-level wrapper around DMPlexInsertBoundaryValues.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[lmoresi]

Adversarial review — PR #500 (rotated copy-back: insert essential boundary values)

Verdict: DEFECTS — the diagnosis and the regression test are right, but the chosen
insertion site (per-field sub-DM) segfaults on a datum that references a MeshVariable,
and the repo already owns a correct full-DM implementation of exactly this operation.

What we verified

• Whole file green on the PR head: tests/test_1018_rotated_freeslip.py 19/19 passed
  (73.9 s), including both enclosed-annulus zero-leakage tests (pressure nullspace /
  gauge handling — no regression from the insertion; pressure sub-DM has no essential
  BCs, so the insertion there is a proven no-op).
• Negative control fires: with the petsc_dm_insert_boundary_values call disabled in
  the installed copy, the new test fails with the exact zero-boundary signature
  (wall u_x off by 0.5 = the full datum). Restored, passes again.
• np2 (rank-safe probe, same setup): max wall error 0.0 across ranks; each rank
  owned ZERO nodes of one wall (rank 0: top=17/bot=0; rank 1: top=0/bot=17) and the
  insertion still completed — no empty-rank or partition-dependent behaviour in the fix.
• Wrapper hygiene: CHKERRQ used, matching the file's pattern; a field with no
  essential BCs is a graceful no-op (observed on the pressure sub-DM).
• _finalize_rotated_solution is called from the ONE exit of solve_rotated_freeslip
  (rotated_bc.py:891), after the loop — converged, unconverged-with-warning, and
  warm-start-already-converged paths all reach it. No skipped exit.

CRITICAL — segfault when the Dirichlet datum references a MeshVariable

The insertion runs on the per-field sub-DM from createFieldDecomposition
(solver._subdict[name][1]). An essential-field BC kernel that reads another
variable dereferences the auxiliary vector (petsc_a) — and the aux vec lives on
solver.dm (setAuxiliaryVec in the solve preamble), not on the sub-DMs.

Repro (scratchpad/aux_probe.py): the regression-test setup with the lid datum changed
from y - 0.5 to w.sym[0] where w is a MeshVariable holding y - 0.5.
Result: SIGSEGV (exit 139) inside the velocity-field insertion, at copy-back,
AFTER a successful solve. Instrumentation shows pressure inserts fine, velocity dies.
Control: with the insertion disabled the same run completes (fields show the original
err = 0.5 gap — the probe is valid). This is not an exotic pattern — a BC datum taken
from a field is how free-surface, remeshing-transfer and fault-datum setups write BCs.

The repo already owns the correct tool. SNES_Stokes_SaddlePt._scatter_global_to_fields
(petsc_generic_snes_solvers.pyx:8805) does this exact job on the FULL dm — where the
aux vec lives — via DMPlexSNESComputeBoundaryFEM, then splits per-field, with the
identical cache-invalidation tail. We ran it on the same rotated solution
(solver._scatter_global_to_fields(info["U"])): wall error 8.9e-16 / 2.2e-16, no
crash. The fix for #497 is one line in _finalize_rotated_solution — replace the
per-field scatter loop + new shim with solver._scatter_global_to_fields(U) — and the
new .pyx wrapper becomes unnecessary. If the shim is kept for other callers it must
document that a sub-DM target is unsafe with field-referencing BCs.

minor

• The time argument is dead as merged. The wrapper takes time=0.0; the caller
  never passes it. Measured context: a mesh.t-dependent lid on the rotated path NaNs
  the solve itself before finalize ever runs (rel |F̂| = nan, zero fields), and the
  NATIVE path with the same BC hits DIVERGED_FUNCTION_NANORINF with the copy-back
  writing ±1e307 garbage into the field (the PETSc PETSC_MIN_REAL insertion — issue
  Time-dependent Dirichlet BCs are broken in solve(): essential values inserted at PETSC_MIN_REAL, mesh.t kernels see -1.8e308 #410; DMPlexSNESComputeBoundaryFEM hardcodes it, dmplexsnes.c:426). So the PR is
  neutral w.r.t. Time-dependent Dirichlet BCs are broken in solve(): essential values inserted at PETSC_MIN_REAL, mesh.t kernels see -1.8e308 #410 — time-dependent essential BCs are broken on both paths
  before this change, and the user-facing SNES_Stokes.solve() doesn't even accept
  time=. But whichever insertion lands must be threaded with the solve time when
  Time-dependent Dirichlet BCs are broken in solve(): essential values inserted at PETSC_MIN_REAL, mesh.t kernels see -1.8e308 #410 is fixed; wants a TODO(BUG) → #410 at the call site.
• The regression test never runs in CI. scripts/test.sh serial batches match
  test_100[0-9]*, test_1010/1011/1050, test_11* — test_1018_* is in no batch
  (pre-existing gap for the whole rotated suite, but this PR's only automated guard is
  currently inert in CI).
• The new test is not parallel-safe. assert mask.sum() > 0 runs on every rank;
  at np2 rank 0 owns no Top nodes and the test fails on the test's own scaffolding,
  not the fix. Moot while 1018 is serial-only, but gather the counts before asserting
  if it ever moves.

nit

• from underworld3.cython.petsc_discretisation import ... inside
  _finalize_rotated_solution runs per solve; module level is fine (no cycle at that
  layer). Irrelevant if the CRITICAL fix removes the call.

Gate

pytest -m "level_1 and tier_a" -q -p no:cacheprovider --ignore=tests/test_0050_utils.py:
555 passed, 17 skipped (need --with-mpi), 1 xfailed, 0 failed in 601.8 s.

Gate infrastructure note (not a PR finding): without the ignore, the gate dies
silently — tests/test_0050_utils.py runs TWO Stokes solves at module scope, so
pytest COLLECTION executes them; we sampled the "collect-only" process 20+ min deep
inside SNESSolve → PCApply_FieldSplit_Schur at 99.7% CPU. scripts/test.sh already
comments this file out ("disable auditor test for now") but any whole-directory
pytest -m ... imports it anyway. Worth its own issue.

[lmoresi]

Response commit abae746 — the CRITICAL is fixed by switching sides, not by patching the shim:

• _finalize_rotated_solution now delegates to _scatter_global_to_fields (full-DM boundary FEM + per-field split + the identical cache tail). The insertion runs where the auxiliary vector lives, so a MeshVariable-referencing essential datum works: the review's segfault repro now completes with wall error 8.9e-16 / 2.2e-16.
• The new cython shim is withdrawn (no callers remain; a sub-DM target is unsafe with field-referencing BCs, so it should not be offered).
• TODO(BUG) → #410 recorded at the call site: DMPlexSNESComputeBoundaryFEM inserts at time=PETSC_MIN_REAL, same defect as the native copy-back — this PR is neutral there.
• Re-verified: test_1018_rotated_freeslip.py 19/19; np2 probe max wall error 0.0 including ranks owning zero wall nodes.

The two review minors that are repo-infrastructure rather than this PR (test_1018 is in no CI batch; test_0050_utils.py solves Stokes at module scope and silently kills whole-directory pytest collection) get their own issues.

Underworld development team with AI support from Claude Code