Return the whole boundary traction from the constrained solver (#607) - #617
Return the whole boundary traction from the constrained solver (#607)#617lmoresi wants to merge 2 commits into
Conversation
The momentum row carries (h + r(n.u - g)).n, so the traction holding the boundary is that sum. multiplier() returns h, and topography() was built on h alone, which is short by the augmented-Lagrangian share: r times the discrete constraint residual. With the viscosity-weighted default r = 1e4.mu that share is a few per cent of the surface traction on a uniform-viscosity annulus, and across SolCx's 1e6 viscosity step it is most of it -- h alone reads a tenth of the exact topography and is ANTI-CORRELATED with it. - traction(boundary) returns h + r(n.u - g), the quantity the CBF back-calculation recovers: at convergence M_Gamma(h + r(n.u - g)) balances the volume residual restricted to the boundary, which is the CBF nodal load. - topography() is built on it, so the public dynamic-topography path is correct without a change at the call site. - multiplier() still returns h and now says what it is not. Why this survived validation: the existing check scored corr(lambda, -n.sigma.n) = 0.9999. A correlation is scale-free and cannot see a systematic amplitude deficit, which is exactly what a missing share of the load is. The new guard scores a relative l2 against the exact SolCx surface topography (0.047 with the share, 1.04 without) and carries the bare multiplier as its negative control, so a regression that quietly reverts to h cannot pass. test_1061's API test asserted topography == h/scale, which pinned the defect. It now asserts the outcome: topography is the traction over the scale, and the traction is not the multiplier. Docs: the design note carries a dated correction and the corrected CBF identity; the user-facing curved-boundary page now writes the penalty recipe against mesh.boundary_normal rather than mesh.Gamma, with the measurement that a facet-normal penalty does not converge (velocity error 0.60 flat, stress error 0.21 -> 0.26 under refinement, while the leak reads 1e-5); CHANGELOG entry added. Underworld development team with AI support from Claude Code
…614) build_rotation skipped any boundary node with a velocity DOF constrained out of the global vector -- which is every node where a rotated wall meets a wall held by an essential condition. The node was dropped from the rotation entirely, so the WALL-NORMAL COMPONENT WAS NEVER CONSTRAINED THERE. The wall leaked at its own end points while every interior node was exact. Measured, unit box, rotated lid, component free slip on the other three walls, uniform viscosity: max|u_y| on the lid was 4.0e-3 against |u|max 2.5e-2 -- 16%, entirely at the two corners -- and the solve differed from the equivalent component-Dirichlet lid by 2.0e-3 globally, with an EXACT LINEAR SOLVE on both sides, so it was never a convergence artefact. After the fix the lid holds u.n to 0.0 and the two solves agree with a direct reference to 8.9e-9. The fix: a partially constrained node keeps the constraint on what is left. With the pinned components at zero, n.v = 0 reduces to n_F.v_F = 0 on the free subspace, so the frame is built there and its normal rows constrained. A node whose normal lies entirely in the pinned subspace is already implied and is still skipped. A node carrying a prescribed v_n datum keeps the old behaviour and says so once -- reducing the affine constraint against the pinned values is not implemented. This also removes the mechanism behind #608: the corner reaction was the essential constraint's because the rotated constraint was not there. boundary_flux() now raises on a multiplier-constrained boundary instead of returning a quiet ~0 (#614). There is no reaction to read there: the constraint enters the same row it constrains, so the assembled residual is balanced at convergence (measured rms 3.9e-13 against a traction of 0.37). The traction is the multiplier, and traction() returns it. My original report of that issue said 1e12; that was a ratio with a near-zero denominator, and the issue is corrected. Guards: tests/test_1066_rotated_meets_essential.py asserts the lid holds u.n to round-off and that both lids reproduce a monolithic direct solve. Existing suites unaffected: 42 rotated/nitsche/constrained tests pass. Underworld development team with AI support from Claude Code
|
Extended to cover the other two paths, so the fix lands as one consistent set The corner was never constrained (#616, and the mechanism behind #608)
Unit box, rotated lid, component free slip on the other three walls, uniform
The second row with an exact linear solve on both sides is the point: it was The fix: a partially constrained node keeps the constraint on what is left. With Guard: boundary_flux on a constrained boundary (#614) — and a correctionI reported that it returns ~1e12. It returns ~0 (rms 3.9e-13 against a It now raises and points at Issue #614 is corrected on the issue itself. |
Closes #607.
Stokes_Constrainedassembles the momentum row's boundary term as(h + r(n·u − g))·n, so the traction holding the boundary is that sum.multiplier()returnsh, andtopography()was built onhalone — short bythe augmented-Lagrangian share, which is
rtimes the discrete constraintresidual.
With the viscosity-weighted default
r = 1e4·μ(x)that share is a few per centof the surface traction on a uniform-viscosity annulus. Across SolCx's
1e6viscosity step, where
rreaches1e10on the stiff half, it is most of it:multiplier()as returnedtraction()(this PR)SolCx.topography_topn·σ·n, same solveThe corrected read is the most accurate surface stress available on that problem
— better than the projection — and the solve itself was never at fault
(velocity error 8.8e-06,
max abs(n·u)1.7e-10).What changed
traction(boundary)— new, returnsh + r(n·u − g)symbolically.topography()— built on it, so the public dynamic-topography path iscorrect with no change at the call site.
multiplier()— unchanged behaviour (it returnsh) and now documents whatit is not.
rnow say the CONSTRAINTand the
traction()read are independent of it, andhalone is not.Why it survived validation
The existing check scored
corr(λ, −n·σ·n) = 0.9999. A correlation is scale-freeand cannot see a systematic amplitude deficit, which is precisely what a missing
share of the load is.
tests/test_1063_constrained_traction.pyscores a relativel2against the exact SolCx topography and carries the bare multiplier as itsnegative control (it must read > 0.5), so a regression that quietly reverts to
hcannot pass.test_1061's API test assertedtopography == h/scale, whichpinned the defect; it now asserts the outcome.
The identity behind it
At convergence
M_Γ(h + r(n·u − g))balances the volume residual restricted tothe boundary, which is the CBF nodal load (Zhong, Gurnis & Hulbert 1993). So the
multiplier route and
boundary_normal_tractionon a rotated constraint are thesame computation; measured across two solves they agree to 3–5%, inside each
route's own error. This is also why the free-surface work was right to reject the
multiplier as returned and keep the rotated lid + CBF.
Documentation
CONSTRAINED_FREESLIP_MULTIPLIER.md: dated correction, corrected CBF identity.docs/advanced/curved-boundary-conditions.md: the penalty free-slip recipe nowuses
mesh.boundary_normalrather thanmesh.Gamma, with the measurement thata facet-normal penalty does not converge — velocity error 0.60 flat and stress
error 0.21 → 0.26 under refinement, while the leak reads 1e-5 throughout. The
projected-normal recipe is marked superseded.
CHANGELOG.mdentry.Tests
tests/test_1061,1062,1063,1064,1024: 26 passed, 2 xfailed.Related, filed separately and NOT addressed here: #608 (the reaction reported at
a corner shared with an essential BC), #614 (
boundary_fluxreturns ~1e12 onStokes_Constrained), #616 (rotated free slip does not reproduce the equivalentcomponent condition on a flat wall, 2e-3, with an exact linear solve on both
sides).
Underworld development team with AI support from Claude Code