underworldcode · lmoresi · Sep 15, 2026 · Sep 14, 2026
@@ -85,10 +85,19 @@ information, not a regression.
 - Do NOT change library code to make one pass. Re-characterise it and say why.
 - If an assertion would break when the code gets better, this is its home.
 
-Examples in the tree: `test_1060_nitsche_freeslip.py`
-(`test_constraint_strength_ordering_characterisation`),
-`test_0773_surface_smoother.py`, `test_0066_integration_point_slcn.py`,
-`test_1070_free_surface_plume.py`.
+Prefer a hard baseline wherever one exists — Charter §8 requires it, and a
+characterisation is the fallback for when none does, not a place to park an
+assertion that was easier to write. Four method-comparison tests were moved here
+in 2026-09 and then removed again, because in every case the baseline was already
+available: the comparisons were less informative than the numbers they were
+computed from. That included the last one: `test_1060_nitsche_freeslip.py` asserted that the
+weak constraints must STAY inaccurate (`leak > 1e-5`), which would have failed if
+Nitsche improved. Its exact half — an essential BC holds `v.n` to machine
+precision — is a hard baseline and is now a tier B contract; the weak-path leaks
+are printed for a reader instead of gated on.
+
+C1 is therefore rare by design. Reach for it only when you have measured
+something worth recording and can show that no baseline exists.
 
 **C2 — Experimental.** Test or code (or both) may be incorrect: written for a
 feature that is not finished, exploring what the behaviour should be, or

@@ -115,17 +115,27 @@ temperature.data[:, 0] = values          # BAD  — compatibility layer in new c
   The number sets a broad sequence, nothing more — selection is by marker, and a
   shared number is not a conflict.
 - Validate a new test's own correctness before changing library code to satisfy it.
-- **Assert against a known answer, not against a rival method.** A test that asserts
-  one method is more accurate than another encodes a preference, not a contract: the
-  result moves with the fixture, the mesh, the forcing and every default the two
-  methods carry. Test the analytic or reference solution with an absolute bound.
-  Convergence ORDER and mathematical exactness are contracts and may be asserted
-  freely; "method A scored better than method B here" may not.
-- **Prefer a relative bound.** An absolute threshold silently tracks whatever sets
-  the scale — a free-slip test asserting `|v_n| < 1e-4` was really asserting
-  5.7e-3 relative, and tracked the buoyancy forcing rather than the method.
-- **If an assertion would break when the code gets better, it belongs at tier C.**
-  That is the test to apply, and it is what tier C is for.
+- **Every test asserts against a HARD BASELINE.** An analytic solution, a published
+  value, a closed-form geometric quantity, a conservation identity, an exactness
+  property. If a failure cannot name what is broken, it is not a test — it is a
+  drift detector, and it belongs in a benchmark rather than the suite.
+- **A baseline must be tight enough to fail the moment a default changes.** This is
+  the point of the rule. A loose test does not fail when behaviour moves: it absorbs
+  the change, drifts inside its own margin, and fails later somewhere else for
+  reasons that are hard to trace back. #692 redefined `mesh.cell_size()`, nothing
+  failed, the Nitsche penalty moved 43%, and a spherical-shell benchmark slid from
+  0.2% to 2.4% — still inside its 5% tolerance — before breaking months later on one
+  platform's triangulation (#734). The test that would have caught it on the day is
+  `cell_size` pinned to its closed form on a known simplex.
+- **Never assert that one method beats another.** A ratio between two errors moves
+  when either moves, so it cannot say which; it is strictly less informative than the
+  numbers it was computed from, and it encodes a preference rather than a contract.
+  Assert each method against the baseline instead. Convergence ORDER and mathematical
+  exactness are contracts and may be asserted freely.
+- **Prefer a relative bound to an absolute one** where the scale is set elsewhere. An
+  absolute threshold silently tracks whatever sets that scale — a free-slip test
+  asserting `|v_n| < 1e-4` was really asserting 5.7e-3 relative, and tracked the
+  buoyancy forcing rather than the method under test.
 
 ### The tiers
 

@@ -6647,8 +6647,16 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
             adaptation-tracking) rather than the single **global** minimum
             cell size (:meth:`Mesh.get_min_radius`). On a non-uniform or
             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.
+            on every facet. Set ``False`` to restore the legacy global-h
+            behaviour exactly.
+
+            The two coincide on **tensor** cells only. On a uniform **simplex**
+            mesh they differ by exactly :math:`\sqrt{2}` — for congruent
+            right-isosceles cells of legs :math:`h`, :meth:`Mesh.cell_size` is
+            :math:`2h/3` while :meth:`Mesh.get_min_radius` is
+            :math:`\sqrt{2}h/3` — so the penalty :math:`\gamma\mu/h` differs
+            between the two settings on the simplex meshes the free-slip and
+            fault models use. See ``tests/test_0010_cell_size_geometry.py``.
         g : sympy expression or float, optional
             Deprecated keyword alias for ``conds`` (one DeprecationWarning).
 
@@ -6745,7 +6753,14 @@ class SNES_Stokes_SaddlePt(SolverBaseClass):
         # or adaptively-refined mesh — the boundary kernel sees the adjacent
         # cell's size. The field tracks mesh deformation/adaptation. Set
         # local_h=False to restore the legacy single global-minimum scalar
-        # (mesh.get_min_radius()); on a uniform mesh the two coincide.
+        # (mesh.get_min_radius()).
+        #
+        # The two coincide on TENSOR cells only. On a uniform SIMPLEX mesh --
+        # which is what the free-slip and fault models are built on -- they
+        # differ by exactly sqrt(2): on congruent right-isosceles cells of legs
+        # h, cell_size is 2h/3 and get_min_radius is sqrt(2)h/3. The penalty
+        # gamma*mu/h moves with that, so the two settings are NOT interchangeable
+        # there (see #734 and tests/test_0010_cell_size_geometry.py).
         if local_h:
             h_sym = mesh.cell_size()
         else:

@@ -64,3 +64,74 @@ def test_regular_square_cell_size_keeps_global_radius():
     assert max(uw.mpi.comm.allgather(error)) < 1e-12
     assert global_radius == pytest.approx(expected, rel=1e-12)
     assert all(uw.mpi.comm.allgather(np.array_equal(mesh._radii, legacy)))
+
+
+@pytest.mark.parametrize("h", [0.5, 0.25])
+def test_regular_simplex_cell_size_is_the_closed_form(h):
+    """cell_size on congruent right-isosceles cells is 2h/3, exactly.
+
+    `test_cell_size_matches_own_vertices_and_tracks_deform` above checks the
+    implementation against an independent reading of the same DEFINITION, so it
+    stays true if the definition itself is changed on both sides. This pins the
+    VALUE against geometry instead: `regular=True` tiles the box with congruent
+    right-isosceles triangles of legs h, whose vertices sit at (0,0), (h,0),
+    (0,h) up to rigid motion, so the RMS distance to the centroid is
+
+        sqrt( ( 2(h/3)^2 + 2[(2h/3)^2 + (h/3)^2] ) / 3 ) = 2h/3
+
+    A redefinition of `cell_size` fails here immediately, at the quantity that
+    changed. #692 changed it for simplices and nothing failed, so the Nitsche
+    penalty (gamma*mu/h) moved 43% unnoticed and surfaced months later as a
+    5.67% miss on a spherical-shell benchmark (#734) — on one platform's
+    triangulation only, which is the hardest kind of failure to read backwards.
+    """
+    mesh = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0),
+        cellSize=h, regular=True, qdegree=2,
+    )
+    mesh.cell_size()
+    radii = np.asarray(mesh._cell_size_variable.array[:, 0, 0])
+
+    error = float(np.abs(radii - 2.0 * h / 3.0).max(initial=0.0))
+    assert max(uw.mpi.comm.allgather(error)) < 1e-12, (
+        f"cell_size on congruent legs-{h} right-isosceles cells is "
+        f"{radii.min():.12g}..{radii.max():.12g}, expected exactly {2.0 * h / 3.0:.12g}"
+    )
+
+
+def test_cell_size_and_min_radius_agree_on_tensor_cells_but_not_simplices():
+    """The two mesh-size measures coincide on TENSOR cells only, by sqrt(2).
+
+    `add_nitsche_bc` offers `local_h=False` to fall back from `mesh.cell_size()`
+    to the global `mesh.get_min_radius()`, and its docstring used to say the two
+    coincide "on a uniform mesh". They do on a regular quad box — which is
+    presumably where that was checked — and they do NOT on a uniform SIMPLEX
+    mesh, which is what every free-slip and fault model is built on.
+
+    Both measures have closed forms on congruent right-isosceles cells of legs h:
+    cell_size is 2h/3 (vertex RMS about the centroid) and get_min_radius is
+    sqrt(2)h/3, so the ratio is exactly sqrt(2). Pinning it means neither measure
+    can be redefined without this saying so, and says which way the Nitsche
+    penalty moves when it is.
+    """
+    quad = uw.meshing.StructuredQuadBox(elementRes=(4, 4), qdegree=2)
+    quad.cell_size()
+    quad_local = np.asarray(quad._cell_size_variable.array[:, 0, 0])
+    quad_ratio = float(quad_local.min()) / quad.get_min_radius()
+    assert quad_ratio == pytest.approx(1.0, rel=1e-12), (
+        f"on tensor cells the two measures must coincide; ratio {quad_ratio:.12g}")
+
+    h = 0.25
+    simplex = uw.meshing.UnstructuredSimplexBox(
+        minCoords=(0.0, 0.0), maxCoords=(1.0, 1.0),
+        cellSize=h, regular=True, qdegree=2,
+    )
+    simplex.cell_size()
+    simplex_local = np.asarray(simplex._cell_size_variable.array[:, 0, 0])
+
+    assert simplex.get_min_radius() == pytest.approx(np.sqrt(2.0) * h / 3.0, rel=1e-12)
+    ratio = float(simplex_local.min()) / simplex.get_min_radius()
+    assert ratio == pytest.approx(np.sqrt(2.0), rel=1e-12), (
+        f"cell_size/get_min_radius on uniform simplices is {ratio:.12g}, "
+        "expected sqrt(2) — the Nitsche penalty gamma*mu/h scales with this"
+    )
@@ -122,43 +122,6 @@ def test_rotating_gaussian_ip_accuracy():
     assert l2_ip < 0.02, f"integration-point trace L2 error {l2_ip:.3e}"
 
 
-# tier_c overrides the module-level tier_a for this test alone: it compares two
-# transport managers, so it can fail because one of them got better.
-@pytest.mark.level_2
-@pytest.mark.tier_c
-def test_rotating_gaussian_ip_against_nodal_characterisation():
-    """Characterisation: the integration-point trace is not worse than nodal.
-
-    This compares two METHODS, so it can fail because the code improved — a
-    better nodal SLCN would break it, and that is good news. Tier C: a failure
-    demands an explanation, not a revert. It is NOT the justification for the
-    integration-point path; `test_rotating_gaussian_ip_accuracy` asserts that
-    against the known solution.
-
-    Measured 2026-09-12 on this fixture (cellSize=0.08, dt=0.1, 16 steps):
-    L2 ip 1.70e-3 against nodal 3.96e-3; peak ip 0.9909 against nodal 0.9696.
-    The relationship is sensitive to the Courant number, the quadrature degree
-    and the element size, so those numbers characterise this fixture rather than
-    making a general claim. Compare
-    `project_integration_point_proxy_pic_lip`, where the bulk diagnostics were
-    identical while the interface answer was not.
-    """
-    mesh = uw.meshing.UnstructuredSimplexBox(
-        minCoords=(-1, -1), maxCoords=(1, 1), cellSize=0.08, qdegree=3
-    )
-    dt, nsteps = 0.1, 16
-    l2_nodal, peak_nodal = _rotating_gaussian(mesh, "nodal", dt, nsteps)
-    l2_ip, peak_ip = _rotating_gaussian(mesh, "ip", dt, nsteps)
-
-    print(f"L2: ip={l2_ip:.4e} nodal={l2_nodal:.4e}; "
-          f"peak: ip={peak_ip:.4f} nodal={peak_nodal:.4f}")
-    explain = ("If the nodal path improved, explain it and re-characterise; "
-               "do not revert to make this pass.")
-    assert l2_ip <= l2_nodal, f"ip {l2_ip:.3e} > nodal {l2_nodal:.3e}. {explain}"
-    assert peak_ip >= peak_nodal, (
-        f"ip peak {peak_ip:.4f} < nodal {peak_nodal:.4f}. {explain}")
-
-
 def _unsteady_uniform_flow_check(kind, vform="var"):
     """Uniform velocity that changes linearly in time, v(t) = a + b t. The
     exact foot for the interval [t1, t1 + dt] is x - dt (a + b (t1 + dt/2)).

@@ -59,44 +59,6 @@ def test_constant_field_preserved_exactly():
     assert np.abs(h.data[:, 0] - 0.37).max() < 1.0e-12
 
 
-# tier_c overrides the module-level tier_a for this test alone: it compares two
-# smoothers, so it can fail because one of them got better.
-@pytest.mark.tier_c
-def test_taubin_against_plain_laplacian_characterisation():
-    """Characterisation: Taubin keeps the passband that plain Laplacian damps.
-
-    This compares two METHODS, so it can fail because the code improved — if
-    plain Laplacian gains a passband, this breaks and that is good news. It is
-    tier C for that reason: a failure demands an explanation, not a revert.
-
-    The contract this rests on is asserted separately and unconditionally in
-    `test_taubin_preserves_low_attenuates_high`: Taubin must preserve the low
-    mode and kill the high one, against fixed bounds and no rival method.
-
-    The absolute contract is asserted separately and unconditionally in
-    `test_taubin_preserves_low_attenuates_high`, against fixed bounds and no
-    rival method. Nothing absolute is asserted here, so nothing gating is lost
-    by this test being tier C.
-
-    Measured 2026-09-12 at n_iters=40, alpha=0.6: Taubin keeps 0.996 of the low
-    mode, plain Laplacian 0.917. The 0.03 margin characterises this fixture and
-    is not a specification.
-    """
-    surf_t, h_t, th = _surface_with_modes()
-    b_low = _amp(h_t.data[:, 0], th, 2)
-    uw.meshing.smooth_surface_field(h_t, n_iters=40, alpha=0.6, taubin=True)
-    taubin_low_kept = _amp(h_t.data[:, 0], th, 2) / b_low
-
-    surf_l, h_l, _ = _surface_with_modes()
-    uw.meshing.smooth_surface_field(h_l, n_iters=40, alpha=0.6, taubin=False)
-    laplacian_low_kept = _amp(h_l.data[:, 0], th, 2) / b_low
-
-    print(f"low mode kept: taubin={taubin_low_kept:.3f} "
-          f"laplacian={laplacian_low_kept:.3f}")
-    assert taubin_low_kept > laplacian_low_kept + 0.03, (
-        f"the passband gap closed: taubin={taubin_low_kept:.3f} vs "
-        f"laplacian={laplacian_low_kept:.3f}. If plain Laplacian improved, "
-        "explain it and re-characterise; do not revert to make this pass.")
 @pytest.mark.tier_a
 
 

@@ -147,27 +147,21 @@ def test_nitsche_constrains_the_wall_normal_velocity(self, solutions):
                 "scale — the free-slip constraint is not being applied"
             )
 
-    @pytest.mark.tier_c
-    def test_constraint_strength_ordering_characterisation(self, solutions):
-        """Characterisation: strong is exact, both weak paths leak; tier C.
-
-        This test CAN fail because the code got better, and that is why it is
-        tier C: a failure here demands an explanation, not a revert. It is not
-        a contract, and nothing should be reverted to make it pass.
-
-        It exists because the ordering is what the free-slip rulings rest on —
-        an essential BC (and a rotated strong free-slip) holds v.n to machine
-        precision, while Nitsche and penalty are weak constraints that leave a
-        finite leak. If a weak path starts coming out exact, or the strong path
-        stops being exact, the documented reasoning in
-        docs/developer/subsystems/rotated-freeslip.md needs revisiting and
-        somebody should say why.
-
-        The numbers are a characterisation of THIS fixture, measured
-        2026-09-12 at res=8: essential 0.0, penalty 1.5e-3, nitsche 5.7e-3
-        relative to the velocity scale. They are deliberately not a ranking —
-        which of the two weak methods leaks less is problem-dependent, moves
-        with gamma and the penalty coefficient, and is not asserted here.
+    @pytest.mark.tier_b
+    def test_essential_bc_holds_the_wall_normal_velocity_exactly(self, solutions):
+        """Contract: an essential BC holds v.n to machine precision.
+
+        Exactness is a hard baseline — a strong Dirichlet constraint eliminates
+        the degree of freedom, so the wall-normal velocity is zero to round-off
+        and any departure names a specific defect in how the BC is applied.
+
+        This is the half of the free-slip reasoning that can be asserted. The
+        other half — that the WEAK paths leave a finite leak — used to be
+        asserted alongside it as `leaks[method] > 1e-5`, i.e. that Nitsche and
+        penalty must STAY inaccurate. That assertion would fail if either method
+        improved, and a failure could not say which of the fixture, the gamma,
+        the penalty coefficient or the method had moved. The measured leaks are
+        printed instead: read them, do not gate on them.
         """
         leaks = {}
         for method in ("essential", "penalty", "nitsche"):
@@ -181,13 +175,6 @@ def test_constraint_strength_ordering_characterisation(self, solutions):
               + ", ".join(f"{k}={v:.3e}" for k, v in leaks.items()))
 
         assert leaks["essential"] < 1.0e-12, (
-            f"an essential BC is expected to hold v.n to machine precision; got "
-            f"{leaks['essential']:.3e}. This is the contract half of this test."
+            f"an essential BC must hold v.n to machine precision; got "
+            f"{leaks['essential']:.3e}"
         )
-        for method in ("penalty", "nitsche"):
-            assert leaks[method] > 1.0e-5, (
-                f"{method} leaked only {leaks[method]:.3e} — a weak constraint "
-                "reaching machine precision is GOOD NEWS and a change in "
-                "behaviour. Explain it and re-characterise; do not revert to "
-                "make this pass."
-            )