165 lines
9.3 KiB
Python
165 lines
9.3 KiB
Python
"""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())
|