791 lines
30 KiB
Python
791 lines
30 KiB
Python
#!/usr/bin/env python3
|
||
"""Frozen v0.2 shear interface diagnostic (test-side) against unchanged CPUREF-0.3.
|
||
|
||
v0.2 change (root/TRON checker fix only): t0 convention mismatches and init_mismatch
|
||
are now explicit Boolean hard checks (INIT_TOL=CONS_TOL for O(1) rho/j restriction,
|
||
TRANS_TOL=1e-12*U0 for t0 velocity matches); per-cell RMS of e_IF patch/outside
|
||
printed as reported diagnostics. No algorithm, threshold or reference change.
|
||
|
||
v0.1 changes (Claude review B1-B4): include A_T in hard checks; print mass/mom
|
||
drifts with relative tol wording; clarify coarse-step vs CPUREF-tick units and
|
||
unnormalised L2 norms; split e_IF patch/outside + print (uA-uB)(y) profiles;
|
||
report Pi_xy mean; drop made-up ±10% bracket boolean (print s alone).
|
||
|
||
TRON freeze (thread 1791864801.121109) + implementation ask (1791865457.667669):
|
||
7 runs, tau_c=0.8 / tau_f=1.1, 32 coarse steps, U0=1e-3, phi in {0, pi/2},
|
||
A/B/C plus A transpose. Equilibrium-only from exact cell-average sine velocity;
|
||
restrict rho/j then divide; common 16^2 projection; startup 0-3 reported, fit 4-32.
|
||
Hard checks: finite active states, conservation, conventions, transpose.
|
||
Rate brackets / acoustic / interface errors are reported diagnostics only.
|
||
NumPy+stdlib only; reference file unchanged; no project integration.
|
||
|
||
Label: init is equilibrium of *averaged velocity*, NOT the exact cell-average of
|
||
the nonlinear equilibrium population — population equality at t=0 is NOT asserted.
|
||
"""
|
||
from __future__ import annotations
|
||
|
||
import hashlib
|
||
import importlib.util
|
||
import os
|
||
import sys
|
||
import time
|
||
from dataclasses import dataclass
|
||
|
||
import numpy as np
|
||
|
||
HERE = os.path.dirname(os.path.abspath(__file__))
|
||
REF_PATH = os.environ.get(
|
||
"CPU_REFERENCE_PATH", os.path.join(HERE, "cpu_reference_v0.3.py")
|
||
)
|
||
|
||
U0 = 1e-3
|
||
N_COARSE = 16
|
||
N_FINE_DOM = 32 # uniform-fine domain in fine cells over physical [0,16)
|
||
PATCH_LO, PATCH_HI = 4, 12
|
||
TAU_C, TAU_F = 0.8, 1.1
|
||
N_STEPS = 32
|
||
FIT_LO, FIT_HI = 4, 32 # inclusive *coarse steps*; CPUREF ticks = 2*step (fine hops)
|
||
CONS_TOL = 1e-12
|
||
TRANS_TOL = 1e-12 * U0
|
||
INIT_TOL = CONS_TOL # round-off tol for restricted rho,j vs coarse formula (O(1) quantities)
|
||
T0_TOL = TRANS_TOL # round-off tol for t0 velocity matches (O(U0) quantities)
|
||
K_PHYS = 2.0 * np.pi / 16.0 # physical wave number on [0,16)
|
||
CS2 = 1.0 / 3.0
|
||
NU = CS2 * (TAU_C - 0.5) # = 0.1 in physical (coarse) units; matches fine scaling
|
||
|
||
CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.int64)
|
||
CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.int64)
|
||
W = np.array([4 / 9] + [1 / 9] * 4 + [1 / 36] * 4, dtype=np.float64)
|
||
NQ = 9
|
||
|
||
|
||
def _load_ref():
|
||
spec = importlib.util.spec_from_file_location("cpu_ref_v03", REF_PATH)
|
||
mod = importlib.util.module_from_spec(spec)
|
||
sys.modules[spec.name] = mod
|
||
spec.loader.exec_module(mod)
|
||
return mod
|
||
|
||
|
||
def sha256_file(path: str) -> str:
|
||
h = hashlib.sha256()
|
||
with open(path, "rb") as f:
|
||
for chunk in iter(lambda: f.read(1 << 16), b""):
|
||
h.update(chunk)
|
||
return h.hexdigest()
|
||
|
||
|
||
def sinc(x: np.ndarray | float) -> np.ndarray | float:
|
||
x = np.asarray(x, dtype=np.float64)
|
||
out = np.empty_like(x, dtype=np.float64)
|
||
small = np.abs(x) < 1e-14
|
||
out[small] = 1.0
|
||
xs = x[~small]
|
||
out[~small] = np.sin(xs) / xs
|
||
return float(out) if out.ndim == 0 else out
|
||
|
||
|
||
def feq(rho, ux, uy):
|
||
"""Second-order D2Q9 equilibrium; rho,ux,uy broadcastable to (ny,nx)."""
|
||
rho = np.asarray(rho, dtype=np.float64)
|
||
ux = np.asarray(ux, dtype=np.float64)
|
||
uy = np.asarray(uy, dtype=np.float64)
|
||
cu = CX[:, None, None] * ux + CY[:, None, None] * uy
|
||
usq = ux * ux + uy * uy
|
||
return W[:, None, None] * rho * (1.0 + 3.0 * cu + 4.5 * cu * cu - 1.5 * usq)
|
||
|
||
|
||
def srt_collide(f, tau, mask=None):
|
||
g = f.copy()
|
||
if mask is None:
|
||
mask = np.ones(f.shape[1:], dtype=bool)
|
||
fm = f[:, mask]
|
||
rho = fm.sum(axis=0)
|
||
ux = (CX[:, None] * fm).sum(axis=0) / rho
|
||
uy = (CY[:, None] * fm).sum(axis=0) / rho
|
||
eq = feq(rho, ux, uy) # shape (9, ny, nx) but we need only masked — rebuild
|
||
# rebuild eq on full grid then mask
|
||
rho_g = np.zeros(f.shape[1:])
|
||
ux_g = np.zeros(f.shape[1:])
|
||
uy_g = np.zeros(f.shape[1:])
|
||
rho_g[mask] = rho
|
||
ux_g[mask] = ux
|
||
uy_g[mask] = uy
|
||
eq_full = feq(rho_g, ux_g, uy_g)
|
||
g[:, mask] = fm - (fm - eq_full[:, mask]) / tau
|
||
return g
|
||
|
||
|
||
def stream_periodic(g):
|
||
out = np.empty_like(g)
|
||
for i in range(NQ):
|
||
out[i] = np.roll(g[i], (int(CY[i]), int(CX[i])), axis=(0, 1))
|
||
return out
|
||
|
||
|
||
def moments(f):
|
||
rho = f.sum(axis=0)
|
||
jx = (CX[:, None, None] * f).sum(axis=0)
|
||
jy = (CY[:, None, None] * f).sum(axis=0)
|
||
return rho, jx, jy
|
||
|
||
|
||
def cell_avg_ux(y_centers, dx, phi, transpose=False):
|
||
"""Exact cell-average of U0 sin(k y + phi) [or sin(k x + phi) if transpose]."""
|
||
factor = float(sinc(K_PHYS * dx / 2.0))
|
||
if transpose:
|
||
# shear u_y(x): need x centers; caller passes x as y_centers arg
|
||
return U0 * factor * np.sin(K_PHYS * y_centers + phi)
|
||
return U0 * factor * np.sin(K_PHYS * y_centers + phi)
|
||
|
||
|
||
def mode_template(phi):
|
||
"""m on 16^2 coarse centers: sinc(k/2) sin(k yc + phi)."""
|
||
yc = np.arange(N_COARSE, dtype=np.float64) + 0.5
|
||
m1d = float(sinc(K_PHYS / 2.0)) * np.sin(K_PHYS * yc + phi)
|
||
return np.broadcast_to(m1d[:, None], (N_COARSE, N_COARSE)).copy() # [y,x]
|
||
|
||
|
||
def project_amp(ux, phi):
|
||
m = mode_template(phi)
|
||
return float(np.sum(ux * m) / np.sum(m * m))
|
||
|
||
|
||
def restrict_fine_to_coarse(rho_f, jx_f, jy_f):
|
||
"""2x2 arithmetic mean of conserved rho,j on 32^2 -> 16^2, then u=j/rho."""
|
||
assert rho_f.shape == (N_FINE_DOM, N_FINE_DOM)
|
||
rho = rho_f.reshape(N_COARSE, 2, N_COARSE, 2).mean(axis=(1, 3))
|
||
jx = jx_f.reshape(N_COARSE, 2, N_COARSE, 2).mean(axis=(1, 3))
|
||
jy = jy_f.reshape(N_COARSE, 2, N_COARSE, 2).mean(axis=(1, 3))
|
||
ux = jx / rho
|
||
uy = jy / rho
|
||
return rho, jx, jy, ux, uy
|
||
|
||
|
||
def pi_xy(f, tau):
|
||
"""Pre-collision Π_xy neq proxy: sum c_x c_y f (lattice units of that grid)."""
|
||
return (CX[:, None, None] * CY[:, None, None] * f).sum(axis=0)
|
||
|
||
|
||
@dataclass
|
||
class Series:
|
||
name: str
|
||
phi: float
|
||
A: np.ndarray # length 33 (coarse steps 0..32)
|
||
e_if: np.ndarray # ||u - u_B|| unnormalised L2 over 256 cells (filled later for A)
|
||
e_cb: np.ndarray
|
||
e_if_patch: np.ndarray # L2 on patch cells only (filled for A)
|
||
e_if_out: np.ndarray # L2 on outside cells only (filled for A)
|
||
nonmodal: np.ndarray
|
||
max_uy: np.ndarray
|
||
max_drho: np.ndarray
|
||
mass: np.ndarray
|
||
momx: np.ndarray
|
||
momy: np.ndarray
|
||
pi_xy_mean: np.ndarray # mean pre-collision Pi_xy proxy (uniform grids); NaN for composite
|
||
finite_ok: bool
|
||
cons_ok: bool
|
||
mass_drift: float = 0.0
|
||
momx_drift: float = 0.0
|
||
momy_drift: float = 0.0
|
||
notes: str = ""
|
||
|
||
|
||
def init_uniform(n, dx, tau, phi, transpose=False):
|
||
"""Initialize n×n periodic SRT lattice; dx = cell size in physical units."""
|
||
centers = (np.arange(n, dtype=np.float64) + 0.5) * dx
|
||
# physical coords
|
||
Y, X = np.meshgrid(centers, centers, indexing="ij") # [y,x]
|
||
rho = np.ones((n, n), dtype=np.float64)
|
||
if transpose:
|
||
ux = np.zeros_like(rho)
|
||
uy = cell_avg_ux(X, dx, phi) # u_y = U0 sinc sin(k x + phi)
|
||
else:
|
||
ux = cell_avg_ux(Y, dx, phi)
|
||
uy = np.zeros_like(rho)
|
||
f = feq(rho, ux, uy)
|
||
return f
|
||
|
||
|
||
def advance_uniform(f, tau, steps):
|
||
for _ in range(steps):
|
||
g = srt_collide(f, tau)
|
||
f = stream_periodic(g)
|
||
return f
|
||
|
||
|
||
def observe_uniform(f, dx, phi, transpose=False):
|
||
rho, jx, jy = moments(f)
|
||
if dx < 1.0 - 1e-12: # fine 32^2 -> restrict
|
||
rho_c, jx_c, jy_c, ux, uy = restrict_fine_to_coarse(rho, jx, jy)
|
||
else:
|
||
rho_c, jx_c, jy_c = rho, jx, jy
|
||
ux, uy = jx / rho, jy / rho
|
||
if transpose:
|
||
# project onto transposed mode using uy vs m(x)
|
||
m = mode_template(phi) # m[y,x] = sinc sin(k yc+phi) — for transpose need m_x
|
||
xc = np.arange(N_COARSE, dtype=np.float64) + 0.5
|
||
mx = float(sinc(K_PHYS / 2.0)) * np.sin(K_PHYS * xc + phi)
|
||
mT = np.broadcast_to(mx[None, :], (N_COARSE, N_COARSE)).copy()
|
||
amp = float(np.sum(uy * mT) / np.sum(mT * mT))
|
||
modal = amp * mT
|
||
nonmodal = float(np.linalg.norm(uy - modal) / U0)
|
||
max_cross = float(np.max(np.abs(ux)) / U0)
|
||
else:
|
||
amp = project_amp(ux, phi)
|
||
m = mode_template(phi)
|
||
modal = amp * m
|
||
nonmodal = float(np.linalg.norm(ux - modal) / U0)
|
||
max_cross = float(np.max(np.abs(uy)) / U0)
|
||
max_drho = float(np.max(np.abs(rho_c - 1.0)))
|
||
return {
|
||
"amp": amp,
|
||
"ux": ux,
|
||
"uy": uy,
|
||
"rho": rho_c,
|
||
"jx": jx_c,
|
||
"jy": jy_c,
|
||
"nonmodal": nonmodal,
|
||
"max_cross": max_cross,
|
||
"max_drho": max_drho,
|
||
"mass": float(rho_c.sum()),
|
||
"momx": float(jx_c.sum()),
|
||
"momy": float(jy_c.sum()),
|
||
"pi_xy_mean": float(np.nanmean(pi_xy(f, None))),
|
||
}
|
||
|
||
|
||
def run_uniform(kind, phi, transpose=False):
|
||
"""kind: 'B' fine 32^2 or 'C' coarse 16^2. Returns Series + per-tick ux fields."""
|
||
if kind == "B":
|
||
n, dx, tau, sub = N_FINE_DOM, 0.5, TAU_F, 2
|
||
else:
|
||
n, dx, tau, sub = N_COARSE, 1.0, TAU_C, 1
|
||
f = init_uniform(n, dx, tau, phi, transpose=transpose)
|
||
amps = []
|
||
nonmodal = []
|
||
max_cross = []
|
||
max_drho = []
|
||
mass = []
|
||
momx = []
|
||
momy = []
|
||
pi_means = []
|
||
fields = []
|
||
finite_ok = True
|
||
for tick in range(N_STEPS + 1):
|
||
if not np.isfinite(f).all():
|
||
finite_ok = False
|
||
obs = observe_uniform(f, dx, phi, transpose=transpose)
|
||
amps.append(obs["amp"])
|
||
nonmodal.append(obs["nonmodal"])
|
||
max_cross.append(obs["max_cross"])
|
||
max_drho.append(obs["max_drho"])
|
||
mass.append(obs["mass"])
|
||
momx.append(obs["momx"])
|
||
momy.append(obs["momy"])
|
||
pi_means.append(obs["pi_xy_mean"])
|
||
fields.append((obs["ux"].copy(), obs["uy"].copy(), obs["rho"].copy(),
|
||
obs["jx"].copy(), obs["jy"].copy()))
|
||
if tick < N_STEPS:
|
||
f = advance_uniform(f, tau, sub)
|
||
mass = np.asarray(mass)
|
||
momx = np.asarray(momx)
|
||
momy = np.asarray(momy)
|
||
mass_drift = float(np.max(np.abs(mass - mass[0])))
|
||
momx_drift = float(np.max(np.abs(momx - momx[0])))
|
||
momy_drift = float(np.max(np.abs(momy - momy[0])))
|
||
tol_abs = CONS_TOL * abs(mass[0])
|
||
cons_ok = (
|
||
mass_drift <= tol_abs
|
||
and momx_drift <= tol_abs
|
||
and momy_drift <= tol_abs
|
||
)
|
||
ser = Series(
|
||
name=f"{kind}_phi{phi:.4f}{'_T' if transpose else ''}",
|
||
phi=phi,
|
||
A=np.asarray(amps),
|
||
e_if=np.full(N_STEPS + 1, np.nan),
|
||
e_cb=np.full(N_STEPS + 1, np.nan),
|
||
e_if_patch=np.full(N_STEPS + 1, np.nan),
|
||
e_if_out=np.full(N_STEPS + 1, np.nan),
|
||
nonmodal=np.asarray(nonmodal),
|
||
max_uy=np.asarray(max_cross),
|
||
max_drho=np.asarray(max_drho),
|
||
mass=mass,
|
||
momx=momx,
|
||
momy=momy,
|
||
pi_xy_mean=np.asarray(pi_means),
|
||
finite_ok=finite_ok,
|
||
cons_ok=bool(cons_ok),
|
||
mass_drift=mass_drift,
|
||
momx_drift=momx_drift,
|
||
momy_drift=momy_drift,
|
||
)
|
||
return ser, fields
|
||
|
||
|
||
def init_composite(ref, spec, maps, phi, transpose=False):
|
||
"""Equilibrium-only init on composite from cell-average sine velocity."""
|
||
n, nf = spec.n_coarse, spec.n_fine
|
||
coarse = np.full((NQ, n, n), np.nan)
|
||
fine = np.full((NQ, nf, nf), np.nan)
|
||
# coarse cells (including covered): physical centers
|
||
yc = np.arange(n, dtype=np.float64) + 0.5
|
||
xc = np.arange(n, dtype=np.float64) + 0.5
|
||
Yc, Xc = np.meshgrid(yc, xc, indexing="ij")
|
||
if transpose:
|
||
ux_c = np.zeros((n, n))
|
||
uy_c = cell_avg_ux(Xc, 1.0, phi)
|
||
else:
|
||
ux_c = cell_avg_ux(Yc, 1.0, phi)
|
||
uy_c = np.zeros((n, n))
|
||
eq_c = feq(np.ones((n, n)), ux_c, uy_c)
|
||
coarse[:, :, :] = eq_c # covered also filled; advance will NaN them as inactive reads
|
||
# but make_state only fills active — we must match API: covered should be NaN for advance
|
||
coarse[:, maps.coarse_covered] = np.nan
|
||
|
||
# fine active: physical centres from maps.fine_centre (includes patch offset)
|
||
fy = maps.fine_centre # length nf
|
||
fx = maps.fine_centre
|
||
Yf, Xf = np.meshgrid(fy, fx, indexing="ij")
|
||
if transpose:
|
||
ux_f = np.zeros((nf, nf))
|
||
uy_f = cell_avg_ux(Xf, 0.5, phi)
|
||
else:
|
||
ux_f = cell_avg_ux(Yf, 0.5, phi)
|
||
uy_f = np.zeros((nf, nf))
|
||
eq_f = feq(np.ones((nf, nf)), ux_f, uy_f)
|
||
fine[:, maps.fine_active] = eq_f[:, maps.fine_active]
|
||
return ref.State(coarse_f=coarse, fine_f=fine, phase="arrival_f", tick=0)
|
||
|
||
|
||
def observe_composite(ref, maps, state, phi, transpose=False):
|
||
fc, ff = state.coarse_f, state.fine_f
|
||
# Build full 32^2 physical fine field for restriction: outside patch use
|
||
# exploded-from-coarse equivalent = coarse cell values duplicated 2x2.
|
||
# Observables: restrict rho/j then u=j/rho on common 16^2.
|
||
# Outside patch: use coarse as stored. Inside patch: 2x2 mean of fine.
|
||
n = N_COARSE
|
||
rho_c = np.full((n, n), np.nan)
|
||
jx_c = np.full((n, n), np.nan)
|
||
jy_c = np.full((n, n), np.nan)
|
||
|
||
# uncovered coarse
|
||
rho_a, jx_a, jy_a = moments(np.nan_to_num(fc, nan=0.0))
|
||
# only trust active
|
||
for y in range(n):
|
||
for x in range(n):
|
||
if maps.coarse_active[y, x]:
|
||
sl = fc[:, y, x]
|
||
rho_c[y, x] = sl.sum()
|
||
jx_c[y, x] = (CX * sl).sum()
|
||
jy_c[y, x] = (CY * sl).sum()
|
||
|
||
# fine patch -> restrict into covered coarse cells
|
||
# fine active indices a=4..19 map to physical coarse [4,12)
|
||
for ky in range(PATCH_LO, PATCH_HI):
|
||
for kx in range(PATCH_LO, PATCH_HI):
|
||
y0 = spec_children_start(ky)
|
||
x0 = spec_children_start(kx)
|
||
kids = []
|
||
for dy in (0, 1):
|
||
for dx in (0, 1):
|
||
fy, fx = y0 + dy, x0 + dx
|
||
sl = ff[:, fy, fx]
|
||
kids.append((sl.sum(), (CX * sl).sum(), (CY * sl).sum()))
|
||
rho_c[ky, kx] = np.mean([k[0] for k in kids])
|
||
jx_c[ky, kx] = np.mean([k[1] for k in kids])
|
||
jy_c[ky, kx] = np.mean([k[2] for k in kids])
|
||
|
||
ux = jx_c / rho_c
|
||
uy = jy_c / rho_c
|
||
if transpose:
|
||
xc = np.arange(n, dtype=np.float64) + 0.5
|
||
mx = float(sinc(K_PHYS / 2.0)) * np.sin(K_PHYS * xc + phi)
|
||
mT = np.broadcast_to(mx[None, :], (n, n)).copy()
|
||
amp = float(np.sum(uy * mT) / np.sum(mT * mT))
|
||
modal = amp * mT
|
||
nonmodal = float(np.linalg.norm(uy - modal) / U0)
|
||
max_cross = float(np.max(np.abs(ux)) / U0)
|
||
else:
|
||
amp = project_amp(ux, phi)
|
||
m = mode_template(phi)
|
||
modal = amp * m
|
||
nonmodal = float(np.linalg.norm(ux - modal) / U0)
|
||
max_cross = float(np.max(np.abs(uy)) / U0)
|
||
|
||
# leaf totals for conservation (volume-weighted)
|
||
leaf = ref._leaf_totals(fc, ff, maps)
|
||
|
||
return {
|
||
"amp": amp,
|
||
"ux": ux,
|
||
"uy": uy,
|
||
"rho": rho_c,
|
||
"jx": jx_c,
|
||
"jy": jy_c,
|
||
"nonmodal": nonmodal,
|
||
"max_cross": max_cross,
|
||
"max_drho": float(np.max(np.abs(rho_c - 1.0))),
|
||
"mass": float(leaf[0]),
|
||
"momx": float(leaf[1]),
|
||
"momy": float(leaf[2]),
|
||
"finite_active": bool(
|
||
np.isfinite(fc[:, maps.coarse_active]).all()
|
||
and np.isfinite(ff[:, maps.fine_active]).all()
|
||
),
|
||
}
|
||
|
||
|
||
def spec_children_start(k):
|
||
return 4 + 2 * (k - PATCH_LO) # halo=4
|
||
|
||
|
||
def run_composite(ref, phi, transpose=False):
|
||
spec = ref.make_spec(collision="srt")
|
||
# force tau (already default)
|
||
assert abs(spec.tau_coarse - TAU_C) < 1e-15 and abs(spec.tau_fine - TAU_F) < 1e-15
|
||
maps = ref.build_maps(spec)
|
||
st = init_composite(ref, spec, maps, phi, transpose=transpose)
|
||
|
||
# hard check: initial restricted rho,j under patch match coarse formula to round-off
|
||
# Compare fine-restricted vs what coarse cell-average would be under patch
|
||
yc = np.arange(N_COARSE) + 0.5
|
||
xc = np.arange(N_COARSE) + 0.5
|
||
init_mismatch = 0.0
|
||
for ky in range(PATCH_LO, PATCH_HI):
|
||
for kx in range(PATCH_LO, PATCH_HI):
|
||
y0 = spec_children_start(ky)
|
||
x0 = spec_children_start(kx)
|
||
rhos, jxs, jys = [], [], []
|
||
for dy in (0, 1):
|
||
for dx in (0, 1):
|
||
sl = st.fine_f[:, y0 + dy, x0 + dx]
|
||
rhos.append(sl.sum())
|
||
jxs.append((CX * sl).sum())
|
||
jys.append((CY * sl).sum())
|
||
rho_r, jx_r, jy_r = np.mean(rhos), np.mean(jxs), np.mean(jys)
|
||
# coarse formula at (kx,ky)
|
||
if transpose:
|
||
ux = 0.0
|
||
uy = float(cell_avg_ux(xc[kx], 1.0, phi))
|
||
else:
|
||
ux = float(cell_avg_ux(yc[ky], 1.0, phi))
|
||
uy = 0.0
|
||
init_mismatch = max(
|
||
init_mismatch,
|
||
abs(rho_r - 1.0),
|
||
abs(jx_r - ux),
|
||
abs(jy_r - uy),
|
||
)
|
||
|
||
amps, nonmodal, max_cross, max_drho = [], [], [], []
|
||
mass, momx, momy = [], [], []
|
||
fields = []
|
||
finite_ok = True
|
||
for tick in range(N_STEPS + 1):
|
||
obs = observe_composite(ref, maps, st, phi, transpose=transpose)
|
||
if not obs["finite_active"]:
|
||
finite_ok = False
|
||
amps.append(obs["amp"])
|
||
nonmodal.append(obs["nonmodal"])
|
||
max_cross.append(obs["max_cross"])
|
||
max_drho.append(obs["max_drho"])
|
||
mass.append(obs["mass"])
|
||
momx.append(obs["momx"])
|
||
momy.append(obs["momy"])
|
||
fields.append((obs["ux"].copy(), obs["uy"].copy(), obs["rho"].copy(),
|
||
obs["jx"].copy(), obs["jy"].copy()))
|
||
if tick < N_STEPS:
|
||
st, _ = ref.advance_coarse_once(spec, maps, st)
|
||
|
||
mass = np.asarray(mass)
|
||
momx = np.asarray(momx)
|
||
momy = np.asarray(momy)
|
||
mass_drift = float(np.max(np.abs(mass - mass[0])))
|
||
momx_drift = float(np.max(np.abs(momx - momx[0])))
|
||
momy_drift = float(np.max(np.abs(momy - momy[0])))
|
||
tol_abs = CONS_TOL * max(abs(mass[0]), 1.0)
|
||
cons_ok = (
|
||
mass_drift <= tol_abs
|
||
and momx_drift <= tol_abs
|
||
and momy_drift <= tol_abs
|
||
)
|
||
ser = Series(
|
||
name=f"A_phi{phi:.4f}{'_T' if transpose else ''}",
|
||
phi=phi,
|
||
A=np.asarray(amps),
|
||
e_if=np.full(N_STEPS + 1, np.nan),
|
||
e_cb=np.full(N_STEPS + 1, np.nan),
|
||
e_if_patch=np.full(N_STEPS + 1, np.nan),
|
||
e_if_out=np.full(N_STEPS + 1, np.nan),
|
||
nonmodal=np.asarray(nonmodal),
|
||
max_uy=np.asarray(max_cross),
|
||
max_drho=np.asarray(max_drho),
|
||
mass=mass,
|
||
momx=momx,
|
||
momy=momy,
|
||
pi_xy_mean=np.full(N_STEPS + 1, np.nan),
|
||
finite_ok=finite_ok,
|
||
cons_ok=bool(cons_ok),
|
||
mass_drift=mass_drift,
|
||
momx_drift=momx_drift,
|
||
momy_drift=momy_drift,
|
||
notes=f"init_rho_j_mismatch_vs_coarse_formula={init_mismatch:.3e}",
|
||
)
|
||
return ser, fields, init_mismatch
|
||
|
||
|
||
def fit_lambda(amps):
|
||
"""Log-linear fit of A(t) over ticks FIT_LO..FIT_HI; return lambda, stderr, R2."""
|
||
t = np.arange(FIT_LO, FIT_HI + 1, dtype=np.float64)
|
||
y = np.log(np.abs(amps[FIT_LO : FIT_HI + 1]))
|
||
# A(t) ~ A0 exp(-lambda t) => log|A| = log|A0| - lambda t
|
||
n = len(t)
|
||
tbar, ybar = t.mean(), y.mean()
|
||
sxx = np.sum((t - tbar) ** 2)
|
||
slope = np.sum((t - tbar) * (y - ybar)) / sxx # = -lambda
|
||
intercept = ybar - slope * tbar
|
||
pred = intercept + slope * t
|
||
resid = y - pred
|
||
dof = max(n - 2, 1)
|
||
s_err = float(np.sqrt(np.sum(resid**2) / dof / sxx))
|
||
ss_tot = np.sum((y - ybar) ** 2)
|
||
r2 = float(1.0 - np.sum(resid**2) / ss_tot) if ss_tot > 0 else 1.0
|
||
return float(-slope), s_err, r2
|
||
|
||
|
||
def field_l2(a, b):
|
||
return float(np.linalg.norm(a - b))
|
||
|
||
|
||
def main():
|
||
t0 = time.time()
|
||
print("=== shear_diagnostic_v0.2 against CPUREF-0.3 ===")
|
||
print(f"REF_PATH={REF_PATH}")
|
||
print(f"REF_SHA256={sha256_file(REF_PATH)}")
|
||
ref = _load_ref()
|
||
print(f"API_VERSION={getattr(ref, 'API_VERSION', '?')}")
|
||
assert getattr(ref, "API_VERSION", None) == "CPUREF-0.3"
|
||
print(f"U0={U0} tau_c={TAU_C} tau_f={TAU_F} steps={N_STEPS} nu={NU} k={K_PHYS}")
|
||
print(f"nu*k^2={NU * K_PHYS**2:.8f} (analytic decay rate)")
|
||
print("INIT: equilibrium of averaged velocity (NOT exact cell-avg of nonlinear feq)")
|
||
print("UNITS: fit/report 'steps' are coarse steps; CPUREF tick = fine step count")
|
||
print(f" fit window coarse steps {FIT_LO}..{FIT_HI} = CPUREF ticks {2*FIT_LO}..{2*FIT_HI}")
|
||
print(f" startup coarse steps 0..3 = CPUREF ticks 0..6")
|
||
print(" nonmodal/e_IF/e_CB are unnormalised L2 sums over 256 coarse cells")
|
||
print(" (per-cell RMS = value/16); nonmodal is already /U0")
|
||
print(f"CONS_TOL={CONS_TOL} relative to mass (~256) => abs ~{CONS_TOL*256:.3e}")
|
||
print()
|
||
|
||
phis = [0.0, 0.5 * np.pi]
|
||
results = {}
|
||
fields = {}
|
||
init_mm = {}
|
||
|
||
# C and B for both phases
|
||
for phi in phis:
|
||
tag = f"phi{phi:.4f}"
|
||
print(f"-- running C (uniform-coarse) {tag}")
|
||
ser_c, fld_c = run_uniform("C", phi)
|
||
results[f"C_{tag}"] = ser_c
|
||
fields[f"C_{tag}"] = fld_c
|
||
print(f"-- running B (uniform-fine) {tag}")
|
||
ser_b, fld_b = run_uniform("B", phi)
|
||
results[f"B_{tag}"] = ser_b
|
||
fields[f"B_{tag}"] = fld_b
|
||
print(f"-- running A (composite CPUREF) {tag}")
|
||
ser_a, fld_a, mism = run_composite(ref, phi)
|
||
init_mm[tag] = mism
|
||
results[f"A_{tag}"] = ser_a
|
||
fields[f"A_{tag}"] = fld_a
|
||
# fill e_IF = ||A-B||, e_CB = ||C-B|| on restricted ux (unnormalised L2)
|
||
patch = np.zeros((N_COARSE, N_COARSE), dtype=bool)
|
||
patch[PATCH_LO:PATCH_HI, PATCH_LO:PATCH_HI] = True
|
||
out = ~patch
|
||
for t in range(N_STEPS + 1):
|
||
du = fld_a[t][0] - fld_b[t][0]
|
||
ser_a.e_if[t] = float(np.linalg.norm(du))
|
||
ser_a.e_if_patch[t] = float(np.linalg.norm(du[patch]))
|
||
ser_a.e_if_out[t] = float(np.linalg.norm(du[out]))
|
||
ser_a.e_cb[t] = field_l2(fld_c[t][0], fld_b[t][0])
|
||
ser_c.e_cb[t] = ser_a.e_cb[t]
|
||
ser_b.e_if[t] = 0.0
|
||
|
||
# A transpose at phi=pi/2
|
||
phi_t = 0.5 * np.pi
|
||
print(f"-- running A transpose phi={phi_t}")
|
||
ser_at, fld_at, mism_t = run_composite(ref, phi_t, transpose=True)
|
||
init_mm["A_T"] = mism_t
|
||
results["A_T_phi1.5708"] = ser_at
|
||
fields["A_T_phi1.5708"] = fld_at
|
||
|
||
# Also need B transpose? For e_IF on transpose — freeze says A transposed only.
|
||
# Transpose hard check: A_T fields equal transpose of A at same phi.
|
||
ser_a_ref = results["A_phi1.5708"]
|
||
fld_a_ref = fields["A_phi1.5708"]
|
||
trans_err = 0.0
|
||
for t in range(N_STEPS + 1):
|
||
ux, uy, rho, jx, jy = fld_a_ref[t]
|
||
ux_t, uy_t, rho_t, jx_t, jy_t = fld_at[t]
|
||
# original: ux(y), uy~0; transpose: uy(x), ux~0
|
||
# transpose of original velocity: u'_x(x,y)=uy(y,x), u'_y(x,y)=ux(y,x)
|
||
err = max(
|
||
float(np.max(np.abs(ux_t - uy.T))),
|
||
float(np.max(np.abs(uy_t - ux.T))),
|
||
float(np.max(np.abs(rho_t - rho.T))),
|
||
)
|
||
trans_err = max(trans_err, err)
|
||
transpose_ok = trans_err <= TRANS_TOL
|
||
|
||
# Convention match A/B/C at t=0: restricted rho,j under patch
|
||
print("\n=== HARD CHECKS ===")
|
||
hard = {}
|
||
hard_keys = []
|
||
for phi in phis:
|
||
tag = f"phi{phi:.4f}"
|
||
hard_keys.extend([f"A_{tag}", f"B_{tag}", f"C_{tag}"])
|
||
hard_keys.append("A_T_phi1.5708")
|
||
for key in hard_keys:
|
||
s = results[key]
|
||
hard[f"finite_{key}"] = s.finite_ok
|
||
hard[f"cons_{key}"] = s.cons_ok
|
||
tol_abs = CONS_TOL * max(abs(s.mass[0]), 1.0)
|
||
print(
|
||
f" finite[{key}]={s.finite_ok} cons={s.cons_ok} "
|
||
f"mass_drift={s.mass_drift:.3e} momx_drift={s.momx_drift:.3e} "
|
||
f"momy_drift={s.momy_drift:.3e} tol_abs={tol_abs:.3e} "
|
||
f"(CONS_TOL*mass={CONS_TOL}*|{s.mass[0]:.6g}|) {s.notes}"
|
||
)
|
||
hard["transpose_ok"] = transpose_ok
|
||
hard["transpose_err"] = trans_err
|
||
print(f" transpose_ok={transpose_ok} max_err={trans_err:.3e} (tol={TRANS_TOL:.3e})")
|
||
|
||
# t=0 convention: A vs B vs C amplitudes and restricted fields
|
||
for phi in phis:
|
||
tag = f"phi{phi:.4f}"
|
||
a0 = results[f"A_{tag}"].A[0]
|
||
b0 = results[f"B_{tag}"].A[0]
|
||
c0 = results[f"C_{tag}"].A[0]
|
||
print(f" t0 amp A/B/C {tag}: {a0:.8e} {b0:.8e} {c0:.8e}")
|
||
# restricted j under patch: A vs theoretical
|
||
ux_a = fields[f"A_{tag}"][0][0]
|
||
ux_b = fields[f"B_{tag}"][0][0]
|
||
ux_c = fields[f"C_{tag}"][0][0]
|
||
# under patch and outside
|
||
patch = np.zeros((16, 16), dtype=bool)
|
||
patch[PATCH_LO:PATCH_HI, PATCH_LO:PATCH_HI] = True
|
||
d_ab_p = float(np.max(np.abs(ux_a[patch] - ux_b[patch])))
|
||
d_ac_o = float(np.max(np.abs(ux_a[~patch] - ux_c[~patch])))
|
||
d_bc = float(np.max(np.abs(ux_b - ux_c)))
|
||
print(f" t0 max|ux A-B| patch={d_ab_p:.3e} |A-C| outside={d_ac_o:.3e} |B-C| all={d_bc:.3e}")
|
||
hard[f"t0_AB_patch_{tag}"] = d_ab_p
|
||
hard[f"t0_AC_out_{tag}"] = d_ac_o
|
||
hard[f"t0_AB_patch_ok_{tag}"] = bool(d_ab_p <= T0_TOL)
|
||
hard[f"t0_AC_out_ok_{tag}"] = bool(d_ac_o <= T0_TOL)
|
||
hard[f"init_mismatch_ok_{tag}"] = bool(init_mm[tag] <= INIT_TOL)
|
||
print(f" t0_AB_patch_ok={hard[f't0_AB_patch_ok_{tag}']} t0_AC_out_ok={hard[f't0_AC_out_ok_{tag}']} "
|
||
f"(T0_TOL={T0_TOL:.1e}) init_mismatch={init_mm[tag]:.3e} ok={hard[f'init_mismatch_ok_{tag}']} "
|
||
f"(INIT_TOL={INIT_TOL:.1e})")
|
||
hard["init_mismatch_ok_A_T_phi1.5708"] = bool(init_mm["A_T"] <= INIT_TOL)
|
||
print(f" init_mismatch[A_T]={init_mm['A_T']:.3e} ok={hard['init_mismatch_ok_A_T_phi1.5708']} (INIT_TOL={INIT_TOL:.1e})")
|
||
_bool_keys = sorted(k for k, v in hard.items() if isinstance(v, bool))
|
||
print(f" N_BOOLEAN_HARD_CHECKS={len(_bool_keys)}: " + ",".join(_bool_keys))
|
||
|
||
print("\n=== DECAY FITS (coarse steps 4..32 = CPUREF ticks 8..64) ===")
|
||
lam_analytic = NU * K_PHYS**2
|
||
print(f" analytic lambda = nu*k^2 = {lam_analytic:.8f}")
|
||
fit_table = {}
|
||
patch = np.zeros((N_COARSE, N_COARSE), dtype=bool)
|
||
patch[PATCH_LO:PATCH_HI, PATCH_LO:PATCH_HI] = True
|
||
for phi in phis:
|
||
tag = f"phi{phi:.4f}"
|
||
row = {}
|
||
for kind in ("A", "B", "C"):
|
||
s = results[f"{kind}_{tag}"]
|
||
lam, se, r2 = fit_lambda(s.A)
|
||
row[kind] = (lam, se, r2, s.A[0], s.A[-1])
|
||
print(
|
||
f" {kind} {tag}: lambda={lam:.8f} stderr={se:.3e} R2={r2:.6f} "
|
||
f"A0={s.A[0]:.6e} A32={s.A[-1]:.6e}"
|
||
)
|
||
fit_table[tag] = row
|
||
la, lb, lc = row["A"][0], row["B"][0], row["C"][0]
|
||
s_flag = (la - lb) / (lc - lb) if abs(lc - lb) > 1e-16 else float("nan")
|
||
print(f" DIAG s=(lambda_A-lambda_B)/(lambda_C-lambda_B) {tag}: s={s_flag:.4f} "
|
||
f"(lambda_C-lambda_B)={lc-lb:.6e}")
|
||
sa = results[f"A_{tag}"]
|
||
print(
|
||
f" DIAG e_IF/e_CB step32 (unnorm L2/256cells): e_IF={sa.e_if[-1]:.6e} "
|
||
f"e_IF_patch={sa.e_if_patch[-1]:.6e} e_IF_out={sa.e_if_out[-1]:.6e} "
|
||
f"e_CB={sa.e_cb[-1]:.6e} ratio_IF/CB={sa.e_if[-1]/max(sa.e_cb[-1],1e-30):.4f}"
|
||
)
|
||
n_p, n_o = int(patch.sum()), int((~patch).sum())
|
||
rms_p = sa.e_if_patch[-1] / np.sqrt(n_p)
|
||
rms_o = sa.e_if_out[-1] / np.sqrt(n_o)
|
||
frac_p = sa.e_if_patch[-1] ** 2 / max(sa.e_if[-1] ** 2, 1e-300)
|
||
print(
|
||
f" DIAG per-cell RMS (uA-uB) step32 {tag}: patch={rms_p:.4e} ({n_p} cells) "
|
||
f"outside={rms_o:.4e} ({n_o} cells) ratio={rms_p/rms_o:.3f} "
|
||
f"patch_energy_frac={frac_p:.4f} (area frac {n_p/(n_p+n_o):.2f})"
|
||
)
|
||
print(
|
||
f" DIAG sensors step32: nonmodal(L2/U0)={sa.nonmodal[-1]:.4e} "
|
||
f"max|uy|/U0={sa.max_uy[-1]:.4e} max|drho|={sa.max_drho[-1]:.4e}"
|
||
)
|
||
# x-averaged (u_A - u_B)(y) at step 32
|
||
du = fields[f"A_{tag}"][32][0] - fields[f"B_{tag}"][32][0]
|
||
profile = du.mean(axis=1) # average over x
|
||
print(f" DIAG (uA-uB)_xavg(y) step32 {tag}: " +
|
||
" ".join(f"{v:.3e}" for v in profile))
|
||
# Pi_xy at step 32 for uniforms (arrival_f / post-stream)
|
||
sb = results[f"B_{tag}"]
|
||
sc = results[f"C_{tag}"]
|
||
print(
|
||
f" DIAG Pi_xy_mean step32 {tag}: B={sb.pi_xy_mean[32]:.6e} "
|
||
f"C={sc.pi_xy_mean[32]:.6e} (composite=NaN; arrival_f proxy on uniforms)"
|
||
)
|
||
|
||
print("\n=== PER-TICK AMPLITUDE (selected) ===")
|
||
hdr = "tick " + " ".join(
|
||
f"{k:>14}" for k in ["A_p0", "B_p0", "C_p0", "A_p90", "B_p90", "C_p90", "exact"]
|
||
)
|
||
print(hdr)
|
||
for t in list(range(0, 4)) + list(range(4, 33, 4)):
|
||
exact = results["B_phi0.0000"].A[0] * np.exp(-lam_analytic * t)
|
||
vals = [
|
||
results["A_phi0.0000"].A[t],
|
||
results["B_phi0.0000"].A[t],
|
||
results["C_phi0.0000"].A[t],
|
||
results["A_phi1.5708"].A[t],
|
||
results["B_phi1.5708"].A[t],
|
||
results["C_phi1.5708"].A[t],
|
||
exact,
|
||
]
|
||
print(f"{t:4d} " + " ".join(f"{v:14.6e}" for v in vals))
|
||
|
||
print("\n=== NUMERIC SUMMARY ===")
|
||
all_hard_pass = all(
|
||
v for k, v in hard.items() if isinstance(v, bool)
|
||
)
|
||
print(f"ALL_HARD_CHECKS_PASS={all_hard_pass}")
|
||
print(f"runtime_sec={time.time()-t0:.2f}")
|
||
# compact JSON-ish summary lines for Drive
|
||
print("SUMMARY_LINES_BEGIN")
|
||
for phi in phis:
|
||
tag = f"phi{phi:.4f}"
|
||
for kind in ("A", "B", "C"):
|
||
lam, se, r2, a0, a32 = fit_table[tag][kind]
|
||
print(
|
||
f"FIT {kind} {tag} lambda={lam:.10f} stderr={se:.6e} R2={r2:.8f} "
|
||
f"A0={a0:.10e} A32={a32:.10e}"
|
||
)
|
||
print(f"TRANSPOSE err={trans_err:.6e} ok={transpose_ok}")
|
||
print(f"ANALYTIC lambda={lam_analytic:.10f}")
|
||
print("SUMMARY_LINES_END")
|
||
return 0 if all_hard_pass else 1
|
||
|
||
|
||
if __name__ == "__main__":
|
||
sys.exit(main())
|