Issue quick-fix batch: units/pint crash, pressure gauge under units, two stale parallel tests, scalar-field default, penalty migration note

docs/advanced/troubleshooting.md

@@ -1,13 +1,52 @@
 ---
-title: "trouuleshooting"
+title: "Troubleshooting"
 ---
 
-# trouuleshooting
+# Troubleshooting
 
-```{note} Documentation In Progress
-Content being developed for advanced users.
-```
+Notes on behaviour changes and common pitfalls that show up as "my old script
+runs but gives different answers / different performance".
 
-## Overview
+## Stokes penalty is now viscosity-scaled (June 2026)
 
-Advanced topic documentation coming soon.
+The Stokes augmented-Lagrangian (grad-div) incompressibility penalty changed
+from a bare constant to a **viscosity-scaled** term. The operator changed
+from
+
+$$
+\sigma + \lambda \, (\nabla \cdot \mathbf{u}) \, I
+$$
+
+to
+
+$$
+\sigma + \lambda \, \mu \, (\nabla \cdot \mathbf{u}) \, I
+$$
+
+where $\mu$ is the local viscosity (`constitutive_model.K`) and $\lambda$ is
+the value you set with `stokes.penalty`.
+
+Why: with spatially variable viscosity, a constant penalty is enormous
+relative to the stress in low-viscosity regions (over-stiffening them into
+velocity locking) and negligible in high-viscosity regions. Scaling by the
+local viscosity keeps the penalty proportional to the local stress scale
+everywhere.
+
+**What to change in older scripts:**
+
+1. `stokes.penalty` is now a *dimensionless* number of order 1. Scripts that
+   tuned a large constant against the viscosity magnitude (say `penalty=1e6`
+   for a model with $\mu \sim 1$, or a value chosen to sit above the largest
+   viscosity in a contrast model) should drop that tuning — use `1.0` (or a
+   modest multiple) instead. Keeping the old large value multiplies it by
+   the viscosity again and can lock the velocity field or destroy the
+   conditioning of the velocity block.
+2. A penalty is still **off by default** (`penalty = 0`) and usually
+   unnecessary: the pressure Schur complement is already preconditioned with
+   the local $1/\mu$ scaling.
+3. The penalty term is part of the operator, so converged pressures include
+   the grad-div contribution; diagnostics that compared pressure against a
+   constant-penalty run will differ at the level of the divergence residual.
+
+See the `Stokes.penalty` property docstring for the full description, and the
+`CONSTRAINED_FREESLIP_MULTIPLIER` design note for the derivation.

docs/developer/CHANGELOG.md

@@ -199,6 +199,11 @@ in `Stokes_Constrained` (#224), then made parallel-correct.
 - `selfp` Schur preconditioner default, viscosity-scaled penalty, and
   nullspace re-setup fix (#229); over-conservative serial guard removed
   (#240); gauge, convergence, knockout, and rotation-gauge fixes (#265).
+- The main `Stokes.penalty` (augmented-Lagrangian grad-div) is likewise
+  viscosity-scaled since June 2026: the parameter is now a dimensionless
+  O(1) number, not a large constant tuned against the viscosity magnitude.
+  Migration note for older scripts in `docs/advanced/troubleshooting.md`
+  (#292).
 
 ### Boundary Conditions: Local-h Nitsche and Boundary-Slip Surfaces (June 2026)
 

src/underworld3/cython/petsc_generic_snes_solvers.pyx

@@ -2406,9 +2406,11 @@ class SNES_Scalar(SolverBaseClass):
 
         ## Todo: some validity checking on the size / type of u_Field supplied
         if u_Field is None:
-            # TODO(BUG): num_components=mesh.dim for a SCALAR unknown looks wrong
-            # (a scalar has one component); kept as-is pending review (READ-22).
-            self.Unknowns.u = uw.discretisation.MeshVariable( mesh=mesh, num_components=mesh.dim,
+            # A scalar unknown has one component. (MeshVariable already forced
+            # this: vtype=SCALAR overrides num_components, so the old
+            # num_components=mesh.dim here was ignored, never over-allocated —
+            # issue #367.)
+            self.Unknowns.u = uw.discretisation.MeshVariable( mesh=mesh, num_components=1,
                                                       varname="Us{}".format(SNES_Scalar._obj_count),
                                                       vtype=uw.VarType.SCALAR, degree=degree, )
         else:

src/underworld3/systems/solvers.py

@@ -1896,18 +1896,35 @@ def set_pressure_gauge(self, boundary, reference=0.0):
         boundary : str
             Mesh boundary label over which the mean pressure is fixed.
         reference : float, default 0.0
-            Target mean pressure on ``boundary``.
+            Target mean pressure on ``boundary``, in the model (non-dimensional)
+            frame.
 
         Returns
         -------
         the registered callback (so it can be identified/removed later).
+
+        Notes
+        -----
+        For **pressure-dependent plasticity** (e.g. Drucker-Prager,
+        :math:`\sigma_y = C + p\sin\varphi`) a pressure gauge is
+        counter-productive: shifting the pressure moves the yield surface at
+        every iteration and the nonlinear solve oscillates instead of
+        converging. The pressure level there is physical, not a free gauge —
+        set the level through the boundary conditions instead.
         """
         p = self.Unknowns.p
         area = uw.maths.BdIntegral(self.mesh, 1.0, boundary).evaluate()
         p_surface_integral = uw.maths.BdIntegral(self.mesh, p.sym[0, 0], boundary)
 
         def _pressure_gauge(solver, iteration):
             mean = p_surface_integral.evaluate() / area
+            # This callback fires inside the SNES non-dimensionalisation
+            # cordon, where p.data holds non-dimensional values. Under an
+            # active units model the boundary integral is a dimensional
+            # quantity, so take its non-dimensional value before applying
+            # the shift (issue #271).
+            if isinstance(mean, uw.function.quantities.UWQuantity):
+                mean = mean.data
             p.data[...] -= (mean - reference)
 
         return self.add_update_callback(_pressure_gauge)

src/underworld3/units.py

@@ -876,15 +876,16 @@ def non_dimensionalise(expression, model=None) -> Any:
 
             try:
                 scale = model.get_scale_for_dimensionality(dimensionality)
-                result_qty = expression / scale
-                # TODO(BUG): UWQuantity.__init__ has no `dimensionality`
-                # kwarg, so this raises TypeError whenever a raw
-                # pint.Quantity is non-dimensionalised with active scaling
-                # and a resolvable scale. Found while verifying BF-12
-                # (2026-07 audit). The uw.quantity (UWQuantity) path works
-                # and is the workaround.
-                # Return UWQuantity to preserve dimensionality
-                return UWQuantity(float(result_qty.magnitude), units="dimensionless", dimensionality=dimensionality)
+                # Convert to SI base units BEFORE dividing, matching the
+                # UnitAwareArray path above: the scale is computed in SI, so
+                # a non-SI input (km, cm/year, ...) must cancel through SI
+                # or the magnitude keeps the unit prefactor.
+                result_qty = expression.to_base_units() / scale
+                # A non-dimensionalised scalar is a plain float, matching the
+                # .magnitude returns above. (A dimensionless UWQuantity cannot
+                # carry the source dimensionality; use dimensionalise() with
+                # an explicit target_dimensionality to round-trip.)
+                return float(result_qty.magnitude)
             except ValueError as e:
                 raise ValueError(f"Cannot non-dimensionalise Pint quantity: {e}")
     except ImportError:

tests/parallel/test_0700_basic_parallel_operations.py

@@ -242,8 +242,10 @@ def test_nodal_swarm_advection_basic():
     """
     Test basic nodal swarm advection in parallel.
 
-    This is a simplified version of ptest_004 that verifies the advection
-    mechanism works without running a full simulation.
+    This carries the coverage of the retired ptest_004 script (issue #340):
+    it verifies the advection mechanism works without running a full
+    simulation. (The particle-count assertion ptest_004 made is obsolete —
+    the SL trace-back no longer allocates a nodal swarm.)
 
     Requires 4 ranks to test cross-processor particle movement.
     """
@@ -261,12 +263,14 @@ def test_nodal_swarm_advection_basic():
     # Vector being advected
     vec_tst = uw.discretisation.MeshVariable("Vn", mesh, mesh.dim, degree=2)
 
-    # Create semi-Lagrangian advection
+    # Create semi-Lagrangian advection (degree/continuous match vec_tst)
     DuDt = uw.systems.ddt.SemiLagrangian(
         mesh,
         vec_tst.sym,
         v.sym,
         vtype=uw.VarType.VECTOR,
+        degree=2,
+        continuous=vec_tst.continuous,
         order=2,
         smoothing=0.0,
     )

tests/parallel/test_0765_internal_boundary_integral_mpi.py

@@ -87,11 +87,10 @@ def test_deformed_spherical_shell_boundary_area_parallel():
     a = math.log(2.0)
     mapped = (np.exp(a * t) - 1.0) / (math.exp(a) - 1.0)
     new_radii = 0.5 + thickness * mapped
-    # TODO(BUG): fails since the _deform_mesh live-mesh guard landed (PR #326):
-    # a direct _deform_mesh call on a mesh that already carries a variable is
-    # rejected. Needs migration to mesh.deform() or a sanctioned mutation
-    # scope. Pre-existing on development; observed while gating the #324 fix.
-    mesh._deform_mesh(coords * (new_radii / radii)[:, None])
+    # Public deform() (pure geometric move, no dt): the raw _deform_mesh
+    # primitive is rejected by the live-mesh coordinate-mutation guard when
+    # a variable exists (issue #331).
+    mesh.deform(coords * (new_radii / radii)[:, None])
 
     lower = float(uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Lower").evaluate())
     upper = float(uw.maths.BdIntegral(mesh=mesh, fn=1.0, boundary="Upper").evaluate())

tests/test_0819_units_nondim_regressions.py

@@ -0,0 +1,79 @@
+"""Regressions at the units <-> non-dimensional boundary.
+
+Covers two crashes found with an active scaling model:
+
+- issue #328: ``uw.non_dimensionalise(pint.Quantity)`` raised ``TypeError``
+  (the protocol-5 branch constructed UWQuantity with an invalid
+  ``dimensionality=`` keyword whenever a scale was resolvable);
+- issue #271: ``set_pressure_gauge`` computed the surface-mean pressure as a
+  dimensional quantity but applied it to the non-dimensional ``p.data``
+  inside the SNES update callback, crashing the solve
+  (``_UFuncOutputCastingError``).
+"""
+
+import numpy as np
+import pytest
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_b]
+
+
+@pytest.fixture()
+def scaling_model():
+    """Active reference scaling covering length, time and mass."""
+    uw.reset_default_model()
+    orchestration_model = uw.get_default_model()
+    orchestration_model.set_reference_quantities(
+        length=uw.quantity(500, "m"),
+        velocity=uw.quantity(1, "cm/year"),
+        viscosity=uw.quantity(1e21, "Pa*s"),
+    )
+    yield orchestration_model
+    uw.reset_default_model()
+
+
+def test_non_dimensionalise_plain_pint_quantity(scaling_model):
+    # issue #328: a raw pint quantity with a resolvable scale must
+    # non-dimensionalise to a plain float, not raise TypeError.
+    result = uw.non_dimensionalise(250 * uw.units("m"))
+    assert isinstance(result, float)
+    assert result == pytest.approx(0.5)  # length scale is 500 m
+
+
+def test_non_dimensionalise_pint_quantity_non_si_units(scaling_model):
+    # The scale is computed in SI, so a non-SI input must reduce through SI:
+    # 0.25 km == 250 m == 0.5 length-scales.
+    assert uw.non_dimensionalise(0.25 * uw.units("km")) == pytest.approx(0.5)
+
+
+def test_pressure_gauge_solves_under_units_model(scaling_model):
+    # issue #271: the gauge callback fires inside the SNES non-dimensional
+    # cordon; the dimensional surface-mean must be non-dimensionalised
+    # before it is subtracted from p.data. Before the fix this solve
+    # crashed with an object-dtype ufunc casting error.
+    mesh = uw.meshing.StructuredQuadBox(
+        elementRes=(8, 8), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0),
+        units="metre", qdegree=3)
+    v = uw.discretisation.MeshVariable("U", mesh, mesh.dim, degree=2, units="cm/year")
+    p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1, units="Pa")
+
+    stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
+    stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    stokes.constitutive_model.Parameters.shear_viscosity_0 = uw.quantity(1e21, "Pa*s")
+    stokes.add_dirichlet_bc((uw.quantity(1.0, "cm/year"), 0.0), "Top")
+    stokes.add_dirichlet_bc((0.0, 0.0), "Bottom")
+    stokes.add_dirichlet_bc((0.0, 0.0), "Left")
+    stokes.add_dirichlet_bc((0.0, 0.0), "Right")
+    stokes.petsc_use_pressure_nullspace = True
+    stokes.tolerance = 1.0e-8
+
+    stokes.set_pressure_gauge("Top")
+    stokes.solve()
+    assert stokes.snes.getConvergedReason() > 0
+
+    # The gauge held: the non-dimensional mean pressure on Top is zero.
+    area = uw.maths.BdIntegral(mesh, 1.0, "Top").evaluate()
+    mean_top = uw.maths.BdIntegral(mesh, p.sym[0, 0], "Top").evaluate() / area
+    if isinstance(mean_top, uw.function.quantities.UWQuantity):
+        mean_top = mean_top.data
+    assert np.isclose(float(mean_top), 0.0, atol=1.0e-6), f"top-mean pressure = {mean_top}"