diff --git a/docs/developer/TESTING-RELIABILITY-SYSTEM.md b/docs/developer/TESTING-RELIABILITY-SYSTEM.md index 6f261aaf..4e8d072c 100644 --- a/docs/developer/TESTING-RELIABILITY-SYSTEM.md +++ b/docs/developer/TESTING-RELIABILITY-SYSTEM.md @@ -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 diff --git a/docs/developer/UW3_STYLE_CHARTER.md b/docs/developer/UW3_STYLE_CHARTER.md index 81a33671..8b574da4 100644 --- a/docs/developer/UW3_STYLE_CHARTER.md +++ b/docs/developer/UW3_STYLE_CHARTER.md @@ -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 diff --git a/src/underworld3/cython/petsc_generic_snes_solvers.pyx b/src/underworld3/cython/petsc_generic_snes_solvers.pyx index c7d53b68..8cbe26c5 100644 --- a/src/underworld3/cython/petsc_generic_snes_solvers.pyx +++ b/src/underworld3/cython/petsc_generic_snes_solvers.pyx @@ -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: diff --git a/tests/test_0010_cell_size_geometry.py b/tests/test_0010_cell_size_geometry.py index 9e955328..cd2de0d0 100644 --- a/tests/test_0010_cell_size_geometry.py +++ b/tests/test_0010_cell_size_geometry.py @@ -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" + ) diff --git a/tests/test_0066_integration_point_slcn.py b/tests/test_0066_integration_point_slcn.py index 14c2e070..9b16d09e 100644 --- a/tests/test_0066_integration_point_slcn.py +++ b/tests/test_0066_integration_point_slcn.py @@ -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)). diff --git a/tests/test_0773_surface_smoother.py b/tests/test_0773_surface_smoother.py index 09f3d071..9553d5a1 100644 --- a/tests/test_0773_surface_smoother.py +++ b/tests/test_0773_surface_smoother.py @@ -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 diff --git a/tests/test_1060_nitsche_freeslip.py b/tests/test_1060_nitsche_freeslip.py index f0de21ec..2f849d63 100644 --- a/tests/test_1060_nitsche_freeslip.py +++ b/tests/test_1060_nitsche_freeslip.py @@ -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." - ) diff --git a/tests/test_1070_free_surface_plume.py b/tests/test_1070_free_surface_plume.py index 44034a40..21225e50 100644 --- a/tests/test_1070_free_surface_plume.py +++ b/tests/test_1070_free_surface_plume.py @@ -304,30 +304,34 @@ def test_freesurface_strong_constraint_passes_no_net_flux(constraint_measurement f"strong constraint leaks net volume flux {leaks['strong']:.2e}" -# The tier goes on the test, not the module: pytest MERGES marks, so a module -# tier would remain on this item and `tier_a or tier_b` would still select it. @pytest.mark.level_2 -@pytest.mark.tier_c -def test_freesurface_strong_constraint_against_penalty(constraint_measurements): - r"""Characterisation: the strong constraint tracks the prescribed rate more - closely than the weak penalty. - - This compares two METHODS, so it can fail because the penalty path improved — - which would be good news. Tier C: a failure demands an explanation, not a - revert. It is NOT the justification for ``consistent_constraint="strong"``; - the contract that justifies it is the no-net-flux test above. - - Measured 2026-09-12 on the annulus fixture, 4 solve/advance steps: datum - error strong 1.06e-2 against penalty 3.01e-2. The 0.5 factor characterises - this fixture and is not a specification. +@pytest.mark.tier_b +def test_freesurface_strong_constraint_tracks_the_prescribed_rate(constraint_measurements): + """Contract: the strong constraint realises the prescribed wall-normal rate. + + The quantity here is already an absolute error against a known answer — the + datum `fs._un_target` that the constraint is asked to reproduce — so it is + asserted directly. + + It replaces a comparison against the penalty path + (`errors["strong"] < 0.5 * errors["penalty"]`), which was strictly less + informative than the number it was computed from: a ratio moves when either + side moves, so a failure could not say which constraint had changed. It also + could not fail when a default changed, only drift within its own margin and + break later somewhere harder to read — which is how the Nitsche penalty + regression in #734 stayed hidden for two months. + + The bound is loose enough to hold across platform triangulations (measured + 2026-09-14: 1.06e-2 on macOS, 1.73e-2 on the Linux CI runner, where the + annulus triangulates differently) and tight enough to fail if the rotated + constraint stops tracking the datum. """ errors, _ = constraint_measurements - print(f"datum error: strong={errors['strong']:.2e} " - f"penalty={errors['penalty']:.2e}") - assert errors["strong"] < 0.5 * errors["penalty"], ( - f"strong constraint not better: {errors['strong']:.2e} vs " - f"penalty {errors['penalty']:.2e}. If the penalty path improved, " - "explain it and re-characterise; do not revert to make this pass.") + print(f"datum error: strong={errors['strong']:.2e} penalty={errors['penalty']:.2e}") + assert errors["strong"] < 5.0e-2, ( + f"the strong constraint misses the prescribed wall-normal rate by " + f"{errors['strong']:.2e}" + ) @pytest.mark.level_1