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-dimensionalWith 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:
- 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.
- 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).
Summary
SNES_Stokes.set_pressure_gauge(boundary)crashes when auw.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_callbackfired bySNESSetUpdate), wherep.datais non-dimensional. But it computes the surface mean withBdIntegral.evaluate(), which returns a dimensional PintQuantity:With units active this raises:
Stripping the Pint wrapper (
.magnitude) avoids the crash but is still wrong: the magnitude is in dimensional (Pa) units whilep.datais 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 itsuw.Modelscaling.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 ofp.data:wprecomputed once → the gauge is identical every iteration (no flapping) and stays in-cordon. Two requirements:w_i = ∫_Γ φ_i dS) to match the documented(1/|Γ|)∫_Γ p dSand keeptest_1016_snes_update_callbacks::test_pressure_gauge_zero_mean_on_boundaryexact (a uniform nodal mean fails it on a non-uniform boundary). For P1 on straight facets the lumped edge-length weights are exact._boundary_dof_maskusing 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 theset_pressure_gaugedocstring.Found via Claude Code while testing the gauge on the Spiegelman benchmark (yield-homotopy / FMG work).