From b983c1d095bc1444aaab2acfcf7c34c3484307e4 Mon Sep 17 00:00:00 2001 From: Codex Date: Sun, 18 Oct 2026 08:27:52 -0700 Subject: [PATCH] test: prototype conservative space-time rectangle coupling --- .../multigrid/spacetime_reference.py | 325 ++++++++++++++++++ .../multigrid/test_spacetime_reference.py | 121 +++++++ 2 files changed, 446 insertions(+) create mode 100644 tests/validation/multigrid/spacetime_reference.py create mode 100644 tests/validation/multigrid/test_spacetime_reference.py diff --git a/tests/validation/multigrid/spacetime_reference.py b/tests/validation/multigrid/spacetime_reference.py new file mode 100644 index 0000000..879aba6 --- /dev/null +++ b/tests/validation/multigrid/spacetime_reference.py @@ -0,0 +1,325 @@ +"""Experimental FP64 BVP coupling: one periodic rectangle, D2Q9, ratio two. + +Guzik et al., https://www.osti.gov/servlets/purl/1226215, sections 3.2--3.5. +Coordinate-generated stream events; CPUREF collision only, no production imports. +This is a mathematical prototype, not a certified accuracy implementation. +""" +from dataclasses import dataclass +from functools import lru_cache +import argparse +import hashlib +import json +from pathlib import Path +import resource +import sys +import subprocess +import time + +import numpy as np +import cpu_reference as cpu + + +class Spec(cpu.Spec): + @property + def n_fine(self): + return 2 * self.n_coarse + + +def geometry(n, bounds): + x0, y0, x1, y1 = bounds + if not (2 <= x0 < x1 <= n - 2 and 2 <= y0 < y1 <= n - 2): + raise ValueError("one interior rectangle with two coarse cells of clearance required") + y, x = np.indices((n, n)) + covered = (x >= x0) & (x < x1) & (y >= y0) & (y < y1) + y, x = np.indices((2*n, 2*n)) + active = (x >= 2*x0) & (x < 2*x1) & (y >= 2*y0) & (y < 2*y1) + inner = (x >= 2*x0-2) & (x < 2*x1+2) & (y >= 2*y0-2) & (y < 2*y1+2) & ~active + rows, links = [], [] + for ky, kx in zip(*np.nonzero(~covered)): + for q, (cy, cx) in enumerate(zip(cpu.CY, cpu.CX)): + if any(active[2*ky+dy-cy, 2*kx+dx-cx] + for dy in (0, 1) for dx in (0, 1) + if 0 <= 2*ky+dy-cy < 2*n and 0 <= 2*kx+dx-cx < 2*n): + rows.append((q, ky, kx)) + sy, sx = (ky-cy) % n, (kx-cx) % n + if not covered[sy, sx]: + links.append(cpu.CornerLink(f"{kx}_{ky}_{q}", q, (sx, sy), (kx, ky))) + kind = np.where(active, 0, np.where(inner, 1, 2)) + return cpu.Maps(~covered, covered, active, inner, kind == 2, kind, + np.arange(2*n)//2, (np.arange(2*n)+.5)/2, + np.array(rows, dtype=int), tuple(links)) + + +@lru_cache(None) +def interpolation(kind, points): + """Return data coordinates, evaluation matrix, and constrained coefficient map.""" + if kind == "face": + samples = np.array([(s, t) for t in (0., 1.) for s in (-1., 0., 1.)]) + a = np.array([[1, s, t, s*t] for s, t in samples]) + e = np.array([[1, s, t, s*t] for s, t in points]) + b = e.mean(axis=0, keepdims=True) + d = np.zeros((1, 6)); d[0, 1] = 1 + else: + cross = ((0, 0), (-1, 0), (1, 0), (0, -1), (0, 1)) + samples = np.array([(x, y, t) for t in (0., 1.) for x, y in cross]) + a = np.array([[1, x, y, t] for x, y, t in samples]) + e = np.array([[1, x, y, t] for x, y, t in points]) + b = np.stack((e.mean(axis=0), [0, 0, 0, 1])) + d = np.zeros((2, 10)); d[0, 0] = 1; d[1, 0] = -1; d[1, 5] = 1 + kkt = np.block([[a.T @ a, b.T], [b, np.zeros((len(b), len(b)))]]) + coeff = np.linalg.solve(kkt, np.vstack((a.T, d)))[:4] + if not np.allclose(b @ coeff, d, rtol=0, atol=2e-14): + raise AssertionError("invalid constrained interpolation matrix") + return samples, a, e, coeff + + +@dataclass +class Group: + q: int + source: int + receiver: int + kind: str + points: tuple + samples: np.ndarray + weights: np.ndarray + events: tuple # (substep, fine ghost slot), or None for a leftover + + +class Plan: + def __init__(self, n, bounds): + self.n, self.nf = n, 2*n + self.maps = geometry(n, bounds) + active = self.maps.fine_active + lo, hi = np.array(bounds[1::-1]), np.array(bounds[3:1:-1]) + self.groups = [] + for iy, ix in zip(*np.nonzero(self.maps.coarse_active)): + origin = np.array((iy, ix)) + for q in range(1, 9): + c = np.array((cpu.CY[q], cpu.CX[q])) + starts = [2*origin + (dy, dx) for dy in (0, 1) for dx in (0, 1)] + events = [] + for p in starts: + event = None + for m in (0, 1): + old, new = p+m*c, p+(m+1)*c + inside = lambda z: bool(np.all(z >= 0) and np.all(z < 2*n) and active[tuple(z)]) + if not inside(old) and inside(new): + event = (m, int(old[0]*2*n+old[1])) + break + events.append(event) + if all(event is None for event in events): + continue + axes = [j for j in range(2) if origin[j] < lo[j] or origin[j] >= hi[j]] + kind = "corner" if len(axes) == 2 else "face" + points = [] + for p in starts: + x = (p+.5)/2 + times = [(float(lo[j] if c[j] > 0 else hi[j])-x[j])/c[j] for j in axes] + time = max(times) + offset = x+c*time-(origin+.5+c*.5) + points.append((float(offset[1]), float(offset[0]), float(time-.5)) + if kind == "corner" else (float(offset[1-axes[0]]), float(time-.5))) + points = tuple(points) + samples, _, e, coeff = interpolation(kind, points) + indices = [] + for row in samples: + shift = np.array((row[1], row[0]), dtype=int) if kind == "corner" else np.eye(2, dtype=int)[1-axes[0]]*int(row[0]) + pos = (origin+shift) % n + indices.append(int(pos[0]*n+pos[1])) + target = (origin+c) % n + self.groups.append(Group(q, int(iy*n+ix), int(target[0]*n+target[1]), kind, + points, np.array(indices), e @ coeff, tuple(events))) + self.outgoing = [] + # m is zero-based: a substep-0 exit travels two fine spacings by + # synchronization; a substep-1 exit travels one. + for m in (0, 1): + records = [] + for q in range(1, 9): + c = np.array((cpu.CY[q], cpu.CX[q])) + for y, x in zip(*np.nonzero(active)): + p = np.array((y, x)); endpoint = p+c + if not active[tuple(endpoint)]: + target = (p+(2-m)*c)//2 + assert self.maps.coarse_active[tuple(target)] + records.append((q, int(y*2*n+x), int(target[0]*n+target[1]))) + self.outgoing.append(np.array(records, dtype=int)) + incoming = [(m, g.q, slot) for g in self.groups for event in g.events + if event is not None for m, slot in [event]] + assert len(incoming) == len(set(incoming)), "duplicate incoming writer" + expected = {(m, int(q), int((y-cpu.CY[q])*2*n+x-cpu.CX[q])) + for m in (0, 1) for q in range(1, 9) + for y, x in zip(*np.nonzero(active)) if not active[y-cpu.CY[q], x-cpu.CX[q]]} + assert set(incoming) == expected, "incoming coverage mismatch" + receivers = {(int(q), int(target)) for out in self.outgoing for q, _, target in out} + receivers |= {(g.q, g.receiver) for g in self.groups if None in g.events} + assert receivers == {(int(q), int(y*n+x)) for q, y, x in self.maps.coalesce} + self.receivers = np.array(sorted(receivers), dtype=int) + + +def directional(state, maps): + return state.coarse_f[:, maps.coarse_active].sum(axis=1) + .25*state.fine_f[:, maps.fine_active].sum(axis=1) + + +def advance(spec, plan, state, mode="spacetime", mutation=None): + """One arrival-to-arrival interval. 'frozen' disables space AND time variation.""" + if mode not in ("frozen", "spacetime"): + raise ValueError(mode) + n, nf, maps = plan.n, plan.nf, plan.maps + if state.phase != "arrival_f" or state.tick % 2: + raise ValueError("synchronized arrival state required") + fc, ff = state.coarse_f.copy(), state.fine_f.copy() + covered_mean = ff.reshape(9, n, 2, n, 2).mean(axis=(2, 4)) + fc[:, maps.coarse_covered] = covered_mean[:, maps.coarse_covered] + all_cells = np.ones((n, n), dtype=bool) + def collide(f, mask, tau): + return cpu._collide(spec, f, np.isfinite(f), tau, mask, np.zeros(mask.shape, dtype=int)) + gc = collide(fc, all_cells, spec.tau_coarse) + collision = (gc[:, maps.coarse_active]-fc[:, maps.coarse_active]).sum(axis=1) + predicted = cpu._stream_periodic(gc) + gp = collide(predicted, all_cells, spec.tau_coarse) if mode == "spacetime" else gc + registers = np.zeros((9, n*n)) + q, target = plan.receivers.T + registers[q, target] = -predicted.reshape(9, -1)[q, target] + incoming = [[], []] + group_error = 0. + min_boundary = float("inf") + negative_boundary = nonfinite_boundary = 0 + dropped = False + for group in plan.groups: + q = group.q + count = len(group.samples)//2 + data = np.concatenate((gc.reshape(9, -1)[q, group.samples[:count]], gp.reshape(9, -1)[q, group.samples[count:]])) + value = (group.weights @ data if mode == "spacetime" + else np.full(4, gc.reshape(9, -1)[q, group.source])) + min_boundary = min(min_boundary, float(value.min())) + negative_boundary += int(np.count_nonzero(value < 0)) + nonfinite_boundary += int(np.count_nonzero(~np.isfinite(value))) + group_error = max(group_error, abs(value.mean()-gc.reshape(9, -1)[q, group.source])) + for k, event in enumerate(group.events): + if event is None: + if mutation == "omit_leftover" and not dropped: + dropped = True + else: + registers[q, group.receiver] += .25*value[k] + else: + m, slot = event + incoming[m].append((q, slot, value[k])) + for m in (0, 1): + gf = collide(ff, maps.fine_active, spec.tau_fine) + collision += .25*(gf[:, maps.fine_active]-ff[:, maps.fine_active]).sum(axis=1) + out = plan.outgoing[m] + if mutation != f"omit_outgoing_{m}": + np.add.at(registers, (out[:, 0], out[:, 2]), .25*gf.reshape(9, -1)[out[:, 0], out[:, 1]]) + for q, slot, value in incoming[m]: + gf.reshape(9, -1)[q, slot] = value + fresh = np.full_like(ff, np.nan) + for q in range(9): + moved = cpu._shift(gf[q], int(cpu.CY[q]), int(cpu.CX[q]), np.nan) + fresh[q, maps.fine_active] = moved[maps.fine_active] + if not np.isfinite(fresh[:, maps.fine_active]).all(): + raise AssertionError("invalid fine incoming population") + ff = fresh + final_c = predicted + registers.reshape(9, n, n) + final_c[:, maps.coarse_covered] = np.nan + final = cpu.State(final_c, ff, tick=state.tick+2) + residual = directional(final, maps)-directional(state, maps)-collision + minimum = min(float(final_c[:, maps.coarse_active].min()), float(ff[:, maps.fine_active].min())) + return final, dict(directional_residual=residual, collision=collision, + group_error=group_error, min_boundary=min_boundary, min_population=minimum, + negative_boundary_count=negative_boundary, nonfinite_boundary_count=nonfinite_boundary) + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--n64", action="store_true", help="explicitly run the frozen 1024-interval experiment") + parser.add_argument("--output", type=Path, required=True) + parser.add_argument("--max-seconds",type=float,default=180.,help="wall-time cap, including plan generation") + parser.add_argument("--max-rss-mib",type=float,default=256.,help="Linux peak-RSS screening cap") + args = parser.parse_args() + if not args.n64: + parser.error("full trajectory requires --n64") + if args.output.exists(): + parser.error("output exists; preserve prior evidence and choose a new path") + if not (np.isfinite(args.max_seconds) and args.max_seconds>0 and np.isfinite(args.max_rss_mib) and args.max_rss_mib>0): + parser.error("positive finite resource caps required") + started = time.monotonic() + sys.path.insert(0, str(Path(__file__).parents[1]/"taylor_green")) + import oracle as tg + plan = Plan(64, (8, 8, 56, 48)) + spec = Spec(n_coarse=64, patch_lo=8, halo=16, tau_coarse=.5576, tau_fine=.6152, collision="srt") + arrays = [tg.initial((64,64), 0, 64, .0192), tg.initial((128,128), 1, 64, .0192)] + for a, mask in zip(arrays, (plan.maps.coarse_active, plan.maps.fine_active)): + a[:, ~mask] = np.nan + initial = cpu.State(*arrays) + mass0 = tg.metrics(arrays, (plan.maps.coarse_active, plan.maps.fine_active), 0, 64, .0192)["mass"] + metrics = ("velocity_l2", "pressure_l2", "energy_relative") + frozen_baseline = dict(zip(metrics, (.01131134574353476, .03592961319244149, .012514764591008386))) + budget_bound = float(64*np.finfo(float).eps*mass0) + root = Path(__file__).resolve().parents[3] + source = dict(commit=subprocess.check_output(["git", "rev-parse", "HEAD"],cwd=root,text=True).strip(), + dirty=bool(subprocess.check_output(["git","status","--porcelain"],cwd=root,text=True)), + hashes={p.name:hashlib.sha256(p.read_bytes()).hexdigest() + for p in (Path(__file__),Path(cpu.__file__),Path(tg.__file__))},numpy=np.__version__) + result = dict(scope="experimental CPU FP64; no GPU/accuracy certification", source=source, + screen=dict(velocity_factor=.8, pressure_factor=1., energy_factor=1., + directional_step_bound=budget_bound, mass_relative=1e-11, + momentum_absolute=1e-9, baseline_metric_atol=1e-10), + limits=dict(seconds=args.max_seconds,rss_mib=args.max_rss_mib),runs=[],completed=False) + def save(): + result["resources"] = dict(seconds=time.monotonic()-started,peak_rss_mib=resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024) + args.output.write_text(json.dumps(result,indent=2,allow_nan=False)+"\n") + save() + try: + for mode in ("constant_ivp", "spacetime"): + measured = mode == "spacetime" + state = initial; record = dict(mode=mode,rows=[], + boundary_budget_status="measured" if measured else "not_measured", + max_directional_residual=0. if measured else None, + minimum_boundary=None,negative_boundary_count=0 if measured else None, + nonfinite_boundary_count=0 if measured else None) + result["runs"].append(record) + for k in range(1, 1025): + if time.monotonic()-started>args.max_seconds or resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024>args.max_rss_mib: + raise RuntimeError("experiment resource cap exceeded") + if mode == "constant_ivp": + state,_ = cpu.advance_coarse_once(spec,plan.maps,state) + else: + state,diagnostic = advance(spec,plan,state) + error = float(np.max(np.abs(diagnostic["directional_residual"]))) + record["max_directional_residual"] = max(record["max_directional_residual"],error) + if np.isfinite(diagnostic["min_boundary"]): + record["minimum_boundary"] = (diagnostic["min_boundary"] if record["minimum_boundary"] is None else min(record["minimum_boundary"],diagnostic["min_boundary"])) + for name in ("negative_boundary_count","nonfinite_boundary_count"): + record[name] += diagnostic[name] + if (not np.isfinite(error) or error>budget_bound or diagnostic["group_error"]>2e-14 + or diagnostic["negative_boundary_count"] or diagnostic["nonfinite_boundary_count"] + or not np.isfinite(diagnostic["min_boundary"]) or diagnostic["min_boundary"]<=0 + or not np.isfinite(diagnostic["min_population"]) or diagnostic["min_population"]<=0): + record["failed_diagnostic"] = {name: repr(value) for name,value in diagnostic.items()} + raise AssertionError("candidate directional/health invariant failed") + if k in (256,512,768,1024): + row = tg.metrics((state.coarse_f,state.fine_f),(plan.maps.coarse_active,plan.maps.fine_active),k,64,.0192,mass0) + row["momentum_change"] = (cpu.MOMENTS.T@(directional(state,plan.maps)-directional(initial,plan.maps)))[1:].tolist() + record["rows"].append(row) + save(); print(mode,k,json.dumps(row),flush=True) + record["worst"] = {key:max(row[key] for row in record["rows"]) for key in metrics} + if mode == "constant_ivp": + for key,value in frozen_baseline.items(): + if abs(record["worst"][key]-value)>1e-10: + raise AssertionError(f"frozen baseline mismatch: {key}") + baseline,candidate = result["runs"] + result["checks"] = {key:candidate["worst"][key]<=factor*baseline["worst"][key] + for key,factor in zip(metrics,(.8,1.,1.))} + result["checks"]["mass_momentum"] = all(row["mass_relative"]<=1e-11 and max(map(abs,row["momentum_change"]))<=1e-9 for row in candidate["rows"]) + result["completed"] = True + result["experimental_pass"] = all(result["checks"].values()) + save() + return 0 if result["experimental_pass"] else 1 + except Exception as exc: + result["failure"] = dict(mode=mode,step=k,error=repr(exc)) + save() + raise + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/validation/multigrid/test_spacetime_reference.py b/tests/validation/multigrid/test_spacetime_reference.py new file mode 100644 index 0000000..dd3387b --- /dev/null +++ b/tests/validation/multigrid/test_spacetime_reference.py @@ -0,0 +1,121 @@ +"""Tiny CPU algebra/event tests; no full trajectory or GPU execution.""" +from dataclasses import replace +import unittest +import numpy as np +import spacetime_reference as st + + +class SpaceTime(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.plan = st.Plan(12, (3, 3, 9, 8)) + cls.spec = st.Spec(n_coarse=12, patch_lo=1, halo=2, + tau_coarse=.5576, tau_fine=.6152, collision="srt") + + def state(self, tagged=False): + arrays = [] + for n, mask in ((12,self.plan.maps.coarse_active),(24,self.plan.maps.fine_active)): + q, y, x = np.indices((9,n,n)) + a = np.broadcast_to(st.cpu.W[:,None,None], (9,n,n)).copy() + if tagged: + a *= 1+1e-3*np.random.default_rng(n).normal(size=a.shape) + a[:,~mask] = np.nan + arrays.append(a) + return st.cpu.State(*arrays) + + def test_events_and_corner_leftovers(self): + groups = [g for g in self.plan.groups if None in g.events] + self.assertEqual(len(groups), 8) + self.assertTrue(all(g.events.count(None)==3 for g in groups)) + for m in (0,1): + events = list(map(tuple,self.plan.outgoing[m])) + self.assertEqual(len(events),len(set(events))) + st.Plan(16,(4,5,11,12)) + + def test_interpolation(self): + for g in self.plan.groups: + samples, a, e, coeff = st.interpolation(g.kind,g.points) + np.testing.assert_allclose(g.weights@np.ones(len(samples)),1,rtol=0,atol=2e-14) + np.testing.assert_allclose(g.weights.mean(axis=0), np.eye(len(samples))[1 if g.kind=='face' else 0],rtol=0,atol=2e-14) + if g.kind=='face': + beta=np.array([1.,.03,.02,0.]) + else: + cx,cy=st.cpu.CX[g.q],st.cpu.CY[g.q] + beta=np.array([1.,.03,.02,-cx*.03-cy*.02]) + steady=a@np.array([1.,.03,.02,0.]) + self.assertLess(abs((coeff@steady)[3]),2e-14) + np.testing.assert_allclose(g.weights@(a@beta),e@beta,rtol=0,atol=2e-14) + + def test_rest(self): + initial=self.state() + final,d=st.advance(self.spec,self.plan,initial) + for a,b in zip((initial.coarse_f,initial.fine_f),(final.coarse_f,final.fine_f)): + np.testing.assert_allclose(a,b,rtol=0,atol=2e-14,equal_nan=True) + + def test_frozen_recovers_ivp(self): + for collision in ('identity','srt'): + spec=replace(self.spec,collision=collision) + a=b=self.state(True) + for _ in range(2): + a,_=st.advance(spec,self.plan,a,'frozen') + b,_=st.cpu.advance_coarse_once(spec,self.plan.maps,b) + for aa,bb in zip((a.coarse_f,a.fine_f),(b.coarse_f,b.fine_f)): + np.testing.assert_allclose(aa,bb,rtol=0,atol=2e-14,equal_nan=True) + + def test_second_substep_uses_characteristic_origin(self): + # Review counterexample in (x,y): ghost (7,10), m=1 traces back + # to (6,9), parent (3,4), not its current parent (3,5). + plan=st.Plan(16,(4,4,12,12)) + owners=[g for g in plan.groups if g.q==5 and (1,10*32+7) in g.events] + self.assertEqual(len(owners),1) + self.assertEqual(owners[0].source,4*16+3) + spec=st.Spec(n_coarse=16,patch_lo=2,halo=4,collision='identity') + arrays=[] + for n,mask in ((16,plan.maps.coarse_active),(32,plan.maps.fine_active)): + a=np.broadcast_to(st.cpu.W[:,None,None],(9,n,n)).copy() + a[:,~mask]=np.nan + arrays.append(a) + arrays[0][5,4,3]=1.; arrays[0][5,5,3]=2. + final,_=st.advance(spec,plan,st.cpu.State(*arrays),'frozen') + self.assertEqual(final.fine_f[5,11,8],1.) + self.assertNotEqual(final.fine_f[5,11,8],arrays[0][5,5,3]) + + def test_corner_affine_tradeoff_is_explicit(self): + group=next(g for g in self.plan.groups if g.kind=='corner') + _,a,e,coeff=st.interpolation(group.kind,group.points) + linear_x=np.array([0.,1.,0.,0.]) + # The two paper constraints are consistent. Arbitrary affine exactness + # is the incompatible additional demand. + self.assertEqual(abs(float((e@linear_x).mean())),.125) + fitted=coeff@(a@linear_x) + self.assertLess(abs(fitted[3]),2e-14) + self.assertLess(abs(float((e@fitted).mean())),2e-14) + self.assertGreater(np.max(np.abs(e@fitted-e@linear_x)),.01) + + def test_directional_budgets(self): + for collision in ('identity','srt'): + spec=replace(self.spec,collision=collision) + a=self.state(True) + for _ in range(2): + a,d=st.advance(spec,self.plan,a) + np.testing.assert_allclose(d['directional_residual'],0,rtol=0,atol=2e-12) + self.assertLess(d['group_error'],2e-14) + + def test_missing_packets_fail(self): + for mutation in ('omit_leftover','omit_outgoing_0','omit_outgoing_1'): + _,d=st.advance(self.spec,self.plan,self.state(True),mutation=mutation) + self.assertGreater(np.max(np.abs(d['directional_residual'])),1e-5) + + def test_negative_interpolation_is_reported_not_clipped(self): + initial=self.state() + group=self.plan.groups[0] + initial.coarse_f.reshape(9,-1)[group.q,group.source]=-.01 + spec=replace(self.spec,collision='identity') + _,d=st.advance(spec,self.plan,initial,'frozen') + self.assertGreaterEqual(d['negative_boundary_count'],4) + self.assertEqual(d['min_boundary'],-.01) + self.assertEqual(d['nonfinite_boundary_count'],0) + + +if __name__=='__main__': + unittest.main(verbosity=2)