set_pressure_gauge crashes under units (computes the gauge in user space, not inside the SNES non-dim cordon)

[lmoresi]

Summary

SNES_Stokes.set_pressure_gauge(boundary) crashes when a uw.Model / units scaling is active, and is architecturally mis-placed across the non-dimensionalisation boundary.

Bug

The gauge callback runs inside the SNES non-dimensionalisation cordon (it is an add_update_callback fired by SNESSetUpdate), where p.data is non-dimensional. But it computes the surface mean with BdIntegral.evaluate(), which returns a dimensional Pint Quantity:

def _pressure_gauge(solver, iteration):
    mean = p_surface_integral.evaluate() / area   # Pint Quantity (dimensional)
    p.data[...] -= (mean - reference)             # p.data is non-dimensional

With units active this raises:

numpy.core._exceptions._UFuncOutputCastingError: Cannot cast ufunc 'subtract'
output from dtype('O') to dtype('float64') with casting rule 'same_kind'

Stripping the Pint wrapper (.magnitude) avoids the crash but is still wrong: the magnitude is in dimensional (Pa) units while p.data is non-dimensional, so the correction is off by the pressure scaling factor. Reproduced with the Spiegelman benchmark (docs/examples/WIP/Benchmark/Ex_VP_Spiegelman_Benchmark.py) under its uw.Model scaling.

Fix (design)

The gauge is a solver-internal operation and must be computed in non-dimensional space — no evaluate() round-trip. Make it a fixed linear functional of p.data:

mean = (w · p.data) / Σw            # w_i = ∫_Γ φ_i dS  (boundary-mass weights, built once)
p.data[...] -= (mean - reference)

w precomputed once → the gauge is identical every iteration (no flapping) and stays in-cordon. Two requirements:

1. Area-weighting (w_i = ∫_Γ φ_i dS) to match the documented (1/|Γ|)∫_Γ p dS and keep test_1016_snes_update_callbacks::test_pressure_gauge_zero_mean_on_boundary exact (a uniform nodal mean fails it on a non-uniform boundary). For P1 on straight facets the lumped edge-length weights are exact.
2. Discontinuous pressure: the boundary label carries no DOFs for a DG pressure, so the DOF set must come from the facet support (adjacent cells), not just its closure. (A working _boundary_dof_mask using closure-then-support was prototyped — selects strict top vertices for P1, top-cell DOFs for DG0/DG1.)

Related finding (not a bug)

For pressure-dependent plasticity (Drucker–Prager, σ_y = C + sinφ·p) a pressure gauge is counter-productive: gauging the pressure shifts the yield surface every iteration, so the solve oscillates (~50 iters / DIVERGED_MAX_IT vs ~5 clean without it). The pressure there is physical, not a free gauge — worth a note in the set_pressure_gauge docstring.

---

Found via Claude Code while testing the gauge on the Spiegelman benchmark (yield-homotopy / FMG work).