diff --git a/tests/validation/multigrid/test_exchange_window.py b/tests/validation/multigrid/test_exchange_window.py index 8776e01..84d19a4 100644 --- a/tests/validation/multigrid/test_exchange_window.py +++ b/tests/validation/multigrid/test_exchange_window.py @@ -1,10 +1,12 @@ """Optional observation-window exchange ledger: host contract (mocked driver) and GPU gate. Run CPU and GPU classes in separate processes; GPU callers select the device. -GPU tolerance: window rows are summed per link over time, the per-step reads over links -then over time. Reordering a sum of m FP64 terms changes it by at most -2*(m-1)*eps*sum|terms| (each order is within (m-1)*eps*sum|terms| of the exact sum), so the -bound is scaled by the absolute sum of the per-step row integrals and is robust to cancellation. +GPU tolerance (conservative engineering tolerance, not a universal proof): the window ledger +sums h*dt*raw per link over fine substeps, while the reference sums per-base interval ledgers +over links and then over base steps. Both are FP64 reorderings of the same FP32-exact substep +contributions, so the scale is the sum of the ACTUAL absolute substep contributions (read from +last_raw after every fine hop), not of per-base totals that could hide substep cancellation: +|window - reference| <= 4 * m * eps * sum|substep terms| + 1e-300, m = links * substeps. """ from contextlib import ExitStack import unittest @@ -26,6 +28,16 @@ def build(m, window=True, raw=True): observe_raw_exchange=raw, observe_window_exchange=window) +def same_bits(test, a, b): + """Identical NaN masks, no infinities, identical finite FP32 bits (incl. signed zero).""" + a, b = np.asarray(a), np.asarray(b) + test.assertEqual(a.dtype, b.dtype) + test.assertFalse(np.isinf(a).any() or np.isinf(b).any()) + np.testing.assert_array_equal(np.isnan(a), np.isnan(b)) + finite = np.isfinite(a) + np.testing.assert_array_equal(a[finite].view(np.uint32), b[finite].view(np.uint32)) + + def accumulations(m): return [c for c in m["calls"] if c[0] == "BlockAccumulateExchange"] @@ -107,9 +119,22 @@ class ExchangeWindowGPU(unittest.TestCase): windowed.start_exchange_window() b = windowed._level_boundaries[1] count = len(b.tables["cell"]) + factors = np.array([.25, .25, .125, .25]) # h*dt, h*h*dt, h*h for level-1 h=.5 + scale = np.zeros(4) + substeps = 0 + kernel = windowed._levels[1].kernel + + def audited(src, dst, *args, **kw): + nonlocal scale, substeps + windowed._stream.synchronize() # BFL and both accumulations are complete + raw = np.empty((count, 4), np.float32) + cuda.memcpy_dtoh(raw, b.last_raw) + scale = scale + (np.abs(raw.astype(np.float64)) * factors).sum(axis=0) + substeps += 1 + kernel(src, dst, *args, **kw) + windowed._levels[1].kernel = audited keys = ("impulse", "angular_impulse", "mass_into_fluid") summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")} - scale = np.zeros(4) steps = 0 for speed in (0., .002, .002, -.001, 0., 0.): # amplitude updates stay in-window windowed.set_jet_amplitude(0, speed) @@ -123,16 +148,13 @@ class ExchangeWindowGPU(unittest.TestCase): summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1], summed[part][2] + o["angular_impulse"], summed[part][3] + o["mass_into_fluid"]] - rows = np.empty((count, 4), np.float64) - cuda.memcpy_dtoh(rows, b.interval_integral) - scale += np.abs(rows).sum(axis=0) for a, e in zip(windowed._download_raw().values, plain._download_raw().values): - np.testing.assert_array_equal(np.isnan(a), np.isnan(e)) - np.testing.assert_array_equal(a[np.isfinite(a)], e[np.isfinite(e)]) + same_bits(self, a, e) out = windowed.read_exchange_window(0) self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]), (steps, 2 * steps)) - bound = 2 * (count * steps) * EPS * scale + np.finfo(np.float64).tiny + self.assertEqual(substeps, 2 * steps) + bound = 4 * (count * substeps) * EPS * scale + 1e-300 for part in summed: o = out[part] got = np.array([o["impulse"][0], o["impulse"][1], o["angular_impulse"],