boundary_flux() returns ~1e12 on Stokes_Constrained

[lmoresi]

solver.boundary_flux(boundary, normal=...) — the consistent-boundary-flux
back-calculation — reads correctly on an ordinary Stokes solve and returns
values of order 1e12 on Stokes_Constrained.

Measured

SolCx, StructuredQuadBox(elementRes=(32, 32)), free slip on all four walls,
normal = [[0, 1]] on Top, compared against SolCx.topography_top (mean
removed, corners excluded):

solver	Top held by	max abs(CBF)	vs exact
Stokes	add_dirichlet_bc((None, 0.0), "Top")	0.381	rel 0.082, corr +0.997
Stokes	add_rotated_freeslip_bc(0.0, "Top")	0.497	rel 0.087, corr +0.996
Stokes_Constrained	add_constraint_bc(0.0, "Top")	~1e12	uncorrelated (corr 0.03)

Exact peak is 0.379. The constrained solve itself is fine — velocity error
8.8e-06 at eta_B = 1e6, and its recovered n·σ·n matches the exact to 0.075 —
so this is the flux primitive on that solver, not the solve.

Likely _assemble_volume_reaction() picking up the multiplier field's rows or
the augmented-Lagrangian boundary stiffness (r = 1e4·μ, which the SolCx step
makes 1e10 on the stiff half); the magnitude is consistent with the latter.

Why it matters

It is the natural cross-check for #607: at convergence the momentum row's
boundary term satisfies

M_Γ (h + r(n·u − g)) = −(A·u − b)|_Γ

so h + r(n·u − g) IS the CBF traction de-smeared by the boundary mass. Being
able to read both off one solve is how that identity gets verified, and at the
moment it cannot be. Across two separate solves the corrected multiplier and the
rotated CBF reaction agree to 3–5%, which is inside each route's own error.

Reproduction: examples/solcx.py in UWTN 2026-016.

Underworld development team with AI support from Claude Code

[lmoresi]

Correcting my own report: the value is ~0, not ~1e12. The 1e12 I quoted was a
relative measure whose denominator was the near-zero CBF itself. Measured
directly, SolCx eta_B=10, 16×16, Top held by add_constraint_bc:

solver	CBF max abs	CBF rms	exact topography max
Stokes + component Dirichlet lid	3.70e-01	2.49e-01	3.68e-01
Stokes_Constrained	8.9e-02 (one corner node)	3.9e-13	3.68e-01

And it is not a defect. A multiplier-constrained boundary has no reaction left to
read: the constraint enters the SAME row it constrains, so at convergence the
assembled residual there is balanced and the back-calculation returns zero. That
is the identity from #607 seen from the other side —
M_Γ(h + r(n·u − g)) IS the CBF nodal load, so once it is subtracted there is
nothing in the residual.

Re-scoped to the real problem: boundary_flux returning a quiet ~0 is a trap. PR
#617 now raises on a multiplier-constrained boundary and points at
traction() / topography().

Apologies for the noise in the original report — the ratio was doing the lying.

Underworld development team with AI support from Claude Code