"""Bounded static-slot validation; frozen gates in research/j2-slot-2026-10-12/REVIEW.md.""" import argparse import json import os from pathlib import Path import sys from tempfile import TemporaryDirectory from unittest.mock import patch import numpy as np import CelerisLab from CelerisLab.cuda import compiler_v2 as compiler from run_j1 import ROOT, make_sim, populations, observables, population_budget, off_path_check NX, WIDTH, CHUNK = 48, 16, 519 GATES = dict(velocity_change=5e-5, mass_change=2e-6, mass_window_span=6e-6, mass_offset=1e-3, rho_error=.005, speed=.01, zero_speed=2e-6, flux_error=.005, flux_span=.002, boundary_error=2e-6, rho_symmetry=2e-5, velocity_symmetry=.002, budget_residual=2e-7) def targets(amplitude): # Independent discrete quadrature, including the zero-target flanks. x = np.arange(NX)-(NX-1)/2 bottom = amplitude*np.clip(1-(x/(WIDTH/2))**2, 0, None)**2 return bottom, bottom.mean(), np.abs(x) < WIDTH/2 def full_width_check(baseline): snapshots = [] for source, width in ((baseline, 0), (Path(compiler.kernel_path('')), NX)): with patch.object(compiler, '_KERNEL_DIR', str(source)): with make_sim(.002, nx=NX, top_u=0., slot_width=width) as sim: sim.initialize() initial = sim.get_ddf().copy() sim.run(64) snapshots.append((initial, sim.get_ddf().copy(), sim.get_flags().copy())) return bool(all(np.array_equal(a, b) for a, b in zip(*snapshots))) def restart_check(sim, amplitude, top_mode='velocity'): with TemporaryDirectory(prefix='celeris-j2-checkpoint-') as folder: path = sim.save_checkpoint(str(Path(folder)/'slot.h5')) sim.run(7) expected, flags, count = sim.get_ddf().copy(), sim.get_flags().copy(), sim.stepper.step_count with make_sim(amplitude, nx=NX, top_u=0., slot_width=WIDTH, top_mode=top_mode) as target: target.load_checkpoint(path) target.run(7) equal = (np.array_equal(expected, target.get_ddf()) and np.array_equal(flags, target.get_flags()) and count == target.stepper.step_count) sim.load_checkpoint(path) return bool(equal) def run_case(amplitude, top_mode='velocity', check_restart=True): bottom, top, slot = targets(amplitude) scale = max(abs(amplitude), .002) with make_sim(amplitude, nx=NX, top_u=0., slot_width=WIDTH, top_mode=top_mode) as sim: sim.initialize() rho, ux, uy = observables(populations(sim)) initial_mass = float(rho.sum()) initial_error = float(max(np.abs(uy[0]-bottom).max(), np.abs(uy[-1]-top).max(), np.abs(ux).max())) old_u, old_v, old_mass = ux.copy(), uy.copy(), initial_mass history = [dict(step=0, mass=initial_mass)] steady = False for sample in range(1, 41): sim.run(CHUNK) f = populations(sim) assert np.isfinite(f).all(), 'Nonfinite slot populations' rho, ux, uy = observables(f) mass = float(rho.sum()) history.append(dict(step=sim.stepper.step_count, mass=mass, u_change=float(np.sqrt(np.mean((ux-old_u)**2))/scale), v_change=float(np.sqrt(np.mean((uy-old_v)**2))/scale), mass_change=abs(mass-old_mass)/initial_mass)) old_u, old_v, old_mass = ux.copy(), uy.copy(), mass span = float(np.ptp([r['mass'] for r in history[-5:]])/initial_mass) steady = bool(sample >= 8 and span <= GATES['mass_window_span'] and all( r['u_change'] <= GATES['velocity_change'] and r['v_change'] <= GATES['velocity_change'] and r['mass_change'] <= GATES['mass_change'] for r in history[-4:])) if steady: break sections = [8, 16, 24] actual_mass_flux = (rho*uy).sum(axis=1) actual_volume_flux = uy.sum(axis=1) request_q = float(bottom.sum()) reference = float(rho.mean()*request_q) flux_error = float(np.max(np.abs(actual_mass_flux[sections]-reference))/abs(reference)) if amplitude else 0. flux_span = float(np.ptp(actual_mass_flux[sections])/abs(reference)) if amplitude else 0. symmetry = dict(rho=float(np.max(np.abs(rho-rho[:, ::-1]))), u=float(np.max(np.abs(ux+ux[:, ::-1]))/scale), v=float(np.max(np.abs(uy-uy[:, ::-1]))/scale)) boundary_error = float(max(np.abs(uy[0]-bottom).max(), np.abs(uy[-1]-top).max() if top_mode == 'velocity' else 0., np.abs(ux[[0, -1]]).max())) top_density_error = float(np.abs(rho[-1]-sim.lbm_cfg.rho).max()/sim.lbm_cfg.rho) offset = abs(mass-initial_mass)/initial_mass speed = float(np.hypot(ux, uy).max()) density_error = float(np.abs(rho-1).max()) nonuniform = bool(uy[8, NX//2]*amplitude > 0 and np.abs(ux).max() > 1e-5) if amplitude else True restart = restart_check(sim, amplitude, top_mode) if check_restart and amplitude > 0 else None budget = population_budget(sim, bottom_slot=slot) checks = dict(steady=steady, initialization=initial_error <= GATES['boundary_error'], boundary=boundary_error <= GATES['boundary_error'], mass=offset <= GATES['mass_offset'] and span <= GATES['mass_window_span'], bounded=density_error <= GATES['rho_error'] and speed <= GATES['speed'], flux=flux_error <= GATES['flux_error'] and flux_span <= GATES['flux_span'], symmetry=symmetry['rho'] <= GATES['rho_symmetry'] and max(symmetry['u'], symmetry['v']) <= GATES['velocity_symmetry'], audit=budget['instrumented_matches_plain'] and budget['max_relative_residual'] <= GATES['budget_residual']) if amplitude: checks['nonuniform'] = nonuniform else: checks.pop('flux') # Relative flux error has no nonzero reference here. checks['zero_speed'] = speed <= GATES['zero_speed'] if restart is not None: checks['restart'] = restart if top_mode == 'reservoir': checks['top_density'] = top_density_error <= 5e-7 checks = {name: bool(passed) for name, passed in checks.items()} return dict(amplitude=amplitude, width=WIDTH, nx=NX, height=32, viscosity=.1, top_u=0., steps=sim.stepper.step_count-16, max_steps=40*CHUNK, checks=checks, passed=all(checks.values()), requested_volume=request_q, top_initial_guess=top, top_target=top if top_mode == 'velocity' else None, bottom_targets=bottom.tolist(), velocity_scale=scale, rho_scale=1., initial_boundary_error=initial_error, top_mode=top_mode, top_density_error=top_density_error, top_velocity_range=[float(uy[-1].min()), float(uy[-1].max())], boundary_error=boundary_error, mass_offset=offset, mass_window_span=span, flux_error=flux_error, flux_span=flux_span, section_y=sections, mass_flux_by_y=actual_mass_flux.tolist(), volume_flux_by_y=actual_volume_flux.tolist(), symmetry=symmetry, max_speed=speed, rho_error=density_error, rho=rho.tolist(), u=ux.tolist(), v=uy.tolist(), history=history, checkpoint_bitwise_equal=restart, population_budget=budget) def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--output', type=Path, required=True) parser.add_argument('--baseline-kernels', type=Path, required=True, help='Pre-J2 source for the required full-width and default-mode comparisons') args = parser.parse_args() assert Path.cwd() == ROOT == Path('/home/frank14f/CelerisLab') assert sys.executable == '/home/frank14f/anaconda3/envs/pycuda_3_10/bin/python' assert Path(CelerisLab.__file__).resolve() == ROOT/'src/CelerisLab/__init__.py' assert Path(compiler.kernel_path('kernel_v2.cu')).resolve() == ROOT/'src/CelerisLab/lbm/kernels/kernel_v2.cu' assert os.environ.get('CUDA_VISIBLE_DEVICES') in ('0', '1') report = dict(python=sys.executable, package=CelerisLab.__file__, cuda_source=compiler.kernel_path('kernel_v2.cu'), physical_gpu=os.environ['CUDA_VISIBLE_DEVICES'], gates=GATES, full_width_bitwise_equal=None, off_path_bitwise_equal=None, cases=[]) report['full_width_bitwise_equal'] = full_width_check(args.baseline_kernels) report['off_path_bitwise_equal'] = off_path_check(args.baseline_kernels) for amplitude in (0., .002, -.002): case = run_case(amplitude) report['cases'].append(case) args.output.write_text(json.dumps(report, indent=2)+'\n') print(json.dumps({k: case[k] for k in ('amplitude', 'steps', 'checks', 'mass_offset', 'mass_window_span', 'flux_error', 'symmetry')}), flush=True) return 0 if (report['full_width_bitwise_equal'] is True and report['off_path_bitwise_equal'] is True and all(c['passed'] for c in report['cases'])) else 1 if __name__ == '__main__': raise SystemExit(main())