Files
CelerisLab/tests/validation/multigrid/run_spacetime_gpu_n64.py

393 lines
24 KiB
Python

"""Opt-in, frozen N64 space-time GPU trajectory; no convergence certification."""
import argparse
import hashlib
import json
import os
from pathlib import Path
import resource
import subprocess
import sys
import time
import numpy as np
ROOT = Path(__file__).resolve().parents[3]
sys.path.insert(0, str(ROOT/'src'))
sys.path.insert(0, str(Path(__file__).parents[1]/'taylor_green'))
from CelerisLab.config import LBMConfig
import oracle as tg
import spacetime_gpu as gpu
import spacetime_table_guard as table_guard
ref = gpu.ref
CPU_SHA = '180eef90a1c3f43113fe18cc54981ee4723cfc0cc2724334907ff188022fa770'
UNIFORM_SHA = '01e7415ae39456f2de7f0c3d331afb971c49d19353b43038650e32e80737f1e7'
CUDA_SHA = 'afea11ea4a1c2c3800b2f6a7a8e81f5568989fc9e1cc02fb52d7e1ea3864e003'
AUDIT_SHA = 'fba340b1533bec308b666ee84666df0ac9d08e2872be31ee379e2d1447565158'
STEPS = (256, 512, 768, 1024)
TIMES = (.12, .24, .36, .48)
FACTORS = dict(velocity_l2=.8, pressure_l2=1., energy_relative=1.)
LIMITS = dict(velocity_l2=.05, pressure_l2=.15, energy_relative=.05, mass_relative=1e-5)
HEALTH = ('boundary_pre_nonfinite', 'boundary_pre_nonpositive', 'boundary_stored_nonfinite',
'boundary_stored_nonpositive', 'physical_population_nonfinite', 'physical_population_nonpositive')
def sha(path):
return hashlib.sha256(Path(path).read_bytes()).hexdigest()
def load_references(cpu_path, uniform_path):
for path, expected in ((cpu_path, CPU_SHA), (uniform_path, UNIFORM_SHA)):
if sha(path) != expected:
raise ValueError(f'pinned reference receipt differs: {path}')
cpu, uniform = [json.loads(Path(p).read_text()) for p in (cpu_path, uniform_path)]
assert cpu['completed'] and cpu['experimental_pass'] and not cpu['source']['dirty']
assert cpu['source']['commit'] == '07879bc52bfef8d95503146c55d3823c8b6f8d0a'
assert [r['mode'] for r in cpu['runs']] == ['constant_ivp', 'spacetime']
assert all(tuple(row['tstar'] for row in r['rows']) == TIMES for r in cpu['runs'])
assert uniform['completed'] and not uniform['source']['dirty']
assert uniform['source']['revision'] == '9805383f9c50d867e8555471baf5dbda8333a652'
assert (uniform['n'], uniform['depth'], uniform['viscosity']) == (64, 1, .0192)
assert uniform['limits'] == tg.LIMITS == LIMITS
return cpu, uniform
def total(state, maps):
arrays = [a.astype(np.float64) for a in (state.coarse_f, state.fine_f)]
return ref.cpu.MOMENTS.T @ ref.directional(ref.cpu.State(*arrays), maps)
def differing_fields(actual, expected):
"""Compare serialized sample numerics exactly, including the sign of zero."""
if actual.keys() != expected.keys():
return ['field_set']
return [key for key in actual if np.asarray(actual[key]).dtype != np.asarray(expected[key]).dtype
or np.asarray(actual[key]).shape != np.asarray(expected[key]).shape
or np.asarray(actual[key]).tobytes() != np.asarray(expected[key]).tobytes()]
def sample_difference(actual, expected):
"""Descriptive comparison only; never an equivalence gate for different codecs."""
return {key:dict(reference=value,signed=(np.asarray(actual[key])-value).tolist(),
absolute=np.abs(np.asarray(actual[key])-value).tolist()) for key,value in expected.items()}
class LedgerTotals:
"""Compensated directional sums; retain absolute step sums to expose cancellation."""
names = ('net_change','physical_collision','coarse_collision','fine0_collision','fine1_collision',
'child_rounding','child_interpolation','child_cast','register_rounding','receiver_rounding','unexplained')
def __init__(self):
self.sums = np.zeros((len(self.names),9)); self.correction = self.sums.copy()
self.absolute = self.sums.copy(); self.absolute_moments = np.zeros((len(self.names),3))
self.maxima = self.sums.copy(); self.max_steps = np.zeros(self.sums.shape,np.int64)
self.cells = np.zeros((3,4)); self.cell_correction = self.cells.copy()
self.steps = 0; self.bound = 0.; self.maximum = 0.
def add(self, ledger):
values = np.array([ledger[name] for name in self.names],np.float64)
cells = np.asarray(ledger['collision_cell_mass'],np.float64)
if (values.shape != self.sums.shape or not np.isfinite(values).all()
or cells.shape != (3,4) or not np.isfinite(cells).all()
or not np.isfinite(ledger['unexplained_bound']) or ledger['unexplained_bound'] <= 0):
raise AssertionError('invalid diagnostic ledger')
delta = values-self.correction
updated = self.sums+delta
self.correction = (updated-self.sums)-delta
self.sums = updated
self.absolute += abs(values)
self.absolute_moments += abs(values @ ref.cpu.MOMENTS)
larger = abs(values) > self.maxima
self.max_steps[larger] = self.steps+1; self.maxima = np.maximum(self.maxima,abs(values))
delta = cells-self.cell_correction
updated = self.cells+delta
self.cell_correction = (updated-self.cells)-delta; self.cells = updated
self.steps += 1; self.bound += ledger['unexplained_bound']
self.maximum = max(self.maximum,float(abs(values[-1]).max()))
def snapshot(self, endpoint_change=None):
values = dict(zip(self.names,self.sums))
accounted = sum(values[key] for key in ('physical_collision','child_rounding','register_rounding','receiver_rounding'))
closure = None if endpoint_change is None else endpoint_change-accounted
telescoping = None if endpoint_change is None else endpoint_change-values['net_change']
closure_pass = None if closure is None else max(float(abs(closure).max()),float(abs(telescoping).max())) <= self.bound
collision_split = values['physical_collision']-sum(values[key] for key in ('coarse_collision','fine0_collision','fine1_collision'))
child_split = values['child_rounding']-values['child_interpolation']-values['child_cast']
split_pass = max(float(abs(collision_split).max()),float(abs(child_split).max())) <= self.bound
return dict(intervals=self.steps, summed_per_step_bound=self.bound, max_step_unexplained=self.maximum,
closure_pass=closure_pass,component_split_pass=split_pass,
endpoint_status='not_measured' if endpoint_change is None else 'measured',
endpoint_directional_change=None if endpoint_change is None else endpoint_change.tolist(),
closure_residual=None if closure is None else closure.tolist(),
telescoping_residual=None if telescoping is None else telescoping.tolist(),
closure_mass_momentum=None if closure is None else (ref.cpu.MOMENTS.T@closure).tolist(),
collision_split_residual=collision_split.tolist(),child_split_residual=child_split.tolist(),
collision_cell_mass=dict(levels=['coarse','fine0','fine1'],columns=['net','positive','negative','l1'],
applied_volume_weights=[1.,.25,.25],values=self.cells.tolist()),
components={key:dict(directional=self.sums[i].tolist(),mass_momentum=(self.sums[i]@ref.cpu.MOMENTS).tolist(),
absolute_step_directional_sum=self.absolute[i].tolist(),
absolute_step_mass_momentum_sum=self.absolute_moments[i].tolist(),
max_absolute_interval=self.maxima[i].tolist(),max_absolute_interval_step=self.max_steps[i].tolist())
for i,key in enumerate(self.names)})
def warm_and_reset(run, initial, check_cap, *, audit=False):
"""Time 32 synchronized intervals, then restore exact uploaded arrival F."""
run.initialize(initial)
before = run.download()
started = time.monotonic()
for _ in range(32):
check_cap()
run.step(audit=audit)
seconds = time.monotonic()-started
warm_health = run.health_counts.tolist()
run.initialize(initial)
after = run.download()
for a, b in zip((before.coarse_f, before.fine_f), (after.coarse_f, after.fine_f)):
np.testing.assert_array_equal(a, b)
if run.step_count != 0 or after.tick != 0 or run.health_counts.any():
raise AssertionError('warm-up reset did not restore tick/counters')
return after, dict(intervals=32, seconds=seconds, estimated_1024_seconds=32*seconds,
health_counts=warm_health, reset_verified=True, ledger_audit=audit)
def assess(rows, cpu, uniform, health_counts):
if tuple(row['tstar'] for row in rows) != TIMES:
raise ValueError('four frozen sample times required')
worst = {key: max(row[key] for row in rows) for key in LIMITS}
campaign = {key: all(np.isfinite(row[key]) and 0 <= row[key] <= limit for row in rows)
for key, limit in LIMITS.items()}
campaign['population_health'] = not any(health_counts)
campaign['sampled_positive_finite'] = all(
np.isfinite(row[key]) and row[key] > 0 for row in rows
for key in ('minimum_population', 'minimum_boundary_pre', 'minimum_boundary_stored'))
baseline, candidate = cpu['runs']
thresholds = {key: factor*baseline['worst'][key] for key, factor in FACTORS.items()}
experimental = {key: all(np.isfinite(row[key]) and row[key] >= 0 for row in rows)
and worst[key] <= threshold for key, threshold in thresholds.items()}
differences = []
for actual, expected in zip(rows, candidate['rows']):
delta = {}
for key, value in expected.items():
d = np.asarray(actual[key])-np.asarray(value)
delta[key] = dict(cpu=value, signed=d.tolist(), absolute=np.abs(d).tolist())
differences.append(dict(tstar=actual['tstar'], fields=delta))
historical_worst = max(row['pressure_l2'] for row in uniform['rows'] if row['tstar'] > 0)
historical_limit = max(2*historical_worst, .02)
assert historical_limit == .024506347612699403
return dict(worst=worst, campaign_checks=campaign, campaign_pass=all(campaign.values()),
experimental_thresholds=thresholds, experimental_checks=experimental,
experimental_pass=all(experimental.values()),
gpu_minus_cpu=dict(status='descriptive_only_no_equivalence_tolerance', rows=differences),
historical_pressure_gate=dict(status='historical_shifted_FP32_uniform_vs_unshifted_FP32_candidate',
matched_storage=False, uniform_worst=historical_worst, threshold=historical_limit,
candidate_worst=worst['pressure_l2'], passed=worst['pressure_l2'] <= historical_limit))
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument('--gpu-admitted', action='store_true')
parser.add_argument('--cpu-receipt', type=Path, required=True)
parser.add_argument('--uniform-receipt', type=Path, required=True)
parser.add_argument('--audit-against', type=Path, help='audit every interval and compare exact saved 1779cf2 samples')
parser.add_argument('--shifted', action='store_true', help='opt-in shifted storage with all interval ledgers')
parser.add_argument('--unshifted-receipt', type=Path, help='pinned1779 comparison, descriptive only')
parser.add_argument('--output', type=Path, required=True)
args = parser.parse_args()
if not args.gpu_admitted or not os.environ.get('CUDA_VISIBLE_DEVICES'):
parser.error('explicit --gpu-admitted and CUDA_VISIBLE_DEVICES required')
if args.output.exists():
parser.error('output exists; preserve evidence and choose a new path')
if args.shifted != (args.unshifted_receipt is not None) or (args.shifted and args.audit_against):
parser.error('--shifted requires --unshifted-receipt and cannot use exact --audit-against')
exact_repeat = args.audit_against is not None
cpu, uniform = load_references(args.cpu_receipt, args.uniform_receipt)
previous = None
comparison_path = args.audit_against or args.unshifted_receipt
if comparison_path is not None:
if sha(comparison_path) != AUDIT_SHA:
raise ValueError('pinned GPU trajectory receipt differs')
previous = json.loads(comparison_path.read_text())
if previous['source']['numpy'] != np.__version__:
raise RuntimeError('audited repeat requires unchanged NumPy environment')
gpu.check_sources()
if sha(Path(gpu.__file__).with_suffix('.cu')) != CUDA_SHA:
raise RuntimeError('reviewed CUDA source changed')
source = dict(commit=subprocess.check_output(['git', 'rev-parse', 'HEAD'], cwd=ROOT, text=True).strip(),
dirty=bool(subprocess.check_output(['git', 'status', '--porcelain'], cwd=ROOT, text=True)),
hashes={p.name: sha(p) for p in (Path(__file__), Path(gpu.__file__), Path(gpu.__file__).with_suffix('.cu'),
Path(__file__).with_name('test_spacetime_gpu.py'), Path(table_guard.__file__))},
reference_hashes=gpu.PINS, numpy=np.__version__, checkout=str(ROOT))
if source['dirty']:
raise RuntimeError('clean committed checkout required')
started = time.monotonic()
report = dict(scope='N64 experimental space-time GPU; no grid-convergence qualification',
completed=False, passed=False, source=source, rows=[], health_counts=[0]*6, health_labels=HEALTH,
setup=dict(n=64, bounds=[8,8,56,48], re_L=100, U=.03, viscosity=.0192, tau=[.5576,.6152],
phase='arrival_f', storage='shifted_FP32' if args.shifted else 'unshifted_FP32', collision='SRT', interval_count=1024,
initialization='common analytic FP64 field, encoded using production FP64 weights before FP32 upload' if args.shifted
else 'same physical FP64 equilibrium as CPU, rounded to FP32 on upload',
cpu_storage='FP64', sample_steps=STEPS, sample_times=TIMES),
references=dict(cpu_receipt_sha256=CPU_SHA, cpu_source=cpu['source'],
cpu_archive_drive_id='197F9nwDOGf_YROK1a46YpVvNrSldnH_w',
cpu_archive_sha256='850ddce3f5203ef38e76e5153e718a9b92e137d1fc388ba7cd12f5d47eab83b4',
uniform_receipt_sha256=UNIFORM_SHA, uniform_source=uniform['source'],
uniform_storage='historical_shifted_FP32'),
frozen_screen=dict(campaign_limits=LIMITS, experimental_factors=FACTORS,
gpu_cpu_equivalence_tolerance=None, historical_pressure_threshold=.024506347612699403),
boundary_budget_status='pending_measurement' if previous else 'not_measured', max_directional_residual=None,
minima_scope='initial and four sampled physical arrival states; boundary minima at four samples',
limits=dict(seconds=180., peak_rss_mib=512.), cuda_visible_devices=os.environ['CUDA_VISIBLE_DEVICES'])
run = context = None
ledger_totals = LedgerTotals() if previous else None
if previous:
report['diagnostic'] = dict(completed=False, passed=False, reference_sha256=AUDIT_SHA,
reference_source=previous['source']['commit'], sample_checks=[], prefixes=[],
scope='net drift attribution; exact saved diagnostics only, not full-field trajectory certification' if exact_repeat
else 'raw stored-value ledger; unshifted and CPU comparisons descriptive, initial codecs differ')
phase = 'setup'
def save():
report['resources'] = dict(seconds=time.monotonic()-started,
peak_rss_mib=resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024,
device_allocations_bytes=None if run is None else run.device_bytes,
device_memory_scope='owner allocations; excludes compiler, driver and context')
if run is not None:
report['health_counts'] = run.health_counts.tolist()
report['completed_intervals'] = run.step_count
if ledger_totals is not None:
report['audited_intervals'] = ledger_totals.steps
if ledger_totals.steps:
report['boundary_budget_status'] = 'measured' if ledger_totals.steps==1024 else 'measured_prefix'
report['diagnostic']['accumulated'] = ledger_totals.snapshot()
ref.write_receipt(args.output, report)
def check_cap():
if time.monotonic()-started > 180 or resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024 > 512:
raise RuntimeError('frozen resource cap exceeded')
save()
try:
from CelerisLab.cuda.context import CudaContext
context = CudaContext(0)
config = LBMConfig(nx=64, ny=64, viscosity=.0192, velocity=0., ddf_shifting=args.shifted)
run = gpu.SpaceTimeGPU(config, (8,8,56,48), context=context)
report.update(device=context.name, sm_arch=context.sm_arch, build_metadata=run.build_metadata,
table_sha256=gpu.table_hash(run.tables))
report['table_validation'] = table_guard.validate_tables(run.plan,run.tables)
if previous:
for key in ('device','sm_arch','table_sha256') + (('build_metadata',) if exact_repeat else ()):
if report[key] != previous[key]:
raise AssertionError(f'audited repeat identity differs: {key}')
analytic = ref.cpu.State(tg.initial((64,64),0,64,.0192), tg.initial((128,128),1,64,.0192))
initial = analytic
phase = 'warm_reset'
initial, report['warmup'] = warm_and_reset(run, initial, check_cap, audit=previous is not None)
maps = run.plan.maps
masks = (maps.coarse_active, maps.fine_active)
arrays = (initial.coarse_f, initial.fine_f)
mass0 = tg.metrics(arrays, masks, 0, 64, .0192)['mass']
moment0 = total(initial, maps)
report['initial'] = dict(mass=mass0, momentum=moment0[1:].tolist(),
minimum_population=min(float(a[:,m].min()) for a,m in zip(arrays,masks)))
if exact_repeat and differing_fields(report['initial'],previous['initial']):
raise AssertionError('audited repeat initial state diagnostics differ')
raw_initial = run._download_stored() if args.shifted else initial
direction0 = ref.directional(ref.cpu.State(*(a.astype(np.float64) for a in
(raw_initial.coarse_f,raw_initial.fine_f))),maps)
if args.shifted:
errors=[]; unshifted_errors=[]
for raw,actual,wanted,mask in zip((raw_initial.coarse_f,raw_initial.fine_f),arrays,
(analytic.coarse_f,analytic.fine_f),masks):
encoded=gpu._encode_physical(wanted,True)
np.testing.assert_array_equal(raw[:,mask].view(np.uint32),encoded[:,mask].view(np.uint32))
errors.append(float(abs(actual[:,mask]-wanted[:,mask]).max()))
unshifted_errors.append(float(abs(gpu._encode_physical(wanted,False)[:,mask]-wanted[:,mask]).max()))
report['warmup']['raw_reset_verified']=True
report['initial_quantization']=dict(maximum_population_error_by_level=errors,
unshifted_codec_error_by_level=unshifted_errors,
moment_change_vs_analytic=(moment0-total(analytic,maps)).tolist(),
cpu_start='pinned analytic FP64 trajectory, not reinitialized to decoded GPU upload',
matched_cpu_initial_bits=False)
report['unshifted_comparison']=dict(status='descriptive_only_codecs_have_different_initial_bits',
reference_sha256=AUDIT_SHA,initial=sample_difference(report['initial'],previous['initial']),rows=[])
save()
check_cap()
phase = 'trajectory'
trajectory_started = time.monotonic()
for step in range(1,1025):
check_cap()
ledger = run.step(audit=previous is not None)
if ledger_totals is not None:
ledger_totals.add(ledger)
report['max_directional_residual'] = ledger_totals.maximum
if step in STEPS:
state = run.download()
if state.tick != 2*step:
raise AssertionError('sample is not synchronized arrival state')
arrays = (state.coarse_f, state.fine_f)
row = tg.metrics(arrays, masks, step, 64, .0192, mass0)
row.update(step=step, momentum_change=(total(state,maps)-moment0)[1:].tolist(),
mass_change=row['mass']-mass0,
minimum_population=min(float(a[:,m].min()) for a,m in zip(arrays,masks)))
if ledger_totals is not None:
row['minimum_boundary_pre'] = ledger['minimum_pre']
row['minimum_boundary_stored'] = ledger['minimum_stored']
else:
with context.activate():
row['minimum_boundary_pre'] = float(run.read(run.pre,(len(run.plan.groups),4),np.float64).min())
row['minimum_boundary_stored'] = float(run.read(run.h,(len(run.plan.groups),4)).min())
if not all(np.isfinite(value).all() for value in row.values()):
report['failed_sample'] = {key: repr(value) for key,value in row.items()}
raise FloatingPointError('nonfinite sample diagnostic')
report['rows'].append(row)
if previous:
expected=previous['rows'][len(report['rows'])-1]
changed = differing_fields(row,expected) if exact_repeat else []
report['diagnostic']['sample_checks'].append(dict(step=step,
status='exact_repeat' if exact_repeat else 'descriptive_only',differing_fields=changed if exact_repeat else None))
if args.shifted:
report['unshifted_comparison']['rows'].append(dict(step=step,fields=sample_difference(row,expected)))
raw_state=run._download_stored() if args.shifted else state
endpoint = ref.directional(ref.cpu.State(*(a.astype(np.float64) for a in
(raw_state.coarse_f,raw_state.fine_f))),maps)-direction0
prefix = ledger_totals.snapshot(endpoint)
if args.shifted:
prefix['physical_minus_raw_endpoint_moments'] = (
total(state,maps)-moment0-ref.cpu.MOMENTS.T@endpoint).tolist()
report['diagnostic']['prefixes'].append(prefix)
if not prefix['closure_pass'] or not prefix['component_split_pass']:
raise AssertionError('cumulative ledger does not close within summed existing per-step bounds')
if changed:
raise AssertionError(f'audited repeat sample differs: {changed}')
save()
print(json.dumps(row, allow_nan=False), flush=True)
report['trajectory_seconds'] = time.monotonic()-trajectory_started
check_cap()
report.update(assess(report['rows'], cpu, uniform, run.health_counts))
if args.shifted:
report['historical_pressure_gate'].update(status='historical_shifted_FP32_uniform_vs_shifted_FP32_candidate',
matched_storage=True,matched_control=False)
report['completed'] = True
report['passed'] = report['campaign_pass'] and report['experimental_pass']
if previous:
report['diagnostic'].update(completed=True,passed=ledger_totals.steps==1024 and len(report['rows'])==4)
except Exception as exc:
report['failure'] = dict(phase=phase, error=repr(exc))
raise
finally:
# Preserve the primary failure and every completed sample even if cleanup fails.
for owner in (run, context):
if owner is not None:
try:
owner.close()
except Exception as exc:
report.setdefault('cleanup_errors', []).append(repr(exc))
report['passed'] = False
save()
return 0 if report['passed'] else 1
if __name__ == '__main__':
raise SystemExit(main())