From 6d2529592403e4d23c419792edc46b547399dc74 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 19 Aug 2026 19:43:38 +1000 Subject: [PATCH 1/2] Return the whole boundary traction from the constrained solver (#607) 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 --- docs/advanced/curved-boundary-conditions.md | 35 ++++-- docs/developer/CHANGELOG.md | 34 +++++ .../design/CONSTRAINED_FREESLIP_MULTIPLIER.md | 57 ++++++--- .../cython/petsc_generic_snes_solvers.pyx | 6 +- src/underworld3/systems/solvers.py | 116 +++++++++++++----- tests/test_1061_constrained_freeslip.py | 18 ++- tests/test_1063_constrained_traction.py | 92 ++++++++++++++ 7 files changed, 305 insertions(+), 53 deletions(-) create mode 100644 tests/test_1063_constrained_traction.py diff --git a/docs/advanced/curved-boundary-conditions.md b/docs/advanced/curved-boundary-conditions.md index 7201d5ef4..ebdbebe10 100644 --- a/docs/advanced/curved-boundary-conditions.md +++ b/docs/advanced/curved-boundary-conditions.md @@ -71,16 +71,31 @@ stokes.add_nitsche_bc(0.0, "Fault", direction=fault_normal, gamma=10) at 1e4 gives 0.15%). -### 2. Penalty Free-Slip (Simple but Fragile) - -Use the mesh-derived normals directly with a penalty parameter: +### 2. Penalty Free-Slip (Simple, and Only With the Node Normal) ```python -Gamma = mesh.Gamma +n = mesh.boundary_normal("Boundary") # measure-weighted node normal penalty = 10000 -stokes.add_natural_bc(penalty * Gamma.dot(v.sym) * Gamma, "Boundary") +stokes.add_natural_bc(penalty * n.dot(v.sym) * n, "Boundary") ``` +**Use `mesh.boundary_normal(boundary)`, not `mesh.Gamma`.** A penalty written +against the per-facet normal asks a node shared by two facets to satisfy two +different constraints, which on a two-component velocity leaves nothing: push +the coefficient up and the boundary freezes. Measured on an annulus with an +exact solution (Kramer et al. 2021), coefficient `1e6`, cell 0.15 → 0.035: + +| normal | leak `u·n` | velocity error | surface stress error | +|---|---|---|---| +| `mesh.Gamma` (facet) | 1e-5 | 0.60, flat under refinement | 0.21 → 0.26, growing | +| `mesh.boundary_normal` (node) | 3e-5 | 1.0e-2 → 4.9e-4 | 2.4e-2 → 1.4e-3 | + +The facet-normal row does not converge, and the leak cannot see it: at `1e3` that +penalty leaks 3e-2 and gets the surface stress right to 2e-3, while at `1e8` it +leaks 1e-7 and is 26% wrong. This is the classical over-constraint that the +*consistent* normal was introduced to avoid (Engelman, Sani & Gresho 1982), and +`mesh.boundary_normal` is that normal. + **When to use:** - Quick prototyping where high accuracy isn't critical - When Nitsche is not yet available for your solver type @@ -88,12 +103,16 @@ stokes.add_natural_bc(penalty * Gamma.dot(v.sym) * Gamma, "Boundary") **Limitations:** - Penalty must be tuned: too small → loose constraint, too large → ill-conditioning - On spherical shells, penalty can become unstable at moderate resolution -- ~25-30% error on elliptical boundaries when using raw facet normals +- Check the answer, not just the leak: a constraint that is satisfied is not + evidence that the solution is right -### 3. Projected Normals (For Curved Boundaries with Penalty) +### 3. Projected Normals (Superseded by `mesh.boundary_normal`) -Project `mesh.Gamma` onto a continuous mesh variable, which interpolates and smooths the normals: +`mesh.boundary_normal(boundary)` assembles the measure-weighted node normal +directly, tracks mesh deformation, and is correct in parallel, so the recipe +below is kept for reference rather than recommended. Project `mesh.Gamma` onto a +continuous mesh variable, which interpolates and smooths the normals: ```python import sympy diff --git a/docs/developer/CHANGELOG.md b/docs/developer/CHANGELOG.md index 53c0126f8..a3cae6699 100644 --- a/docs/developer/CHANGELOG.md +++ b/docs/developer/CHANGELOG.md @@ -6,6 +6,40 @@ This log tracks significant development work at a conceptual level, suitable for ## 2026 Q3 (July – September) +### The Multiplier Was Not the Whole Traction (August 2026) + +**`Stokes_Constrained.topography()` now returns the traction the boundary is +actually held with**, and a new `traction()` exposes it directly. The momentum +row carries `λ + r(n·u − g)`, so the bare multiplier is short by the +augmented-Lagrangian share — `r` times the discrete constraint residual. With the +viscosity-weighted default `r = 1e4·μ(x)` that share is a few per cent of the +surface traction on a uniform-viscosity annulus and most of it across a `1e6` +viscosity step, where `λ` alone reads a tenth of the exact SolCx topography and +is anti-correlated with it. `multiplier()` still returns `λ` and now says what it +is not. + +The defect survived because the validation scored a **correlation** (0.9999) +between the multiplier and the recovered normal stress. A correlation is +scale-free and cannot see a systematic amplitude deficit, which is precisely what +a missing share of the load is. The new guard, +`tests/test_1063_constrained_traction.py`, scores a relative `l2` against the +exact SolCx surface topography and carries the bare multiplier as its negative +control. + +The corrected quantity is the consistent boundary flux: at convergence +`M_Γ(λ + r(n·u − g))` balances the volume residual restricted to the boundary, +which is the CBF nodal load (Zhong, Gurnis & Hulbert 1993). So the multiplier +route and the rotated constraint's `boundary_normal_traction` are the same +computation, and they agree to 3–5% — inside each route's own error against the +exact answer. + +Documentation: `docs/advanced/curved-boundary-conditions.md` now writes the +penalty free-slip recipe against `mesh.boundary_normal` rather than `mesh.Gamma`. +A penalty against the per-facet normal over-constrains the shared nodes and does +not converge — measured on an annulus at coefficient `1e6`, the velocity error +stays at 0.60 and the surface-stress error grows from 0.21 to 0.26 as the mesh is +refined, while the leak reads 1e-5 throughout. (underworld3#607, #608, #614) + ### The Free Surface Reaches the Spherical Shell (July 2026) **`uw.systems.FreeSurface` now runs in 3D on a spherical shell** — the same diff --git a/docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md b/docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md index ce5192e16..982d8f68b 100644 --- a/docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md +++ b/docs/developer/design/CONSTRAINED_FREESLIP_MULTIPLIER.md @@ -2,8 +2,24 @@ **Status**: shipped as `uw.systems.Stokes_Constrained` (serial). The constraint is enforced by a multiplier carried **inside** the saddle point (one coupled -solve); the converged boundary multiplier is the normal traction = dynamic -topography. An earlier augmented-Lagrangian **outer-loop** variant was removed in +solve); the converged boundary traction is `λ + r(n·u − g)`, which is the dynamic +topography, and it is returned by `traction()` / `topography()`. + +> **Correction, 2026-08-19.** This document said in several places that `λ` alone +> is the normal traction and that the augmentation `r` is a pure speed knob. Both +> are true only in the exact limit. Discretely the constraint row is satisfied to +> the solver's tolerance and `r` multiplies that residual back into the traction, +> so `λ` is short by `r(n·u − g)`. With the viscosity-weighted default +> `r = 1e4·μ(x)` the omitted 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 — `λ` alone reads a tenth of the exact topography and is *anti-correlated* +> with it. `traction()` and `topography()` now return the sum; `multiplier()` +> still returns `λ`. See underworld3#607 and +> `tests/test_1063_constrained_traction.py`. +> +> Why the existing validation did not catch it: it scored `corr(λ, −n·σ·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. An earlier augmented-Lagrangian **outer-loop** variant was removed in favour of this in-saddle formulation (it is straightforward to reproduce in Python if needed). Validated against the exact SolCx analytic solution (`tests/test_1062_constrained_solcx.py`). @@ -19,9 +35,9 @@ is ~100× too weak at Ra=1e6); too strong and the system ill-conditions and the Stokes solve diverges in line search. This feature enforces `u·n = g` on a curved boundary with a **true Lagrange -multiplier** `λ` instead of a penalty. Because the converged multiplier *is* the -normal traction holding the boundary, it is simultaneously a direct estimate of -**dynamic surface topography**, `h = λ / (Δρ g)`. The equilibrium `λ` is also the +multiplier** `λ` instead of a penalty. Because the converged boundary term *is* +the normal traction holding the boundary, it is simultaneously a direct estimate +of **dynamic surface topography**, `h = (λ + r(n·u − g)) / (Δρ g)`. The equilibrium `λ` is also the target end-state toward which a free surface can be integrated over a time interval (connecting to the ETD free-surface work on `feature/exp-integrator-freesurface`). @@ -104,9 +120,10 @@ topo = stokes.topography("Upper", buoyancy_scale=delta_rho_g) # h = lambda/(dr ``` `solve()` does **one coupled solve** — no outer iteration or constraint tuning. -The augmentation defaults to `1e4·μ(x)` (local-viscosity-weighted); accuracy is -independent of it (the λ-row carries the exact constraint), so no per-problem -tuning is needed. +The augmentation defaults to `1e4·μ(x)` (local-viscosity-weighted). The +CONSTRAINT is independent of it (the λ-row carries the exact constraint), so no +per-problem tuning is needed, and so is the traction read through `traction()` / +`topography()`. The bare `λ` is NOT: it is short by `r(n·u − g)`. Key design points: @@ -139,15 +156,27 @@ Two regression tests cover the shipped solver: **exact analytic** solution: velocity `rel ≈ 8.7e-6` (== the Dirichlet baseline), constraint `RMS(u·n) ≈ 1.6e-10`. -The consistent-boundary-flux identity `λ = −n·σ·n|_Γ` is the independent -cross-check: the multiplier's boundary trace equals the recovered normal Cauchy -stress (negative sign = the reaction traction holding the boundary), confirming -`λ` is the dynamic topography signal. +The consistent-boundary-flux identity is the independent cross-check, and it is +an identity for the WHOLE boundary term rather than for `λ`: + +``` +M_Γ (λ + r(n·u − g)) = −(A·u − b)|_Γ (the CBF nodal load) +``` + +so `λ + r(n·u − g)` is the CBF traction de-smeared by the boundary mass — the +same computation `boundary_normal_traction` performs on a rotated constraint, +arrived at by carrying the traction as an unknown instead of reading it out of +the residual. Measured on SolCx at a `1e6` contrast, the corrected multiplier +gives relative `l2` 0.047 against the exact topography where `λ` alone gives 1.04. + +A third test, `tests/test_1063_constrained_traction.py`, guards this with the +bare multiplier as its negative control. ## The augmentation parameter `r`: true-work trade-off -`r` is a *speed* knob, not an *accuracy* knob — this is the key advantage over a -pure penalty, and it carries over to the in-saddle solver (accuracy is +`r` is a *speed* knob for the CONSTRAINT, not an accuracy knob — this is the key +advantage over a pure penalty, and it carries over to the in-saddle solver +(accuracy is `r`-independent; `r` only sets the iteration count). The sweep table below is from the **historical outer-loop** variant (its "outer iterations" have no analogue in the one-shot coupled solve), but the shape and the conclusion stand. For the diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c42bbb474..b0400b517 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -8110,8 +8110,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass): # u-row residual: fn_f = h·n + r(n·u − g)·n # The r-term is the augmented-Lagrangian penalty: it adds a uu # boundary stiffness r·(n⊗n) that conditions the Schur complement - # but does NOT bias the multiplier (the h-row stays the exact - # constraint, so h still converges to the true normal traction). + # but does not change what the constraint ENFORCES (the h-row stays + # the exact constraint). It does change what h IS: the traction is + # h + r(n.u - g), and only the sum is r-independent. Stokes_Constrained + # .traction() / .topography() return that sum; .multiplier() returns h. fn_f = sympy.Matrix( [(hsym + r_sym * (u_dot_n - g_sym)) * n_row[i] for i in range(dim)] ).as_immutable() diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py index 98d0c1b7f..1365b2d1e 100644 --- a/src/underworld3/systems/solvers.py +++ b/src/underworld3/systems/solvers.py @@ -2327,7 +2327,7 @@ class SNES_Stokes_Constrained(SNES_Stokes): The constraint is enforced in one coupled solve (no outer iteration). The augmented-Lagrangian term conditions the :math:`[p,h]` Schur complement - without biasing the multiplier, and the interior multiplier DOFs are reduced + without changing what the constraint enforces, and the interior multiplier DOFs are reduced away so the solved block is boundary-sized. Runs in parallel: the interior-multiplier reduction is rank-local section @@ -2612,9 +2612,11 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No Adds a scalar multiplier field ``h`` coupled into the saddle-point system so that :math:`\mathbf{u}\cdot\mathbf{n}=\mathrm{conds}` is enforced on - ``boundary`` in the coupled solve; at convergence ``h`` on the boundary - is the normal traction (dynamic topography), recoverable via - :meth:`multiplier` / :meth:`topography`. + ``boundary`` in the coupled solve. At convergence the boundary traction — + the dynamic topography before scaling — is + :math:`h + r(\mathbf{u}\cdot\mathbf{n} - g)`, returned by :meth:`traction` + and :meth:`topography`. :meth:`multiplier` returns the field :math:`h`, + which is that sum less the augmented-Lagrangian share. Parameters ---------- @@ -2635,14 +2637,20 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No :math:`r(\mathbf{n}\cdot\mathbf{u}-g)\,\mathbf{n}` to the u-row, giving a ``uu`` boundary stiffness :math:`r\,(\mathbf{n}\otimes \mathbf{n})` that conditions the :math:`[p,h]` Schur complement - **without biasing the multiplier** (the h-row is still the exact - constraint). Defaults to ``augmentation_base · μ(x)`` (viscosity- + without changing what the constraint enforces (the h-row is still the + exact constraint). Defaults to ``augmentation_base · μ(x)`` (viscosity- weighted, mesh-independent). Pass ``0`` for the bare KKT system. + + It DOES change what :math:`h` is: the traction is + :math:`h + r(\mathbf{u}\cdot\mathbf{n} - g)` and only the sum is + :math:`r`-independent. Read it through :meth:`traction` or + :meth:`topography`, never off :meth:`multiplier` alone. augmentation_base : float, default 1e4 - Base multiple used when ``augmentation`` is not given. Accuracy is - independent of this value (the multiplier carries the exact - constraint); larger values reduce the iteration count up to a broad - plateau, well below the roundoff limit. + Base multiple used when ``augmentation`` is not given. The CONSTRAINT + is independent of this value and so is the traction read through + :meth:`traction`; larger values reduce the iteration count up to a + broad plateau. The bare multiplier :math:`h` is not independent of it + — see :meth:`multiplier`. g : float or sympy expression, optional Deprecated keyword alias for ``conds`` (one DeprecationWarning). @@ -2742,19 +2750,69 @@ def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=No return h def multiplier(self, boundary): - """Return the multiplier field for ``boundary`` (None if not constrained). - - After :meth:`solve`, the multiplier's boundary trace is the normal - traction holding the constraint. Divide by :math:`\\Delta\\rho\\,g` for - dynamic topography (see :meth:`topography`). + """Return the multiplier FIELD for ``boundary`` (None if not constrained). + + .. warning:: + This is :math:`h` itself, which is **not** the whole boundary traction + whenever an augmented-Lagrangian term is in place (it is by default). + The momentum row carries + :math:`(h + r(\\mathbf{u}\\cdot\\mathbf{n} - g))\\mathbf{n}`, so the + traction holding the boundary is that sum — see :meth:`traction`, and + :meth:`topography` which is built on it. Measured on SolCx at a + viscosity contrast of :math:`10^6`, where the viscosity-weighted default + :math:`r = 10^4\\mu` reaches :math:`10^{10}` on the stiff half, ``h`` + alone carries about a tenth of the surface traction and is + anti-correlated with it. """ for cbc in self._block_constraint_bcs: if cbc.boundary == boundary: return cbc.lam return None + def _constraint_bc(self, boundary): + """The registered constraint record for ``boundary``, or raise.""" + for cbc in self._block_constraint_bcs: + if cbc.boundary == boundary: + return cbc + raise ValueError(f"No constraint registered on boundary '{boundary}'.") + + def traction(self, boundary): + r"""Boundary normal traction on ``boundary``, as a symbolic expression. + + The momentum row's boundary term is + :math:`(h + r(\mathbf{u}\cdot\mathbf{n} - g))\,\mathbf{n}`, so the traction + holding the constraint is + + .. math:: \sigma_{nn} = h + r(\mathbf{u}\cdot\mathbf{n} - g), + + with :math:`r` the augmented-Lagrangian parameter and :math:`g` the + prescribed normal velocity. The second term vanishes only where the + constraint row is satisfied exactly; discretely it is satisfied to the + solver's tolerance, and :math:`r` multiplies that residual straight back + into the traction. With the viscosity-weighted default + :math:`r = 10^4\mu(x)` the omitted share is a few per cent of the surface + traction on a uniform-viscosity annulus and most of it across a + :math:`10^6` viscosity step. + + This is the same quantity the consistent boundary flux back-calculation + recovers (Zhong, Gurnis & Hulbert 1993): at convergence the assembled + boundary load :math:`M_\Gamma(h + r(\mathbf{u}\cdot\mathbf{n}-g))` balances + the volume residual restricted to the boundary, which is the CBF nodal + load, so the two differ only by the mass de-smear. + + Valid ON ``boundary``; the expression involves the multiplier field, whose + interior degrees of freedom are constrained out of the solve. + """ + cbc = self._constraint_bc(boundary) + u_dot_n = sum(cbc.normal[i] * self.u.sym[i] for i in range(self.mesh.dim)) + return cbc.lam.sym[0] + cbc.augmentation * (u_dot_n - cbc.g) + def topography(self, boundary, buoyancy_scale=1.0, reference=None): - r"""Dynamic topography expression :math:`h / (\Delta\rho\, g)` on ``boundary``. + r"""Dynamic topography on ``boundary``, as a symbolic expression. + + :math:`\sigma_{nn} / (\Delta\rho\, g)` with :math:`\sigma_{nn}` the full + boundary traction :math:`h + r(\mathbf{u}\cdot\mathbf{n} - g)` — see + :meth:`traction`, and :meth:`multiplier` for why :math:`h` alone is not it. For an **enclosed** problem (no net normal flow through any boundary) the multiplier :math:`h` is determined only up to the :math:`[p,\lambda]` gauge @@ -2763,14 +2821,14 @@ def topography(self, boundary, buoyancy_scale=1.0, reference=None): the absolute level of :math:`h` is not reproducible across ranks. For such problems pass ``reference="mean"`` to subtract the boundary mean and obtain a gauge-fixed, partition-independent topography. The default - (``reference=None``) returns the raw multiplier — correct for problems with - **no** gauge freedom (e.g. an open boundary), where the mean of :math:`h` is - the physical mean traction and must NOT be removed. + (``reference=None``) returns the traction unshifted — correct for problems + with **no** gauge freedom (e.g. an open boundary), where its mean is the + physical mean traction and must NOT be removed. Note that the automatic pressure gauge (:attr:`auto_pressure_gauge`) fixes the raw *pressure* level but NOT the raw *multiplier* level (the constant multiplier is an independent gauge freedom). So on an enclosed problem the - raw multiplier (``reference=None``) is still partition-dependent — + unshifted traction (``reference=None``) is still partition-dependent — ``reference="mean"`` is the gauge-invariant, partition-reproducible read for dynamic topography and is the recommended path. @@ -2781,8 +2839,8 @@ def topography(self, boundary, buoyancy_scale=1.0, reference=None): buoyancy_scale : float, default 1.0 Divide by :math:`\Delta\rho\,g` to convert traction to length. reference : {None, "mean"}, default None - ``None`` returns the raw multiplier (correct when there is no gauge - freedom). ``"mean"`` subtracts the boundary mean (gauge-fixed, + ``None`` returns the traction unshifted (correct when there is no + gauge freedom). ``"mean"`` subtracts the boundary mean (gauge-fixed, reproducible) — use for enclosed problems. Notes @@ -2793,17 +2851,19 @@ def topography(self, boundary, buoyancy_scale=1.0, reference=None): single-rank branch). ``reference=None`` is a pure symbolic accessor with no reduction. """ - lam = self.multiplier(boundary) - if lam is None: - raise ValueError(f"No constraint registered on boundary '{boundary}'.") - expr = lam.sym[0] + # The WHOLE traction, not the bare multiplier: with an augmented-Lagrangian + # term in place the momentum row carries h + r(u.n - g), and r multiplies + # the discrete constraint residual back into the traction. Reading h alone + # was wrong by a few per cent on a uniform-viscosity annulus and by an + # order of magnitude (and a sign) across a 1e6 viscosity step. + expr = self.traction(boundary) if reference == "mean": - # Subtract the boundary mean of h via parallel-safe surface integrals + # Subtract the boundary mean via parallel-safe surface integrals # (BdIntegral handles the cross-rank reduction); this fixes the gauge. blen = uw.maths.BdIntegral( mesh=self.mesh, fn=sympy.Integer(1), boundary=boundary).evaluate() hbar = uw.maths.BdIntegral( - mesh=self.mesh, fn=lam.sym[0], boundary=boundary).evaluate() / blen + mesh=self.mesh, fn=expr, boundary=boundary).evaluate() / blen expr = expr - hbar elif reference is not None: raise ValueError( diff --git a/tests/test_1061_constrained_freeslip.py b/tests/test_1061_constrained_freeslip.py index 3d0b01a73..c81387ee1 100644 --- a/tests/test_1061_constrained_freeslip.py +++ b/tests/test_1061_constrained_freeslip.py @@ -94,10 +94,26 @@ def test_box_topography_is_normal_traction(box): def test_multiplier_and_topography_api(box): + """The accessors return what they say, and topography is the WHOLE traction. + + This asserted `topography == h / scale` until 2026-08-19, which pinned the + defect in underworld3#607: the momentum row carries `h + r(n.u - g)`, so the + bare multiplier is short by the augmented-Lagrangian share. The outcome to + assert is that topography is the traction divided by the scale, and that the + traction is not just the multiplier. + """ blk, hL = box["blk"], box["hL"] assert blk.multiplier("Left") is hL assert blk.multiplier("Nonexistent") is None - assert blk.topography("Left", buoyancy_scale=2.0) == hL.sym[0] / 2.0 + + traction = blk.traction("Left") + assert sympy.simplify(blk.topography("Left", buoyancy_scale=2.0) + - traction / 2.0) == 0 + + # the share the multiplier alone leaves out is really there + share = sympy.simplify(traction - hL.sym[0]) + assert share != 0 + assert blk.u.sym[0] in share.atoms(sympy.Function) def test_rejects_unknown_boundary(box): diff --git a/tests/test_1063_constrained_traction.py b/tests/test_1063_constrained_traction.py new file mode 100644 index 000000000..003e747bb --- /dev/null +++ b/tests/test_1063_constrained_traction.py @@ -0,0 +1,92 @@ +"""The multiplier is not the whole traction — the augmented-Lagrangian share is. + +`Stokes_Constrained` assembles the momentum row's boundary term as +`(h + r(n.u - g)) n`, so the traction holding the boundary is that sum. `h` alone +is short by `r` times the discrete constraint residual, and the default `r` is +viscosity-weighted (`1e4 * mu`), so a lateral viscosity contrast makes the +omitted share most of the answer. + +SolCx is the case that shows it: `mu` steps from 1 to 1e6 at x = 0.5, so +`r = 1e10` on the stiff half. `uw.analytic.SolCx.topography_top` is the exact +surface topography, which is what makes this an oracle rather than a comparison. + +Guards `Stokes_Constrained.traction()` and `topography()`, and carries the +NEGATIVE CONTROL in the same test: the bare multiplier fails the same threshold +by a wide margin, so a regression that quietly reverts to `h` cannot pass. + +Run: pixi run python -m pytest tests/test_1063_constrained_traction.py -v +""" + +import pytest + +pytestmark = [pytest.mark.level_2] + +import numpy as np +import sympy +import underworld3 as uw +from underworld3 import analytic as A + +ETA_B, RES = 1.0e6, 32 + + +def _solve(): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(RES, RES), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3 + ) + sol = A.SolCx(mesh, eta_A=1.0, eta_B=ETA_B, x_c=0.5, n=1) + + s = uw.systems.Stokes_Constrained(mesh) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = sol.fn_viscosity + s.bodyforce = sol.fn_bodyforce + # Three walls by the ordinary component condition, so the multiplier under + # test is the only constraint on the wall being measured. + s.add_dirichlet_bc((0.0, None), "Left") + s.add_dirichlet_bc((0.0, None), "Right") + s.add_dirichlet_bc((None, 0.0), "Bottom") + s.add_constraint_bc(0.0, "Top") + s.petsc_use_pressure_nullspace = True + s.tolerance = 1.0e-9 + s.solve() + assert s.snes.getConvergedReason() > 0 + return mesh, s, sol + + +def _relative_error(values, coords, sol): + """Relative l2 against the exact topography, both mean-removed. + + The box is enclosed, so the level is a gauge and only the deviation is + determined — which is what topography is anyway. + """ + exact = np.asarray(sol.topography_top(coords)).reshape(-1) + got = np.asarray(values).reshape(-1) - np.mean(values) + exact = exact - exact.mean() + return float(np.linalg.norm(got - exact) / np.linalg.norm(exact)) + + +def test_constrained_topography_carries_the_augmentation_share(): + mesh, s, sol = _solve() + + top = np.abs(s.u.coords[:, 1] - 1.0) < 1.0e-9 + coords = s.u.coords[top] + + traction = uw.function.evaluate(s.traction("Top"), coords) + bare = uw.function.evaluate(s.multiplier("Top").sym[0], coords) + + whole = _relative_error(traction, coords, sol) + without = _relative_error(bare, coords, sol) + + # The traction is the surface topography, to the accuracy of the discretisation. + assert whole < 0.15, f"traction wrong by {whole:.3f}" + + # NEGATIVE CONTROL. Without the augmentation share the same read is useless + # here: measured 1.04 (it is anti-correlated with the right answer). If this + # ever passes, `traction()` has quietly become `multiplier()` again. + assert without > 0.5, ( + f"the bare multiplier read {without:.3f}: the negative control no longer " + "fires, so this test is not guarding anything") + + # topography() is the traction divided by the buoyancy scale, so it inherits + # the fix; check the public path a user actually calls. + height = uw.function.evaluate(s.topography("Top", reference="mean"), coords) + assert _relative_error(height, coords, sol) < 0.15 From 12ddbba1714a05a9fd44e76791a6c25b4e934c77 Mon Sep 17 00:00:00 2001 From: lmoresi Date: Wed, 19 Aug 2026 20:15:13 +1000 Subject: [PATCH 2/2] Hold the corner: a rotated wall meeting an essential one (#616, #608, #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 --- src/underworld3/utilities/boundary_flux.py | 14 ++++ src/underworld3/utilities/rotated_bc.py | 54 ++++++++++++++- tests/test_1066_rotated_meets_essential.py | 79 ++++++++++++++++++++++ 3 files changed, 146 insertions(+), 1 deletion(-) create mode 100644 tests/test_1066_rotated_meets_essential.py diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index c0b8e5b2e..20aef5bd3 100644 --- a/src/underworld3/utilities/boundary_flux.py +++ b/src/underworld3/utilities/boundary_flux.py @@ -637,6 +637,20 @@ 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).""" + # A multiplier-constrained boundary has no reaction left to read. The constraint + # is imposed by a term in the SAME row, so the assembled residual there is + # balanced at convergence and this back-calculation returns ~0 (measured: rms + # 4e-13 against a traction of 0.37). The traction is carried by the multiplier — + # M_Gamma(h + r(n.u - g)) IS the CBF nodal load — so send the caller there rather + # than hand back a quiet zero. + for cbc in getattr(solver, "_block_constraint_bcs", ()): + if cbc.boundary == boundary: + raise RuntimeError( + f"'{boundary}' is held by a multiplier constraint, so the consistent " + "boundary flux reads ~0 there: the constraint term balances the row it " + "sits in. Use solver.traction(boundary) — h + r(n.u - g) — which is the " + "same quantity, or solver.topography(boundary) for the scaled version.") + 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/src/underworld3/utilities/rotated_bc.py b/src/underworld3/utilities/rotated_bc.py index 8400a887b..b62c48695 100644 --- a/src/underworld3/utilities/rotated_bc.py +++ b/src/underworld3/utilities/rotated_bc.py @@ -425,7 +425,44 @@ def build_rotation(solver, boundaries, datum_specs=None): for q, nrms in node_normals.items(): lo = lsec.getFieldOffset(q, _VELOCITY_FIELD) grows = [int(l2g.apply([lo + c])[0]) for c in range(dim)] - if any(g < 0 for g in grows): + free = [c for c in range(dim) if grows[c] >= 0] + if not free: # every component pinned already + continue + if len(free) < dim: + # PARTIALLY CONSTRAINED NODE — a rotated wall meeting an essential one. + # Some components are constrained out of the global vector (g < 0) and + # the rest are free. This used to `continue`, which left the wall-normal + # component UNCONSTRAINED at those nodes: the wall leaked at its own end + # points while every interior node was exact. Measured on a unit box with + # a rotated lid and component free slip on the other three walls, + # 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 2e-3 globally, with an exact linear solve on both + # sides (issue #616; the corner reaction of #608 is the same node). + # + # The constraint is still imposable on what is left: with the pinned + # components held at zero, n̂·v = 0 reduces to n̂_F·v_F = 0 on the free + # subspace F. Build the frame there and constrain its normal rows. + if q in node_dspec: + # A prescribed v_n datum at such a node needs the pinned components' + # values to reduce the affine constraint, which are not read here. + # Preserve the previous behaviour rather than impose the wrong datum. + _warn_once_partial_datum() + continue + Mf = np.array(nrms, dtype=float)[:, free] + scale = float(np.linalg.norm(np.array(nrms, dtype=float))) + if float(np.linalg.norm(Mf)) <= 1e-12 * max(scale, 1.0): + # the normal lies entirely in the pinned subspace: already implied + continue + rows = [grows[c] for c in free] + if not (rstart <= rows[0] < rend): # not owned by this rank → skip + continue + _, svf, Vtf = np.linalg.svd(Mf) + rf = int((svf > 1e-8 * (svf[0] if svf.size else 1.0)).sum()) + for i in range(len(free)): + for j in range(len(free)): + Q.setValue(rows[i], rows[j], float(Vtf[i, j])) + normal_rows.extend(rows[:rf]) continue if not (rstart <= grows[0] < rend): # not owned by this rank → skip continue @@ -544,6 +581,21 @@ def build_rotation(solver, boundaries, datum_specs=None): return Q, Qt, sorted(set(normal_rows)), datum_map +_PARTIAL_DATUM_WARNED = [False] + + +def _warn_once_partial_datum(): + """A prescribed wall-normal datum at a node shared with an essential BC is not + reduced here, so that node keeps the pre-#616 behaviour (unconstrained). Say so + once rather than silently.""" + if not _PARTIAL_DATUM_WARNED[0]: + _PARTIAL_DATUM_WARNED[0] = True + print("[rotated_bc] WARNING: a prescribed v_n datum sits on a node shared " + "with an essential BC; the wall-normal component is left free there " + "(the affine reduction against the pinned components is not " + "implemented). Free-slip nodes are unaffected.") + + def _zero_rows_local(vec, normal_rows): """Zero ``vec`` at the global rows ``normal_rows`` using ownership-relative local indices (indexing the local slice with global rows overflows on any rank diff --git a/tests/test_1066_rotated_meets_essential.py b/tests/test_1066_rotated_meets_essential.py new file mode 100644 index 000000000..f50c74c71 --- /dev/null +++ b/tests/test_1066_rotated_meets_essential.py @@ -0,0 +1,79 @@ +"""A rotated free-slip wall meeting an essential wall must still hold its corner. + +On a flat, axis-aligned lid the rotated constraint and `add_dirichlet_bc((None, +0.0), "Top")` are the same discrete constraint: the measure-weighted node normal +is exactly (0,1), so striking the rotated normal row is striking u_y. + +They were not the same. `build_rotation` skipped any node with a velocity DOF +constrained out of the global vector — which is every node where the lid meets a +wall held by an essential condition — so the wall-normal component was left FREE +at those nodes. The lid leaked at its own end points (max|u_y| = 4.0e-3 against +|u|max 2.5e-2, entirely at the two corners) and the solve differed from the +component-Dirichlet one by 2e-3 globally, with an exact linear solve on both +sides. Issue #616; the corner reaction of #608 is the same node. + +Run: pixi run python -m pytest tests/test_1066_rotated_meets_essential.py -v +""" + +import pytest + +pytestmark = [pytest.mark.level_1] + +import numpy as np +import sympy +import underworld3 as uw + +RES = 16 + + +def _solve(lid, direct=False): + mesh = uw.meshing.StructuredQuadBox( + elementRes=(RES, RES), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), qdegree=3 + ) + x, z = mesh.X + v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2) + p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1) + s = uw.systems.Stokes(mesh, velocityField=v, pressureField=p) + s.constitutive_model = uw.constitutive_models.ViscousFlowModel + s.constitutive_model.Parameters.shear_viscosity_0 = 1.0 + s.saddle_preconditioner = 1.0 + s.bodyforce = sympy.Matrix([[0, sympy.cos(sympy.pi * x) * sympy.sin(sympy.pi * z)]]) + s.tolerance = 1.0e-9 + s.petsc_use_pressure_nullspace = True + s.add_dirichlet_bc((0.0, None), "Left") + s.add_dirichlet_bc((0.0, None), "Right") + s.add_dirichlet_bc((None, 0.0), "Bottom") + if lid == "dirichlet": + s.add_dirichlet_bc((None, 0.0), "Top") + if direct: + s.petsc_options["ksp_type"] = "preonly" + s.petsc_options["pc_type"] = "lu" + else: + s.add_rotated_freeslip_bc(0.0, "Top") + if direct: + s._rotated_use_lu = True + s.solve() + assert s.snes.getConvergedReason() > 0 + return np.squeeze(np.asarray(v.array)).copy(), v.coords.copy() + + +def test_rotated_lid_holds_its_corners(): + """u.n = 0 on the whole lid, corners included, to round-off.""" + u, coords = _solve("rotated") + top = np.abs(coords[:, 1] - 1.0) < 1.0e-9 + leak = np.abs(u[top, 1]).max() / np.abs(u).max() + # It was 1.6e-1, at the two corner nodes and nowhere else. + assert leak < 1.0e-12, f"the lid leaks at u.n = {leak:.2e} of |u|max" + + +def test_rotated_lid_reproduces_the_component_condition(): + """The two are the same discrete constraint on a flat wall, so the same answer. + + Judged against a monolithic DIRECT solve of the component-Dirichlet problem, so + a difference cannot be blamed on either iterative solve. + """ + reference, _c = _solve("dirichlet", direct=True) + for lid in ("dirichlet", "rotated"): + u, _c = _solve(lid) + rel = np.linalg.norm(u - reference) / np.linalg.norm(reference) + assert rel < 1.0e-7, f"{lid} lid differs from the direct reference by {rel:.2e}"