Draft UWTN 2026-016: Boundary conditions on non-planar boundaries - #26
Open
lmoresi wants to merge 13 commits into
Open
Draft UWTN 2026-016: Boundary conditions on non-planar boundaries#26lmoresi wants to merge 13 commits into
lmoresi wants to merge 13 commits into
Conversation
R1 in the writing plan. 'No flow through this wall' is a velocity component on a box and is not a component of anything on a sphere, a deformed mesh, or a surface with topography, and the note is about what to do instead. Three approaches in the order in which each answers the previous one's weakness: a direct penalty, which is consistent only in the limit so closing the leak means conditioning the operator worse; Nitsche, which carries the boundary traction terms and is consistent above a stability threshold; and rotating the degrees of freedom, which imposes the constraint exactly. The weak forms are written out, because the differences between them are the argument and 'add two more terms and it becomes consistent' is not believable unseen. The rotated approach is presented as the classical answer it is -- Engelman, Sani and Gresho were reviewing the alternatives in 1982 -- so the question the note actually answers is why the exact method is the least used of the three. The answer is structural: a per-node rotation leaves the problem in a mixed basis, and the figure shows where that obligation falls. It is contained. The velocity solve is rotated and carries its multigrid; the Schur complement and pressure solve beside it never handle a rotated vector. The figure moves here from the underworld3 repository, where it was parked because this article did not exist yet. Its typst source, generator and JSON geometry ship in examples/. Two things are deliberately visible rather than quietly resolved. The section on when the choice matters carries a warning that it is not yet measured -- the argument is structural and we have not run the experiment. And the note asks whether Behr (2004) reaches the same measure-weighted normal we did, because if it does, our #560 was a rediscovery and should say so. Carries the allocator fix from the exact-solutions branch. Without it new_article.py issued 2026-012, which is already gthyagi's -- the same collision the fix exists to prevent, because the fix is not yet merged. Underworld development team with AI support from Claude Code
Contributor
|
Preview Or the whole site. Built from |
The note asserted that a weakly imposed constraint holds to the accuracy of the discretisation and a strong one to the accuracy of the arithmetic. Annulus, no slip inside, the treatment under test outside, degree-four radial forcing; largest normal velocity on the outer boundary against the true radial direction, over the flow speed: cell Nitsche rotated 0.150 4.6e-03 7.3e-11 0.100 2.2e-03 6.5e-11 0.075 1.7e-03 1.0e-10 0.050 5.8e-04 8.4e-11 Refinement is what separates them and one resolution would not have. Nitsche falls as about h^1.9 -- discretisation-limited, which is what consistency buys and all it buys. The rotated constraint does not move with the mesh at all. Three things the measurement cost, all worth recording: The first control was worthless. Two free-slip circles leave the rigid rotation unconstrained, and that nullspace is purely TANGENTIAL, so a radial leak metric read 2e-14 on a solve that had diverged with |u| = 2.7e5. The inner boundary is now no-slip and the control -- outer boundary left natural -- reads 0.98. A vector MeshVariable's .array is (N, 1, dim). It broadcasts silently against an (N, dim) array of normals and returns projections of 1e-16 for a velocity of 1e-2, without raising. Every number in the first run was that artefact. Underworld exposes no separate direct-penalty condition and building one from a natural BC would not converge, so that row has no measured counterpart. Said in the note rather than quietly dropped, alongside the other limit: this measures the constraint, not the traction recovered from it. Underworld development team with AI support from Claude Code
The previous commit said Underworld exposes no direct-penalty condition. It
does, it is documented in docs/advanced/curved-boundary-conditions.md, and it
is in a good many working scripts:
G = mesh.Gamma
stokes.add_natural_bc(1.0e4 * G.dot(v.sym) * G, "Upper")
What I had written used the nodal boundary_normal and a NEGATIVE coefficient.
A negative one is anti-damping, so the linear solve failed on the first
iteration, and I read that as the method being unavailable rather than as my
sign being wrong.
With all three measured, the resolution table gains its first column, and the
two weak methods look alike in it -- both leak parts in a thousand and both
improve with the mesh. What separates them is their own parameter:
penalty 1e2 2.6e-01 | 1e3 2.3e-02 | 1e4 8.5e-04 | 1e5 9.6e-04 | 1e6 diverged
gamma 1 diverged | 10 1.7e-03 | 100 2.7e-04 | 1000 3.0e-05
The penalty improves in proportion to how hard it pushes until the conditioning
catches up: no gain from 1e4 to 1e5, failure at 1e6. Nitsche fails below its
stability threshold and improves steadily above it, with no wall in this range,
which is why gamma=10 can be a documented default instead of a number to fit
per problem.
Two runs in the first parameter sweep had failed the line search and were about
to be tabulated as measurements. The sweep now records the converged reason and
prints 'diverged' instead of a number.
Underworld development team with AI support from Claude Code
The note said Nitsche 'keeps improving without the conditioning wall'. That was
written from a sweep that stopped at gamma = 1000, which is one value short of
where it fails.
gamma 1e3 3.0e-05 converged
1e4 -- DIVERGED_LINE_SEARCH
1e5 -- DIVERGED_LINE_SEARCH
It does not stop converging in the metric, it stops converging: reason -6, and
|u| collapses toward zero, so the leak figures above the threshold are from
solves that never finished.
So both weak methods are bounded above by conditioning and neither escapes
tuning. What differs is the floor each reaches before failing -- the penalty
bottoms out near 1e-3, Nitsche at 3e-5 -- which is what consistency actually
buys: a lower floor, not an unbounded parameter.
Nitsche is bounded below as well, at gamma = 1, where the form stops being
coercive. Its usable range has a threshold at each end, and gamma = 10 sits in
the middle of that window on any mesh because gamma is dimensionless and the
term it scales already carries mu/h. The penalty coefficient carries no such
scaling, which is why its working value is a property of the problem.
The shipped sweep now runs into the failing range rather than stopping just
below it.
Underworld development team with AI support from Claude Code
The note had three approaches and there are four. add_constraint_bc carries a scalar multiplier field h coupled into the saddle point, so the constraint is a ROW of the system rather than a term added to one -- and at convergence h on the boundary is the normal traction, which makes it the second of the four that returns the stress rather than recovering it. Measured alongside the others, same annulus, same forcing: cell penalty nitsche multiplier rotated 0.150 4.5e-03 4.6e-03 8.3e-04 7.3e-11 0.100 2.5e-03 2.2e-03 1.6e-04 6.5e-11 0.075 8.5e-04 1.7e-03 5.6e-05 1.0e-10 0.050 9.5e-04 5.8e-04 1.2e-05 8.4e-11 It is an approximation, as expected, but a much better behaved one: it starts an order of magnitude below the weak forms and falls near h^3.9 against their h^1.9, because the only thing approximated is the screening of the interior multiplier degrees of freedom rather than the enforcement itself. The note says where that screening enters, and that the augmented-Lagrangian term conditions the [p,h] block without biasing the multiplier, so accuracy does not depend on it the way a penalty's does. Also started the surface-traction comparison, which is the test that matters, and it is NOT in the note. Two exact routes to sigma_nn -- the multiplier and the rotated reaction -- agree to about 1% at one resolution, but under refinement the vertex agreement converges while the midpoint agreement gets steadily worse, which is backwards: the solver documents midpoints as superconvergent on a curved boundary. The fault is most likely this script selecting the multiplier's boundary trace by a radius band that narrows with the mesh; it needs to come from the boundary label instead. The symptom and the suspicion are recorded in examples/traction.py rather than the numbers being quoted anywhere. Underworld development team with AI support from Claude Code
The traction comparison in traction.py has no oracle -- it measures two methods against each other, so its diverging midpoint column cannot be interpreted. uw.analytic.Zhong2008 fixes that: .response() returns the surface topography kernel directly, 0.4191904 for the default degree-2 case, referenced to Zhong et al. (2008) GGG 9 Q10017 after Hager and O'Connell. The quantity it publishes is the one the note cares about rather than a velocity field to infer a stress from. One consequence worth stating before anyone starts: Zhong is SPHERICAL. The harmonic degree is a spherical-harmonic degree and the solution carries a planet radius, so the matching model is a 3-D shell rather than the 2-D annulus the leak tables use. That sets the cost and is why this is separate work rather than another column. Also records why SolCx stays useful even though it cannot test a curved boundary: on a box all four treatments reduce to the same component constraint, but it is still a fair test of whether a penalty or Nitsche form constrains the normal degrees of freedom it is handed, with the geometry factored out and an exact answer to hand. Underworld development team with AI support from Claude Code
Checked what the analytic suites actually vary. assess's cylindrical and spherical solutions (Kramer et al. 2021) take nu as a SCALAR -- isoviscous. Zhong2008 layers viscosity by radius through viscosity_interfaces, also not laterally. SolCx has the lateral jump and a Cartesian box. So no analytic solution has both a curved boundary and a lateral viscosity variation, and the test splits in two rather than one being a weaker version of the other: Kramer or Zhong for the geometry with the rheology trivial, SolCx for the rheology with the geometry trivial. Lateral contrast is the condition under which these methods are known to give trouble and it is the only one of the two for which an oracle exists at all. Also noted that a shell with laterally varying viscosity has no exact solution, so that case needs a different kind of evidence -- convergence or agreement between methods -- and the note should say which rather than implying a measurement it cannot make. Underworld development team with AI support from Claude Code
The leak tables measure the constraint. This measures what the constraint is wanted for, against `uw.analytic.CylindricalStokes` (Kramer et al. 2021), whose `assess` package publishes the radial stress rather than leaving it to be recovered. On the outer boundary the exact answer is a single harmonic, 0.1506696 cos 2t, fitted to a residual of 1e-16, so the metric is one amplitude. Three findings, two of which correct what the note said. The direct penalty written against the FACET normal does not converge. Its leak falls to 1e-5 while the velocity is 60% wrong and the surface stress 26% wrong, and refining from cell 0.15 to 0.035 does not move either. The same coefficient against the measure-weighted node normal converges cleanly. This is the classic over-constraint the consistent normal was introduced to avoid, and it means the note's penalty column was measuring the onset of locking: at a coefficient of 1e3 that penalty leaks 3e-2 and gets the stress right to 2e-3, and at 1e8 it leaks 1e-7 and is 26% wrong. Both penalty columns are now tabulated. Once the constraint is imposed against the node normal, which treatment imposed it stops mattering to the recovered stress -- rotated, multiplier, Nitsche and penalty all read 6.3e-3 at cell 0.075. The recovery sets the floor, not the boundary condition, so the note's reasoning that a weakly imposed constraint poisons the traction was too strong. What the rotated constraint buys is the reaction, which is three times better on the same solve and differentiates nothing. The multiplier is not quite the whole traction: the momentum row carries h + r(u.n - g), and with the default r = 1e4.mu the missing share is a few per cent of the surface stress. underworld3#607. Also measured: the SolCx half, where Nitsche at gamma = 10 gets the surface stress 277% wrong at a viscosity contrast of 1e6 while holding the constraint to 1e-7, and where the rotated constraint disagrees with a component condition at the corner they share (underworld3#608). Behr (2004) is read and cited for what it says: the consistent normal is standard and preferred for conservation, and it does not on its own guarantee a good discrete slip condition. The derivation in the "Which normal" section is credited to Engelman, Sani and Gresho (1982) rather than presented as ours. `traction.py` is removed -- `stress.py` supersedes it, with an oracle and with the boundary trace selected by label rather than by a radius band. The mesh cache the examples write beside themselves is untracked and gitignored. Underworld development team with AI support from Claude Code
Contributor
Author
|
The stress test is in. It measures the surface normal stress against an exact The facet-normal penalty does not converge. Its leak falls to 1e-5 while the The recovery is the floor, not the boundary condition. At cell 0.075 the Two defects fell out, both filed against underworld3:
The SolCx half is in too. It cannot rank the treatments — the component Behr (2004) is now read and cited for what it says — the consistent normal is Also: a banner, Still open: the spherical case (Zhong 2008, Underworld development team with AI support from Claude Code |
Predicted against computed, along the top wall, at a viscosity contrast of 100 and of a million. The exact answer is published by `uw.analytic.SolCx`, so this comparison can be a figure rather than a column of numbers. Drawing it turned up a defect the tables had missed. At a contrast of 1e6 the multiplier field as `multiplier()` returns it peaks at 0.042 against an exact 0.383 and is anti-correlated with it (-0.53), on a solve whose velocity error is 8.8e-06. The default augmentation is r = 1e4.mu(x), which the viscosity step makes 1e10 on the stiff half, so almost the whole traction sits in the augmented-Lagrangian term rather than in h. Adding it back gives 0.047, the most accurate surface stress in the comparison -- better than the projected recovery at 0.075. Turning r down is not available: at this contrast r = 0 and r = 1e2 both fail to solve. Recorded on underworld3#607, which the annulus had already found at the few-per-cent level. The figure carries both multiplier curves, as returned and corrected, so the size of the missing term is visible rather than asserted. The script caches its curves beside itself: the solves take about ten minutes and the figure is redrawn more often than it is recomputed. A dash pattern is a linestyle and not a format string, and a sign convention is asserted rather than assumed -- the first draft of this figure stopped on that assertion, which is how the collapse was noticed. Underworld development team with AI support from Claude Code
… are masked Two corrections and one identity, all from Louis reading the draft against the free-surface work. The note said the multiplier's interior degrees of freedom are screened and that the screening is where the approximation enters. They are not screened, they are CONSTRAINED OUT of the global system in the PetscSection before the solve sees them (_constrain_interior_multipliers_in_section), which is documented as lossless: the interior rows are the screening block alone, so the interior multiplier is determined by the boundary trace and pinning it moves a converged solve by ~1e-8 in velocity. They neither cost anything nor limit the accuracy. The identity: the momentum row's boundary load is M_Gamma (h + r(u.n - g)), and at convergence that balances the volume residual restricted to the boundary, which is the nodal load the consistent boundary flux back-calculation reads (Zhong, Gurnis & Hulbert 1993). So the corrected multiplier is the CBF traction de-smeared by the boundary mass -- the same computation, arrived at by carrying the traction as an unknown rather than reading it out of the residual. Measured across two solves they agree to 3.2% and 4.9%, inside each route's own error. That is why the free-surface work was right on both counts: it rejected the multiplier AS RETURNED, which is missing the augmentation share, and kept the rotated lid with the CBF reaction, which is the same quantity with nothing missing. The corner spike is re-framed. It is the documented behaviour of the reaction recovery at a node where a rotated wall meets an essential one -- the value reported there is the essential constraint's reaction, not the wall's traction -- and the practice, which the free-surface SolCx work already followed, is to leave the corner out of the rotation. Control: the same CBF read on a lid held by the component condition alone peaks at 0.381 against an exact 0.379 with no spike, where the rotated run peaks at 0.497. underworld3#608 is re-scoped to match. New: underworld3#614, boundary_flux() returns ~1e12 on Stokes_Constrained while reading the exact topography to 8% on an ordinary solve. It is the cross-check that would verify the identity within a single solve, and it is unavailable. The Zhong citation is pinned in a references.bib: the DOI contains parentheses and MyST's inline DOI form truncates at the first one. An at-sign in a bibtex COMMENT fails the whole file. Underworld development team with AI support from Claude Code
…taminated table The three defects this note turned up are fixed in underworld3 PR #617, so every table here is re-measured against a solver that behaves. The SolCx half is the one that changes. The three exact treatments -- component Dirichlet, multiplier, rotated -- now agree to three figures at every viscosity contrast from 10 to 1e6, whole wall and trimmed alike. Before, the rotated column read 0.322 at a contrast of 10 against the reference's 0.048, all of it two corner nodes, and the multiplier was an order of magnitude out at 1e6. Getting those three to agree was the point. A bare penalty coefficient now has numbers instead of "diverged": 0.045 at a contrast of 10, and 0.992 at 1e6, which is to say none of the signal. Scaling it by the local viscosity still does not solve. WITHDRAWN: the earlier SolCx numbers, including the Nitsche column. They were taken while concurrent runs shared one mesh-cache file -- StructuredQuadBox keys its cache on the box corners and NOT on the element resolution, so runs at different resolutions silently swap meshes (underworld3#618, filed). A marginal solve then flips between converged and diverged for reasons that look like the method, which is exactly how it fooled me. Everything in that section was re-measured sequentially on a cleared cache. Nitsche has no column now: our configuration of it on this box converges at a contrast of 1e6 and fails the line search at 10, at both resolutions and at gamma = 10, 100 and 1000, and we are not confident enough in that setup to publish numbers from it. New: what each treatment COSTS. The solve is the same for all four within the spread; the recovery is 15 to 25 times cheaper for the reaction and the multiplier than for a projection, because they read the state the solve already returned. And a "which one to use" section, which the accuracy and the cost now point at together. The multiplier section is rewritten around the fix rather than the defect, with the clean SolCx numbers: h alone reads 0.042 against an exact 0.381 and is anti-correlated with it; h + r(u.n - g) reads 0.382. Underworld development team with AI support from Claude Code
The preview build failed on CI with "Citation data from doi.org was not available or malformed" for the Behr and Kramer DOIs, while resolving fine locally. That is the failure the MMPDE note's references.bib was written to prevent: a build that depends on doi.org answering will eventually publish a note with a broken citation, and a deposited PDF cannot be repaired afterwards. All five are now pinned and cited by key. Underworld development team with AI support from Claude Code
Louis's note: the interesting thing here is two penalty approaches against two ways of imposing the constraint exactly -- by construction, or by a Lagrange multiplier -- and that pairing is now stated where the four are introduced. Removed: the section that told the multiplier's missing augmentation share as a defect story, the withdrawn-measurement note, the "how it survived validation" paragraph, and the issue numbers. What survives is the method. The traction a multiplier constraint holds the boundary with is h + r(u.n - g), which is stated where the method is described, and the two exact routes are shown to be the same computation, which is the point worth keeping. The timings are re-measured at sizes where they mean something: 71 000 velocity nodes rather than 2 000, where the numbers were hundredths of a second and inside the run-to-run spread. At that size the multiplier's solve costs 50% more than the others, and the recovery separates by a factor of 25 to 40 -- the projection costs MORE than the Stokes solve (0.63 s against 0.45 s) while the two exact routes read their traction off the state the solve already returned (27 ms and 16 ms). The smaller meshes are kept only to show the scaling, and the table says so. Underworld development team with AI support from Claude Code
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
R1 in
WRITING-PLAN.md, rescoped from "rotated boundary conditions" — the question is the better frame and it is the one a reader arrives with. Atreview."No flow through this wall" is a velocity component on a box, and is not a component of anything on a sphere, a deformed mesh, or a surface with topography.
What it says
Three approaches, in the order in which each answers the previous one's weakness, with the weak forms written out:
The maths is there because the differences are the argument — "add two more terms and it becomes consistent" is not believable unseen.
The rotated approach is presented as the classical answer it is. Engelman, Sani & Gresho were reviewing the alternatives in 1982, so the question the note actually answers is why the exact method is the least used of the three. That answer is structural rather than numerical: a per-node rotation leaves the problem in a mixed basis, and everything downstream has to agree about which nodes are which.
The figure shows where that obligation falls, and the point is that it is contained — the velocity solve is rotated and carries its multigrid with it, while the Schur complement and pressure solve abutting it never handle a rotated vector. One un-rotation on the boundary between them feeds both.
Also: a coordinate-system formulation is treated as the same rotation applied globally rather than as a fourth approach, which explains both its appeal and why it stops working the moment the boundary leaves a coordinate surface.
Two things left deliberately visible
Notes
.typ, generator and JSON geometry ship inexamples/, with build instructions in the source header.scripts/rebuild_figures.normalise_svg— Typst emits 652 random glyph ids, whichtest-unitrejects.new_article.pyissued 2026-012, which is already gthyagi's — the very collision that fix exists to prevent, because it is not yet merged.validateclean,test-unit130 passed, PDF builds to 5 pages.Underworld development team with AI support from Claude Code