From a8a393b4068b9cdac4e3b28618bf260259da43cf Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:07:03 +0100
Subject: [PATCH 01/11] refactor(api): WC-01/WC-02 value-first BC order and
 conds datum for the newer BC methods

Migrate add_nitsche_bc (SNES_Vector and Stokes), add_rotated_freeslip_bc and
add_constraint_bc to the canonical value-first signature (conds, boundary, ...)
per maintainer decisions D2 (2026-07-04) and D3 (conds, 2026-07-06). Legacy
boundary-first calls and the g= spelling keep working through a conservative
shim (first positional is a string while the second is not) that emits exactly
one DeprecationWarning; the shared detection lives in
SolverBaseClass._value_first_bc_args. The legacy trio already conforms and is
untouched apart from bringing the components= deprecation warning up to the
Charter form (stacklevel, replacement named). The Charter's pending BC-value-
parameter placeholder is replaced with the confirmed conds ruling.

Contract tests land in tests/test_0641_wave_c_api_shims.py at the end of the
Wave C series.

Underworld development team with AI support from Claude Code
---
 docs/developer/UW3_STYLE_CHARTER.md           |   4 +-
 .../cython/petsc_generic_snes_solvers.pyx     | 154 +++++++++++++++---
 src/underworld3/systems/solvers.py            |  32 +++-
 3 files changed, 159 insertions(+), 31 deletions(-)

diff --git a/docs/developer/UW3_STYLE_CHARTER.md b/docs/developer/UW3_STYLE_CHARTER.md
index 65996e05..f607d520 100644
--- a/docs/developer/UW3_STYLE_CHARTER.md
+++ b/docs/developer/UW3_STYLE_CHARTER.md
@@ -82,8 +82,8 @@ These settle the June-2026 drift (see `docs/reviews/2026-07/API-CONSISTENCY-REVI
 
 | Topic | Rule |
 |---|---|
-| BC argument order | `add_<kind>_bc(value, boundary, ...)` — value first, matching the original trio, the published examples/benchmarks, and previous major versions (maintainer decision 2026-07-04). The newer boundary-first methods (`add_nitsche_bc`, `add_rotated_freeslip_bc`, `add_constraint_bc`) migrate to it; their current signatures are shimmed. |
-| BC value parameter | ONE name across all BC methods ⟨pending maintainer confirmation: `conds` (implied by the original-order decision) vs `value`⟩. Until decided, do not introduce a third spelling. |
+| BC argument order | `add_<kind>_bc(conds, boundary, ...)` — value first, matching the original trio, the published examples/benchmarks, and previous major versions (maintainer decision 2026-07-04). The newer boundary-first methods (`add_nitsche_bc`, `add_rotated_freeslip_bc`, `add_constraint_bc`) migrate to it; their current signatures are shimmed. |
+| BC value parameter | ONE name across all BC methods: **`conds`** (maintainer decision 2026-07-06, recorded in `docs/reviews/2026-07/REMEDIATION-WORKLIST.md`). Other spellings (`g`, `value`, ...) exist only as deprecated keyword aliases where they existed before the decision; never introduce a new spelling. |
 | Direction selection | Component masking via `None`/`sympy.oo` entries in the value vector; `direction=` only for a scalar constraint along a vector (defaults to outward normal); `normal=` strictly overrides the geometric surface-normal source. `components=` is deprecated — never in new code. |
 | Solver capabilities | Anything that configures or reads one solver is a METHOD on that solver, lazily importing its `utilities/*` implementation (the `boundary_flux` pattern) — never a free function as the documented entry point. |
 | Namespaces | Every user-facing module is exported from its subpackage `__init__`/`__all__` in the PR that creates it. No deep-import-only features. |
diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 048b52a6..e131279e 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -1447,12 +1447,12 @@ class SolverBaseClass(uw_object):
                   "sympy.Matrix, i.e. conds = sympy.Matrix([sympy.oo, 5, 1.2])\n")
 
         if isinstance(components, (tuple, list, int)):
-            # TODO: DECPRECATE
             import warnings
-            warnings.warn(category=DeprecationWarning,
-                          message="Using the 'components' argument is being DEPRECATED in the next release\n" +
-                                  "The same functionality can be setup with the 'conds' argument and using\n" +
-                                  "'sympy.oo' or 'None', see docstring")
+            warnings.warn(
+                "The 'components' argument is deprecated; select components "
+                "with None / sympy.oo entries in 'conds' instead, e.g. "
+                "conds=(None, 5, 1.2)",
+                DeprecationWarning, stacklevel=3)
             components = np.array(components, dtype=np.int32, ndmin=1)
 
         elif components is None:
@@ -1525,6 +1525,49 @@ class SolverBaseClass(uw_object):
             self.essential_bcs.append(BC(f_id, components,sympy_fn, label, -1,  'essential', -1))
 
 
+    def _value_first_bc_args(self, method, conds, boundary, alias=None,
+                             alias_name="g"):
+        """Normalize BC arguments to the canonical value-first order.
+
+        The canonical BC signature is ``method(conds, boundary, ...)`` with the
+        prescribed datum named ``conds`` (Style Charter, API conventions;
+        maintainer decisions D2/D3, 2026-07). Two legacy spellings are shimmed
+        here, each with exactly one DeprecationWarning:
+
+        * **boundary-first order** — detected conservatively: the first
+          positional argument is a string (a boundary label; a BC datum is
+          never a string — see :meth:`add_condition`) while the second is not.
+          The two arguments are swapped.
+        * **the** ``g=`` **keyword alias** for the datum — forwarded to
+          ``conds``. Supplying both ``conds`` and ``g`` is an error.
+
+        Returns the normalized ``(conds, boundary)`` pair; ``boundary`` is
+        required to be a string on exit.
+        """
+        legacy = False
+        if isinstance(conds, str) and not isinstance(boundary, str):
+            conds, boundary = boundary, conds
+            legacy = True
+        if alias is not None:
+            if conds is not None:
+                raise TypeError(
+                    f"{method}() received the boundary datum twice "
+                    f"(positionally and as '{alias_name}='); pass it once, "
+                    f"as 'conds'")
+            conds = alias
+            legacy = True
+        if legacy:
+            import warnings
+            warnings.warn(
+                f"{method}(boundary, {alias_name}=...) is deprecated; "
+                f"use {method}(conds, boundary, ...)",
+                DeprecationWarning, stacklevel=3)
+        if not isinstance(boundary, str):
+            raise TypeError(
+                f"{method}() requires a boundary label string; "
+                f"got {type(boundary).__name__}")
+        return conds, boundary
+
     # Use FE terminology note f_id is 0.
     @timing.routine_timer_decorator
     def add_essential_bc(self, conds, boundary, components=None):
@@ -3314,19 +3357,21 @@ class SNES_Vector(SolverBaseClass):
         self.petsc_options["ksp_atol"]  = self._tolerance * 1.0e-6
 
 
-    def add_nitsche_bc(self, boundary, g=None, direction=None, gamma=10.0, theta=1, local_h=True):
+    def add_nitsche_bc(self, conds=None, boundary=None, direction=None,
+                       gamma=10.0, theta=1, local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
         For vector solvers (no pressure field), this constrains
-        :math:`\mathbf{u} \cdot \mathbf{d} = g` on the boundary using
-        Nitsche's method with penalty, consistency, and symmetry terms.
+        :math:`\mathbf{u} \cdot \mathbf{d} = \mathrm{conds}` on the boundary
+        using Nitsche's method with penalty, consistency, and symmetry terms.
 
         Parameters
         ----------
+        conds : sympy expression or float, optional
+            Prescribed velocity along the constraint direction. Default zero
+            (free-slip when the direction is the surface normal).
         boundary : str
             Boundary label.
-        g : sympy expression or float, optional
-            Prescribed velocity along constraint direction. Default zero.
         direction : sympy.Matrix or list, optional
             Constraint direction. Default ``None`` uses surface normal.
         gamma : float, default=10.0
@@ -3338,12 +3383,21 @@ class SNES_Vector(SolverBaseClass):
             (:meth:`Mesh.cell_size`) rather than the global minimum
             (:meth:`Mesh.get_min_radius`). See
             ``SNES_Stokes_SaddlePt.add_nitsche_bc`` for details.
+        g : sympy expression or float, optional
+            Deprecated keyword alias for ``conds`` (one DeprecationWarning).
 
         Warnings
         --------
         Exterior boundaries only. See ``SNES_Stokes_SaddlePt.add_nitsche_bc``
         for details on why internal boundaries are not supported.
 
+        Notes
+        -----
+        The legacy boundary-first call ``add_nitsche_bc(boundary, g=...)`` is
+        detected conservatively (first positional argument a string while the
+        second is not — a BC datum is never a string) and shimmed with one
+        DeprecationWarning; see :meth:`_value_first_bc_args`.
+
         See Also
         --------
         SNES_Stokes_SaddlePt.add_nitsche_bc : Stokes version with pressure coupling.
@@ -3351,6 +3405,10 @@ class SNES_Vector(SolverBaseClass):
         import sympy
         from collections import namedtuple
 
+        conds, boundary = self._value_first_bc_args(
+            "add_nitsche_bc", conds, boundary, alias=g)
+        g = conds
+
         self.is_setup = False
 
         mesh = self.mesh
@@ -5167,7 +5225,7 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
     #     BC = namedtuple('EssentialBC', ['components', 'fn', 'boundary', 'boundary_label_val', 'type', 'PETScID'])
     #     self.essential_p_bcs.append(BC(components, sympy_fn, boundary, -1,  'essential', -1))
 
-    def add_rotated_freeslip_bc(self, boundary, normal=None):
+    def add_rotated_freeslip_bc(self, conds=None, boundary=None, normal=None):
         r"""Add STRONG free-slip (:math:`\mathbf{u}\cdot\hat{\mathbf n}=0`) by rotating
         the boundary velocity DOFs into a per-node (normal, tangential) frame and
         imposing the rotated normal component as an exact Dirichlet constraint.
@@ -5181,6 +5239,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
 
         Parameters
         ----------
+        conds : float or None, optional
+            Prescribed wall-normal velocity datum, in the canonical value-first
+            BC order (Style Charter, API conventions). Only the homogeneous
+            free-slip constraint :math:`\mathbf{u}\cdot\hat{\mathbf n}=0` is
+            implemented, so this must be zero (or ``None``, meaning zero); a
+            non-zero datum raises ``NotImplementedError`` (use
+            :meth:`add_nitsche_bc` or ``add_constraint_bc`` for prescribed
+            normal in/outflow).
         boundary : str
             Boundary label to constrain.
         normal : None or sympy 1×dim Matrix or array, optional
@@ -5197,7 +5263,38 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         face frees two tangential directions, a 3D edge frees one (the edge
         tangent), a corner is fully pinned. Registering delegates the solve to
         :mod:`underworld3.utilities.rotated_bc`.
+
+        The legacy boundary-first call ``add_rotated_freeslip_bc(boundary,
+        normal)`` is detected conservatively (the first positional argument is
+        a string — the datum is never a string) and shimmed with one
+        DeprecationWarning: the string becomes ``boundary`` and a second
+        positional argument, if present, becomes ``normal``.
         """
+        if isinstance(conds, str):
+            # legacy boundary-first call: (boundary[, normal])
+            if boundary is not None:
+                if normal is not None:
+                    raise TypeError(
+                        "add_rotated_freeslip_bc() received 'normal' twice "
+                        "(positionally, legacy order, and as a keyword)")
+                normal = boundary
+            boundary = conds
+            conds = None
+            import warnings
+            warnings.warn(
+                "add_rotated_freeslip_bc(boundary, normal) is deprecated; "
+                "use add_rotated_freeslip_bc(conds, boundary, normal=...) "
+                "with conds=0 (or boundary=... by keyword)",
+                DeprecationWarning, stacklevel=2)
+        if not isinstance(boundary, str):
+            raise TypeError(
+                f"add_rotated_freeslip_bc() requires a boundary label string; "
+                f"got {type(boundary).__name__}")
+        if conds is not None and sympy.sympify(conds) != 0:
+            raise NotImplementedError(
+                "add_rotated_freeslip_bc imposes u.n = 0 only; a non-zero "
+                "wall-normal datum is not implemented (use add_nitsche_bc or "
+                "add_constraint_bc for prescribed normal in/outflow)")
         self._rotated_freeslip_bcs.append((boundary, normal))
         self.is_setup = False
         return
@@ -5290,7 +5387,8 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         return _dtf(self, boundary, self._rotated_freeslip_info, field,
                     buoyancy_scale=buoyancy_scale, mass=mass)
 
-    def add_nitsche_bc(self, boundary, g=None, direction=None, normal=None, gamma=10.0, theta=1, mask=None, local_h=True):
+    def add_nitsche_bc(self, conds=None, boundary=None, direction=None, normal=None,
+                       gamma=10.0, theta=1, mask=None, local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
         Nitsche's method provides a variationally consistent alternative to
@@ -5298,9 +5396,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         and gives optimal convergence rates.
 
         By default, constrains the normal velocity component
-        :math:`\mathbf{u} \cdot \mathbf{n} = g` (free-slip when *g* = 0).
-        When *direction* is provided, constrains
-        :math:`\mathbf{u} \cdot \mathbf{d} = g` along that direction instead.
+        :math:`\mathbf{u} \cdot \mathbf{n} = \mathrm{conds}` (free-slip when
+        ``conds`` is zero). When *direction* is provided, constrains
+        :math:`\mathbf{u} \cdot \mathbf{d} = \mathrm{conds}` along that
+        direction instead.
 
         The method constructs boundary residuals and Jacobians for:
 
@@ -5313,11 +5412,11 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
 
         Parameters
         ----------
-        boundary : str
-            Boundary label (e.g., ``"Upper"``, ``"Lower"``).
-        g : sympy expression or float, optional
+        conds : sympy expression or float, optional
             Prescribed velocity along the constraint direction. Default
             ``None`` means zero (:math:`\mathbf{u} \cdot \mathbf{d} = 0`).
+        boundary : str
+            Boundary label (e.g., ``"Upper"``, ``"Lower"``).
         direction : sympy.Matrix or list, optional
             Constraint direction vector. Default ``None`` uses the boundary
             surface normal (free-slip). Can be spatially varying (e.g.,
@@ -5347,18 +5446,27 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             adaptive mesh the local size scales the stabilisation correctly
             on every facet; on a uniform mesh the two coincide. Set ``False``
             to restore the legacy global-h behaviour exactly.
+        g : sympy expression or float, optional
+            Deprecated keyword alias for ``conds`` (one DeprecationWarning).
 
         Examples
         --------
         >>> # Free-slip (u.n = 0)
-        >>> stokes.add_nitsche_bc("Upper", gamma=10)
+        >>> stokes.add_nitsche_bc(0.0, "Upper", gamma=10)
 
         >>> # Prescribed normal inflow
-        >>> stokes.add_nitsche_bc("Left", g=1.0, gamma=10)
+        >>> stokes.add_nitsche_bc(1.0, "Left", gamma=10)
 
         >>> # Constrain along a specific direction (e.g. fault normal)
         >>> fault_normal = sympy.Matrix([0.6, 0.8])
-        >>> stokes.add_nitsche_bc("Fault", direction=fault_normal, gamma=10)
+        >>> stokes.add_nitsche_bc(0.0, "Fault", direction=fault_normal, gamma=10)
+
+        Notes
+        -----
+        The legacy boundary-first call ``add_nitsche_bc(boundary, g=...)`` is
+        detected conservatively (first positional argument a string while the
+        second is not — a BC datum is never a string) and shimmed with one
+        DeprecationWarning; see :meth:`_value_first_bc_args`.
 
         Warnings
         --------
@@ -5376,6 +5484,10 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         import sympy
         from collections import namedtuple
 
+        conds, boundary = self._value_first_bc_args(
+            "add_nitsche_bc", conds, boundary, alias=g)
+        g = conds
+
         self.is_setup = False
 
         mesh = self.mesh
diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 2b098464..89245b97 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -2332,22 +2332,25 @@ def _viscosity_scale(self):
         except (TypeError, ValueError, AttributeError):
             return 1.0
 
-    def add_constraint_bc(self, boundary, g=0.0, normal=None, screening=None,
-                          augmentation=None, augmentation_base=1.0e4, degree=None):
+    def add_constraint_bc(self, conds=None, boundary=None, normal=None, screening=None,
+                          augmentation=None, augmentation_base=1.0e4, degree=None,
+                          g=None):
         r"""Register a multiplier-enforced normal-velocity constraint on ``boundary``.
 
         Adds a scalar multiplier field ``h`` coupled into the saddle-point system
-        so that :math:`\mathbf{u}\cdot\mathbf{n}=g` 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`.
+        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`.
 
         Parameters
         ----------
+        conds : float or sympy expression, optional
+            Prescribed normal velocity :math:`\mathbf{u}\cdot\mathbf{n}`,
+            in the canonical value-first BC order. Default ``None`` means zero
+            (free-slip).
         boundary : str
             Mesh boundary label (e.g. ``"Upper"``).
-        g : float or sympy expression, default 0.0
-            Prescribed normal velocity :math:`\mathbf{u}\cdot\mathbf{n} = g`.
         normal : sympy matrix, optional
             Row-vector constraint normal. Defaults to ``mesh.Gamma_P1``.
         screening : float or sympy expression, optional
@@ -2366,12 +2369,25 @@ def add_constraint_bc(self, boundary, g=0.0, normal=None, screening=None,
             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.
+        g : float or sympy expression, optional
+            Deprecated keyword alias for ``conds`` (one DeprecationWarning).
 
         Returns
         -------
         h : MeshVariable
             The scalar multiplier field.
+
+        Notes
+        -----
+        The legacy boundary-first call ``add_constraint_bc(boundary, g=...)``
+        is detected conservatively (first positional argument a string while
+        the second is not — a BC datum is never a string) and shimmed with one
+        DeprecationWarning; see
+        ``SolverBaseClass._value_first_bc_args``.
         """
+        conds, boundary = self._value_first_bc_args(
+            "add_constraint_bc", conds, boundary, alias=g)
+        g = conds if conds is not None else 0.0
         # Parallel-safe: the interior-multiplier reduction
         # (_constrain_interior_multipliers_in_section) is rank-local section
         # surgery (it uses the distributed boundary label IS and iterates the

From 99dc4147b32d0d67c182b97e9c503f8e30800198 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:07:49 +0100
Subject: [PATCH 02/11] refactor(api): WC-08 align SNES_Vector.add_nitsche_bc
 with the Stokes variant

Accept normal= with the same geometric-normal-override semantics as
SNES_Stokes_SaddlePt.add_nitsche_bc, so code written against the Stokes
variant no longer raises TypeError on a vector solver. Accept mask= for
signature parity but raise a clear NotImplementedError (one-sided masking is
implemented on the Stokes variant only). Signature and docs only; the default
path (normal=None, mask=None) is unchanged.

Underworld development team with AI support from Claude Code
---
 .../cython/petsc_generic_snes_solvers.pyx     | 31 ++++++++++++++++---
 1 file changed, 27 insertions(+), 4 deletions(-)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index e131279e..79bdf80d 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -3358,7 +3358,8 @@ class SNES_Vector(SolverBaseClass):
 
 
     def add_nitsche_bc(self, conds=None, boundary=None, direction=None,
-                       gamma=10.0, theta=1, local_h=True, g=None):
+                       normal=None, gamma=10.0, theta=1, mask=None,
+                       local_h=True, g=None):
         r"""Add Nitsche weak enforcement of a velocity constraint along a direction.
 
         For vector solvers (no pressure field), this constrains
@@ -3374,10 +3375,19 @@ class SNES_Vector(SolverBaseClass):
             Boundary label.
         direction : sympy.Matrix or list, optional
             Constraint direction. Default ``None`` uses surface normal.
+        normal : sympy.Matrix or list, optional
+            Boundary unit normal used in the Nitsche consistency and symmetry
+            terms — the same geometric-normal override as on the Stokes
+            variant. Default ``None`` uses the per-boundary,
+            deformation-tracking ``mesh.boundary_normal(boundary)``.
         gamma : float, default=10.0
             Dimensionless stabilisation parameter.
         theta : {-1, 0, 1}, default=1
             Symmetry parameter (1=symmetric, -1=skew-symmetric).
+        mask : sympy expression, optional
+            Accepted for signature parity with the Stokes variant, but
+            one-sided masking is **not implemented** on vector solvers:
+            passing a mask raises ``NotImplementedError``.
         local_h : bool, default=True
             Scale the penalty by a local per-cell mesh size
             (:meth:`Mesh.cell_size`) rather than the global minimum
@@ -3409,6 +3419,12 @@ class SNES_Vector(SolverBaseClass):
             "add_nitsche_bc", conds, boundary, alias=g)
         g = conds
 
+        if mask is not None:
+            raise NotImplementedError(
+                "mask= (one-sided internal-boundary application) is not "
+                "implemented on SNES_Vector.add_nitsche_bc; it is supported "
+                "on SNES_Stokes_SaddlePt.add_nitsche_bc only")
+
         self.is_setup = False
 
         mesh = self.mesh
@@ -3419,10 +3435,17 @@ class SNES_Vector(SolverBaseClass):
         dim = mesh.cdim
 
         # Surface normal components — use this boundary's own deformation-
-        # tracking facet normal (see Mesh.boundary_normal); the legacy global
+        # tracking facet normal (see Mesh.boundary_normal) unless the caller
+        # overrides the geometric-normal source; the legacy global
         # mesh.Gamma_P1 stays radial on a deformed surface.
-        bnorm = mesh.boundary_normal(boundary)
-        n = [bnorm[i] for i in range(dim)]
+        if normal is not None:
+            if isinstance(normal, sympy.MatrixBase):
+                n = [normal[i] for i in range(dim)]
+            else:
+                n = list(normal)
+        else:
+            bnorm = mesh.boundary_normal(boundary)
+            n = [bnorm[i] for i in range(dim)]
 
         # Constraint direction: defaults to surface normal
         if direction is not None:

From c6584639b31ea18aa8e0508342f404adcd4c04c3 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:08:45 +0100
Subject: [PATCH 03/11] fix(api): WC-03 consistent_jacobian becomes a validated
 property

The tri-state mode switch was a bare attribute: any unrecognized truthy value
('picard', 'Continuation', 1) silently selected the full-Newton tangent in
_jacobian_source, and a typo silently skipped the continuation solve dispatch.
The setter now accepts exactly {False, True, 'continuation'}, normalizes falsy
values (None, 0) to False (they already meant Picard), and raises ValueError
otherwise. The explanatory comment block moves into the property's NumPy
docstring. Behaviour is bit-identical for the three valid values.

Underworld development team with AI support from Claude Code
---
 .../cython/petsc_generic_snes_solvers.pyx     | 77 ++++++++++++++-----
 1 file changed, 56 insertions(+), 21 deletions(-)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 79bdf80d..9ab955f1 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -68,27 +68,8 @@ class SolverBaseClass(uw_object):
         self.compiled_extensions = None
         self.constants_manifest = []
 
-        # Jacobian tangent selection: False | True | "continuation".
-        #
-        #   False (default): differentiate the residual flux *as wrapped* — the
-        #     effective viscosity is frozen, giving a Picard / defect-correction
-        #     tangent. BIT-IDENTICAL to the long-standing behaviour. Globally
-        #     robust; load-bearing for the tuned hard-yield viscoplastic paths.
-        #   True: unwrap the flux before differentiation so the tangent captures
-        #     d(eta)/d(grad v) (full Newton). Fast near the solution; its yield
-        #     kink can stall the line search far from it.
-        #   "continuation": Picard -> Newton. Blend J(alpha) = J_picard +
-        #     alpha*(J_newton - J_picard) with alpha a constants[] parameter
-        #     ramped 0 -> 1 by a SNES monitor as the residual drops. Picard
-        #     locates the basin, Newton gives quadratic convergence inside it
-        #     (cf. Spiegelman et al. 2016; ASPECT defect-correction-then-Newton).
-        #     alpha=0 is bit-identical to Picard, so no recompile to switch.
-        #
-        # The Newton flux for a model whose flux has a non-smooth yield kink is
-        # the model's own smooth law (constitutive_model.flux_jacobian) when it
-        # provides one; otherwise the exact unwrapped flux.
-        #
-        # See docs/developer/design/jacobian-unwrap-constants-bug.md.
+        # Jacobian tangent selection — validated property, see the
+        # consistent_jacobian docstring below for the mode semantics.
         self.consistent_jacobian = False
         # Picard->Newton continuation parameter (constants[]-routed so it can be
         # ramped at solve time without a JIT recompile). 0 = Picard, 1 = Newton.
@@ -157,6 +138,60 @@ class SolverBaseClass(uw_object):
         # utilities.custom_mg). None => standard FMG/GAMG path, unchanged.
         self._custom_mg = None
 
+    @property
+    def consistent_jacobian(self):
+        r"""Jacobian tangent selection: ``False`` | ``True`` | ``"continuation"``.
+
+        Selects the tangent used by :meth:`_jacobian_source` and the solve
+        dispatch; the residual is never affected, so the converged solution
+        always satisfies the exact constitutive law.
+
+        ``False`` (default)
+            Differentiate the residual flux *as wrapped* — the effective
+            viscosity is frozen, giving a Picard / defect-correction tangent.
+            Bit-identical to the long-standing behaviour. Globally robust;
+            load-bearing for the tuned hard-yield viscoplastic paths.
+        ``True``
+            Unwrap the flux before differentiation so the tangent captures
+            :math:`\partial\eta/\partial(\nabla v)` (full Newton). Fast near
+            the solution; its yield kink can stall the line search far from it.
+        ``"continuation"``
+            Picard :math:`\rightarrow` Newton. Blend
+            :math:`J(\alpha) = J_{\mathrm{picard}} + \alpha\,(J_{\mathrm{newton}}
+            - J_{\mathrm{picard}})` with :math:`\alpha` a ``constants[]``
+            parameter ramped :math:`0 \rightarrow 1` by a SNES monitor as the
+            residual drops. Picard locates the basin, Newton gives quadratic
+            convergence inside it (cf. Spiegelman et al. 2016; ASPECT
+            defect-correction-then-Newton). :math:`\alpha = 0` is bit-identical
+            to Picard, so no recompile is needed to switch.
+
+        The Newton flux for a model whose flux has a non-smooth yield kink is
+        the model's own smooth law (``constitutive_model.flux_jacobian``) when
+        it provides one; otherwise the exact unwrapped flux. See
+        ``docs/developer/design/jacobian-unwrap-constants-bug.md``.
+
+        Raises
+        ------
+        ValueError
+            On assignment of anything other than ``False``, ``True`` or
+            ``"continuation"``. Falsy values (``None``, ``0``, ``""``)
+            normalize to ``False`` (they already selected the Picard tangent).
+            Before validation, any other truthy value (``1``, ``"picard"``,
+            ``"Continuation"``) silently selected the full-Newton tangent.
+        """
+        return self._consistent_jacobian
+
+    @consistent_jacobian.setter
+    def consistent_jacobian(self, mode):
+        if not mode:
+            self._consistent_jacobian = False
+        elif mode is True or mode == "continuation":
+            self._consistent_jacobian = mode
+        else:
+            raise ValueError(
+                f"consistent_jacobian must be False, True or 'continuation'; "
+                f"got {mode!r}")
+
     def _jacobian_source(self, expr, newton_expr=None):
         """Prepare a residual flux for Jacobian differentiation.
 

From 2f5af80d46ef0638ae090014f982c2e3965b0109 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:09:30 +0100
Subject: [PATCH 04/11] refactor(api): WC-04 SolverBaseClass.set_custom_fmg
 method; deprecate set_custom_mg

The discoverable solver method for custom multigrid was the legacy serial-
only, finest-only-reduction, single-field path (set_custom_mg) while the
correct parallel-capable BC-per-level path (set_custom_fmg) was reachable only
by deep import. Add a thin set_custom_fmg method delegating via lazy import to
utilities.custom_mg (the boundary_flux pattern) and emit a DeprecationWarning
from set_custom_mg naming the replacement. No change to custom_mg.py itself
(coordination with the active bugfix/custom-mg-parallel branch).

Underworld development team with AI support from Claude Code
---
 .../cython/petsc_generic_snes_solvers.pyx     | 47 +++++++++++++++++++
 1 file changed, 47 insertions(+)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 9ab955f1..615784d6 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -324,7 +324,19 @@ class SolverBaseClass(uw_object):
         ``setFromOptions`` / nullspace attach) via
         :func:`underworld3.utilities.custom_mg.inject_custom_mg`. See
         :mod:`underworld3.utilities.custom_mg`.
+
+        .. deprecated:: 2026-07
+            This is the legacy **serial-only, finest-only-reduction,
+            single-field** path; the Stokes velocity-block support described
+            above is delivered by :meth:`set_custom_fmg`, which is the
+            canonical entry point (parallel-capable, BC-per-level reduction).
         """
+        import warnings
+        warnings.warn(
+            "set_custom_mg is deprecated (legacy serial-only, single-field "
+            "custom-MG path); use set_custom_fmg(coarse_meshes, "
+            "builder=..., field_id=...) instead",
+            DeprecationWarning, stacklevel=2)
         if kind not in ("barycentric", "rbf"):
             raise ValueError("kind must be 'barycentric' or 'rbf'")
         if not coarse_meshes:
@@ -333,6 +345,41 @@ class SolverBaseClass(uw_object):
                            "kind": kind, "verbose": verbose}
         self.is_setup = False
 
+    def set_custom_fmg(self, coarse_meshes, *, builder="barycentric",
+                       field_id=None, verbose=False):
+        r"""Drive geometric multigrid with a prolongation built from
+        ``coarse_meshes`` — the canonical custom-MG entry point.
+
+        Registers a BC-per-level reduced hierarchy on the solver so the next
+        :meth:`solve` builds and installs it (build-time injection). Works in
+        parallel and on the **Stokes velocity block** (pass ``field_id=0`` on a
+        saddle-point solver). This supersedes the legacy
+        :meth:`set_custom_mg` path (serial-only, finest-only reduction,
+        single-field).
+
+        Parameters
+        ----------
+        coarse_meshes : list of Mesh
+            Coarsest-first list of coarse meshes (the finest level is the
+            solver's own mesh). They need only carry the same boundary labels
+            as the solver's mesh.
+        builder : {"barycentric", "rbf"}, default "barycentric"
+            Prolongation builder. ``barycentric`` is FE-exact; ``rbf`` is a
+            polyharmonic RBF (Shepard-normalised).
+        field_id : int, optional
+            Target field for a multi-field (saddle-point) solver; ``0`` is the
+            Stokes velocity block. ``None`` (default) for single-field solvers.
+        verbose : bool, default False
+            Print the per-level DOF counts at injection.
+
+        See Also
+        --------
+        underworld3.utilities.custom_mg.set_custom_fmg : the implementation.
+        """
+        from underworld3.utilities.custom_mg import set_custom_fmg as _set_custom_fmg
+        _set_custom_fmg(self, coarse_meshes, builder=builder,
+                        field_id=field_id, verbose=verbose)
+
     def add_update_callback(self, callback):
         r"""Register a callback fired at the start of every nonlinear (SNES) iteration.
 

From e8bd45021bcfb153bda20b1e1211cde0db4175f8 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:10:59 +0100
Subject: [PATCH 05/11] feat(api): WC-05 export
 rotated_bc/boundary_flux/custom_mg and BoundingSurface from their packages

utilities/__init__ now imports the three solver-capability implementation
modules (and re-exports set_custom_fmg), so docstring cross-references resolve
without deep imports. meshing/__init__ imports and lists BoundingSurface and
the three register_*_surfaces helpers, which Mesh.register_tangent_slip_provider
requires but which previously raised AttributeError as uw.meshing attributes.
Pure additions (Charter API conventions: every user-facing module is exported
from its subpackage __init__/__all__).

Underworld development team with AI support from Claude Code
---
 src/underworld3/meshing/__init__.py   | 12 ++++++++++++
 src/underworld3/utilities/__init__.py | 10 ++++++++++
 2 files changed, 22 insertions(+)

diff --git a/src/underworld3/meshing/__init__.py b/src/underworld3/meshing/__init__.py
index 3bbd541c..896e8481 100644
--- a/src/underworld3/meshing/__init__.py
+++ b/src/underworld3/meshing/__init__.py
@@ -64,6 +64,13 @@
     ADAPT_STRATEGIES,
 )
 
+from .bounding_surface import (
+    BoundingSurface,
+    register_radial_surfaces,
+    register_plane_surfaces,
+    register_box_face_surfaces,
+)
+
 # Make all functions available at module level for backward compatibility
 __all__ = [
     # Cartesian meshes
@@ -108,4 +115,9 @@
     "mesh_metric_mismatch",
     "follow_metric",
     "ADAPT_STRATEGIES",
+    # Bounding surfaces (tangent-slip providers for deforming boundaries)
+    "BoundingSurface",
+    "register_radial_surfaces",
+    "register_plane_surfaces",
+    "register_box_face_surfaces",
 ]
diff --git a/src/underworld3/utilities/__init__.py b/src/underworld3/utilities/__init__.py
index 1731805b..0d7c627e 100644
--- a/src/underworld3/utilities/__init__.py
+++ b/src/underworld3/utilities/__init__.py
@@ -84,3 +84,13 @@ def _append_petsc_path():
 
 from . import retention_curves
 from . import memprobe
+
+# Solver-scoped capability modules (rotated strong free-slip, consistent
+# boundary-flux recovery, custom multigrid). Their documented entry points are
+# methods on the solvers, which lazy-import these implementations; exporting
+# the modules here makes the docstring cross-references resolvable without a
+# deep import (Style Charter, API conventions: namespaces).
+from . import rotated_bc
+from . import boundary_flux
+from . import custom_mg
+from .custom_mg import set_custom_fmg

From 8e400eaa3c72d06f026e1b48486be2de0c852a1c Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:12:33 +0100
Subject: [PATCH 06/11] refactor(api): WC-06 uw.quantity is THE quantity
 factory (D14)

uw.quantity (returns UWQuantity) and uw.create_quantity (returns a raw Pint
Quantity) duplicated the same concept with silently different return types.
Per maintainer decision D14: create_quantity keeps its exact behaviour for one
deprecation cycle but emits a DeprecationWarning naming uw.quantity; UWQuantity
is exposed at top level for isinstance checks. test_0640, which pinned
create_quantity as public, is updated in the same commit.

Underworld development team with AI support from Claude Code
---
 src/underworld3/__init__.py                   |  4 +++-
 src/underworld3/units.py                      | 13 +++++++++++-
 tests/test_0640_api_consistency_regression.py | 21 +++++++++++++++----
 3 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/src/underworld3/__init__.py b/src/underworld3/__init__.py
index 4a316cc7..801e1d86 100644
--- a/src/underworld3/__init__.py
+++ b/src/underworld3/__init__.py
@@ -219,7 +219,9 @@ def view():
 from .parameters import ParameterRegistry, ParameterType
 from .materials import MaterialRegistry, MaterialProperty
 from .constitutive_models import MultiMaterialConstitutiveModel
-from .function import quantity, expression, with_units, expand, unwrap
+# uw.quantity is THE quantity factory (returns UWQuantity, exposed alongside
+# for isinstance checks); uw.create_quantity is deprecated (see units.py).
+from .function import quantity, UWQuantity, expression, with_units, expand, unwrap
 from .coordinates import uwdiff  # Differentiation helper for UWCoordinates
 from .utilities import retention_curves
 
diff --git a/src/underworld3/units.py b/src/underworld3/units.py
index c837e5a5..28352ddb 100644
--- a/src/underworld3/units.py
+++ b/src/underworld3/units.py
@@ -1122,18 +1122,29 @@ def create_quantity(value, units: Union[str, Any], backend: Optional[str] = None
     """
     Create a dimensional quantity from a value and units.
 
+    .. deprecated:: 2026-07
+        ``uw.quantity`` is THE quantity factory (maintainer decision D14):
+        it returns a :class:`~underworld3.function.quantities.UWQuantity`,
+        which composes with expressions and non-dimensionalisation. This
+        function returns a raw Pint ``Quantity`` and will be removed after one
+        deprecation cycle.
+
     Args:
         value: Numeric value or array
         units: Units specification (string or units object)
         backend: Backend to use ('pint' or 'sympy'), defaults to 'pint'
 
     Returns:
-        Dimensional quantity
+        Dimensional quantity (raw Pint ``Quantity``)
 
     Examples:
         >>> velocity_qty = create_quantity([1.0, 2.0], "m/s")
         >>> pressure_qty = create_quantity(101325, "Pa")
     """
+    warnings.warn(
+        "create_quantity is deprecated; use uw.quantity(value, units) — note "
+        "it returns a UWQuantity rather than a raw Pint Quantity",
+        DeprecationWarning, stacklevel=2)
     # backend parameter is deprecated - Pint is the only supported backend
     if backend is not None and backend.lower() != "pint":
         raise ValueError(f"Unknown backend: {backend}. Only 'pint' is supported.")
diff --git a/tests/test_0640_api_consistency_regression.py b/tests/test_0640_api_consistency_regression.py
index 934c3dd1..2c7438d4 100644
--- a/tests/test_0640_api_consistency_regression.py
+++ b/tests/test_0640_api_consistency_regression.py
@@ -137,7 +137,7 @@ def test_create_functions_exist(self):
 
         factory_functions = [
             "create_enhanced_mesh_variable",
-            "create_quantity",  # From units system
+            "create_quantity",  # deprecated alias for one cycle (D14); see uw.quantity
         ]
 
         for func_name in factory_functions:
@@ -146,6 +146,12 @@ def test_create_functions_exist(self):
             func = getattr(uw, func_name)
             assert callable(func), f"{func_name} is not callable"
 
+    def test_quantity_is_the_factory(self):
+        """uw.quantity is THE quantity factory (maintainer decision D14)."""
+
+        q = uw.quantity(1.0, "m")
+        assert isinstance(q, uw.UWQuantity)  # type exposed at top level
+
     def test_factory_function_documentation(self):
         """Test that factory functions have proper documentation."""
 
@@ -249,14 +255,19 @@ class TestUnitsSystemAPIConsistency:
     def test_units_creation_patterns(self):
         """Test consistent patterns for creating unit-aware objects."""
 
-        # Pattern 1: Direct units specification
+        # Pattern 1: Direct units specification (THE factory, decision D14)
         try:
-            quantity1 = uw.create_quantity(10.0, "m/s")
+            quantity1 = uw.quantity(10.0, "m/s")
             assert quantity1 is not None
 
         except Exception as e:
             pytest.fail(f"Direct units creation failed: {e}")
 
+        # Deprecated spelling still works for one cycle, with a warning
+        with pytest.warns(DeprecationWarning, match="uw.quantity"):
+            legacy = uw.create_quantity(10.0, "m/s")
+        assert legacy is not None
+
         # Pattern 2: uw.units multiplication
         try:
             quantity2 = 10.0 * uw.units.m / uw.units.s
@@ -271,7 +282,9 @@ def test_units_api_exposure(self):
         # Core units functionality should be accessible
         units_api = [
             "units",  # Units registry
-            "create_quantity",  # Factory function
+            "quantity",  # THE factory function (D14)
+            "UWQuantity",  # the factory's return type
+            "create_quantity",  # deprecated alias, kept one cycle
         ]
 
         for api_item in units_api:

From 09a7e42caacbbb59614b81f453733bc79405941e Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:13:03 +0100
Subject: [PATCH 07/11] refactor(api): WC-07 SNES_Poisson constructor order
 (mesh, u_Field, degree, verbose)

SNES_Poisson was uniquely inverted against its base class and every sibling
solver: (mesh, u_Field, verbose, degree). A positional SNES_Poisson(mesh,
None, 3) intended as degree=3 silently set verbose=3 and auto-created a
degree-2 unknown. Reorder to the family-wide (mesh, u_Field, degree, verbose)
with a legacy shim keyed on type(degree) is bool (not isinstance -
isinstance(True, int) is True), emitting one DeprecationWarning.

Underworld development team with AI support from Claude Code
---
 src/underworld3/systems/solvers.py | 27 +++++++++++++++++++++++----
 1 file changed, 23 insertions(+), 4 deletions(-)

diff --git a/src/underworld3/systems/solvers.py b/src/underworld3/systems/solvers.py
index 89245b97..01bbb144 100644
--- a/src/underworld3/systems/solvers.py
+++ b/src/underworld3/systems/solvers.py
@@ -175,15 +175,24 @@ class SNES_Poisson(SNES_Scalar):
         The computational mesh.
     u_Field : MeshVariable, optional
         Pre-existing mesh variable for the solution. If None, one is created.
-    verbose : bool, optional
-        Enable verbose output during solve.
     degree : int, optional
         Polynomial degree for the solution field (default: 2).
+    verbose : bool, optional
+        Enable verbose output during solve.
     DuDt : SemiLagrangian_DDt or Lagrangian_DDt, optional
         Time derivative operator for time-dependent problems.
     DFDt : SemiLagrangian_DDt or Lagrangian_DDt, optional
         Time derivative operator for the flux.
 
+    Notes
+    -----
+    The constructor follows the family-wide parameter order
+    ``(mesh, u_Field, degree, verbose, ...)`` (Style Charter, API
+    conventions); SNES_Poisson historically inverted ``(verbose, degree)``.
+    A legacy positional call is detected by ``type(degree) is bool`` (a
+    polynomial degree is never a bool) and shimmed with one
+    DeprecationWarning.
+
     """
 
     @timing.routine_timer_decorator
@@ -191,12 +200,22 @@ def __init__(
         self,
         mesh: uw.discretisation.Mesh,
         u_Field: uw.discretisation.MeshVariable = None,
-        verbose=False,
         degree=2,
+        verbose=False,
         DuDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None,
         DFDt: Union[SemiLagrangian_DDt, Lagrangian_DDt] = None,
     ):
-        ## Keep track
+        if type(degree) is bool:
+            # Legacy positional order (mesh, u_Field, verbose, degree): the
+            # bool in the degree slot is the legacy verbose flag; a legacy
+            # positional degree, if present, landed in the verbose slot.
+            import warnings
+            warnings.warn(
+                "SNES_Poisson(mesh, u_Field, verbose, degree) is deprecated; "
+                "use SNES_Poisson(mesh, u_Field, degree, verbose)",
+                DeprecationWarning, stacklevel=2)
+            degree, verbose = (
+                verbose if type(verbose) is not bool else 2), degree
 
         ## Parent class will set up default values etc
         super().__init__(

From 915ef4428488937801c42d44831398fcbfd5f580 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:14:15 +0100
Subject: [PATCH 08/11] refactor(api): WC-09 rename boundary_flux_to_field to
 boundary_flux_field; document scale vs buoyancy_scale

The free function is the implementation behind the solver method
boundary_flux_field but was spelled differently. Rename it to match; the old
name remains for one cycle as a warning alias. Document (free function and
method) that scale is a generic multiplier with scale = -1/buoyancy_scale for
dynamic topography - the parameters are deliberately NOT aliased (sign and
reciprocal differ; see API-09 in the 2026-07 API consistency review).

Underworld development team with AI support from Claude Code
---
 .../cython/petsc_generic_snes_solvers.pyx     |  9 ++++--
 src/underworld3/utilities/boundary_flux.py    | 29 ++++++++++++++++---
 2 files changed, 32 insertions(+), 6 deletions(-)

diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
index 615784d6..04d0b269 100644
--- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx
+++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx
@@ -2311,8 +2311,13 @@ class SolverBaseClass(uw_object):
         MeshVariable ``field`` at the boundary nodes (interior untouched), multiplied by
         ``scale``. This is the field hand-off for downstream machinery (surface heat
         flux for coupling, or — with ``remove_mean=True`` and ``scale=-1/(\Delta\rho g)``
-        — dynamic topography). Returns ``field``."""
-        from underworld3.utilities.boundary_flux import boundary_flux_to_field as _bff
+        — dynamic topography). Returns ``field``.
+
+        Note that ``scale`` is a generic multiplier, NOT the ``buoyancy_scale``
+        taken by :meth:`dynamic_topography` / ``topography``: for topography the
+        relationship is ``scale = -1 / buoyancy_scale`` (there the division by
+        :math:`\Delta\rho\,g` and the minus sign are internal)."""
+        from underworld3.utilities.boundary_flux import boundary_flux_field as _bff
         return _bff(self, boundary, field, mass=mass, remove_mean=remove_mean,
                     scale=scale, normal=normal)
 
diff --git a/src/underworld3/utilities/boundary_flux.py b/src/underworld3/utilities/boundary_flux.py
index a06ab808..a36c5466 100644
--- a/src/underworld3/utilities/boundary_flux.py
+++ b/src/underworld3/utilities/boundary_flux.py
@@ -253,10 +253,21 @@ def write_boundary_scalar_field(solver, field, value_by_key, dim):
     return field
 
 
-def boundary_flux_to_field(solver, boundary, field, mass="lumped",
-                           remove_mean=False, scale=1.0, normal=None):
-    """See ``SolverBaseClass.boundary_flux_field``. Writes ``scale * flux`` onto the
-    scalar MeshVariable ``field`` at the boundary nodes (interior untouched)."""
+def boundary_flux_field(solver, boundary, field, mass="lumped",
+                        remove_mean=False, scale=1.0, normal=None):
+    r"""See ``SolverBaseClass.boundary_flux_field`` (the documented entry point;
+    this free function is its implementation and shares its name). Writes
+    ``scale * flux`` onto the scalar MeshVariable ``field`` at the boundary
+    nodes (interior untouched).
+
+    ``scale`` is a generic multiplier on the recovered flux. For dynamic
+    topography it is the **negated reciprocal** of the buoyancy scale used by
+    the expression-return paths: ``scale = -1 / buoyancy_scale``, where
+    ``buoyancy_scale`` is :math:`\Delta\rho\, g` as taken by
+    ``dynamic_topography`` / ``topography`` (there the division and the minus
+    sign are internal). The two parameters are deliberately NOT aliased —
+    treating them as one factor invites sign/reciprocal errors.
+    """
     dim = solver.mesh.dim
     xs, flux = boundary_flux(solver, boundary, mass=mass, remove_mean=remove_mean, normal=normal)
     flux = np.asarray(flux)
@@ -270,3 +281,13 @@ def boundary_flux_to_field(solver, boundary, field, mass="lumped",
             "component, or use boundary_flux() directly for the full vector.")
     fmap = {_key(x, dim): scale * float(f) for x, f in zip(np.asarray(xs), flux.ravel())}
     return write_boundary_scalar_field(solver, field, fmap, dim)
+
+
+def boundary_flux_to_field(*args, **kwargs):
+    """Deprecated alias for :func:`boundary_flux_field` (renamed 2026-07 so the
+    free function matches the solver method it implements; kept one cycle)."""
+    import warnings
+    warnings.warn(
+        "boundary_flux_to_field is renamed; use boundary_flux_field(...)",
+        DeprecationWarning, stacklevel=2)
+    return boundary_flux_field(*args, **kwargs)

From b88e294365b4f0f7cf87ee7245317227c82d3bc7 Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:16:09 +0100
Subject: [PATCH 09/11] refactor(api): WC-12 deprecate the no-op sync= kwarg on
 swarm pack/unpack methods

The sync argument on the four SwarmVariable pack/unpack methods
(pack_uw_data_to_petsc, unpack_uw_data_from_petsc, pack_raw_data_to_petsc,
unpack_raw_data_from_petsc) guarded only an empty TODO stub and never had any
effect. The default becomes a None sentinel; passing sync explicitly emits one
DeprecationWarning; the dead 'if sync: pass' blocks are removed and the six
internal call sites stop passing it. The same-named MeshVariable methods are a
separate copy scheduled for the FO-01 shared array-view refactor and are not
touched here.

Underworld development team with AI support from Claude Code
---
 src/underworld3/swarm.py | 71 +++++++++++++++++++---------------------
 1 file changed, 33 insertions(+), 38 deletions(-)

diff --git a/src/underworld3/swarm.py b/src/underworld3/swarm.py
index 7e631b0a..c648f01d 100644
--- a/src/underworld3/swarm.py
+++ b/src/underworld3/swarm.py
@@ -397,7 +397,7 @@ def _create_variable_array(self, initial_data=None):
             Array object with callback for automatic PETSc synchronization
         """
         if initial_data is None:
-            initial_data = self.unpack_uw_data_from_petsc(squeeze=False, sync=True)
+            initial_data = self.unpack_uw_data_from_petsc(squeeze=False)
 
         # Create NDArray_With_Callback (following swarm._points pattern)
         array_obj = uw.utilities.NDArray_With_Callback(
@@ -432,7 +432,7 @@ def variable_update_callback(array, change_context):
                 return
 
             # Persist changes to PETSc (like swarm callback updates coordinates)
-            var.pack_uw_data_to_petsc(array, sync=True)
+            var.pack_uw_data_to_petsc(array)
 
         # Register the callback (following swarm.points pattern)
         array_obj.add_callback(variable_update_callback)
@@ -457,7 +457,7 @@ def _create_canonical_data_array(self, initial_data=None):
         """
         if initial_data is None:
             # Use unpack_raw to get flat format (-1, num_components)
-            initial_data = self.unpack_raw_data_from_petsc(squeeze=False, sync=True)
+            initial_data = self.unpack_raw_data_from_petsc(squeeze=False)
 
             # Handle case where unpack returns None (swarm not initialized)
             if initial_data is None:
@@ -504,7 +504,7 @@ def canonical_data_callback(array, change_context):
                 canonical_array = canonical_array.reshape(-1, self.num_components)
 
             # STEP 1: Sync to PETSc using established method with correct shape
-            self.pack_raw_data_to_petsc(canonical_array, sync=True)
+            self.pack_raw_data_to_petsc(canonical_array)
 
             # Coordinate writes may strand particles on the wrong rank. Mark
             # the swarm for DEFERRED migration — migrate() itself is
@@ -1271,7 +1271,18 @@ def _rbf_reduce_to_meshVar(self, meshVar, verbose=False):
 
     #     return
 
-    def pack_uw_data_to_petsc(self, data_array, sync=True):
+    @staticmethod
+    def _warn_deprecated_sync(sync):
+        """One-cycle keyword shim for the removed no-op ``sync=`` argument on
+        the four pack/unpack methods (it never had an effect)."""
+        if sync is not None:
+            import warnings
+            warnings.warn(
+                "the 'sync' argument never had an effect and is deprecated; "
+                "remove it",
+                DeprecationWarning, stacklevel=3)
+
+    def pack_uw_data_to_petsc(self, data_array, sync=None):
         """
         Enhanced pack method that directly accesses PETSc field without access() context.
         Designed for the new swarmVariable.array interface.
@@ -1280,9 +1291,10 @@ def pack_uw_data_to_petsc(self, data_array, sync=True):
         ----------
         data_array : numpy.ndarray
             Array data to pack into PETSc field
-        sync : bool
-            Whether to sync parallel operations (default True)
+        sync : deprecated
+            Never had an effect; deprecated (one DeprecationWarning if passed).
         """
+        self._warn_deprecated_sync(sync)
         shape = self.shape
         data_array_3d = data_array.reshape(-1, *self.shape)
 
@@ -1302,11 +1314,6 @@ def pack_uw_data_to_petsc(self, data_array, sync=True):
             # Update the proxy mesh variable if one exists (for integral calculations)
             self._update()
 
-            # Sync parallel operations if requested
-            if sync:
-                # TODO: Add parallel sync logic here if needed
-                pass
-
         finally:
             # Always restore the field
             self.swarm.dm.restoreField(self.clean_name)
@@ -1334,7 +1341,7 @@ def pack_uw_data_to_petsc(self, data_array, sync=True):
     #     else:
     #         return data_array_3d
 
-    def unpack_uw_data_from_petsc(self, squeeze=True, sync=True):
+    def unpack_uw_data_from_petsc(self, squeeze=True, sync=None):
         """
         Enhanced unpack method that directly accesses PETSc field without access() context.
         Designed for the new swarmVariable.array interface.
@@ -1343,20 +1350,16 @@ def unpack_uw_data_from_petsc(self, squeeze=True, sync=True):
         ----------
         squeeze : bool
             Whether to squeeze singleton dimensions (default True)
-        sync : bool
-            Whether to sync parallel operations (default True)
+        sync : deprecated
+            Never had an effect; deprecated (one DeprecationWarning if passed).
         """
+        self._warn_deprecated_sync(sync)
         shape = self.shape
 
         # Direct PETSc field access without context manager
         petsc_data = self.swarm.dm.getField(self.clean_name).reshape((-1, self.num_components))
 
         try:
-            # Sync parallel operations if requested
-            if sync:
-                # TODO: Add parallel sync logic here if needed
-                pass
-
             # Unpack data using same layout as original method
             points = petsc_data.shape[0]
             data_array_3d = np.empty(shape=(points, *shape), dtype=petsc_data.dtype)
@@ -1375,7 +1378,7 @@ def unpack_uw_data_from_petsc(self, squeeze=True, sync=True):
         else:
             return data_array_3d
 
-    def pack_raw_data_to_petsc(self, data_array, sync=True):
+    def pack_raw_data_to_petsc(self, data_array, sync=None):
         """
         Pack data array to PETSc using traditional data shape (-1, num_components).
         Direct PETSc access without access() context for backward compatibility.
@@ -1384,9 +1387,10 @@ def pack_raw_data_to_petsc(self, data_array, sync=True):
         ----------
         data_array : numpy.ndarray
             Array data in traditional flat format (-1, num_components)
-        sync : bool
-            Whether to sync parallel operations (default True)
+        sync : deprecated
+            Never had an effect; deprecated (one DeprecationWarning if passed).
         """
+        self._warn_deprecated_sync(sync)
         import numpy as np
 
         # Convert to expected shape: (-1, num_components)
@@ -1409,18 +1413,13 @@ def pack_raw_data_to_petsc(self, data_array, sync=True):
             # Update the proxy mesh variable if one exists (for integral calculations)
             self._update()
 
-            # Sync parallel operations if requested
-            if sync:
-                # TODO: Add parallel sync logic here if needed
-                pass
-
         finally:
             # Always restore the field
             self.swarm.dm.restoreField(self.clean_name)
 
         return
 
-    def unpack_raw_data_from_petsc(self, squeeze=True, sync=True):
+    def unpack_raw_data_from_petsc(self, squeeze=True, sync=None):
         """
         Unpack data from PETSc in traditional data shape (-1, num_components).
         Direct PETSc access without access() context for backward compatibility.
@@ -1429,14 +1428,15 @@ def unpack_raw_data_from_petsc(self, squeeze=True, sync=True):
         ----------
         squeeze : bool
             Whether to remove singleton dimensions (default True)
-        sync : bool
-            Whether to sync parallel operations (default True)
+        sync : deprecated
+            Never had an effect; deprecated (one DeprecationWarning if passed).
 
         Returns
         -------
         numpy.ndarray
             Array data in traditional flat format (-1, num_components)
         """
+        self._warn_deprecated_sync(sync)
         import numpy as np
 
         # Check if swarm has any particles before accessing field
@@ -1458,11 +1458,6 @@ def unpack_raw_data_from_petsc(self, squeeze=True, sync=True):
             # Return data in traditional flat format
             result = petsc_data.copy()
 
-            # Sync parallel operations if requested
-            if sync:
-                # TODO: Add parallel sync logic here if needed
-                pass
-
         finally:
             # Always restore the field
             self.swarm.dm.restoreField(self.clean_name)
@@ -1522,7 +1517,7 @@ def rbf_interpolate(self, new_coords, verbose=False, nnn=None):
         import numpy as np
 
         # Get data directly from PETSc to avoid circular callback dependencies
-        raw_data = self.unpack_raw_data_from_petsc(squeeze=False, sync=False)
+        raw_data = self.unpack_raw_data_from_petsc(squeeze=False)
         data_size = raw_data.shape
 
         # What to do if there are no particles: never SILENTLY return zeros
@@ -2897,7 +2892,7 @@ def _flush_pending_petsc_sync(self):
                     f"DMSwarm holds {self.dm.getLocalSize()} particles. The "
                     "particle layout changed while migration was disabled."
                 )
-            var.pack_raw_data_to_petsc(arr, sync=True)
+            var.pack_raw_data_to_petsc(arr)
             if self._coord_var is var:
                 self._needs_migration = True
             if hasattr(var, "_on_data_changed"):

From b38cc06662e854d90a9b83e78177823aed66495c Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:16:51 +0100
Subject: [PATCH 10/11] refactor(api): WC-13 _winslow_spring drops dead
 relax/step_frac; n_sweeps renamed max_cg_iters

relax= and step_frac= were unused on the spring-equilibrium path ('kept only
for signature stability') and are removed. n_sweeps misnamed the nonlinear-CG
iteration cap and is renamed max_cg_iters; the old spelling is accepted for
one cycle with a DeprecationWarning (it is reachable from user code via
smooth_mesh_interior(method_kwargs=...)). Signature-only; the D-smoothing
readability group is a separate wave.

Underworld development team with AI support from Claude Code
---
 src/underworld3/meshing/smoothing.py | 25 +++++++++++++++++--------
 1 file changed, 17 insertions(+), 8 deletions(-)

diff --git a/src/underworld3/meshing/smoothing.py b/src/underworld3/meshing/smoothing.py
index 02157f34..7d35a4cc 100644
--- a/src/underworld3/meshing/smoothing.py
+++ b/src/underworld3/meshing/smoothing.py
@@ -504,8 +504,9 @@ def _edge_pairs(dm):
 
 
 def _winslow_spring(mesh, metric, pinned_labels, verbose,
-                    n_sweeps=300, relax=None, step_frac=None,
-                    boundary_slip=False, shape_w=1.0, size_w=8.0):
+                    max_cg_iters=300,
+                    boundary_slip=False, shape_w=1.0, size_w=8.0,
+                    n_sweeps=None):
     r"""Metric-driven mesh grading by elastic-spring equilibrium.
 
     Every mesh edge is a linear spring whose *rest length* is set
@@ -540,12 +541,20 @@ def _winslow_spring(mesh, metric, pinned_labels, verbose,
     the mean edge length ⇒ only a benign mild regularisation toward
     uniform spacing (no grading change).
 
-    ``n_sweeps`` caps the CG iterations (CG converges far faster
-    than the old Jacobi sweep budget). ``relax`` / ``step_frac`` are
-    unused on the equilibrium path (the CG line search controls the
-    step and the inversion guard) and are kept only for signature
-    stability. ``n_iters`` / ``alpha`` do not apply.
+    ``max_cg_iters`` caps the CG iterations (CG converges far faster
+    than the old Jacobi sweep budget); ``n_sweeps`` is its deprecated
+    former name, accepted for one cycle with a DeprecationWarning.
+    The old ``relax`` / ``step_frac`` parameters were unused on the
+    equilibrium path (the CG line search controls the step and the
+    inversion guard) and have been removed. ``n_iters`` / ``alpha``
+    do not apply.
     """
+    if n_sweeps is not None:
+        warnings.warn(
+            "the 'n_sweeps' argument is renamed; use max_cg_iters= "
+            "(it caps the nonlinear-CG iterations, not Jacobi sweeps)",
+            DeprecationWarning, stacklevel=2)
+        max_cg_iters = n_sweeps
     pinned_labels = tuple(pinned_labels)
     dm = mesh.dm
     pStart, pEnd = dm.getDepthStratum(0)
@@ -738,7 +747,7 @@ def _energy_grad(X):
         dmax = uw.mpi.comm.allreduce(dmax, op=_MPI.MAX)
     t0 = 0.5 * L0_mean / dmax
     c_arm = 1.0e-4
-    max_iter = int(n_sweeps)
+    max_iter = int(max_cg_iters)
     for it in range(max_iter):
         gnorm = _allsum((G * G).sum()) ** 0.5
         if gnorm <= 1.0e-8 * g0:

From 692b96a697ada7db4c44fceba5c652d3d3fa2bce Mon Sep 17 00:00:00 2001
From: lmoresi <louis.moresi@anu.edu.au>
Date: Mon, 6 Jul 2026 11:20:03 +0100
Subject: [PATCH 11/11] test(api): WC contract tests for the Wave C shims
 (test_0641)

Two-test pattern for every Wave C shim: the legacy spelling produces the
identical result with exactly one DeprecationWarning naming the replacement;
the canonical spelling runs warning-free under an error filter. Covers
WC-01/02 (value-first BC order + conds datum on the three migrated methods,
trio regression, components= warning), WC-03 (consistent_jacobian
validation), WC-04 (set_custom_fmg / set_custom_mg), WC-05 (namespace
exposure), WC-06 (quantity factory), WC-07 (SNES_Poisson order), WC-08
(vector Nitsche normal=/mask=), WC-09 (boundary_flux_field rename), WC-12
(swarm sync=), WC-13 (max_cg_iters rename). Validated to fail (34/40)
against the pre-shim build. Markers: level_1, tier_a.

Underworld development team with AI support from Claude Code
---
 tests/test_0641_wave_c_api_shims.py | 423 ++++++++++++++++++++++++++++
 1 file changed, 423 insertions(+)
 create mode 100644 tests/test_0641_wave_c_api_shims.py

diff --git a/tests/test_0641_wave_c_api_shims.py b/tests/test_0641_wave_c_api_shims.py
new file mode 100644
index 00000000..5c8015e4
--- /dev/null
+++ b/tests/test_0641_wave_c_api_shims.py
@@ -0,0 +1,423 @@
+"""Wave C API-harmonization contract tests (WC-01..WC-13).
+
+Every deprecation shim introduced by the 2026-07 Wave C API harmonization
+carries the two-test contract here:
+
+1. the OLD signature/spelling produces the identical result and emits exactly
+   one ``DeprecationWarning`` naming the replacement;
+2. the NEW canonical signature emits no ``DeprecationWarning`` at all.
+
+Canonical conventions under test (see ``docs/developer/UW3_STYLE_CHARTER.md``
+paragraph 6 and ``docs/reviews/2026-07/API-CONSISTENCY-REVIEW.md``):
+
+- BC argument order is value-first: ``add_<kind>_bc(conds, boundary, ...)``
+  (maintainer decision D2, 2026-07-04); the datum name is ``conds`` (D3).
+- ``consistent_jacobian`` is a validated property (WC-03).
+- ``set_custom_fmg`` is the solver-method entry point for custom MG (WC-04).
+- ``uw.quantity`` is THE quantity factory (D14, WC-06).
+"""
+
+import warnings
+
+import pytest
+import sympy
+
+import underworld3 as uw
+
+pytestmark = [pytest.mark.level_1, pytest.mark.tier_a]
+
+
+# ---------------------------------------------------------------------------
+# helpers
+# ---------------------------------------------------------------------------
+
+def _one_deprecation(record):
+    """Count the DeprecationWarnings in a pytest.warns record."""
+    return sum(
+        1 for w in record if issubclass(w.category, DeprecationWarning)
+    )
+
+
+class _no_deprecation(warnings.catch_warnings):
+    """Context manager: any DeprecationWarning inside is a test failure."""
+
+    def __enter__(self):
+        log = super().__enter__()
+        warnings.simplefilter("error", DeprecationWarning)
+        return log
+
+
+@pytest.fixture(scope="module")
+def mesh():
+    return uw.meshing.StructuredQuadBox(
+        elementRes=(2, 2), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)
+    )
+
+
+@pytest.fixture()
+def stokes(mesh):
+    solver = uw.systems.Stokes(mesh)
+    solver.constitutive_model = uw.constitutive_models.ViscousFlowModel
+    return solver
+
+
+# ---------------------------------------------------------------------------
+# WC-01 / WC-02 - value-first BC argument order, `conds` datum name
+# ---------------------------------------------------------------------------
+
+class TestStokesNitscheValueFirst:
+    def test_new_order_no_warning(self, mesh, stokes):
+        with _no_deprecation():
+            stokes.add_nitsche_bc(0.5, "Top")
+        assert stokes.natural_bcs[-1].boundary == "Top"
+
+    def test_legacy_positional_order_warns_once(self, mesh, stokes):
+        with pytest.warns(DeprecationWarning, match=r"add_nitsche_bc\(conds, boundary") as rec:
+            stokes.add_nitsche_bc("Top", 0.5)
+        assert _one_deprecation(rec) == 1
+
+    def test_legacy_g_keyword_warns_once(self, mesh, stokes):
+        with pytest.warns(DeprecationWarning, match=r"add_nitsche_bc\(conds, boundary") as rec:
+            stokes.add_nitsche_bc("Top", g=0.5)
+        assert _one_deprecation(rec) == 1
+
+    def test_legacy_and_new_give_identical_bc(self, mesh, stokes):
+        # Register the same condition through both spellings on ONE solver so
+        # every symbol (velocity, pressure, viscosity, normals) is shared -
+        # the stored BC tuples must be identical.
+        with pytest.warns(DeprecationWarning):
+            stokes.add_nitsche_bc("Top", 0.5)
+        with _no_deprecation():
+            stokes.add_nitsche_bc(0.5, "Top")
+
+        bc_old = stokes.natural_bcs[-2]
+        bc_new = stokes.natural_bcs[-1]
+        assert bc_old.boundary == bc_new.boundary
+        assert bc_old.fn_f == bc_new.fn_f
+        assert bc_old.fn_F == bc_new.fn_F
+        assert bc_old.fn_p == bc_new.fn_p
+
+    def test_datum_given_twice_is_an_error(self, mesh, stokes):
+        with pytest.raises(TypeError, match="conds"):
+            stokes.add_nitsche_bc(0.5, "Top", g=0.5)
+
+    def test_missing_boundary_is_an_error(self, mesh, stokes):
+        with pytest.raises(TypeError, match="boundary"):
+            stokes.add_nitsche_bc(0.5)
+
+
+class TestRotatedFreeslipValueFirst:
+    def test_new_order_no_warning(self, mesh, stokes):
+        with _no_deprecation():
+            stokes.add_rotated_freeslip_bc(0, "Top")
+        assert stokes._rotated_freeslip_bcs[-1] == ("Top", None)
+
+    def test_boundary_keyword_no_warning(self, mesh, stokes):
+        with _no_deprecation():
+            stokes.add_rotated_freeslip_bc(boundary="Bottom")
+        assert stokes._rotated_freeslip_bcs[-1] == ("Bottom", None)
+
+    def test_legacy_boundary_first_warns_once(self, mesh, stokes):
+        with pytest.warns(DeprecationWarning, match=r"add_rotated_freeslip_bc\(conds, boundary") as rec:
+            stokes.add_rotated_freeslip_bc("Top")
+        assert _one_deprecation(rec) == 1
+        assert stokes._rotated_freeslip_bcs[-1] == ("Top", None)
+
+    def test_legacy_positional_normal_is_preserved(self, mesh, stokes):
+        n = sympy.Matrix([[0, 1]])
+        with pytest.warns(DeprecationWarning) as rec:
+            stokes.add_rotated_freeslip_bc("Top", n)
+        assert _one_deprecation(rec) == 1
+        assert stokes._rotated_freeslip_bcs[-1] == ("Top", n)
+
+    def test_nonzero_datum_not_implemented(self, mesh, stokes):
+        with pytest.raises(NotImplementedError):
+            stokes.add_rotated_freeslip_bc(1.0, "Top")
+
+
+class TestConstraintBCValueFirst:
+    def test_new_order_no_warning(self, mesh):
+        solver = uw.systems.Stokes_Constrained(mesh)
+        solver.constitutive_model = uw.constitutive_models.ViscousFlowModel
+        with _no_deprecation():
+            h = solver.add_constraint_bc(0.0, "Top")
+        assert solver._block_constraint_bcs[-1].boundary == "Top"
+        assert h is not None
+
+    def test_legacy_order_warns_once(self, mesh):
+        solver = uw.systems.Stokes_Constrained(mesh)
+        solver.constitutive_model = uw.constitutive_models.ViscousFlowModel
+        with pytest.warns(DeprecationWarning, match=r"add_constraint_bc\(conds, boundary") as rec:
+            solver.add_constraint_bc("Top", 0.0)
+        assert _one_deprecation(rec) == 1
+        cbc = solver._block_constraint_bcs[-1]
+        assert cbc.boundary == "Top"
+        assert float(cbc.g) == 0.0
+
+    def test_legacy_g_keyword_warns_once(self, mesh):
+        solver = uw.systems.Stokes_Constrained(mesh)
+        solver.constitutive_model = uw.constitutive_models.ViscousFlowModel
+        with pytest.warns(DeprecationWarning) as rec:
+            solver.add_constraint_bc(boundary="Top", g=0.5)
+        assert _one_deprecation(rec) == 1
+        assert solver._block_constraint_bcs[-1].g == sympy.Float(0.5)
+
+
+class TestLegacyTrioUnchanged:
+    """The original trio is already canonical - it must NOT warn."""
+
+    def test_dirichlet_no_warning(self, mesh):
+        solver = uw.systems.Poisson(mesh)
+        with _no_deprecation():
+            solver.add_dirichlet_bc(0.0, "Top")
+            solver.add_essential_bc(1.0, "Bottom")
+            solver.add_natural_bc(0.0, "Left")
+
+    def test_components_kwarg_warns_once(self, mesh, stokes):
+        with pytest.warns(DeprecationWarning, match="components") as rec:
+            stokes.add_dirichlet_bc((0.0, 0.0), "Top", components=(0,))
+        assert _one_deprecation(rec) == 1
+
+
+# ---------------------------------------------------------------------------
+# WC-03 - consistent_jacobian validated property
+# ---------------------------------------------------------------------------
+
+class TestConsistentJacobianValidation:
+    def test_valid_values_round_trip(self, mesh, stokes):
+        for value in (False, True, "continuation"):
+            stokes.consistent_jacobian = value
+            assert stokes.consistent_jacobian == value
+
+    def test_falsy_normalizes_to_false(self, mesh, stokes):
+        for value in (None, 0, 0.0, ""):
+            stokes.consistent_jacobian = value
+            assert stokes.consistent_jacobian is False
+
+    def test_invalid_values_raise(self, mesh, stokes):
+        for value in ("picard", "Continuation", 1, "newton", 2.0):
+            with pytest.raises(ValueError, match="consistent_jacobian"):
+                stokes.consistent_jacobian = value
+
+    def test_default_is_false(self, mesh):
+        solver = uw.systems.Poisson(mesh)
+        assert solver.consistent_jacobian is False
+
+
+# ---------------------------------------------------------------------------
+# WC-04 - set_custom_fmg solver method; set_custom_mg deprecated
+# ---------------------------------------------------------------------------
+
+class TestCustomMGEntryPoint:
+    def test_set_custom_fmg_method_no_warning(self, mesh):
+        coarse = uw.meshing.StructuredQuadBox(
+            elementRes=(1, 1), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)
+        )
+        solver = uw.systems.Poisson(mesh)
+        with _no_deprecation():
+            solver.set_custom_fmg([coarse])
+        assert solver._custom_mg["mode"] == "hierarchy"
+
+    def test_set_custom_mg_warns_once_and_still_registers(self, mesh):
+        coarse = uw.meshing.StructuredQuadBox(
+            elementRes=(1, 1), minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0)
+        )
+        solver = uw.systems.Poisson(mesh)
+        with pytest.warns(DeprecationWarning, match="set_custom_fmg") as rec:
+            solver.set_custom_mg([coarse])
+        assert _one_deprecation(rec) == 1
+        assert solver._custom_mg["kind"] == "barycentric"
+
+
+# ---------------------------------------------------------------------------
+# WC-05 - namespace exposure (no deep imports needed)
+# ---------------------------------------------------------------------------
+
+class TestNamespaceExposure:
+    def test_utilities_modules_exposed(self):
+        assert hasattr(uw.utilities, "rotated_bc")
+        assert hasattr(uw.utilities, "boundary_flux")
+        assert hasattr(uw.utilities, "custom_mg")
+        assert callable(uw.utilities.set_custom_fmg)
+
+    def test_meshing_bounding_surface_exposed(self):
+        assert hasattr(uw.meshing, "BoundingSurface")
+        for name in (
+            "register_radial_surfaces",
+            "register_plane_surfaces",
+            "register_box_face_surfaces",
+        ):
+            assert callable(getattr(uw.meshing, name))
+            assert name in uw.meshing.__all__
+        assert "BoundingSurface" in uw.meshing.__all__
+
+
+# ---------------------------------------------------------------------------
+# WC-06 - uw.quantity is THE factory (D14)
+# ---------------------------------------------------------------------------
+
+class TestQuantityFactory:
+    def test_quantity_no_warning_and_type_exposed(self):
+        with _no_deprecation():
+            q = uw.quantity(1.0, "m")
+        assert isinstance(q, uw.UWQuantity)
+
+    def test_create_quantity_warns_once(self):
+        with pytest.warns(DeprecationWarning, match="uw.quantity") as rec:
+            q = uw.create_quantity(10.0, "m/s")
+        assert _one_deprecation(rec) == 1
+        # behaviour preserved for one cycle: still the raw Pint quantity
+        assert float(q.magnitude) == 10.0
+        assert str(q.units) in ("meter / second", "m / s")
+
+
+# ---------------------------------------------------------------------------
+# WC-07 - SNES_Poisson constructor order (mesh, u_Field, degree, verbose)
+# ---------------------------------------------------------------------------
+
+class TestPoissonConstructorOrder:
+    def test_positional_degree_no_warning(self, mesh):
+        with _no_deprecation():
+            solver = uw.systems.Poisson(mesh, None, 3)
+        assert solver.u.degree == 3
+        assert solver.verbose is False
+
+    def test_legacy_positional_verbose_warns_once(self, mesh):
+        with pytest.warns(DeprecationWarning, match=r"degree, verbose") as rec:
+            solver = uw.systems.Poisson(mesh, None, True)
+        assert _one_deprecation(rec) == 1
+        assert solver.verbose is True
+        assert solver.u.degree == 2  # legacy default preserved
+
+    def test_legacy_four_positional_args(self, mesh):
+        with pytest.warns(DeprecationWarning) as rec:
+            solver = uw.systems.Poisson(mesh, None, True, 3)
+        assert _one_deprecation(rec) == 1
+        assert solver.verbose is True
+        assert solver.u.degree == 3
+
+    def test_keyword_calls_unchanged(self, mesh):
+        with _no_deprecation():
+            solver = uw.systems.Poisson(mesh, degree=1, verbose=False)
+        assert solver.u.degree == 1
+
+
+# ---------------------------------------------------------------------------
+# WC-08 - SNES_Vector.add_nitsche_bc aligned with the Stokes variant
+# ---------------------------------------------------------------------------
+
+class TestVectorNitscheAlignment:
+    @pytest.fixture()
+    def vector_solver(self, mesh):
+        v = uw.discretisation.MeshVariable(
+            "v_wc08", mesh, mesh.dim, degree=2, vtype=uw.VarType.VECTOR
+        )
+        solver = uw.systems.Vector_Projection(mesh, v)
+        solver.constitutive_model = uw.constitutive_models.ViscousFlowModel
+        return solver
+
+    def test_normal_kwarg_accepted(self, mesh, vector_solver):
+        n = sympy.Matrix([[0, 1]])
+        with _no_deprecation():
+            vector_solver.add_nitsche_bc(0.0, "Top", normal=n)
+        assert vector_solver.natural_bcs[-1].boundary == "Top"
+
+    def test_mask_raises_not_implemented(self, mesh, vector_solver):
+        with pytest.raises(NotImplementedError, match="mask"):
+            vector_solver.add_nitsche_bc(0.0, "Top", mask=sympy.Integer(1))
+
+    def test_legacy_order_warns_once(self, mesh, vector_solver):
+        with pytest.warns(DeprecationWarning, match=r"add_nitsche_bc\(conds, boundary") as rec:
+            vector_solver.add_nitsche_bc("Top", 0.0)
+        assert _one_deprecation(rec) == 1
+
+
+# ---------------------------------------------------------------------------
+# WC-09 - boundary_flux_to_field renamed boundary_flux_field
+# ---------------------------------------------------------------------------
+
+class TestBoundaryFluxRename:
+    def test_new_name_importable_no_warning(self):
+        with _no_deprecation():
+            from underworld3.utilities.boundary_flux import boundary_flux_field
+        assert callable(boundary_flux_field)
+
+    def test_old_name_warns_once(self):
+        from underworld3.utilities.boundary_flux import boundary_flux_to_field
+
+        with pytest.warns(DeprecationWarning, match="boundary_flux_field") as rec:
+            # The warning is emitted before any work; the junk arguments then
+            # fail inside the real implementation, which is fine here.
+            with pytest.raises(Exception):
+                boundary_flux_to_field(None, "Top", None)
+        assert _one_deprecation(rec) == 1
+
+
+# ---------------------------------------------------------------------------
+# WC-12 - no-op `sync=` kwarg deprecated on swarm pack/unpack
+# ---------------------------------------------------------------------------
+
+class TestSwarmSyncDeprecated:
+    @pytest.fixture(scope="class")
+    def swarm_var(self, mesh):
+        swarm = uw.swarm.Swarm(mesh)
+        var = swarm.add_variable("s_wc12", 1)
+        swarm.populate(fill_param=1)
+        # yield keeps the parent swarm alive for the lifetime of the tests
+        # (the variable holds only a weak link to it)
+        yield var
+        del swarm
+
+    def test_unpack_no_sync_no_warning(self, swarm_var):
+        with _no_deprecation():
+            data = swarm_var.unpack_uw_data_from_petsc(squeeze=False)
+        assert data is not None
+
+    def test_unpack_sync_warns_once_identical_result(self, swarm_var):
+        import numpy as np
+
+        clean = swarm_var.unpack_uw_data_from_petsc(squeeze=False)
+        with pytest.warns(DeprecationWarning, match="sync") as rec:
+            legacy = swarm_var.unpack_uw_data_from_petsc(squeeze=False, sync=True)
+        assert _one_deprecation(rec) == 1
+        assert np.array_equal(np.asarray(clean), np.asarray(legacy))
+
+    def test_pack_sync_warns_once(self, swarm_var):
+        data = swarm_var.unpack_uw_data_from_petsc(squeeze=False)
+        with pytest.warns(DeprecationWarning, match="sync") as rec:
+            swarm_var.pack_uw_data_to_petsc(data, sync=True)
+        assert _one_deprecation(rec) == 1
+        with _no_deprecation():
+            swarm_var.pack_uw_data_to_petsc(data)
+
+
+# ---------------------------------------------------------------------------
+# WC-13 - smoothing: dead params dropped, n_sweeps renamed max_cg_iters
+# ---------------------------------------------------------------------------
+
+class TestWinslowSpringSignature:
+    def test_dead_params_gone_new_name_present(self):
+        import inspect
+
+        from underworld3.meshing.smoothing import _winslow_spring
+
+        params = inspect.signature(_winslow_spring).parameters
+        assert "relax" not in params
+        assert "step_frac" not in params
+        assert "max_cg_iters" in params
+        assert "n_sweeps" in params  # deprecated alias, one cycle
+
+    def test_n_sweeps_alias_warns_once(self):
+        import numpy as np
+
+        from underworld3.meshing.smoothing import _winslow_spring
+
+        tri_mesh = uw.meshing.UnstructuredSimplexBox(
+            minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0), cellSize=0.5
+        )
+        pinned = ("Top", "Bottom", "Left", "Right")
+        with pytest.warns(DeprecationWarning, match="max_cg_iters") as rec:
+            _winslow_spring(tri_mesh, sympy.Integer(1), pinned, False, n_sweeps=1)
+        assert _one_deprecation(rec) == 1
+        with _no_deprecation():
+            _winslow_spring(tri_mesh, sympy.Integer(1), pinned, False, max_cg_iters=1)