343 lines
19 KiB
Python
343 lines
19 KiB
Python
"""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())
|