From 707afe6bb21cf6d4bd993f76a3f42f3ab9dc7c8b Mon Sep 17 00:00:00 2001 From: Codex Date: Sun, 18 Oct 2026 12:35:05 -0700 Subject: [PATCH] test: prepare matched shifted GPU spacing ladder --- .../multigrid/run_spacetime_gpu_ladder.py | 263 ++++++++++++++++++ tests/validation/multigrid/spacetime_gpu.py | 131 ++++++++- .../multigrid/spacetime_table_guard.py | 11 +- .../multigrid/test_spacetime_gpu_ladder.py | 205 ++++++++++++++ 4 files changed, 601 insertions(+), 9 deletions(-) create mode 100644 tests/validation/multigrid/run_spacetime_gpu_ladder.py create mode 100644 tests/validation/multigrid/test_spacetime_gpu_ladder.py diff --git a/tests/validation/multigrid/run_spacetime_gpu_ladder.py b/tests/validation/multigrid/run_spacetime_gpu_ladder.py new file mode 100644 index 0000000..84ec9da --- /dev/null +++ b/tests/validation/multigrid/run_spacetime_gpu_ladder.py @@ -0,0 +1,263 @@ +"""Four opt-in shifted GPU controls: U64/U128/U256 and refined base128/fine256.""" +import argparse +import hashlib +import inspect +import json +import os +from pathlib import Path +import resource +import subprocess +import time + +import numpy as np +import run_spacetime_gpu_n64 as old + +gpu,ref,tg,ROOT = old.gpu,old.ref,old.tg,old.ROOT +BASE='26d0d3853c9d4d3cd6ae3f27653dced56f8c5fec' +REFERENCE_SHA='9d154c989b087bf288f86d413db0f3067b73767b618c329d4b9e53c4eee00032' +CASES=(('uniform64',64,False),('uniform128',128,False),('uniform256',256,False),('refined128',128,True)) +LIMITS=old.LIMITS +SECONDS,RSS_MIB=600.,512. + + +def setup(n, refined): + return dict(base_n=n,finest_n=2*n if refined else n,normalized_domain=[0.,1.,0.,1.], + re_L=100,U=.03,mach=.03*np.sqrt(3),viscosity=.0003*n, + tau=[.5+.0009*n]+([.5+.0018*n] if refined else []), + bounds=[n//8,n//8,7*n//8,3*n//4] if refined else None, + normalized_patch=[.125,.125,.875,.75] if refined else None, + intervals=16*n,sample_steps=[4*n,8*n,12*n,16*n],sample_times=list(old.TIMES), + early_steps=list(range(0,4*n+1,n//16)),early_times=(np.arange(65)*.001875).tolist(), + phase='arrival_f',storage='shifted_FP32',collision='SRT', + initialization='common analytic FP64 equilibrium at actual cell centers, production shifted codec') + + +def pressure_projection(arrays,masks,t,n,viscosity): + """One global gauge and a weighted Gram fit, including unequal leaf areas.""" + rows=[] + for level,(a,mask) in enumerate(zip(arrays,masks)): + rho=np.asarray(a[:,mask],np.float64).sum(axis=0) + x,y=tg.centers(mask.shape,level); x,y=x[mask],y[mask] + exact=tg.exact(x,y,t,n,viscosity)[2] + angle=4*np.pi/n + basis=np.column_stack((np.cos(angle*x),np.sin(angle*x),np.cos(angle*y),np.sin(angle*y))) + rows.append(((rho-1)/3,exact,np.full(len(rho),.25**level),basis)) + p,exact,w,b=[np.concatenate(x) for x in zip(*rows)] + volume=w.sum(); p-=w@p/volume; exact-=w@exact/volume; e=p-exact + b-=w@b/volume + gram=b.T@(w[:,None]*b); rhs=b.T@(w*e) + coefficients=np.linalg.solve(gram,rhs); residual=e-b@coefficients + norm=float(w@(e*e)); exact_norm=float(w@(exact*exact)); remaining=float(w@(residual*residual)) + projected=float(coefficients@gram@coefficients) + return dict(step=t,tstar=t*.03/n,absolute_rms=np.sqrt(norm/volume),relative_rms=np.sqrt(norm/exact_norm), + coefficients=coefficients.tolist(),basis=['cos_4pi_x','sin_4pi_x','cos_4pi_y','sin_4pi_y'], + residual_rms=np.sqrt(remaining/volume),residual_relative=np.sqrt(remaining/exact_norm), + gram_condition=float(np.linalg.cond(gram)), + orthogonality=(b.T@(w*residual)/volume).tolist(),energy_closure=(norm-projected-remaining)/volume) + + +def load_reference(path): + if old.sha(path)!=REFERENCE_SHA: + raise ValueError('verified shifted N64 receipt differs') + data=json.loads(Path(path).read_text()) + if (not data['completed'] or not data['passed'] or data['source']['commit']!=BASE + or data['setup']['storage']!='shifted_FP32' or data['setup']['bounds']!=[8,8,56,48] + or data['setup']['interval_count']!=1024 or len(data['rows'])!=4): + raise ValueError('N64 reference identity differs') + if tuple(row['tstar'] for row in data['rows'])!=old.TIMES: + raise ValueError('N64 reference sample times differ') + return data + + +def source_identity(reference): + git=lambda *args:subprocess.check_output(['git',*args],cwd=ROOT,text=True).strip() + if git('status','--porcelain'): + raise RuntimeError('clean committed checkout required') + unchanged=['src','tests/validation/taylor_green/oracle.py']+[f'tests/validation/multigrid/{name}' for name in + ('spacetime_gpu.cu','spacetime_reference.py','cpu_reference.py','run_spacetime_gpu_n64.py','test_spacetime_gpu.py')] + subprocess.check_call(['git','diff','--exit-code',BASE,'--',*unchanged],cwd=ROOT) + previous=subprocess.check_output(['git','show',f'{BASE}:tests/validation/multigrid/spacetime_gpu.py'],cwd=ROOT,text=True) + # The raw interval schedule remains byte-for-byte unchanged despite the new owner. + method=inspect.getsource(gpu.SpaceTimeGPU.step) + if method not in previous or np.__version__!=reference['source']['numpy']: + raise RuntimeError('raw schedule or NumPy identity changed') + files=[Path(__file__),Path(gpu.__file__),Path(old.table_guard.__file__),Path(__file__).with_name('test_spacetime_gpu_ladder.py')] + return dict(commit=git('rev-parse','HEAD'),dirty=False,numpy=np.__version__,production_tree=git('rev-parse','HEAD:src'), + hashes={p.name:old.sha(p) for p in files},cuda_sha256=old.sha(Path(gpu.__file__).with_suffix('.cu')), + raw_schedule_sha256=hashlib.sha256(method.encode()).hexdigest(),unchanged_against=BASE, + reference_hashes=gpu.PINS) + + +def directional(arrays,masks): + return sum(.25**level*np.asarray(a[:,m],np.float64).sum(axis=1) for level,(a,m) in enumerate(zip(arrays,masks))) + + +def populations(run,refined,*,raw=False): + state=run._download_stored() if raw else run.download() + if not refined: + return (state,) + if state.phase!='arrival_f' or state.tick!=2*run.step_count: + raise AssertionError('refined sample is not synchronized arrival F') + return state.coarse_f,state.fine_f + + +def assess(rows,health): + checks={key:all(np.isfinite(row[key]) and 0<=row[key]<=limit for row in rows) for key,limit in LIMITS.items()} + checks['population_health']=not any(health) + checks['sampled_population_positive']=all(row['minimum_population']>0 for row in rows) + return dict(worst={key:max(row[key] for row in rows) for key in LIMITS},campaign_checks=checks, + campaign_pass=all(checks.values()),accuracy_comparison_status='descriptive_only_no_new_threshold') + + +def ratios(a,b,*,trend=False): + out=dict(ratio=a/b if b>0 else None,ratio_status='defined' if b>0 else 'zero_denominator') + if trend: + out.update(two_point_slope=np.log2(a/b) if a>0 and b>0 else None, + slope_status='defined' if a>0 and b>0 else 'zero_error_log_undefined') + return out + + +def comparisons(runs,reference): + data={r['name']:r for r in runs}; data['refined64']=reference + out=[] + for refined,base,fine in (('refined64','uniform64','uniform128'),('refined128','uniform128','uniform256')): + for control in (base,fine): + out.append(dict(candidate=refined,control=control,control_spacing='base' if control==base else 'finest', + samples=[dict(tstar=a['tstar'],**{key:dict(candidate=a[key],control=b[key],**ratios(a[key],b[key])) + for key in ('velocity_l2','pressure_l2','energy_relative')}) + for a,b in zip(data[refined]['rows'],data[control]['rows'])])) + trends=[] + for coarse,fine in (('uniform64','uniform128'),('uniform128','uniform256'),('refined64','refined128')): + trends.append(dict(coarse=coarse,fine=fine,samples=[dict(tstar=a['tstar'],**{ + key:ratios(a[key],b[key],trend=True) + for key in ('velocity_l2','pressure_l2','energy_relative')}) + for a,b in zip(data[coarse]['rows'],data[fine]['rows'])])) + return dict(controls=out,trends=trends,caveat='norm differences are not interface error; fixed Mach and two refined sizes do not certify convergence') + + +def main(): + p=argparse.ArgumentParser(description=__doc__) + p.add_argument('--gpu-admitted',action='store_true'); p.add_argument('--reference-n64',type=Path,required=True) + p.add_argument('--output',type=Path,required=True); args=p.parse_args() + if not args.gpu_admitted or not os.environ.get('CUDA_VISIBLE_DEVICES'): + p.error('explicit GPU admission and CUDA_VISIBLE_DEVICES required') + if args.output.exists() or args.output.with_suffix(args.output.suffix+'.tmp').exists(): p.error('output or partial receipt exists; preserve it') + reference=load_reference(args.reference_n64); source=source_identity(reference); gpu.check_sources() + started=time.monotonic() + report=dict(completed=False,passed=False,source=source,runs=[],limits=dict(seconds=SECONDS,peak_rss_mib=RSS_MIB), + campaign_limits=LIMITS,reference=dict(sha256=REFERENCE_SHA,source=reference['source'], + coverage='verified refined base64/fine128 four-time rows only; early pressure trace unavailable'), + historical_pressure_gate=dict(scope='reused_refined64_only',**reference['historical_pressure_gate']), + n64_experimental_screen=dict(scope='reused_refined64_only',checks=reference['experimental_checks']), + frozen_cases=[setup(n,r) for _,n,r in CASES],early_trace_scope='65 predeclared rows per new run, no field snapshots', + cuda_visible_devices=os.environ['CUDA_VISIBLE_DEVICES']) + run=context=ledger=None + case_started=None + def save(): + report['resources']=dict(seconds=time.monotonic()-started,peak_rss_mib=resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024, + rss_scope='Python process peak; external wrapper monitors process group and wall time') + report['resources']['within_limits']=report['resources']['seconds']<=SECONDS and report['resources']['peak_rss_mib']<=RSS_MIB + if report['completed'] and not report['resources']['within_limits']: report['passed']=False + if run is not None: + report['runs'][-1].update(completed_intervals=run.step_count,health_counts=run.health_counts.tolist(), + seconds=time.monotonic()-case_started) + if ledger is not None: + report['runs'][-1].update(ledger=ledger.snapshot(),audited_intervals=ledger.steps) + temporary=args.output.with_suffix(args.output.suffix+'.tmp') + ref.write_receipt(temporary,report); os.replace(temporary,args.output) + def cap(): + if time.monotonic()-started>SECONDS or resource.getrusage(resource.RUSAGE_SELF).ru_maxrss/1024>RSS_MIB: + raise RuntimeError('frozen aggregate resource cap exceeded') + save() + try: + from CelerisLab.cuda.context import CudaContext + context=CudaContext(0) + if context.name!=reference['device'] or context.sm_arch!=reference['sm_arch']: + raise RuntimeError('matched reference device identity differs') + report.update(device=context.name,sm_arch=context.sm_arch) + for name,n,refined in CASES: + cap(); ledger=None; case_started=time.monotonic(); cfg=old.LBMConfig(nx=n,ny=n,viscosity=.0003*n,velocity=0.,ddf_shifting=True) + case=dict(name=name,setup=setup(n,refined),completed=False,phase='setup',rows=[],early_pressure=[],prefixes=[],health_counts=[0]*6, + ledger_scope='all intervals; fine/transfer components are mathematically absent for uniform runs') + report['runs'].append(case); save() + run=gpu.SpaceTimeGPU(cfg,tuple(case['setup']['bounds']),context=context,allow_n128=True) if refined else gpu.UniformGPU(cfg,context=context) + case.update(build_metadata=run.build_metadata,device_allocations_bytes=run.device_bytes, + device_memory_scope='owner allocations; excludes driver/context/compiler') + expected_modules=None + if not refined and n in (64,128): expected_modules=[reference['build_metadata'][0 if n==64 else 1]] + if refined: expected_modules=[report['runs'][1]['build_metadata'][0],report['runs'][2]['build_metadata'][0]] + if expected_modules is not None: + if run.build_metadata!=expected_modules: raise RuntimeError('matched source/config/PTX metadata differs') + case['matched_build_metadata']=True + if any(m['flags']!=reference['build_metadata'][0]['flags'] for m in run.build_metadata): + raise RuntimeError('compiler flags differ from reference') + if refined: + case['tables']=old.table_guard.validate_tables(run.plan,run.tables,allow_n128=True) + masks=(run.plan.maps.coarse_active,run.plan.maps.fine_active) + else: masks=(np.ones((n,n),bool),) + analytic=tuple(tg.initial(m.shape,l,n,cfg.viscosity) for l,m in enumerate(masks)) + initial=ref.cpu.State(*analytic) if refined else analytic[0] + physical=lambda:populations(run,refined) + stored=lambda:populations(run,refined,raw=True) + run.initialize(initial); raw0=stored(); direction0=directional(raw0,masks) + case['phase']='warm_reset'; warm_started=time.monotonic() + for _ in range(32): cap(); run.step() + case['warmup_seconds']=time.monotonic()-warm_started + run.initialize(initial) + for a,b in zip(raw0,stored()): np.testing.assert_array_equal(a.view(np.uint32),b.view(np.uint32)) + if run.step_count or run.health_counts.any(): raise AssertionError('warm reset failed') + case['raw_reset_verified']=True; arrays=physical() + first=tg.metrics(arrays,masks,0,n,cfg.viscosity); mass0=first['mass']; moment0=ref.cpu.MOMENTS.T@directional(arrays,masks) + case['initial']=dict(mass=mass0,momentum=moment0[1:].tolist(),minimum_population=min(float(a[:,m].min()) for a,m in zip(arrays,masks))) + case['initial_quantization']=[float(abs(a[:,m]-b[:,m]).max()) for a,b,m in zip(arrays,analytic,masks)] + case['early_pressure'].append(pressure_projection(arrays,masks,0,n,cfg.viscosity)) + ledger=old.LedgerTotals(); case['phase']='trajectory'; trajectory_started=time.monotonic() + case['all_interval_mass']=dict(maximum_relative=0.,step=0,signed_change_at_maximum=0.,limit=LIMITS['mass_relative'],passed=True) + for step in range(1,16*n+1): + cap(); measured=run.step(); ledger.add(measured) + if run.step_count!=step: raise AssertionError('interval count differs') + mass_change=float(ledger.sums[0].sum()); drift=abs(mass_change)/mass0 + if drift>case['all_interval_mass']['maximum_relative']: + case['all_interval_mass'].update(maximum_relative=drift,step=step,signed_change_at_maximum=mass_change,passed=drift<=LIMITS['mass_relative']) + if step in case['setup']['early_steps'] or step in case['setup']['sample_steps']: + arrays=physical() + if step in case['setup']['early_steps']: + case['early_pressure'].append(pressure_projection(arrays,masks,step,n,cfg.viscosity)); save() + if step in case['setup']['sample_steps']: + row=tg.metrics(arrays,masks,step,n,cfg.viscosity,mass0) + row.update(step=step,momentum_change=(ref.cpu.MOMENTS.T@directional(arrays,masks)-moment0)[1:].tolist(), + mass_change=row['mass']-mass0,minimum_population=min(float(a[:,m].min()) for a,m in zip(arrays,masks))) + if not all(np.isfinite(v).all() for v in row.values()): raise FloatingPointError('nonfinite sampled metric') + case['rows'].append(row); prefix=ledger.snapshot(directional(stored(),masks)-direction0) + prefix['physical_minus_raw_endpoint_moments']=(ref.cpu.MOMENTS.T@directional(arrays,masks)-moment0-ref.cpu.MOMENTS.T@np.array(prefix['endpoint_directional_change'])).tolist() + case['prefixes'].append(prefix) + if not prefix['closure_pass'] or not prefix['component_split_pass']: raise AssertionError('cumulative budget failed') + save(); print(json.dumps(dict(case=name,**row),allow_nan=False),flush=True) + case.update(assess(case['rows'],run.health_counts)) + cap(); case.update(completed=True,phase='completed',health_counts=run.health_counts.tolist(),audited_intervals=ledger.steps, + ledger=ledger.snapshot(),trajectory_seconds=time.monotonic()-trajectory_started,seconds=time.monotonic()-case_started) + case['passed']=case['campaign_pass'] and case['all_interval_mass']['passed'] and ledger.steps==16*n and len(case['early_pressure'])==65 + run.close(); run=None; save() + cap(); report.update(completed=True,passed=all(c['passed'] for c in report['runs']),comparisons=comparisons(report['runs'],reference)) + except Exception as exc: + report['failure']=repr(exc) + raise + finally: + if run is not None and ledger is not None: + current=report['runs'][-1] + current['failure_endpoint_status']='unavailable_owner_not_ready' + if run.status=='ready': + try: + current['failure_endpoint']=ledger.snapshot(directional(populations(run,refined,raw=True),masks)-direction0) + current['failure_endpoint_status']='measured' + except Exception as exc: current['failure_endpoint_status']=repr(exc) + 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()) diff --git a/tests/validation/multigrid/spacetime_gpu.py b/tests/validation/multigrid/spacetime_gpu.py index 1e8269a..9a5addc 100644 --- a/tests/validation/multigrid/spacetime_gpu.py +++ b/tests/validation/multigrid/spacetime_gpu.py @@ -19,15 +19,19 @@ def check_sources(): raise RuntimeError(f"reference source changed: {name}") -def validate_config(config): +def validate_config(config, *, uniform=False, allow_n128=False): + if type(uniform) is not bool or type(allow_n128) is not bool: + raise ValueError("admission flags must be bool") cfg = replace(config) cfg.validate() - if (cfg.nx != cfg.ny or not 8 <= cfg.nx <= 64 or cfg.nz != 1 + if (cfg.nx != cfg.ny or not 8 <= cfg.nx <= (256 if uniform else (128 if allow_n128 else 64)) or cfg.nz != 1 or cfg.collision != "SRT" or not cfg.is_d2q9 or cfg.streaming != "double_buffer" or cfg.store_precision != "FP32" or cfg.data_type != "FP32" or type(cfg.ddf_shifting) is not bool or cfg.les_enabled or cfg.velocity != 0 or cfg.wall_slot_width or cfg.wall_normal_velocity): - raise ValueError("experiment requires periodic unforced square N8..64, FP32 D2Q9/SRT") + raise ValueError("experiment requires periodic unforced square FP32 D2Q9/SRT within admitted size") + if uniform: + return (cfg,) fine = replace(cfg, nx=2*cfg.nx, ny=2*cfg.ny, viscosity=2*cfg.viscosity) fine.validate() return cfg, fine @@ -136,9 +140,9 @@ def rounding_ledger(plan, tables, before, after, g0, predicted, fine_phases, pre class SpaceTimeGPU: """Caller-owned context, one stream, explicit synchronized arrival state.""" - def __init__(self, config, bounds, *, context, mode='spacetime', identity=False): + def __init__(self, config, bounds, *, context, mode='spacetime', identity=False, allow_n128=False): check_sources() - self.configs = validate_config(config) + self.configs = validate_config(config,allow_n128=allow_n128) if mode not in ('spacetime', 'frozen') or type(identity) is not bool: raise ValueError("unsupported interpolation or debug collision mode") self.bounds = tuple(bounds) @@ -349,3 +353,120 @@ class SpaceTimeGPU: except Exception: if kind is None: raise + + +def uniform_ledger(before, after, collided, *, shifted): + """Periodic streaming permutes each direction; only collision changes its sum.""" + a,b,g = [np.asarray(x,np.float64) for x in (before,after,collided)] + collision=(g-a).sum(axis=(1,2)); net=b.sum(axis=(1,2))-a.sum(axis=(1,2)) + residual=net-collision + bound=float(64*np.finfo(float).eps*(a.sum()+int(shifted)*a.shape[1]*a.shape[2])) + np.testing.assert_allclose(residual,0,rtol=0,atol=bound) + cells=(g-a).sum(axis=0) + out={key:[0.]*9 for key in ('fine0_collision','fine1_collision','child_rounding', + 'child_interpolation','child_cast','register_rounding','receiver_rounding')} + out.update(net_change=net.tolist(),physical_collision=collision.tolist(),coarse_collision=collision.tolist(), + unexplained=residual.tolist(),unexplained_bound=bound, + collision_cell_mass=[[float(cells.sum()),float(cells[cells>0].sum()), + float(cells[cells<0].sum()),float(abs(cells).sum())],[0.]*4,[0.]*4], + boundary_budget_status='not_applicable_uniform',inactive_components='fine, child, register, receiver') + return out + + +class UniformGPU(SpaceTimeGPU): + """One explicit arrival-F level; reuse the reviewed kernels and owner helpers.""" + def __init__(self, config, *, context, identity=False): + check_sources() + self.configs=validate_config(config,uniform=True) + if type(identity) is not bool: + raise ValueError('identity must be bool') + self.ctx,self.identity=context,identity + self.allocations,self.modules,self.build_metadata=[],[],[] + self.status,self.step_count,self.device_bytes='new',0,0 + self.health_counts=np.zeros(6,np.uint64) + import pycuda.driver as cuda + from CelerisLab.cuda import compiler_v2 as compiler + if Path(compiler._KERNEL_DIR).resolve()!=Path(__file__).resolve().parents[3]/'src/CelerisLab/lbm/kernels': + raise RuntimeError('wrong production kernel checkout') + self.cuda=cuda + with context.activate(): + self.stream=cuda.Stream() + try: + source=Path(__file__).with_suffix('.cu') + build=compiler.generate_config(self.configs[0]) + try: + Path(build,source.name).write_bytes(source.read_bytes()) + ptx=compiler.compile_kernel(context.sm_arch,source=source.name,build_dir=build) + self.modules.append(compiler.load_module(ptx)) + import json + metadata=json.loads(Path(ptx+'.json').read_text()) + metadata['test_entry_sha256']=hashlib.sha256(source.read_bytes()).hexdigest() + self.build_metadata.append(metadata) + finally: + compiler.cleanup_build_dir(build) + n=config.nx + self.buffers=[self.allocate(9*n*n*4) for _ in range(3)] + self.active=self.upload(np.ones((n,n),np.uint8)) + self.counts=self.allocate(6*4) + except BaseException: + try: + self.close() + except BaseException: + pass + raise + + def initialize(self, physical): + cfg=self.configs[0] + if self.status not in ('new','ready') or np.shape(physical)!=(9,cfg.ny,cfg.nx): + raise ValueError('ready owner and exact uniform population shape required') + stored=_encode_physical(physical,cfg.ddf_shifting) + decoded=_decode_stored(stored,cfg.ddf_shifting) + if not np.isfinite(decoded).all() or (decoded<=0).any(): + raise ValueError('positive finite physical populations required') + try: + with self.ctx.activate(): + self.cuda.memcpy_htod(self.buffers[0],stored) + self.health_counts[:]=0; self.step_count=0; self.status='ready' + except BaseException: + self.status='failed' + raise + + def _download_stored(self): + if self.status!='ready': + raise RuntimeError('uniform owner is not ready') + try: + with self.ctx.activate(): + n=self.configs[0].nx + return self.read(self.buffers[0],(9,n,n)) + except BaseException: + self.status='failed' + raise + + def download(self): + return _decode_stored(self._download_stored(),self.configs[0].ddf_shifting) + + def step(self, *, audit=True): + if self.status!='ready' or type(audit) is not bool: + raise RuntimeError('ready owner and bool audit required') + cfg=self.configs[0]; shape=(9,cfg.ny,cfg.nx) + before=self._download_stored() if audit else None + with self.ctx.activate(): + try: + self.status='advancing' + self.cuda.memset_d32_async(self.counts,0,6,self.stream) + f,g,out=self.buffers + self.launch(0,'Collide',f,g,self.active,np.float32(0 if self.identity else cfg.omega),self.counts) + self.launch(0,'Pull',g,out,self.active,self.counts) + counts=self.read(self.counts,(6,),np.uint32) + self.health_counts+=counts.astype(np.uint64) + if counts.any(): + raise AssertionError(f'uniform population health failed: {counts.tolist()}') + ledger=uniform_ledger(before,self.read(out,shape),self.read(g,shape),shifted=cfg.ddf_shifting) if audit else dict( + boundary_budget_status='not_measured',unexplained=None) + ledger['health_counts']=counts.tolist() + self.buffers[0],self.buffers[2]=out,f + self.step_count+=1; self.status='ready' + return ledger + except BaseException: + self.status='failed' + raise diff --git a/tests/validation/multigrid/spacetime_table_guard.py b/tests/validation/multigrid/spacetime_table_guard.py index 5334df2..c6f3820 100644 --- a/tests/validation/multigrid/spacetime_table_guard.py +++ b/tests/validation/multigrid/spacetime_table_guard.py @@ -1,4 +1,4 @@ -"""N64 admission check: exact topology and independently rational interpolation. +"""N64/N128 admission check: exact topology and independently rational interpolation. Raw FP64 bytes are diagnostic: different LAPACK builds need not produce them identically. The 2e-14 ceiling is the existing interpolation constraint policy, @@ -11,6 +11,8 @@ import numpy as np from spacetime_gpu import table_hash STRUCTURE_SHA = 'cd1842952c4dd4df82da8015454c25da9be12b983270236cc65b4044fa38e6e4' +STRUCTURES = {64: (STRUCTURE_SHA, 532), + 128: ('ee4ecb645674ae04eb87220a5217851405ccc0f45fbb0fe210499dcea5cd93b1', 1060)} ATOL = 2e-14 @@ -57,13 +59,14 @@ def exact_weights(kind, points): return np.array(_multiply(e,solution[:4]), dtype=np.float64) -def validate_tables(plan, tables): +def validate_tables(plan, tables, *, allow_n128=False): """Reject topology, layout, padding, operator, or constraint changes separately.""" structure = table_hash({key:value for key,value in tables.items() if key != 'weights'}) - if plan.n != 64 or structure != STRUCTURE_SHA: + admitted = (64,128) if allow_n128 is True else (64,) + if plan.n not in admitted or structure != STRUCTURES[plan.n][0]: raise ValueError('frozen integer topology changed') w = tables['weights'] - if (w.dtype != np.dtype(np.float64) or w.shape != (532,4,10) or len(plan.groups) != w.shape[0] + if (w.dtype != np.dtype(np.float64) or w.shape != (STRUCTURES[plan.n][1],4,10) or len(plan.groups) != w.shape[0] or not all(a.flags.c_contiguous for a in tables.values()) or not np.isfinite(w).all()): raise ValueError('weight dtype/shape/layout/finiteness changed') errors = dict(coefficient=0., constant=0., mean=0., temporal=0.) diff --git a/tests/validation/multigrid/test_spacetime_gpu_ladder.py b/tests/validation/multigrid/test_spacetime_gpu_ladder.py new file mode 100644 index 0000000..1cf2c8a --- /dev/null +++ b/tests/validation/multigrid/test_spacetime_gpu_ladder.py @@ -0,0 +1,205 @@ +"""Host checks or separately admitted tiny GPU ladder capability; no trajectories.""" +import argparse +from contextlib import nullcontext +import json +import os +from pathlib import Path +import subprocess +import tempfile +import time +from types import SimpleNamespace +import unittest +from unittest.mock import Mock + +import numpy as np +import run_spacetime_gpu_ladder as ladder + +gpu,ref,tg=ladder.gpu,ladder.ref,ladder.tg +RECORDS=[] +ADMITTED=False + + +def config(n): + return ladder.old.LBMConfig(nx=n,ny=n,viscosity=.0003*n,velocity=0.,ddf_shifting=True) + + +class Host(unittest.TestCase): + def test_scaling_and_admission(self): + for _,n,refined in ladder.CASES: + s=ladder.setup(n,refined) + self.assertEqual(s['finest_n'],n*(2 if refined else 1)) + self.assertEqual(len(s['early_steps']),65) + np.testing.assert_allclose(np.array(s['early_steps'])*.03/n,s['early_times'],rtol=0,atol=3e-17) + self.assertEqual(s['intervals']*.03/n,.48) + gpu.validate_config(config(n),uniform=not refined,allow_n128=refined) + with self.assertRaises(ValueError): gpu.validate_config(config(128)) + with self.assertRaises(ValueError): gpu.validate_config(config(256),allow_n128=True) + with self.assertRaises(ValueError): gpu.validate_config(config(512),uniform=True) + + def test_n128_tables_and_mutations(self): + plan=ref.Plan(128,(16,16,112,96)); tables=gpu.pack_plan(plan) + out=ladder.old.table_guard.validate_tables(plan,tables,allow_n128=True) + self.assertEqual(out['structure_sha256'],'ee4ecb645674ae04eb87220a5217851405ccc0f45fbb0fe210499dcea5cd93b1') + self.assertEqual(len(plan.groups),1060) + with self.assertRaises(ValueError): ladder.old.table_guard.validate_tables(plan,tables) + for kind in ('topology','balanced_weights'): + bad={k:v.copy() for k,v in tables.items()} + if kind=='topology': bad['sources'][0,0]+=1 + else: bad['weights'][0,:2,:2]+=1e-5*np.array([[1,-1],[-1,1]]) + with self.assertRaises(ValueError): ladder.old.table_guard.validate_tables(plan,bad,allow_n128=True) + + def test_projection_manufactured_nonuniform_leaves(self): + n=16; plan=ref.Plan(n,(3,3,13,12)); masks=(plan.maps.coarse_active,plan.maps.fine_active) + coords=[tg.centers(m.shape,l) for l,m in enumerate(masks)] + x=np.concatenate([v[0][m] for v,m in zip(coords,masks)]) + y=np.concatenate([v[1][m] for v,m in zip(coords,masks)]) + w=np.concatenate([np.full(m.sum(),.25**l) for l,m in enumerate(masks)]) + b=np.column_stack((np.cos(4*np.pi*x/n),np.sin(4*np.pi*x/n),np.cos(4*np.pi*y/n),np.sin(4*np.pi*y/n))) + b-=w@b/w.sum(); gram=b.T@(w[:,None]*b) + residual=np.cos(6*np.pi*x/n)+np.sin(6*np.pi*y/n); residual-=w@residual/w.sum() + residual-=b@np.linalg.solve(gram,b.T@(w*residual)); residual*=1e-5 + coefficients=np.array([2,-3,4,-1])*1e-5 + pressure=tg.exact(x,y,0,n,.0003*n)[2]+b@coefficients+residual+.017 + arrays=[]; offset=0 + for m in masks: + a=np.full((9,*m.shape),np.nan); count=m.sum() + a[:,m]=ref.cpu.W[:,None]*(1+3*pressure[offset:offset+count]); arrays.append(a); offset+=count + result=ladder.pressure_projection(arrays,masks,0,n,.0003*n) + np.testing.assert_allclose(result['coefficients'],coefficients,rtol=0,atol=3e-17) + self.assertAlmostEqual(result['residual_rms'],np.sqrt(w@(residual*residual)/w.sum()),places=15) + np.testing.assert_allclose(result['orthogonality'],0,rtol=0,atol=1e-18) + self.assertLess(abs(result['energy_closure']),1e-20) + expected=tg.metrics(arrays,masks,0,n,.0003*n)['pressure_l2'] + self.assertAlmostEqual(result['relative_rms'],expected,places=12) + + def test_uniform_budget_and_negative_control(self): + rng=np.random.default_rng(12); a=rng.normal(0,.001,(9,12,12)).astype(np.float32) + g=(a*.99).astype(np.float32) + b=np.stack([np.roll(g[q],(int(ref.cpu.CY[q]),int(ref.cpu.CX[q])),axis=(0,1)) for q in range(9)]) + ledger=gpu.uniform_ledger(a,b,g,shifted=True) + self.assertLess(max(abs(np.array(ledger['unexplained']))),ledger['unexplained_bound']) + b[1,0,0]+=1e-4 + with self.assertRaises(AssertionError): gpu.uniform_ledger(a,b,g,shifted=True) + + def test_uniform_failed_io_poisoning(self): + run=object.__new__(gpu.UniformGPU); run.configs=(config(12),); run.status='ready' + run.ctx=SimpleNamespace(activate=nullcontext); run.buffers=['f']; run.cuda=SimpleNamespace(memcpy_htod=Mock(side_effect=RuntimeError('copy'))) + with self.assertRaisesRegex(RuntimeError,'copy'): run.initialize(np.broadcast_to(ref.cpu.W[:,None,None],(9,12,12))) + self.assertEqual(run.status,'failed'); run.status='ready'; run.read=Mock(side_effect=RuntimeError('read')) + with self.assertRaisesRegex(RuntimeError,'read'): run._download_stored() + self.assertEqual(run.status,'failed') + + def test_uniform_health_mask_reaches_both_kernels(self): + run=object.__new__(gpu.UniformGPU); run.configs=(config(16),); run.status='ready'; run.identity=False + run.ctx=SimpleNamespace(activate=nullcontext); run.buffers=['f','g','out']; run.active=object(); run.counts='counts' + run.stream=object(); run.cuda=SimpleNamespace(memset_d32_async=Mock()); run.launch=Mock() + run.read=Mock(return_value=np.zeros(6,np.uint32)); run.health_counts=np.zeros(6,np.uint64); run.step_count=0 + run.step(audit=False) + self.assertEqual([call.args[1] for call in run.launch.call_args_list],['Collide','Pull']) + self.assertTrue(all(call.args[4] is run.active for call in run.launch.call_args_list)) + + def test_zero_error_comparisons(self): + rows=[dict(tstar=t,velocity_l2=0.,pressure_l2=0.,energy_relative=0.) for t in ladder.old.TIMES] + result=ladder.comparisons([dict(name=name,rows=rows) for name,_,_ in ladder.CASES],dict(rows=rows)) + self.assertIsNone(result['trends'][0]['samples'][0]['pressure_l2']['ratio']) + self.assertNotIn('two_point_slope',result['controls'][0]['samples'][0]['pressure_l2']) + self.assertIn('two_point_slope',result['trends'][0]['samples'][0]['pressure_l2']) + json.dumps(result,allow_nan=False) + + def test_frozen_screen_and_json_receipt(self): + rows=[dict(tstar=t,minimum_population=.02,**{key:np.float64(limit*.5) for key,limit in ladder.LIMITS.items()}) for t in ladder.old.TIMES] + self.assertTrue(ladder.assess(rows,[0]*6)['campaign_pass']) + rows[0]['mass_relative']=2e-5 + failed=ladder.assess(rows,[0]*6); self.assertFalse(failed['campaign_pass']) + with tempfile.TemporaryDirectory() as tmp: + path=Path(tmp)/'receipt.json'; ref.write_receipt(path,dict(completed=True,passed=np.bool_(False),assessment=failed)) + saved=json.loads(path.read_text()); self.assertFalse(saved['passed']); self.assertTrue(saved['completed']) + + +class Device(unittest.TestCase): + @classmethod + def setUpClass(cls): + if not ADMITTED: raise unittest.SkipTest('explicit tiny GPU admission required') + from CelerisLab.cuda.context import CudaContext + cls.context=CudaContext(0) + + @classmethod + def tearDownClass(cls): cls.context.close() + + def test_uniform_actual_upload_oracle_and_reset(self): + for n in (16,64,128,256): + with gpu.UniformGPU(config(n),context=self.context) as run: + analytic=tg.initial((n,n),0,n,.0003*n); run.initialize(analytic); start=run._download_stored() + actual_start=run.download(); mask=np.ones((n,n),bool) + spec=ref.Spec(n_coarse=n,patch_lo=1,halo=2,tau_coarse=.5+.0009*n,tau_fine=.5+.0018*n,collision='srt') + collided=ref.cpu._collide(spec,actual_start,np.isfinite(actual_start),spec.tau_coarse,mask,np.zeros((n,n),int)) + expected=np.stack([np.roll(collided[q],(int(ref.cpu.CY[q]),int(ref.cpu.CX[q])),axis=(0,1)) for q in range(9)]) + ledger=run.step(); actual=run.download() + np.testing.assert_allclose(actual,expected,rtol=0,atol=3e-6) + audited=run._download_stored(); run.initialize(analytic); run.step(audit=False) + np.testing.assert_array_equal(audited.view(np.uint32),run._download_stored().view(np.uint32)) + run.initialize(analytic); np.testing.assert_array_equal(start.view(np.uint32),run._download_stored().view(np.uint32)) + self.assertEqual(run.step_count,0); self.assertFalse(run.health_counts.any()) + RECORDS.append(dict(case=f'uniform{n}_one_interval',error=float(abs(actual-expected).max()),ledger=ledger, + build_metadata=run.build_metadata,device_bytes=run.device_bytes,raw_reset_verified=True)) + + def test_uniform_tagged_phase_and_health_controls(self): + n=16; rng=np.random.default_rng(n) + analytic=np.broadcast_to(ref.cpu.W[:,None,None],(9,n,n))*(1+1e-3*rng.standard_normal((9,n,n))) + for identity in (True,False): + with gpu.UniformGPU(config(n),context=self.context,identity=identity) as run: + run.initialize(analytic); start=run.download(); raw_start=run._download_stored(); mask=np.ones((n,n),bool) + spec=ref.Spec(tau_coarse=.5+.0009*n,collision='identity' if identity else 'srt') + collided=ref.cpu._collide(spec,start,np.isfinite(start),spec.tau_coarse,mask,np.zeros((n,n),int)) + expected=np.stack([np.roll(collided[q],(int(ref.cpu.CY[q]),int(ref.cpu.CX[q])),axis=(0,1)) for q in range(9)]) + ledger=run.step(); actual=run.download(); np.testing.assert_allclose(actual,expected,rtol=0,atol=3e-6) + if identity: + rolled=np.stack([np.roll(raw_start[q],(int(ref.cpu.CY[q]),int(ref.cpu.CX[q])),axis=(0,1)) for q in range(9)]) + np.testing.assert_array_equal(run._download_stored().view(np.uint32),rolled.view(np.uint32)) + if not identity: + wrong=ref.cpu._collide(spec,expected,np.isfinite(expected),spec.tau_coarse,mask,np.zeros((n,n),int)) + self.assertGreater(float(abs(actual-wrong).max()),3e-6) + RECORDS.append(dict(case='uniform_tagged',identity=identity,error=float(abs(actual-expected).max()),ledger=ledger)) + with gpu.UniformGPU(config(n),context=self.context) as run: + run.initialize(analytic) + with self.context.activate(): + bad=np.broadcast_to(-gpu._WEIGHTS-.01,(9,n,n)).astype(np.float32).copy() + run.cuda.memcpy_htod(run.buffers[0],bad) + with self.assertRaises(AssertionError): run.step() + self.assertTrue(run.health_counts.any()); self.assertEqual(run.status,'failed') + RECORDS.append(dict(case='uniform_health_negative_control',health_counts=run.health_counts.tolist())) + + def test_refined_n128_actual_upload_oracle(self): + with gpu.SpaceTimeGPU(config(128),(16,16,112,96),context=self.context,allow_n128=True) as run: + tables=ladder.old.table_guard.validate_tables(run.plan,run.tables,allow_n128=True) + run.initialize(ref.cpu.State(tg.initial((128,128),0,128,.0384),tg.initial((256,256),1,128,.0384))) + spec=ref.Spec(n_coarse=128,patch_lo=16,halo=32,tau_coarse=.6152,tau_fine=.7304,collision='srt') + expected,_=ref.advance(spec,run.plan,run.download()); ledger=run.step(); actual=run.download(); errors=[] + for a,b in zip((actual.coarse_f,actual.fine_f),(expected.coarse_f,expected.fine_f)): + np.testing.assert_allclose(a,b,rtol=0,atol=3e-6,equal_nan=True); errors.append(float(np.nanmax(abs(a-b)))) + RECORDS.append(dict(case='refined128_one_interval',errors=errors,ledger=ledger,tables=tables, + build_metadata=run.build_metadata,device_bytes=run.device_bytes)) + + +def main(): + global ADMITTED + p=argparse.ArgumentParser(description=__doc__); modes=p.add_mutually_exclusive_group(required=True) + modes.add_argument('--host-only',action='store_true'); modes.add_argument('--gpu-admitted',action='store_true') + p.add_argument('--output',type=Path,required=True); args=p.parse_args(); ADMITTED=args.gpu_admitted + if args.output.exists(): p.error('output exists; preserve evidence') + if ADMITTED and not os.environ.get('CUDA_VISIBLE_DEVICES'): p.error('explicit CUDA_VISIBLE_DEVICES required') + paths=(Path(__file__),Path(ladder.__file__),Path(gpu.__file__),Path(gpu.__file__).with_suffix('.cu'),Path(ladder.old.table_guard.__file__)) + report=dict(completed=False,passed=False,mode='tiny_gpu' if ADMITTED else 'host_only',records=RECORDS, + source=dict(commit=subprocess.check_output(['git','rev-parse','HEAD'],cwd=ladder.ROOT,text=True).strip(), + dirty=bool(subprocess.check_output(['git','status','--porcelain'],cwd=ladder.ROOT,text=True)), + hashes={p.name:ladder.old.sha(p) for p in paths})) + ref.write_receipt(args.output,report); started=time.monotonic() + suite=unittest.TestLoader().loadTestsFromTestCase(Host) + if ADMITTED: suite.addTests(unittest.TestLoader().loadTestsFromTestCase(Device)) + result=unittest.TextTestRunner(verbosity=2).run(suite) + report.update(completed=True,passed=result.wasSuccessful(),tests=result.testsRun,seconds=time.monotonic()-started, + failures=[dict(test=str(t),error=e) for t,e in result.errors+result.failures]) + ref.write_receipt(args.output,report); return 0 if report['passed'] else 1 + + +if __name__=='__main__': raise SystemExit(main())