ClaudeandClaude Opus 5.5
106aa620b9
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>
102 lines
3.9 KiB
Python
102 lines
3.9 KiB
Python
"""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()
|