diff --git a/docs/validation/taylor_green.md b/docs/validation/taylor_green.md index fd2a981..8793328 100644 --- a/docs/validation/taylor_green.md +++ b/docs/validation/taylor_green.md @@ -418,3 +418,39 @@ Labels: - Neither contrast is additive by construction. - The thin-band/acoustic explanation for the nested L1 pressure is a hypothesis, not established. + +### Signed interface-probe traces (opt-in `--probe-trace`, mechanism diagnostic) + +`shear_probes.py` (pinned) defines 24 fixed points on two interface-normal +transects, all at least 8 root cells from any box corner: +- V: y = 64.1, crossing the vertical faces x = 16 (L0/L1) and x = 22 (L1/L2); +- H: x = 90.1, crossing the horizontal faces y = 24 and y = 30; +- offsets from the outer face: −3.9, −1.9, −0.4, 0.1, 0.6, 2.1, 4.1, 5.6, 6.1, + 6.6, 8.1, 10.1. + +Allowed cases: +- nested U128, nested U512, nested projected depth 3; +- twin projected depth 2. + +The option can't be combined with `--final-snapshot`. + +Sampling: +- Each point is read from the ACTIVE leaf cell containing it. The receipt + stores the leaf level, `i`/`j`, and the actual cell center, all in + root-physical units (N=128). +- Grids differ between cases, so centers are *not* collocated. Analytic + fields are evaluated at each actual center (`shear_probes.residuals`), and + comparisons across cases are between residuals at nearby but distinct + centers. +- Times are root-physical steps 0–128 and 1920–2048. U512 is sampled every 4 + lattice steps to match. +- Each sample stores `rho`, `ux`, `uy` from the FP64 physical populations, + POST-collision G at root synchronization only. No substep content is + observed, and grid-scale alternation can alias. + +What's kept and enforced: +- The five qualification samples stay unchanged. +- No field dumps. The trace is validated for exact schedule, finiteness and + a 600 kB cap. If invalid it is dropped and an error is recorded, keeping + the scalar receipt. +- There are no new pass gates; this is a mechanism diagnostic only. diff --git a/tests/validation/taylor_green/run_case.py b/tests/validation/taylor_green/run_case.py index 91d8b9a..60b8ad9 100644 --- a/tests/validation/taylor_green/run_case.py +++ b/tests/validation/taylor_green/run_case.py @@ -25,6 +25,7 @@ NESTED_SHEAR = "nested_shear_n128" # fixed three-level/two-region shear preset SHEAR_NESTED_SHA256 = "fa6cfe87fb4df8c6d732c637c885276833c81973b950c7724d734f8fc350bcb0" TWIN_SHEAR = "twin_shear_n128" # depth-2 twin: nested L0->L1 boxes without L2 SHEAR_TWIN_SHA256 = "1cbf37f18e89aea577581631c4565f51e68d22a3739b698816ef64acbf2b2bb4" +SHEAR_PROBES_SHA256 = "98ea8090ace48f04b83dfb015db10e681e414e57eb748b34c0fbb90fb790d13b" # signed interface-probe mechanism diagnostic def flow_functions(flow="taylor_green"): @@ -94,7 +95,7 @@ def provenance(repo): ) -def verify_source(repo, revision, flow="taylor_green"): +def verify_source(repo, revision, flow="taylor_green", probes=False): """Exact clean harness commit, fixed production tree, and reviewed capture helper.""" import oracle from tests.validation.multigrid import check_projected_pairs as gate @@ -113,6 +114,13 @@ def verify_source(repo, revision, flow="taylor_green"): if result["source_tree_sha256"] != SOURCE_SHA256: raise RuntimeError("Wrong numerical source fingerprint") result["build_capture_sha256"] = BUILD_CAPTURE_SHA256 + if probes: + import shear_probes + module = harness / "shear_probes.py" + if (Path(shear_probes.__file__).resolve() != module + or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_PROBES_SHA256): + raise RuntimeError("Wrong shear_probes module checkout or content") + result["shear_probes_sha256"] = SHEAR_PROBES_SHA256 if flow == "shear_twin": import shear_twin module = harness / "shear_twin.py" @@ -206,6 +214,9 @@ def main(): "(scaled_n256: Taylor-Green; scaled_shear_n256: shear wave) or the " "fixed three-level/two-region shear preset (nested_shear_n128) " "or its depth-2 twin (twin_shear_n128)") + p.add_argument("--probe-trace", action="store_true", + help="signed interface-probe trace (nested U128/U512/projected d3, " + "twin projected d2); mechanism diagnostic, no gates") p.add_argument("--final-snapshot", action="store_true", help="nested_shear_n128 projected depth-3 only: final macro snapshot") p.add_argument("--early-pressure", action="store_true", @@ -216,6 +227,13 @@ def main(): if args.flow == "shear_wave" and (args.re != 100 or (args.depth == 3 and not nested) or args.early_pressure): p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure") + probe_cases = {(NESTED_SHEAR, 128, 1, "constant"), (NESTED_SHEAR, 512, 1, "constant"), + (NESTED_SHEAR, 128, 3, "projected_linear"), + (TWIN_SHEAR, 128, 2, "projected_linear")} + if args.probe_trace and (args.final_snapshot or ( + args.diagnostic, args.n, args.depth, args.explosion) not in probe_cases): + p.error("--probe-trace requires nested U128/U512/projected depth 3 or twin projected " + "depth 2, without --final-snapshot") if args.final_snapshot and not (nested and args.depth == 3 and args.explosion == "projected_linear"): p.error("--final-snapshot requires nested_shear_n128 projected depth 3") @@ -253,7 +271,8 @@ def main(): repo = args.repo.resolve() sys.path[:0] = [str(repo / "src"), str(repo)] source = verify_source(repo, args.expected_revision, - "shear_nested" if nested else "shear_twin" if twin else args.flow) + "shear_nested" if nested else "shear_twin" if twin else args.flow, + probes=args.probe_trace) runtime = runtime_source(repo) if args.dry_run: print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True, @@ -351,8 +370,23 @@ def main(): ] ) initial_mass = None - for step in sorted(set(main_steps + early_steps)): + probe_steps = set() + if args.probe_trace: + import shear_probes + probe_root = dict(zip(shear_probes.case_steps(args.n), shear_probes.STEPS)) + probe_steps = set(probe_root) + probes = shear_probes.select(masks, args.n) + report["probe_trace"] = shear_probes.new_trace( + probes, cfg.viscosity, [args.diagnostic, args.n, args.depth, args.explosion]) + for step in sorted(set(main_steps + early_steps) | probe_steps): run.step(step - run.step_count) + if step in probe_steps: + trace = report["probe_trace"] + trace["steps"].append(probe_root[step]) # root-physical step + trace["values"].append(shear_probes.sample( + (run.download(),) if args.depth == 1 else run.download(), probes)) + if step not in main_steps and step not in early_steps: + continue # probe-only step: no qualification metrics row = sample(run, masks, step, args.n, cfg.viscosity, initial_mass, metrics, args.depth) row["step"] = step @@ -405,6 +439,11 @@ def main(): ), flush=True, ) + if args.probe_trace: + try: + report["probe_trace"]["bytes"] = shear_probes.validate(report["probe_trace"]) + except ValueError as error: # keep the scalar receipt; drop an invalid trace + report["probe_trace"] = dict(error=str(error)) report.update( completed=True, seconds=time.monotonic() - started, diff --git a/tests/validation/taylor_green/shear_probes.py b/tests/validation/taylor_green/shear_probes.py new file mode 100644 index 0000000..47521d0 --- /dev/null +++ b/tests/validation/taylor_green/shear_probes.py @@ -0,0 +1,113 @@ +"""Signed interface-probe traces for the nested/twin shear presets (mechanism diagnostic). + +No pass gates. 24 fixed physical points on two interface-normal transects, away +from corners, crossing L0->L1 and (nested only) L1->L2: + V: y = 64.1, across the vertical faces x = 16 (L0/L1) and x = 22 (L1/L2) of box 1; + H: x = 90.1, across the horizontal faces y = 24 (L0/L1) and y = 30 (L1/L2) of box 2. +Each point is sampled at the ACTIVE leaf cell containing it, so cases with +different grids (U128, U512, nested, twin) give different cell centers: the +actual center and level are stored and analytic fields are evaluated there. +Samples are POST-collision G at root synchronization only (root-step units); +no substep content is observed, and grid-scale alternation can alias. +""" + +import math + +import numpy as np + +import shear_wave as sw +from oracle import CX, CY, U + +ROOT_N = 128 +OFFSETS = (-3.9, -1.9, -0.4, 0.1, 0.6, 2.1, 4.1, 5.6, 6.1, 6.6, 8.1, 10.1) # from outer face +TRANSECTS = {"V": dict(axis=0, face=16.0, other=64.1), # x varies, y fixed + "H": dict(axis=1, face=24.0, other=90.1)} # y varies, x fixed +STEPS = tuple(range(0, 129)) + tuple(range(1920, 2049)) # root steps, both inclusive +TRACE_MAX_BYTES = 600_000 + + +def points(): + out = [] + for name, t in TRANSECTS.items(): + for offset in OFFSETS: + along = t["face"] + offset + xy = (along, t["other"]) if t["axis"] == 0 else (t["other"], along) + out.append((name, offset, xy)) + return out + + +def case_steps(n): + """Case lattice steps at the fixed root-physical sample times (U512: 4 per root step).""" + return [s * (n // ROOT_N) for s in STEPS] + + +def select(masks, n=ROOT_N): + """Leaf cell containing each point: finest level whose ACTIVE mask covers it. + + Points and stored centers are in root-physical units (root cell = 1); a case + with root N = n has cells of size ROOT_N / (n * 2**level) there. + """ + if n % ROOT_N: + raise ValueError("case root N must be a multiple of 128") + probes = [] + for name, offset, (x, y) in points(): + for level in reversed(range(len(masks))): + scale = (n // ROOT_N) * 2 ** level + i, j = int(math.floor(x * scale)), int(math.floor(y * scale)) + if masks[level][j, i]: + h = 1.0 / scale + probes.append(dict(transect=name, offset=offset, point=[x, y], level=level, + i=i, j=j, center=[h * (i + 0.5), h * (j + 0.5)])) + break + else: + raise ValueError(f"point {(x, y)} is not covered by an ACTIVE leaf") + return probes + + +def sample(arrays, probes): + """Flat [rho, ux, uy] per probe from FP64 physical populations.""" + row = [] + for p in probes: + f = np.asarray(arrays[p["level"]][:, p["j"], p["i"]], dtype=float) + rho = float(f.sum()) + row += [rho, float(CX @ f) / rho, float(CY @ f) / rho] + return row + + +def new_trace(probes, case_viscosity, case): + """Times/coordinates are root-physical (N=128); nu_root = U*128/Re for all cases.""" + return dict(schema=1, kind="mechanism diagnostic, no pass gates", case=case, + time_unit="root step (N=128)", sample_phase="postcollision G at root synchronization", + root_n=ROOT_N, viscosity_root=U * ROOT_N / sw.RE, case_viscosity=case_viscosity, + mode=list(sw.MODE), phase=sw.PHI, + amplitude=U, fields=["rho", "ux", "uy"], probes=probes, steps=[], values=[]) + + +def residuals(trace): + """Signed (drho, dux, duy) per step and probe, analytic field at each actual center.""" + out = [] + for step, row in zip(trace["steps"], trace["values"]): + res = [] + for n, p in enumerate(trace["probes"]): + x, y = p["center"] + ux, uy, _ = sw.exact(x, y, step, trace["root_n"], trace["viscosity_root"]) + rho, vx, vy = row[3 * n: 3 * n + 3] + res.append((rho - 1.0, vx - float(ux), vy - float(uy))) + out.append(res) + return out + + +def validate(trace): + """Exact schedule, finite values, fixed probe count and byte cap.""" + import json + if trace.get("steps") != list(STEPS) or len(trace["values"]) != len(STEPS): + raise ValueError("trace steps differ from the fixed schedule") + if len(trace["probes"]) != len(OFFSETS) * len(TRANSECTS): + raise ValueError("unexpected probe count") + width = 3 * len(trace["probes"]) + if any(len(r) != width or not all(math.isfinite(v) for v in r) for r in trace["values"]): + raise ValueError("malformed or nonfinite trace row") + size = len(json.dumps(trace, allow_nan=False).encode()) + if size > TRACE_MAX_BYTES: + raise ValueError(f"trace exceeds {TRACE_MAX_BYTES} bytes ({size})") + return size diff --git a/tests/validation/taylor_green/test_shear_wave.py b/tests/validation/taylor_green/test_shear_wave.py index ba6f7de..9579791 100644 --- a/tests/validation/taylor_green/test_shear_wave.py +++ b/tests/validation/taylor_green/test_shear_wave.py @@ -930,3 +930,189 @@ class TwinPresetTests(unittest.TestCase): self.assertTrue(all(len(r["spurious_pressure_scaled_by_level"]) == 2 for r in rec["rows"])) self.assertEqual(selectors, [("projected_linear", [list(self.twin.REFINEMENTS[0])])]) self.assertEqual(sorted(p.name for p in Path(tmp).iterdir()), ["projected_twin.json"]) + + +class ProbeTraceTests(unittest.TestCase): + """Signed interface-probe mechanism diagnostic (--probe-trace); no pass gates.""" + + def setUp(self): + import shear_nested + import shear_probes + import shear_twin + self.nested, self.probes, self.twin = shear_nested, shear_probes, shear_twin + + def exact_arrays(self, masks, n, t): + """Equilibrium of the exact wave at root-physical time t on a case grid of root N=n.""" + out = [] + for level, m in enumerate(masks): + x, y = centers(m.shape, level) + x, y = x * 128 / n, y * 128 / n + ux, uy, _ = sw.exact(x, y, t, 128, U * 128 / sw.RE) + out.append(sw.equilibrium(np.ones_like(ux), ux, uy)) + return out + + def test_frozen_history_identities(self): + for name, digest in (("shear_wave.py", "81d5c8ee60864372887cddf9239580c82fb6750d956038f1e3ee9de61aa4e877"), + ("shear_nested.py", "fa6cfe87fb4df8c6d732c637c885276833c81973b950c7724d734f8fc350bcb0"), + ("shear_twin.py", "1cbf37f18e89aea577581631c4565f51e68d22a3739b698816ef64acbf2b2bb4"), + ("oracle.py", ORACLE_SHA256_EEC4677)): + self.assertEqual(hashlib.sha256((HERE / name).read_bytes()).hexdigest(), digest, name) + self.assertEqual(run_case.SHEAR_PROBES_SHA256, + hashlib.sha256((HERE / "shear_probes.py").read_bytes()).hexdigest()) + + def test_probe_selection_per_case(self): + sp = self.probes + self.assertEqual(len(sp.points()), 24) + corners = [(x, y) for level in (self.nested.REFINEMENTS[0], + [tuple(v / 2 for v in r) for r in self.nested.REFINEMENTS[1]]) + for x0, y0, x1, y1 in level for x in (x0, x1) for y in (y0, y1)] + for _, _, (x, y) in sp.points(): + self.assertGreaterEqual(min(math.hypot(x - cx, y - cy) for cx, cy in corners), 8.0) + pattern = [0, 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2] + nested = sp.select(self.nested.leaf_masks()) + self.assertEqual([p["level"] for p in nested], pattern * 2) + twin = sp.select(self.twin.leaf_masks()) + self.assertEqual([p["level"] for p in twin], [min(v, 1) for v in pattern] * 2) + u128 = sp.select([np.ones((128, 128), bool)], 128) + u512 = sp.select([np.ones((512, 512), bool)], 512) + self.assertTrue(all(p["level"] == 0 for p in u128 + u512)) + for case in (nested, twin, u128, u512): + for p in case: + # The stored center is the actual cell center containing the point (root units). + h = 1.0 / (2 ** p["level"]) if case is not u512 else 0.25 + self.assertLessEqual(abs(p["center"][0] - p["point"][0]), h / 2) + self.assertLessEqual(abs(p["center"][1] - p["point"][1]), h / 2) + # Different grids give different centers: no pretended collocation. + self.assertNotEqual(u128[3]["center"], u512[3]["center"]) + self.assertEqual(sp.case_steps(512)[:3], [0, 4, 8]) + self.assertEqual(sp.case_steps(128), list(sp.STEPS)) + + def test_residuals_vanish_for_the_exact_field(self): + sp = self.probes + for masks, n in ((self.nested.leaf_masks(), 128), ([np.ones((512, 512), bool)], 512)): + probes = sp.select(masks, n) + trace = sp.new_trace(probes, U * n / sw.RE, ["x", n, len(masks), "c"]) + for t in (0, 100, 2048): + trace["steps"].append(t) + trace["values"].append(sp.sample(self.exact_arrays(masks, n, t), probes)) + for res in sp.residuals(trace): + self.assertLess(max(abs(v) for r in res for v in r), 1e-14) + + def test_validate_schedule_values_and_cap(self): + sp = self.probes + probes = sp.select(self.nested.leaf_masks()) + trace = sp.new_trace(probes, 0.0384, ["n", 128, 3, "projected_linear"]) + trace["steps"] = list(sp.STEPS) + trace["values"] = [[1.0, 0.01, -0.02] * 24 for _ in sp.STEPS] + size = sp.validate(trace) + self.assertLess(size, sp.TRACE_MAX_BYTES) + realistic = dict(trace, values=[[1.0000012345678901, 0.012345678901234567, + -0.023456789012345678] * 24 for _ in sp.STEPS]) + self.assertLess(sp.validate(realistic), sp.TRACE_MAX_BYTES) + for bad in (dict(trace, steps=trace["steps"][:-1]), + dict(trace, values=trace["values"][:-1] + [[float("nan")] * 72]), + dict(trace, probes=trace["probes"][:-1])): + with self.assertRaises(ValueError): + sp.validate(bad) + with patch.object(sp, "TRACE_MAX_BYTES", 1000): + with self.assertRaises(ValueError): + sp.validate(trace) + + def test_probe_runner_options(self): + repo = Path(__file__).resolve().parents[3] + with tempfile.TemporaryDirectory() as tmp: + base = ["run_case.py", "--repo", str(repo), "--expected-revision", "f" * 40, + "--output", str(Path(tmp) / "case.json"), "--dry-run", "--re", "100", + "--flow", "shear_wave", "--probe-trace"] + nest = ["--diagnostic", run_case.NESTED_SHEAR] + tw = ["--diagnostic", run_case.TWIN_SHEAR] + proj = ["--explosion", "projected_linear"] + for extra, valid in ( + (["--n", "128", "--depth", "1"] + nest, True), + (["--n", "512", "--depth", "1"] + nest, True), + (["--n", "128", "--depth", "3"] + proj + nest, True), + (["--n", "128", "--depth", "2"] + proj + tw, True), + (["--n", "128", "--depth", "3"] + nest, False), # constant d3 + (["--n", "128", "--depth", "1"] + tw, False), + (["--n", "128", "--depth", "3"] + proj + nest + ["--final-snapshot"], False), + (["--n", "128", "--depth", "2"], False)): + argv = base + extra + with self.subTest(argv=argv), patch.object(sys, "argv", argv), \ + patch.object(run_case, "verify_source", return_value={}) as verify, \ + patch.object(run_case, "runtime_source", return_value={}), \ + redirect_stdout(io.StringIO()), patch("sys.stderr", io.StringIO()): + if valid: + self.assertEqual(run_case.main(), 0) + self.assertIs(verify.call_args.kwargs["probes"], True) + else: + with self.assertRaises(SystemExit): + run_case.main() + + def run_mocked(self, n, depth, diagnostic, scheme, masks): + from types import SimpleNamespace + from contextlib import nullcontext + arrays = self.exact_arrays(masks, n, 0.0) + counts = [] + + class FakeRun: + def __init__(self, *args, **kwargs): + self.step_count = 0 + + def __enter__(self): + return self + + def __exit__(self, *args): + pass + + def initialize(self, *given): + pass + + def step(self, count): + self.step_count += count + counts.append(self.step_count) + + def download(self): + return arrays[0] if depth == 1 else arrays + + modules = {"CelerisLab.cuda.context": SimpleNamespace(CudaContext=FakeRun), + "CelerisLab.lbm.block_stepper": SimpleNamespace(PeriodicBlockStepper=FakeRun), + "CelerisLab.lbm.refinement_stepper": SimpleNamespace(FixedPatchStepper=FakeRun)} + repo = Path(__file__).resolve().parents[3] + with tempfile.TemporaryDirectory() as tmp: + output = Path(tmp) / "probe.json" + argv = ["run_case.py", "--repo", str(repo), "--re", "100", "--n", str(n), + "--depth", str(depth), "--explosion", scheme, "--flow", "shear_wave", + "--diagnostic", diagnostic, "--probe-trace", "--output", str(output), + "--expected-revision", "f" * 40, "--run-gpu"] + with patch.dict(sys.modules, modules), patch.object(sys, "argv", argv), \ + patch.dict("os.environ", CUDA_VISIBLE_DEVICES="mock"), \ + patch.object(run_case, "verify_source", return_value=dict(revision="f" * 40)), \ + patch.object(run_case, "runtime_source", return_value={}), \ + patch.object(run_case, "runtime_identity", return_value={}), \ + patch.object(run_case, "build_capture", return_value=nullcontext()), \ + redirect_stdout(io.StringIO()): + run_case.main() + return json.loads(output.read_text()), counts, sorted(p.name for p in Path(tmp).iterdir()) + + def test_mocked_probe_traces(self): + sp = self.probes + rec, counts, files = self.run_mocked(128, 3, run_case.NESTED_SHEAR, "projected_linear", + self.nested.leaf_masks()) + trace = rec["probe_trace"] + self.assertEqual(trace["steps"], list(sp.STEPS)) + self.assertEqual(len(trace["values"]), len(sp.STEPS)) + self.assertLessEqual(trace["bytes"], sp.TRACE_MAX_BYTES) + self.assertEqual(trace["case"], [run_case.NESTED_SHEAR, 128, 3, "projected_linear"]) + self.assertEqual(len(rec["rows"]), 5) # qualification samples preserved + self.assertEqual(files, ["probe.json"]) # no field dumps + self.assertEqual(sp.validate(trace) <= sp.TRACE_MAX_BYTES, True) + rec512, counts512, _ = self.run_mocked(512, 1, run_case.NESTED_SHEAR, "constant", + [np.ones((512, 512), bool)]) + t512 = rec512["probe_trace"] + self.assertEqual(t512["steps"], list(sp.STEPS)) # root-physical times + self.assertEqual(t512["viscosity_root"], U * 128 / sw.RE) + self.assertEqual(t512["case_viscosity"], U * 512 / sw.RE) + self.assertIn(4 * 1920, counts512) # U512 advanced to 4x lattice steps + self.assertTrue(all(p["level"] == 0 for p in t512["probes"])) + res = sp.residuals(t512) + self.assertLess(abs(res[0][0][1]), 1e-14) # exact field at t=0 at actual centers