diff --git a/docs/advanced/curved-boundary-conditions.md b/docs/advanced/curved-boundary-conditions.md index 7201d5ef..ebdbebe1 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 53c0126f..a3cae669 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 ce5192e1..982d8f68 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 c42bbb47..b0400b51 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 98d0c1b7..1365b2d1 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/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py index c0b8e5b2..20aef5bd 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 8400a887..b62c4869 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_1061_constrained_freeslip.py b/tests/test_1061_constrained_freeslip.py index 3d0b01a7..c81387ee 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 00000000..003e747b --- /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 diff --git a/tests/test_1066_rotated_meets_essential.py b/tests/test_1066_rotated_meets_essential.py new file mode 100644 index 00000000..f50c74c7 --- /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}"