From 0613a48f168827bb8cdb3e3f1260a2fbcf45e11d Mon Sep 17 00:00:00 2001 From: lmoresi Date: Thu, 27 Aug 2026 21:31:21 +1000 Subject: [PATCH 1/3] Pin the DIRECTION of the spherical surface response (#496, partial) test_1072's surface assertion is symmetric in sign: assert 0.0 < np.abs(shape).max() <= 1.5 * np.abs(h_inf).max() Flip the 3-D recovery and |shape| is unchanged, so it passes while the surface moves the wrong way. #496 records exactly that: CI cannot catch a 3-D sign error. Starting from flat, the displacement IS the shape, so it must correlate POSITIVELY with the equilibrium it is relaxing toward. Measured +0.995; a sign error gives about -0.995, which the magnitude bound cannot see. Asserted at > 0.9. PARTIAL, and the docstring now says so. The other half of #496 -- the spherical relaxation RATE -- is not pinned. The review measured a fitted decay rate of 0.0842 against a half-space Cathles rate rho*g/(2*eta*k), k = sqrt(6)/R, of 0.204: ratio 0.41, an O(1) shell correction below the half-space value and in the physically correct direction. Asserting that band needs an initial-topography relaxation setup, and no test in the suite has one -- every free-surface test starts flat and is driven by a load. Building it hastily would risk a test that passes for the wrong reason, which is the failure #496 exists to fix, so it is left for its own session with the numbers recorded in the docstring. The mean-free assertion above is also weak (it applies _surface_mean to a field _demean'ed with the same weights, so it is close to true by construction), as #496 notes. Left alone here: replacing it properly means checking the gauge against an independently computed surface integral, which belongs with the rate work. Underworld development team with AI support from Claude Code --- tests/test_1072_free_surface_spherical.py | 28 ++++++++++++++++++++--- 1 file changed, 25 insertions(+), 3 deletions(-) diff --git a/tests/test_1072_free_surface_spherical.py b/tests/test_1072_free_surface_spherical.py index 7e83d48a..799e6709 100644 --- a/tests/test_1072_free_surface_spherical.py +++ b/tests/test_1072_free_surface_spherical.py @@ -4,8 +4,18 @@ The quantitative benchmarking (analytic Y_lm shell rate, convergence of the h_inf modal bias, 3D parallel) is the review-team's scope — this test pins the -CAPABILITY: construction, one solve/advance cycle, finite mean-free h_inf, and -the explicit refusal of the 2D-only features. +CAPABILITY: construction, one solve/advance cycle, finite mean-free h_inf, the +DIRECTION of the surface response, and the explicit refusal of the 2D-only +features. + +Still NOT pinned here, and the remaining half of #496: the spherical relaxation +RATE. Measured in the review at cell size 0.35 with Y20 initial topography, the +fitted decay rate was 0.0842 against a half-space Cathles rate rho*g/(2*eta*k), +k = sqrt(6)/R, of 0.204 — a ratio of 0.41, an O(1) shell correction below the +half-space value and in the physically correct direction. Asserting that band +needs an initial-topography relaxation setup, which no test here has: every +free-surface test starts flat and is driven by a load. A rate regression is +therefore still invisible to CI. """ import numpy as np import pytest @@ -44,11 +54,23 @@ def test_freesurface_spherical_shell_end_to_end(): "h_inf datum is not mean-free under the trace-mass gauge" assert np.abs(h_inf).max() > 1.0e-4, "no topographic response to the load" fs.advance(fs.estimate_dt(advect_scale=10.0)) - shape = fs._current_shape() + shape = np.asarray(fs._current_shape()) assert np.isfinite(shape).all() assert 0.0 < np.abs(shape).max() <= 1.5 * np.abs(h_inf).max(), \ "surface did not move toward (or overshot) equilibrium" + # DIRECTION, not just magnitude. The bound above is symmetric in sign: flip + # the 3-D recovery and |shape| is unchanged, so it passes while the surface + # moves the wrong way (#496). Starting from flat, the displacement IS the + # shape, so it must correlate POSITIVELY with the equilibrium it is moving + # toward. Measured +0.995 here; a sign error gives about -0.995. + interesting = np.abs(h_inf) > 1.0e-12 + direction = float(np.corrcoef(shape.ravel()[interesting.ravel()], + np.asarray(h_inf).ravel()[interesting.ravel()])[0, 1]) + assert direction > 0.9, ( + f"surface moved away from equilibrium (corr {direction:+.3f}) — the " + "magnitude bound above cannot see a sign error") + def test_freesurface_spherical_refuses_2d_only_features(): """The 2D-only features fail loudly at construction in 3D, not silently.""" From 6c46c7832521d85a104c741c5799b5a55d42274f Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 28 Aug 2026 14:50:12 +1000 Subject: [PATCH 2/3] Pin the spherical relaxation RATE, and give the datum check teeth (#496) Closes the second half of #496. Every other free-surface test starts flat and is driven by a load, so nothing in the suite exercised the decay of an imposed topography: a rate regression was invisible to CI. This adds a Y20 relaxation on a constant-density shell under radial gravity, where the topographic self-load is the only driver, and asserts the fitted decay rate against the half-space Cathles rate rho*g/(2*eta*k), k = sqrt(6)/R. Measured ratio 0.528 at cell 0.35 and 0.516 at cell 0.25; np=1 and np=2 agree to five digits. Asserted band 0.45 to 0.62. The band catches a sign error, an order-of-magnitude error, and the half-space value itself being returned; it does NOT catch a 2x error in buoyancy_scale (measured 0.476), and the assertion message says so rather than implying more coverage than it has. The decay lands on a floor - the #431-class bias in the recovered h_inf, 16% of the initial amplitude at cell 0.35 and 11% at 0.25 - so the extracted rate depends entirely on how the floor is treated: the same record gives 0.064 with no floor, 0.110 with the floor free, 0.142 with the floor pinned to the last sample. The protocol (floor free) is therefore fixed in code and documented. This also means the docstring's previous 0.0842 / ratio 0.41 was not reproducible without its protocol, and that the earlier "the rate is not yet converged" reading was mostly the naive fit reacting to a changing floor. Also replaces this file's mean-free assertion, which had no content: free_surface.py:901 builds _h_inf as _demean(-_demean(h)) and _demean IS subtraction of _surface_mean, so the asserted quantity is zero by construction - measured 1e-18 against a 1e-8 threshold - and it stays zero even if the trace-mass weights are wrong, because both sides use the same weights. It is now the FE boundary integral of the P1 carrier, an independent implementation of the same integral, plus a negative control asserting the UNWEIGHTED nodal mean is NOT zero (2.5e-4 of max|h_inf|), so the check cannot pass by construction. The rank-local max and shape bounds in the end-to-end test are now reduced globally; they were rank-local diagnostics that could disagree across ranks. Measurements, negative controls, and the surface-array ordering bug that had to be fixed in the probe first (matching _field_rows against the LIVE _ring_coords after a deform re-sorts by deformed x and silently permutes the carrier fill; free_surface.py:211 warns of exactly this) are in ~/+Simulations/spherical_relaxation_rate_496/README.md. That bug is also present in the earlier ylm_relaxation.py probe the workstream-C numbers came from. Serial and np=2: 3 passed, ~55s. Underworld development team with AI support from Claude Code --- tests/test_1072_free_surface_spherical.py | 195 +++++++++++++++++++--- 1 file changed, 175 insertions(+), 20 deletions(-) diff --git a/tests/test_1072_free_surface_spherical.py b/tests/test_1072_free_surface_spherical.py index 799e6709..85a7b13d 100644 --- a/tests/test_1072_free_surface_spherical.py +++ b/tests/test_1072_free_surface_spherical.py @@ -2,22 +2,26 @@ physically sensible topography (guards the dimension-general surface machinery: owned-facet trace-mass gauge, P1-projected sigma_nn recovery, radial deform). -The quantitative benchmarking (analytic Y_lm shell rate, convergence of the -h_inf modal bias, 3D parallel) is the review-team's scope — this test pins the -CAPABILITY: construction, one solve/advance cycle, finite mean-free h_inf, the -DIRECTION of the surface response, and the explicit refusal of the 2D-only -features. - -Still NOT pinned here, and the remaining half of #496: the spherical relaxation -RATE. Measured in the review at cell size 0.35 with Y20 initial topography, the -fitted decay rate was 0.0842 against a half-space Cathles rate rho*g/(2*eta*k), -k = sqrt(6)/R, of 0.204 — a ratio of 0.41, an O(1) shell correction below the -half-space value and in the physically correct direction. Asserting that band -needs an initial-topography relaxation setup, which no test here has: every -free-surface test starts flat and is driven by a load. A rate regression is -therefore still invisible to CI. +Pinned here: construction, one solve/advance cycle, a finite and genuinely +volume-preserving h_inf, the DIRECTION of the surface response, the explicit +refusal of the 2D-only features, and the relaxation RATE of an initial Y20 +topography against the half-space Cathles rate. + +The rate test closes the second half of #496. It is a REGRESSION pin at one +fixed resolution: the ratio measured 0.528 at cell 0.35 and 0.516 at cell 0.25, +so it is close to converged, but two coarse resolutions do not establish that. +The decay lands on a floor -- 16% of the initial amplitude here, 11% at cell +0.25 -- which is the resolution-convergent #431-class bias in the recovered +h_inf, not part of the physics. Because the decay is onto a floor rather than to +zero, the FIT PROTOCOL is part of the specification: the same record yields +rates of 0.064, 0.110 or 0.142 depending on how the floor is treated. It is +fixed in _fit_relaxation and must not be changed casually. + +Measurements, negative controls and the ordering bug that had to be fixed first: +~/+Simulations/spherical_relaxation_rate_496/README.md. """ import numpy as np +from mpi4py import MPI import pytest import sympy import underworld3 as uw @@ -41,22 +45,70 @@ def _shell_stokes(cell=0.35): return mesh, stokes, rhat +def _surface_carrier(fs, mesh, name): + """A P1 field on the surface plus the rows that hold a surface array. + + Ordering trap, and the reason the rows come from ``_surf_coords``: every + surface ARRAY (``_current_shape``, ``_h_inf``, the increment ``_carry_and_deform`` + consumes) is indexed in ``_surf_rows`` order, which is ``_surf_coords`` sorted + by REFERENCE x. ``_field_rows`` re-sorts whatever coordinates it is handed, so + matching against ``_ring_coords`` -- documented at free_surface.py:211 as LIVE, + deformed geometry -- sorts by DEFORMED x instead and silently permutes the fill. + """ + field = uw.discretisation.MeshVariable(name, mesh, 1, degree=1, continuous=True) + rows, _ = fs._field_rows(field, fs._surf_coords) + return field, rows + + def test_freesurface_spherical_shell_end_to_end(): """Construction + one full solve/advance on the shell; h_inf finite and - mean-free; the surface responds toward equilibrium (|h| grows from flat - under the one-sided load and stays bounded by |h_inf|).""" + volume-preserving; the surface responds toward equilibrium (|h| grows from + flat under the one-sided load and stays bounded by |h_inf|).""" mesh, stokes, rhat = _shell_stokes() fs = uw.systems.FreeSurface(stokes, "Upper", buoyancy_scale=50.0, normal=rhat) fs.solve() h_inf = np.asarray(fs._h_inf) assert np.isfinite(h_inf).all(), "3D h_inf recovery produced non-finite values" - assert abs(fs._surface_mean(h_inf)) < 1.0e-8 * (np.abs(h_inf).max() + 1e-30), \ - "h_inf datum is not mean-free under the trace-mass gauge" - assert np.abs(h_inf).max() > 1.0e-4, "no topographic response to the load" + scale = uw.mpi.comm.allreduce( + float(np.abs(h_inf).max()) if h_inf.size else 0.0, op=MPI.MAX) + assert scale > 1.0e-4, "no topographic response to the load" + + # The datum must be VOLUME-preserving, measured by an integrator that does not + # share the gauge's own weight vector. Asserting _surface_mean(h_inf) instead + # is vacuous: free_surface.py:901 builds h_inf as _demean(-_demean(h)) and + # _demean IS subtraction of _surface_mean, so the quantity is zero by + # construction -- measured 1e-18 relative against a 1e-8 threshold, and it + # stays zero even if the trace-mass weights themselves are wrong, because both + # sides use the same weights. The FE boundary integral of the P1 carrier is an + # independent implementation of the same integral and does catch that. + h_field, h_rows = _surface_carrier(fs, mesh, "hinfP1") + h_field.array[...] = 0.0 + h_field.array[h_rows, 0, 0] = h_inf + area = uw.maths.BdIntegral(mesh=mesh, fn=sympy.sympify(1), boundary="Upper") + flux = uw.maths.BdIntegral(mesh=mesh, fn=h_field.sym[0], boundary="Upper") + datum = float(flux.evaluate()) / float(area.evaluate()) + assert abs(datum) < 1.0e-12 * scale, ( + f"h_inf datum is not volume-preserving: = {datum:.3e} over the " + f"boundary, {abs(datum) / scale:.2e} of max|h_inf|") + + # NEGATIVE CONTROL for the line above: the trace-mass weighting must actually + # do work here, or "the weighted mean vanishes" says nothing. The UNWEIGHTED + # nodal mean of the same data is not zero (measured 2.5e-4 of max|h_inf| at + # this resolution) -- surface nodes are not equally spaced, so a node-count + # gauge would not be volume-preserving. + nodal = (uw.mpi.comm.allreduce(float(h_inf.sum()) if h_inf.size else 0.0) + / max(uw.mpi.comm.allreduce(int(h_inf.size)), 1)) + assert abs(nodal) > 1.0e-6 * scale, ( + f"the unweighted nodal mean is also ~zero ({abs(nodal) / scale:.2e} of " + "max|h_inf|), so the volume-preservation check above has no content on " + "this mesh -- it cannot distinguish the trace-mass gauge from node counting") + fs.advance(fs.estimate_dt(advect_scale=10.0)) shape = np.asarray(fs._current_shape()) assert np.isfinite(shape).all() - assert 0.0 < np.abs(shape).max() <= 1.5 * np.abs(h_inf).max(), \ + shape_max = uw.mpi.comm.allreduce( + float(np.abs(shape).max()) if shape.size else 0.0, op=MPI.MAX) + assert 0.0 < shape_max <= 1.5 * scale, \ "surface did not move toward (or overshot) equilibrium" # DIRECTION, not just magnitude. The bound above is symmetric in sign: flip @@ -72,6 +124,109 @@ def test_freesurface_spherical_shell_end_to_end(): "magnitude bound above cannot see a sign error") +def _fit_relaxation(t, A): + r"""Fit :math:`A(t) = A_\infty + A_1 e^{-\lambda t}` and return + ``(lambda, A_inf, residual_rms)``. + + The floor :math:`A_\infty` is a FREE parameter, and that is the whole protocol. + The decay does not run to zero -- it lands on the Y20 bias in the recovered + h_inf -- so the extracted rate depends entirely on how the floor is treated: + on this same record, fitting with no floor gives 0.064, floor free 0.108, and + floor pinned to the last sample 0.142. Any of the three can be called "the + fitted decay rate"; only one of them is a stable number, so the test fixes it. + """ + from scipy.optimize import curve_fit + + def model(tt, A_inf, A_1, lam): + return A_inf + A_1 * np.exp(-lam * tt) + + p, _ = curve_fit(model, t, A, p0=[A[-1] * 0.5, A[0], 0.1], maxfev=200000) + return p[2], p[0], float((A - model(t, *p)).std()) + + +def test_freesurface_spherical_relaxation_rate(): + """An initial Y20 topography must relax at the shell rate (#496, second half). + + Every other free-surface test starts FLAT and is driven by a load, so nothing + in the suite exercises the decay of an imposed topography and a rate + regression is invisible. Here a constant-density shell under radial gravity + carries a degree-2 bump and nothing else: the topographic self-load is the + only driver, so the modal amplitude must decay exponentially at a rate set by + the Stokes solve. + """ + r_out, r_in, eps, dt, nsteps = 1.0, 0.547, 0.01, 1.5, 8 + rho_g = eta = 1.0 + + mesh = uw.meshing.SphericalShell(radiusOuter=r_out, radiusInner=r_in, + cellSize=0.35, qdegree=3) + x, y, z = mesh.X + r = sympy.sqrt(x ** 2 + y ** 2 + z ** 2) + rhat = sympy.Matrix([[x / r, y / r, z / r]]) + stokes = uw.systems.Stokes(mesh) + stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel + stokes.constitutive_model.Parameters.shear_viscosity_0 = eta + stokes.bodyforce = -rho_g * rhat.T + stokes.add_essential_bc((0.0, 0.0, 0.0), "Lower") + stokes.tolerance = 1.0e-6 + fs = uw.systems.FreeSurface(stokes, "Upper", buoyancy_scale=rho_g, + normal=rhat, background_buoyancy="analytic") + + h_field, h_rows = _surface_carrier(fs, mesh, "hY20") + Y20 = 0.5 * (3.0 * (z / r) ** 2 - 1.0) + num = uw.maths.BdIntegral(mesh=mesh, fn=h_field.sym[0] * Y20, boundary="Upper") + den = uw.maths.BdIntegral(mesh=mesh, fn=Y20 ** 2, boundary="Upper") + + def modal_amplitude(): + h_field.array[...] = 0.0 + h_field.array[h_rows, 0, 0] = fs._current_shape() + return float(num.evaluate()) / float(den.evaluate()) + + # impose the bump through the manager's own carrier, so the initial state is + # exactly a valid FS state (smooth interior decay + radial deform) + sc = fs._ring_coords + y20 = 0.5 * (3.0 * (sc[:, 2] / np.linalg.norm(sc, axis=1)) ** 2 - 1.0) + fs._carry_and_deform(eps * y20, dt=0.0) + + ts, As = [0.0], [modal_amplitude()] + assert As[0] > 0.5 * eps, \ + f"the imposed Y20 bump did not register on the surface (A0 = {As[0]:.3e})" + for _ in range(nsteps): + fs.solve() + fs.advance(dt) # FIXED step: estimate_dt is CFL-based and is + ts.append(ts[-1] + dt) # ~10 tau here, which hops the whole decay + As.append(modal_amplitude()) + t, A = np.array(ts), np.array(As) + + rate, floor, resid = _fit_relaxation(t, A) + + # A rate is only meaningful if the decay IS an exponential; a fit rammed + # through a non-exponential record would otherwise report a number and pass. + assert resid < 5.0e-3 * A[0], ( + f"the decay is not a clean exponential (residual rms {resid:.2e} = " + f"{resid / A[0] * 100:.2f}% of A0) — the fitted rate is meaningless") + assert 0.0 < floor < 0.4 * A[0], ( + f"h_inf Y20 bias floor {floor / A[0] * 100:.1f}% of A0 is out of band; " + "this is the #431-class recovery defect and should shrink with resolution") + + # The shell over a no-slip base must relax SLOWER than a half-space (finite + # depth, fixed bottom) but on the same order -- an O(1) correction below 1. + cathles = rho_g / (2.0 * eta * np.sqrt(2.0 * 3.0) / r_out) # k = sqrt(l(l+1))/R + ratio = rate / cathles + assert 0.45 < ratio < 0.62, ( + f"Y20 relaxation rate {rate:.4f} is {ratio:.3f} of the half-space Cathles " + f"rate {cathles:.4f}; expected an O(1) shell correction below 1 (finite " + "depth over a no-slip base relaxes slower than a half-space).\n" + "Measured: 0.528 at cell 0.35 and 0.516 at cell 0.25 (so the ratio is " + "nearly resolution-converged, unlike the floor); 0.518 to 0.552 as the " + "fitted record grows from 5 to 20 steps; np=1 and np=2 agree to 5 digits.\n" + "What this band catches: a sign error (growth, not decay), an order-of-" + "magnitude error, and the half-space value itself (1.0) being returned. " + "What it does NOT catch: a 2x error in buoyancy_scale, measured at 0.476 " + "-- the exponential update is only weakly sensitive to it, and a band " + "tight enough to see it would sit within a few percent of the resolution " + "spread and flake.") + + def test_freesurface_spherical_refuses_2d_only_features(): """The 2D-only features fail loudly at construction in 3D, not silently.""" mesh, stokes, rhat = _shell_stokes() From 713fc40f6ce6943ee865f161459311441c5969bc Mon Sep 17 00:00:00 2001 From: lmoresi Date: Fri, 28 Aug 2026 15:45:30 +1000 Subject: [PATCH 3/3] Actually run test_1072 in CI, and name the mesh when the rate band fails The rate test added in the previous commit does not run in CI. scripts/test.sh excludes test_106*py and test_107*py pending the #504 triage, so test_1072_free_surface_spherical.py has never been executed there - which means BOTH halves of #496 were still invisible to CI, the sign half included, and merging the previous commit would have closed the issue without closing the gap it exists to close. Pulls that one file forward out of the excluded group, the same way tests/parallel/test_1069_boundary_normal_parallel.py already is, and for the same kind of reason. It is level_2/tier_b, ~55s serial, and passes at np=1 and np=2. The rest of the group still awaits its triage decision. Also names the mesh in the band-failure message. gmsh triangulates differently across platforms - serial_reference.py:202 records 1417 vs 1395 cells for a spherical shell between hosts - so the Linux runner will not build the 585-cell mesh the band was measured on. The resolution study makes a failure unlikely (cell 0.35 -> 0.25, a far larger change, moved the ratio only 2.3%), but if the band does trip on another host, the message now carries the cell count, A0 and the floor, which is what distinguishes a platform mesh difference from a real regression. Underworld development team with AI support from Claude Code --- scripts/test.sh | 8 ++++++++ tests/test_1072_free_surface_spherical.py | 9 +++++++++ 2 files changed, 17 insertions(+) diff --git a/scripts/test.sh b/scripts/test.sh index 93a69edf..30c37913 100755 --- a/scripts/test.sh +++ b/scripts/test.sh @@ -116,6 +116,14 @@ if [ $PARALLEL_ONLY -eq 0 ]; then # NOT yet batched (issue #504 audit): test_106*py and test_107*py contain # level_2/level_3 + slow + tier_b/tier_c suites (e.g. test_1064) and need # a triage/deselect decision before being wired into CI. + # + # test_1072 is pulled forward out of that group, the same way test_1069 is + # above, because leaving it there defeats its purpose: it is the ONLY guard on + # the 3-D free-surface sign and relaxation rate, and #496 exists precisely + # because those regressions were invisible to CI. Landing the test into an + # unbatched file would have closed the issue without closing the gap. + # level_2/tier_b, ~55s serial; passes at np=1 and np=2. + $PYTEST tests/test_1072_free_surface_spherical.py || status=1 # Diffusion / Advection tests $PYTEST tests/test_1100*py || status=1 diff --git a/tests/test_1072_free_surface_spherical.py b/tests/test_1072_free_surface_spherical.py index 85a7b13d..f5951aa7 100644 --- a/tests/test_1072_free_surface_spherical.py +++ b/tests/test_1072_free_surface_spherical.py @@ -199,6 +199,14 @@ def modal_amplitude(): rate, floor, resid = _fit_relaxation(t, A) + # gmsh triangulates differently across platforms (serial_reference.py:202 + # records e.g. 1417 vs 1395 cells for a spherical shell), so a band failure + # has to name the mesh it was measured on or it is not diagnosable. + # Reference: 585 cells, macOS/arm64. + ncells = uw.mpi.comm.allreduce(mesh.dm.getStratumSize("depth", mesh.dim)) + where = (f"[{ncells} cells, A0 = {A[0]:.4e}, " + f"floor {floor / A[0] * 100:.1f}% of A0]") + # A rate is only meaningful if the decay IS an exponential; a fit rammed # through a non-exponential record would otherwise report a number and pass. assert resid < 5.0e-3 * A[0], ( @@ -214,6 +222,7 @@ def modal_amplitude(): ratio = rate / cathles assert 0.45 < ratio < 0.62, ( f"Y20 relaxation rate {rate:.4f} is {ratio:.3f} of the half-space Cathles " + f"rate, on this mesh {where}.\n" f"rate {cathles:.4f}; expected an O(1) shell correction below 1 (finite " "depth over a no-slip base relaxes slower than a half-space).\n" "Measured: 0.528 at cell 0.35 and 0.516 at cell 0.25 (so the ratio is "