"""One bounded J3 restart/zero/negative run; frozen plan in research/j3-coverage-2026-10-12.""" from dataclasses import asdict import hashlib import json import os from pathlib import Path import resource import subprocess import sys import traceback import h5py import numpy as np import pycuda.driver as cuda import CelerisLab from CelerisLab.cuda import compiler_v2 as compiler from run_j1 import ROOT, AUDIT, make_sim, populations, observables from run_j2_slot import GATES, targets OUT = ROOT/'research/j3-coverage-2026-10-12' SOURCE = ROOT/'research/j3-lag-2026-10-12/end-21280.h5' SOURCE_SHA = '19a94e9e8db6dd9ac981effe06032c351992250b58768d9c090409193893d588' SCALE, CHUNK, LIMIT = .002, 519, 20760 def sha(path): digest = hashlib.sha256() with Path(path).open('rb') as stream: for chunk in iter(lambda: stream.read(1024*1024), b''): digest.update(chunk) return digest.hexdigest() def create(amplitude): return make_sim(amplitude, nx=48, top_u=0., slot_width=16, top_mode='reservoir') def state(sim): arrays = dict(ddf=np.empty(sim.field.n*9, dtype=np.float32), temp=np.empty(sim.field.n*9, dtype=np.float32), flag=np.empty(sim.field.n, dtype=np.uint16)) with sim.ctx.activate(): for name, array in arrays.items(): cuda.memcpy_dtoh(array, getattr(sim.field, name+'_gpu')) return dict(arrays=arrays, step=sim.stepper.step_count, config=asdict(sim.lbm_cfg), runtime=sim.field._runtime_params.copy()) def compare(a, b): checks = {name: a['arrays'][name].tobytes() == b['arrays'][name].tobytes() for name in a['arrays']} checks.update({name: a[name] == b[name] for name in ('step', 'config', 'runtime')}) return checks def persist(sim, name): path = OUT/name assert not path.exists(), path sim.save_checkpoint(str(path)) return dict(path=str(path), sha256=sha(path), step=sim.stepper.step_count) def measure(f, amplitude, initial_mass, previous=None): rho, u, v = observables(f) bottom, top, _ = targets(amplitude) mass = float(rho.sum()) flux = (rho*v).sum(axis=1) reference = float(rho.mean()*bottom.sum()) row = dict(mass=mass, mass_offset=abs(mass-initial_mass)/initial_mass, rho_range=[float(rho.min()), float(rho.max())], rho_error=float(np.abs(rho-1).max()), max_speed=float(np.hypot(u, v).max()), u_mean=float(u.mean()), v_mean=float(v.mean()), boundary_error=float(max(np.abs(v[0]-bottom).max(), np.abs(u[[0, -1]]).max())), top_density_error=float(np.abs(rho[-1]-1).max()), top_velocity_range=[float(v[-1].min()), float(v[-1].max())], mass_flux_by_y=flux.tolist(), volume_flux_by_y=v.sum(axis=1).tolist(), flux_error=float(np.max(np.abs(flux[[8, 16, 24]]-reference))/abs(reference)) if amplitude else None, flux_span=float(np.ptp(flux[[8, 16, 24]])/abs(reference)) if amplitude else None, symmetry=dict(rho=float(np.max(np.abs(rho-rho[:, ::-1]))), u=float(np.max(np.abs(u+u[:, ::-1]))/SCALE), v=float(np.max(np.abs(v-v[:, ::-1]))/SCALE)), center_v=float(v[8, 24]), max_abs_u=float(np.abs(u).max()), nonuniform=bool(v[8, 24]*amplitude > 0 and np.abs(u).max() > 1e-5) if amplitude else None, finite=bool(np.isfinite(f).all() and np.isfinite(rho).all() and np.isfinite(u).all() and np.isfinite(v).all())) if previous is not None: row.update(u_change=float(np.sqrt(np.mean((u-previous[0])**2))/SCALE), v_change=float(np.sqrt(np.mean((v-previous[1])**2))/SCALE), u_signed_change=float((u-previous[0]).mean()), v_signed_change=float((v-previous[1]).mean())) return row, (rho, u, v) def check_step(sim, mass0, extrema): f = populations(sim) rho, u, v = observables(f) values = dict(rho_error=float(np.abs(rho-1).max()), speed=float(np.hypot(u, v).max()), mass_offset=abs(float(rho.sum())-mass0)/mass0) finite = bool(np.isfinite(f).all() and all(np.isfinite(a).all() for a in (rho, u, v))) for key, value in values.items(): extrema[key] = max(extrema.get(key, 0.), value) if not finite or any(value > GATES[key] for key, value in values.items()): raise RuntimeError(f'Hard stop at step {sim.stepper.step_count}: finite={finite}, {values}') return f, (rho, u, v) def restart(report, save): with h5py.File(SOURCE, 'r') as hf: assert int(hf.attrs['step_count']) == 21280 assert str(hf.attrs['distribution_phase']).startswith('postcollision g(t)') config, runtime = json.loads(hf.attrs['config_json']), json.loads(hf.attrs['runtime_params_json']) saved = {name: hf[name][:] for name in ('ddf', 'temp', 'flag')} row = report['restart'] = dict(status='running', source=str(SOURCE), source_sha256=sha(SOURCE), start=21280, end=21296, a_steps=0, b_steps=0, extrema={}) for label in ('a', 'b'): with create(.002) as sim: sim.load_checkpoint(str(SOURCE if label == 'a' else OUT/'restart-21280.h5')) before = state(sim) assert before['config'] == config and before['runtime'] == runtime checks = compare(dict(arrays=saved, step=21280, config=config, runtime=runtime), before) row[label+'_initial_equal'] = checks assert all(checks.values()), checks mass0 = float(observables(populations(sim))[0].sum()) check_step(sim, mass0, row['extrema']) if label == 'a': row['checkpoint'] = persist(sim, 'restart-21280.h5') for _ in range(16): sim.run(1) row[label+'_steps'] += 1 check_step(sim, mass0, row['extrema']) if label == 'a': expected = state(sim) row['reference_end'] = persist(sim, 'restart-reference-21296.h5') else: row['continued_equal'] = compare(expected, state(sim)) row['status'] = 'passed' if all(row['continued_equal'].values()) else 'failed' save() assert row['status'] == 'passed', row['continued_equal'] print('RESTART', json.dumps(row), flush=True) def zero(report, save): row = report['zero'] = dict(status='running', steps=0, history=[], extrema={}, steady='not_tested_short_case', flux_relative='not_applicable_zero_reference') with create(0.) as sim: sim.initialize() mass0 = float(observables(populations(sim))[0].sum()) for step in range(65): if step: sim.run(1) row['steps'] = step f, fields = check_step(sim, mass0, row['extrema']) sample, _ = measure(f, 0., mass0) sample['step'] = step row['history'].append(sample) if sample['max_speed'] > GATES['zero_speed']: raise RuntimeError(f'Zero speed exceeds gate at {step}') if np.ptp([s['mass'] for s in row['history']])/mass0 > GATES['mass_window_span']: raise RuntimeError(f'Zero mass span exceeds gate at {step}') if step and abs(sample['mass']-row['history'][-2]['mass'])/mass0 > GATES['mass_change']: raise RuntimeError(f'Zero mass change exceeds gate at {step}') row['mass_span'] = float(np.ptp([s['mass'] for s in row['history']])/mass0) row['checks'] = dict(zero_speed=row['extrema']['speed'] <= GATES['zero_speed'], mass_offset=row['extrema']['mass_offset'] <= GATES['mass_offset'], mass_span=row['mass_span'] <= GATES['mass_window_span'], mass_change=max(abs(s['mass']-p['mass'])/mass0 for p,s in zip(row['history'],row['history'][1:])) <= GATES['mass_change'], boundary=max(s['boundary_error'] for s in row['history']) <= GATES['boundary_error'], top_density=max(s['top_density_error'] for s in row['history']) <= 5e-7, finite=all(s['finite'] for s in row['history'])) row['status'] = 'passed' if all(row['checks'].values()) else 'failed' save() assert row['status'] == 'passed', row['checks'] print('ZERO', json.dumps({k:v for k,v in row.items() if k != 'history'}), flush=True) def audit(reference, mass0, report): result = dict(status='running', start=LIMIT-16, end=LIMIT, plain_steps=16, audit_steps=0, steps=[]) report['negative']['population_budget'] = result with create(-.002) as sim: sim.load_checkpoint(str(OUT/'negative-budget-20744.h5')) build = compiler.generate_config(sim.lbm_cfg) original = sim.stepper.step_fn try: source = Path(build)/'step/one_step_double.cu' code = source.read_text() call, marker = 'collide_dispatch(f, rho_n, ux, uy);', '__global__ __launch_bounds__(NT, 2)' assert code.count(call) == code.count(marker) == 1 source.write_text(code.replace(call, 'j1_collide_audit(f, rho_n, ux, uy);').replace(marker, AUDIT+'\n'+marker)) result['audit_source_sha256'] = sha(source) ptx = compiler.compile_kernel(build_dir=build) with sim.ctx.activate(): module = compiler.load_module(ptx) ptr, size = module.get_global('d_params') cuda.memcpy_dtod(ptr, sim.field.module.get_global('d_params')[0], size) iptr, _ = module.get_global('j1_incoming') cptr, size = module.get_global('j1_collision') cuda.memset_d8(cptr, 0, size) incoming, collision = np.empty((2,48)), np.empty((35,48)) sim.stepper.step_fn = module.get_function('OneStep') slot = targets(-.002)[2] for _ in range(16): before = populations(sim) outgoing = np.stack((before[[4,6,7],1].sum(axis=0), before[[3,5,8],-2].sum(axis=0))) mass = float(before[:,1:-1].sum()) sim.run(1) result['audit_steps'] += 1 after, _ = check_step(sim, mass0, report['negative']['extrema']) cuda.memcpy_dtoh(incoming, iptr) cuda.memcpy_dtoh(collision, cptr) delta = float(after[:,1:-1].sum()-mass) net = float(incoming.sum()-outgoing.sum()) rounding = float(collision.sum()) entry = dict(step=sim.stepper.step_count, mass=mass, delta_mass=delta, net_boundary=net, collision_roundoff=rounding, physical_residual=delta-net, budget_residual=delta-net-rounding, exchange={name:dict(incoming=float(i.sum()), outgoing=float(o.sum())) for name,i,o in [('slot',incoming[0,slot],outgoing[0,slot]), ('flanks',incoming[0,~slot],outgoing[0,~slot]), ('top',incoming[1],outgoing[1])]}) result['steps'].append(entry) if abs(entry['budget_residual'])/mass > GATES['budget_residual']: raise RuntimeError('Population budget gate exceeded') result['continued_equal'] = compare(reference, state(sim)) result['max_relative_residual'] = max(abs(s['budget_residual'])/s['mass'] for s in result['steps']) result['max_physical_relative_residual'] = max(abs(s['physical_residual'])/s['mass'] for s in result['steps']) result['status'] = 'passed' if all(result['continued_equal'].values()) else 'failed' finally: sim.stepper.step_fn = original compiler.cleanup_build_dir(build) return result def negative(report, save): row = report['negative'] = dict(status='running', amplitude=-.002, steps=0, main_steps=0, history=[], extrema={}, population_budget=dict(status='not_tested')) with create(-.002) as sim: sim.initialize() initial = populations(sim) mass0 = float(observables(initial)[0].sum()) first, fields = measure(initial, -.002, mass0) first.update(step=0) row['history'].append(first) bottom, top, _ = targets(-.002) row['initial_error'] = float(max(np.abs(fields[2][0]-bottom).max(), np.abs(fields[2][-1]-top).max(), np.abs(fields[1]).max())) previous, mass_previous = fields[1:], mass0 check_step(sim, mass0, row['extrema']) try: for interval in range(1,41): mean = np.zeros((2,33,48)); m2 = np.zeros_like(mean) low, high = np.full_like(mean, np.inf), np.full_like(mean, -np.inf) for n in range(1,CHUNK+1): if sim.stepper.step_count == LIMIT-16: row['budget_start'] = persist(sim, 'negative-budget-20744.h5') sim.run(1) row['steps'] = row['main_steps'] = sim.stepper.step_count f, fields = check_step(sim, mass0, row['extrema']) uv = np.stack(fields[1:]) delta = uv-mean mean += delta/n m2 += delta*(uv-mean) low, high = np.minimum(low,uv), np.maximum(high,uv) sample, fields = measure(f, -.002, mass0, previous) sample.update(step=sim.stepper.step_count, mass_change=abs(sample['mass']-mass_previous)/mass0, fluctuation=dict(samples=CHUNK, first_step=sim.stepper.step_count-CHUNK+1, rms=np.sqrt(np.mean(m2/CHUNK,axis=(1,2))).tolist(), max_node_peak_to_peak=np.max(high-low,axis=(1,2)).tolist())) row['history'].append(sample) sample['mass_window_span'] = float(np.ptp([s['mass'] for s in row['history'][-5:]])/mass0) previous, mass_previous = fields[1:], sample['mass'] save() print('NEGATIVE', json.dumps({k:sample[k] for k in ('step','mass_offset','mass_change','mass_window_span','u_change','v_change','flux_error','boundary_error')}), flush=True) if sample['mass_change'] > GATES['mass_change'] or sample['mass_window_span'] > GATES['mass_window_span']: raise RuntimeError('Negative sampled mass hard gate exceeded') assert resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024 < 512 reference = state(sim) row['end_checkpoint'] = persist(sim, 'negative-end-20760.h5') row['final_boundary'] = dict(bottom_rho=fields[0][0].tolist(), bottom_u=fields[1][0].tolist(), bottom_v=fields[2][0].tolist(), top_rho=fields[0][-1].tolist(), top_u=fields[1][-1].tolist(), top_v=fields[2][-1].tolist()) except Exception: row['failure_checkpoint'] = persist(sim, 'negative-failure.h5') raise budget = audit(reference, mass0, report) last, tail = row['history'][-1], row['history'][-4:] checks = dict(steady=len(row['history']) >= 9 and all(s['u_change'] <= GATES['velocity_change'] and s['v_change'] <= GATES['velocity_change'] and s['mass_change'] <= GATES['mass_change'] for s in tail) and last['mass_window_span'] <= GATES['mass_window_span'], initialization=row['initial_error'] <= GATES['boundary_error'], boundary=last['boundary_error'] <= GATES['boundary_error'], top_density=last['top_density_error'] <= 5e-7, mass=last['mass_offset'] <= GATES['mass_offset'] and last['mass_window_span'] <= GATES['mass_window_span'], bounded=last['rho_error'] <= GATES['rho_error'] and last['max_speed'] <= GATES['speed'], flux=last['flux_error'] <= GATES['flux_error'] and last['flux_span'] <= GATES['flux_span'], symmetry=last['symmetry']['rho'] <= GATES['rho_symmetry'] and max(last['symmetry']['u'],last['symmetry']['v']) <= GATES['velocity_symmetry'], nonuniform=last['nonuniform'], audit=budget['status'] == 'passed' and budget['max_relative_residual'] <= GATES['budget_residual']) row['checks'] = {k:bool(v) for k,v in checks.items()} row['status'] = 'passed' if all(checks.values()) else 'failed' save() print('NEGATIVE RESULT',json.dumps(row['checks']),flush=True) def main(): 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') == '0' assert not (OUT/'results.json').exists(), 'Do not repeat this run' baseline = json.loads((OUT/'baseline.json').read_text()) assert all(sha(ROOT/name) == digest for name,digest in baseline['hashes'].items()) assert sha(SOURCE) == SOURCE_SHA gpu = subprocess.check_output(['nvidia-smi','-i','0','--query-gpu=memory.used,utilization.gpu','--format=csv,noheader,nounits'],text=True) used, utilization = map(int,gpu.strip().split(',')) apps = subprocess.check_output(['nvidia-smi','-i','0','--query-compute-apps=pid','--format=csv,noheader'],text=True) assert not apps.strip() and used < 100 and utilization == 0, (gpu,apps) report = dict(python=sys.executable, package=CelerisLab.__file__, cuda_source=compiler.kernel_path('kernel_v2.cu'), script_sha256=sha(__file__), gpu='0', gpu_before=gpu, gates=dict(GATES,top_density_error=5e-7), restart=dict(status='not_tested'),zero=dict(status='not_tested'),negative=dict(status='not_tested')) def save(): report['peak_rss_mib'] = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024 (OUT/'results.json').write_text(json.dumps(report,indent=2)+'\n') print('SOURCES',report['python'],report['package'],report['cuda_source'],flush=True) save() try: restart(report,save) zero(report,save) negative(report,save) except Exception as error: report['error'] = repr(error) for name in ('restart','zero','negative'): if report[name]['status'] == 'running': report[name]['status'] = 'failed' budget = report['negative'].get('population_budget', {}) if budget.get('status') == 'running': budget['status'] = 'failed' traceback.print_exc() finally: report['old_files_unchanged'] = all(sha(ROOT/name) == digest for name,digest in baseline['hashes'].items()) save() print('FINAL',json.dumps({name:report[name]['status'] for name in ('restart','zero','negative')}),flush=True) return 0 if report['old_files_unchanged'] and all(report[n]['status']=='passed' for n in ('restart','zero','negative')) else 1 if __name__ == '__main__': raise SystemExit(main())