test: annulus mass ledger definition with synthetic seam/edge checks

Host-only stdlib definition of the per-hop identity dM + Q_exact = C + H with Q_exact from pre-staging owner populations and H = sum over every annulus cell/direction of (pulled - ideal owner input), split into seam (internal) and edge parts. Synthetic tests: ideal pull, collision source, interior-seam and edge perturbations, the double count of staged flux + full H, a late owner snapshot hiding staging, and a wrong pull direction. No solver change, no GPU. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-20 18:32:42 +00:00
2 changed files with 176 additions and 0 deletions
@@ -0,0 +1,101 @@
+"""Synthetic checks of the annulus mass ledger. Host-only, stdlib; not solver evidence."""
+
+import random
+import sys
+from pathlib import Path
+import unittest
+
+sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
+import annulus_ledger as al  # noqa: E402
+
+NX = NY = 24
+# Annulus 4 <= r <= 8 around (12, 12); dyadic values keep every sum exact.
+A = {(x, y) for x in range(NX) for y in range(NY) if 16 <= (x - 12) ** 2 + (y - 12) ** 2 <= 64}
+DELTA = 0.375
+
+
+def owner(seed=1):
+    rng = random.Random(seed)
+    return [[[rng.randint(1, 1024) / 1024 for _ in range(NX)] for _ in range(NY)] for _ in range(9)]
+
+
+def collide(P, source=None):
+    """Mass-conserving stand-in: rotate populations within each cell (+ optional source)."""
+    D = [[[P[(q + 1) % 9][y][x] for x in range(NX)] for y in range(NY)] for q in range(9)]
+    if source:
+        (x, y), s = source
+        D[0][y][x] += s
+    return D
+
+
+def link(interior):
+    """(x, y, q) with x in A and x - c_q inside (seam) or outside (edge) A."""
+    for (x, y) in sorted(A):
+        for q in range(1, 9):
+            src = ((x - al.CX[q]) % NX, (y - al.CY[q]) % NY)
+            if (src in A) == interior:
+                return x, y, q
+    raise AssertionError
+
+
+class Ledger(unittest.TestCase):
+    def setUp(self):
+        self.O = owner()
+        self.P = al.pull(self.O, NX, NY)
+
+    def run_ledger(self, P, D):
+        return al.ledger(self.O, P, D, A, NX, NY)
+
+    def test_ideal_pull_has_no_h_and_closes(self):
+        r = self.run_ledger(self.P, collide(self.P))
+        self.assertEqual((r["H_seam"], r["H_edge"], r["C"], r["residual"]), (0, 0, 0, 0))
+        self.assertEqual(r["dM"], -r["Q_exact"])
+        self.assertNotEqual(r["Q_exact"], 0)
+
+    def test_collision_source_is_c(self):
+        x, y = sorted(A)[3]
+        r = self.run_ledger(self.P, collide(self.P, ((x, y), DELTA)))
+        self.assertEqual((r["C"], r["H"], r["residual"]), (DELTA, 0, 0))
+
+    def test_interior_seam_perturbation(self):
+        x, y, q = link(interior=True)
+        self.P[q][y][x] += DELTA
+        r = self.run_ledger(self.P, collide(self.P))
+        self.assertEqual((r["H_seam"], r["H_edge"], r["residual"]), (DELTA, 0, 0))
+        # Invisible to any flux through the edge of A; only H explains the mass change.
+        self.assertEqual(al.staged_flux(self.O, self.P, A, NX, NY), r["Q_exact"])
+        self.assertEqual(r["dM"] + r["Q_exact"], DELTA)
+
+    def test_edge_perturbation_and_staged_flux_double_count(self):
+        x, y, q = link(interior=False)
+        self.P[q][y][x] += DELTA
+        r = self.run_ledger(self.P, collide(self.P))
+        self.assertEqual((r["H_seam"], r["H_edge"], r["residual"]), (0, DELTA, 0))
+        staged = al.staged_flux(self.O, self.P, A, NX, NY)
+        self.assertEqual(staged, r["Q_exact"] - DELTA)
+        # Staged flux + full H counts the staging change twice.
+        self.assertEqual(r["dM"] + staged - r["C"] - r["H"], -DELTA)
+        # Staged flux alone closes only without the edge part of H.
+        self.assertEqual(r["dM"] + staged - r["C"] - r["H_seam"], 0)
+
+    def test_owner_snapshot_after_staging_hides_the_change(self):
+        # Snapshot taken after staging wrote the same slot into the owner image.
+        x, y, q = link(interior=False)
+        sx, sy = (x - al.CX[q]) % NX, (y - al.CY[q]) % NY
+        self.P[q][y][x] += DELTA
+        late = [[row[:] for row in plane] for plane in self.O]
+        late[q][sy][sx] += DELTA
+        r = al.ledger(late, self.P, collide(self.P), A, NX, NY)
+        self.assertEqual((r["H"], r["residual"]), (0, 0))
+
+    def test_wrong_pull_direction_is_flagged(self):
+        # Indexing error (push instead of pull): H is nonzero with no physical cause.
+        bad = [[[self.O[q][(y + al.CY[q]) % NY][(x + al.CX[q]) % NX] for x in range(NX)]
+                for y in range(NY)] for q in range(9)]
+        r = self.run_ledger(bad, collide(bad))
+        self.assertNotEqual(r["H"], 0)
+        self.assertEqual(r["residual"], 0)
+
+
+if __name__ == "__main__":
+    unittest.main()
@@ -0,0 +1,75 @@
+"""Per-hop mass ledger of a fluid cell set A (host-only, stdlib, D2Q9 pull scheme).
+
+Diagnostic definition only; no solver code, no GPU. Arrays are indexed [q][y][x] on a
+periodic nx x ny lattice; A is a set of (x, y) ACTIVE cells.
+
+  O  owner populations before any staging (canonical owner slots, start of the hop)
+  P  populations each cell actually pulls: P[q](x) is what x reads for direction q
+     (after halo copy, BFL/coalesce staging, ...)
+  D  populations written by the cell after collision
+
+With M0 = sum_{A,q} O and M1 = sum_{A,q} D (both owner snapshots of A):
+
+  M1 - M0 + Q_exact = C + H                                            (exact)
+
+  C       = sum_{x in A, q} (D_q(x) - P_q(x))              collision source
+  H       = sum_{x in A, q} (P_q(x) - O_q(x - c_q))        actual minus ideal owner input,
+            split into H_seam (x - c_q in A: internal pulls, e.g. tile seams) and
+            H_edge (x - c_q not in A: inputs crossing the edge of A)
+  Q_exact = sum_{y in A, y+c_q not in A} O_q(y) - sum_{y not in A, y+c_q in A} O_q(y)
+            owner-based net outward flux (outgoing minus incoming)
+
+Proof: sum_{A,q} O_q(x - c_q) - M0 = -Q_exact (relabel y = x - c_q per direction), and
+sum_{A,q} D = C + H + sum_{A,q} O_q(x - c_q). Q_exact must use O, not staged values: a
+flux built from staged inputs already contains H_edge, so adding the full H counts it
+twice (see staged_flux). The edge of A must not cross cut links, otherwise O at the solid
+end enters Q_exact.
+"""
+
+CX = (0, 1, -1, 0, 0, 1, -1, 1, -1)
+CY = (0, 0, 0, 1, -1, 1, -1, -1, 1)
+
+
+def _src(x, y, q, nx, ny):
+    return (x - CX[q]) % nx, (y - CY[q]) % ny
+
+
+def pull(O, nx, ny):
+    """Ideal pulled input: P[q](x) = O[q](x - c_q)."""
+    return [[[O[q][_src(x, y, q, nx, ny)[1]][_src(x, y, q, nx, ny)[0]]
+              for x in range(nx)] for y in range(ny)] for q in range(9)]
+
+
+def ledger(O, P, D, A, nx, ny):
+    A = set(A)
+    out = dict(dM=0.0, C=0.0, H_seam=0.0, H_edge=0.0, Q_exact=0.0)
+    for (x, y) in A:
+        for q in range(9):
+            out["dM"] += D[q][y][x] - O[q][y][x]
+            out["C"] += D[q][y][x] - P[q][y][x]
+            sx, sy = _src(x, y, q, nx, ny)
+            key = "H_seam" if (sx, sy) in A else "H_edge"
+            out[key] += P[q][y][x] - O[q][sy][sx]
+            if (sx, sy) not in A:
+                out["Q_exact"] -= O[q][sy][sx]          # incoming owner population
+            tx, ty = (x + CX[q]) % nx, (y + CY[q]) % ny
+            if (tx, ty) not in A:
+                out["Q_exact"] += O[q][y][x]            # outgoing owner population
+    out["H"] = out["H_seam"] + out["H_edge"]
+    out["residual"] = out["dM"] + out["Q_exact"] - out["C"] - out["H"]
+    return out
+
+
+def staged_flux(O, P, A, nx, ny):
+    """Outward flux with incoming terms from the staged inputs P (the wrong combination)."""
+    A = set(A)
+    q_out = 0.0
+    for (x, y) in A:
+        for q in range(9):
+            sx, sy = _src(x, y, q, nx, ny)
+            if (sx, sy) not in A:
+                q_out -= P[q][y][x]
+            tx, ty = (x + CX[q]) % nx, (y + CY[q]) % ny
+            if (tx, ty) not in A:
+                q_out += O[q][y][x]
+    return q_out