"""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()