Refuse CBF recovery on a multiplier-constrained boundary (#614)

[lmoresi]

Stops the silent wrong answer in #614. Root-caused first — see the analysis on the issue.

What is actually wrong

CBF 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. The raw reaction along the boundary takes exactly two magnitudes, in ratio

2.2971e-03 : 1.1486e-03 : 5.7428e-04   =   4 : 2 : 1

which is the P2 line-trace lumped mass at midpoint / interior vertex / end vertex. So R = c·m_i to the digit, and de-smearing divides m_i back out and returns c. Recovered std across 63 nodes: 2e-15, against 0.256 and correlation −0.997 for the equivalent Stokes solve.

Nothing is corrupted. The reaction assembly and the de-smear both do exactly what they are told; the primitive reads a place that no longer holds the quantity on this solver.

Why this is more urgent than the issue implies

#614 reported ~1e12. That is no longer what happens — it now returns a plausible O(0.1) number carrying no spatial information. 1e12 announces itself; 0.11 does not. Anyone taking a mean or a peak off it gets something entirely reasonable-looking.

What this does

Refuses, naming multiplier() / topography() — where the traction actually is.

The guard is on the boundary, not the solver. The regression asserts that the Dirichlet boundary of the same constrained solve still recovers, and that what it recovers has spatial content rather than being constant too — otherwise this would be a guard on the wrong thing.

Also retracted on the issue: the augmented-Lagrangian hypothesis (r = 1e4·μ) rested on the 1e12 magnitude, and there is no r-scaling left — the peak moves only 0.110 → 0.143 across three decades of viscosity contrast.

The half this deliberately does not do

The real fix is for the traction accessor to dispatch on how each boundary is held — rotated reaction, multiplier, or CBF — so callers do not have to know which primitive matches their boundary condition. Today there are three entry points with three different contracts: boundary_normal_traction raises unless the boundary is rotated, boundary_flux silently misread a constrained one (until this PR), and multiplier()/topography() exist only on the constrained class.

The solver already has everything needed to decide (_block_constraint_bcs, _rotated_freeslip_info). What is missing is a common return shape: multiplier() returns a field where the others return nodal arrays. That is exactly what #617 is building, so the dispatch belongs there rather than here.

Tests: test_1019 15 passed; test_1061 + test_1018 39 passed.

Underworld development team with AI support from Claude Code

[Copilot]

Copilot was unable to review this pull request because the user who requested the review has reached their quota limit.