diff --git a/tests/validation/run_j3_coverage.py b/tests/validation/run_j3_coverage.py new file mode 100644 index 0000000..d66c299 --- /dev/null +++ b/tests/validation/run_j3_coverage.py @@ -0,0 +1,342 @@ +"""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())