From 84638d4607c71b4361e9601dbc68eb938e619288 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 26 Aug 2026 22:00:54 +1000 Subject: [PATCH] Refuse CBF recovery on a multiplier-constrained boundary (#614) The consistent-boundary-flux back-calculation reads the VELOCITY residual. On a boundary held by `add_constraint_bc` the traction has been moved into the multiplier block, so the velocity rows retain only a CONSTANT traction, and de-smearing hands that constant straight back. That is not an inaccuracy, it is a different quantity. Measured on SolCx, the raw reaction along the boundary took exactly two magnitudes in the ratio 4 : 2 : 1 at midpoint / interior vertex / end vertex -- the P2 line-trace lumped mass, i.e. R = c*m_i to the digit -- so the recovered flux had a standard deviation of 2e-15 across 63 nodes where the equivalent `Stokes` solve varied correctly and correlated -0.997 with the exact topography. The issue reported this as ~1e12, which is no longer what happens; today it returns a plausible O(0.1) number that happens to carry no spatial information at all. That is worse, not better: 1e12 announces itself and 0.11 does not, and a caller taking a mean or a peak off it gets something entirely reasonable looking. Nothing here is corrupted -- the reaction assembly and the de-smear both do exactly what they are told. The primitive reads a place that, on this solver, no longer holds the quantity. So this refuses instead, naming multiplier() and topography(), which is where the traction actually is. The guard is on the BOUNDARY, not the solver: the regression checks that the Dirichlet boundary of the same constrained solve still recovers, and that what it recovers has spatial content rather than being constant too. This is the cheap half. The full fix is for the traction accessor to DISPATCH on how each boundary is held -- rotated reaction, multiplier, or CBF -- so a caller does not have to know. The solver already has what it needs to decide (`_block_constraint_bcs`, `_rotated_freeslip_info`); what is missing is a common return shape, since multiplier() returns a field where the others return nodal arrays. That belongs with #607/#617, which is already making the constrained route return the whole traction. Underworld development team with AI support from Claude Code --- src/underworld3/utilities/boundary_flux.py | 28 ++++++++++++++ tests/test_1019_boundary_flux.py | 43 ++++++++++++++++++++++ 2 files changed, 71 insertions(+) diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index c0b8e5b2..9ac32c70 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -633,10 +633,38 @@ def vcoord(q): return cvec[csec.getOffset(q) // dim] return np.array([sig[gi[_key(x, dim)]] for x in xs]) +def _refuse_a_constrained_boundary(solver, boundary): + """CBF recovery reads the VELOCITY residual. On a boundary held by + ``add_constraint_bc`` the traction has been moved into the multiplier block, so + what remains in the velocity rows is a CONSTANT traction — the reaction comes back + exactly proportional to the nodal boundary mass, and de-smearing hands that constant + straight back (#614). + + That is not an approximation, it is the wrong quantity: measured on SolCx the + recovered flux had a standard deviation of 2e-15 across the boundary while the + equivalent ``Stokes`` solve varied correctly and correlated -0.997 with the exact + topography. It returns a plausible-looking number, which is why it went unnoticed. + + Refuse rather than mislead. The traction on such a boundary is the multiplier + (plus its augmented-Lagrangian share — #607).""" + constraints = getattr(solver, "_block_constraint_bcs", None) + if not constraints: + return + if any(getattr(c, "boundary", None) == boundary for c in constraints): + raise NotImplementedError( + f"boundary_flux cannot recover the traction on {boundary!r}: it is held " + "by add_constraint_bc, so the traction is carried by the multiplier and " + "the velocity residual this reads holds only a constant (#614). Use " + f"solver.multiplier({boundary!r}) or solver.topography({boundary!r}) " + "instead." + ) + + def boundary_flux(solver, boundary, mass="auto", remove_mean=False, normal=None): """See ``SolverBaseClass.boundary_flux``. Returns ``(xs, flux)`` for this rank's boundary nodes; scalar solver → normal flux, vector solver → traction (or its normal component if ``normal`` is given).""" + _refuse_a_constrained_boundary(solver, boundary) dm = solver.dm; dim = solver.mesh.dim ra = np.asarray(solver._assemble_volume_reaction()).ravel() nodes, lsec, csec, cvec, v0, v1, edge_nodes = _boundary_field_nodes( diff --git a/tests/test_1019_boundary_flux.py b/tests/test_1019_boundary_flux.py index 86e67faa..805f8675 100644 --- a/tests/test_1019_boundary_flux.py +++ b/tests/test_1019_boundary_flux.py @@ -11,6 +11,7 @@ import sympy import pytest import underworld3 as uw +from underworld3 import analytic as A pytestmark = [pytest.mark.level_1, pytest.mark.tier_a] @@ -316,3 +317,45 @@ def test_volume_residual_fields_insert_essential_values(): f"{name}: compute_volume_residual_fields diverges from " f"_assemble_volume_reaction on g != 0 walls — essential values not inserted" ) + + +@pytest.mark.level_2 +def test_cbf_refuses_a_multiplier_constrained_boundary(): + """CBF must not silently return a constant on a constrained boundary (#614). + + The back-calculation reads the VELOCITY residual. Where `add_constraint_bc` + holds the boundary, the traction has been moved into the multiplier block and + the velocity rows retain only a constant one — the reaction comes back exactly + proportional to the nodal boundary mass (ratio 4:2:1 at midpoint / interior + vertex / end vertex, the P2 line lumping), so de-smearing returns that constant. + Measured std across the boundary was 2e-15 where the equivalent `Stokes` solve + varied correctly. The danger is that the value LOOKS reasonable. + + The unconstrained boundary of the same solver must keep working, or this is a + guard on the solver rather than on the boundary. + """ + mesh = uw.meshing.StructuredQuadBox( + elementRes=(16, 16), minCoords=(0, 0), maxCoords=(1, 1), qdegree=3) + sol = A.SolCx(mesh, eta_A=1.0, eta_B=1.0e3, x_c=0.5, n=1) + solver = uw.systems.Stokes_Constrained(mesh) + solver.constitutive_model = uw.constitutive_models.ViscousFlowModel + solver.constitutive_model.Parameters.shear_viscosity_0 = sol.fn_viscosity + solver.bodyforce = sol.fn_bodyforce + solver.tolerance = 1.0e-9 + solver.petsc_use_pressure_nullspace = True + solver.add_constraint_bc(0.0, "Top") + solver.add_dirichlet_bc((None, 0.0), "Bottom") + solver.add_dirichlet_bc((0.0, None), "Left") + solver.add_dirichlet_bc((0.0, None), "Right") + solver.solve() + + with pytest.raises(NotImplementedError, match="add_constraint_bc"): + solver.boundary_flux("Top", normal=[[0.0, 1.0]]) + + # negative control: the Dirichlet boundary of the SAME solver still recovers, + # and recovers something with spatial content rather than a constant. + xs, flux = solver.boundary_flux("Bottom", normal=[[0.0, -1.0]]) + assert len(xs) > 0 + assert np.asarray(flux).std() > 1.0e-3, ( + "the unconstrained boundary came back constant too — the guard is on the " + "solver, not on the boundary")