Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

test: prototype conservative space-time rectangle coupling

Browse Source
This commit is contained in:
Codex
2026-10-18 08:55:41 -07:00
parent fc0e544c0a
commit b983c1d095
2 changed files with 446 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/validation/multigrid/spacetime_reference.py
+325
View File
@@ -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())