diff --git a/docs/projected-explosion-experiment.md b/docs/projected-explosion-experiment.md index 823c76d..7ac2ced 100644 --- a/docs/projected-explosion-experiment.md +++ b/docs/projected-explosion-experiment.md @@ -27,8 +27,36 @@ Use the same candidate checkout for paired Re100 N64 depth2 runs, differing only `python tests/validation/taylor_green/run_case.py --repo "$PWD" --re 100 --n 64 --depth 2 --explosion projected_linear --output ` -Select the physical CUDA device explicitly. Never overwrite historical reports. Compare frozen TG error gates and full trajectories, including mass, energy, pressure and velocity. Amplitude/nonmodal shape decomposition can inform mechanism but cannot substitute for the gates. First establish CUDA compilation and CPU/device parity, then paired end-to-end budgets. No accuracy improvement, production suitability or GPU pass is claimed by these CPU results. +Select the physical CUDA device explicitly. Never overwrite historical reports. Compare frozen TG error gates and full trajectories, including mass, energy, pressure and velocity. Amplitude/nonmodal shape decomposition can inform mechanism but cannot substitute for the gates. First establish CUDA compilation and CPU/device parity, then paired end-to-end budgets. The CPU diagnostic below establishes improvement only in its FP64 reference trajectories; it does not establish production suitability or a candidate GPU pass. ## Local verification limits -Focused candidate/refinement/block-layout/build-contract checks pass (62 unittest tests). Full unittest discovery is not green: the attempted 200-test run reports 18 missing-dependency/import errors and three historical source-pin refusals (jet-circle and LES guards), plus one skipped test. Those pin refusals correctly reject the changed `blocks.cu`; their immutable fixtures were not weakened or rewritten. Python syntax checks and `git diff --check` pass. No CUDA build, device execution, GPU performance or TG improvement has been established. +Focused candidate/refinement/block-layout/build-contract checks pass (62 unittest tests). Full unittest discovery is not green: the attempted 200-test run reports 18 missing-dependency/import errors and three historical source-pin refusals (jet-circle and LES guards), plus one skipped test. Those pin refusals correctly reject the changed `blocks.cu`; their immutable fixtures were not weakened or rewritten. Python syntax checks and `git diff --check` pass. No candidate CUDA build, device execution, GPU performance or GPU TG improvement has been established. + + +## Donor-support audit + +The current support is exactly the original two-coarse-cell-wide explosion ring. Centered normal derivatives therefore vanish along straight faces: the inner donor lacks an uncovered inward neighbor, and the outer donor lacks a third-ring outward neighbor. This is an explicit limitation, not evidence of a missing donor layer. Independent inspection of the authors' [LinearExplosion implementation](https://github.com/lssfau/walberla/blob/master/src/lbm/refinement/LinearExplosion.h) found the same two-cell directional slabs and both-neighbors-valid condition; their union equals all 368 distinct donors in the N64 TG geometry. No upstream code was copied. This does not claim identical block-local tangential fallback at every fine/fine block seam. + +A separate read-only dependency trace compared the ring to all uncovered periodic parent neighbors for four rectangles. Wider support changed some unused outer ghost populations but none consumed by final fine ACTIVE outputs or coarse restriction through the two hops, including full collision dependencies. Widening the support is not part of this candidate. + +## Reproducible CPU Taylor–Green A/B + +Run from this checkout, with NumPy and Python installed: + +`python tests/validation/taylor_green/run_cpu_explosion.py --repo "$PWD" --n 32 --output-dir /tmp/explosion-ab` + +Repeat with `--n 64` for the second retained case. These are the only supported grids. The script reuses unchanged CPUREF `advance_coarse_once`, existing TG initialization/metrics, actual singleton-rectangle transfer tables, and the tested projected reconstruction. The rectangular adapter verifies its coalescence children against the production table and includes all eight grazing-corner links. Only the explosion operation differs within each pair. + +The reference stores physical FP64 arrival populations; production stores shifted FP32 postcollision populations. Collision preserves the macroscopic moments compared here, and initialization is equilibrium. This supports a macro-trajectory comparison, not population-by-population/device equivalence. Outputs record actual checkout/harness provenance and hashes of both imported diagnostic helpers. The compact retained `tests/validation/taylor_green/cpu_explosion_summary.json` records original N32/N64 results and historical GPU report hashes, without importing those reports into this branch. + +At Re100, depth 2, t*=0.48: + +| Root N | Velocity L2, constant → projected | Pressure L2 | Orthogonal velocity residual | +|---|---|---|---| +| 32 | 3.24942% → 1.52264% | 8.98696% → 4.85601% | 2.28814% → 0.98615% | +| 64 | 1.13113% → 0.49304% | 3.55050% → 1.65801% | 0.93699% → 0.40289% | + +Across sampled times, constant CPU versus historical constant GPU differs by at most 6.93e-8 in velocity L2 and 3.58e-6 in pressure L2 for N32; 1.74e-8 and 8.64e-7 for N64. Relative mass drift remains below 1.42e-13 and absolute momentum change below 2.11e-13 in these CPU pairs. Independent N32 rerunning reproduced trajectory/error metrics exactly; recomputed quadrature-factor metadata agrees to roundoff. + +Amplitude is the volume-weighted projection onto the exact velocity field; the orthogonal residual isolates the remaining shape error. The energy decomposition includes the discrete composite-grid exact-energy quadrature factor recorded per row, not simply amplitude squared plus residual squared. These are bounded diagnostic improvements, not proof of complete method fidelity, uniform-grid accuracy equivalence or a validated CUDA implementation. Frozen production validation gates are unchanged. diff --git a/tests/validation/taylor_green/cpu_explosion_summary.json b/tests/validation/taylor_green/cpu_explosion_summary.json new file mode 100644 index 0000000..0bc327a --- /dev/null +++ b/tests/validation/taylor_green/cpu_explosion_summary.json @@ -0,0 +1,101 @@ +{ + "schema": 1, + "implementation_revision": "43c527b", + "historical_baseline_revision": "9805383f9c50d867e8555471baf5dbda8333a652", + "scope": "Re100 depth2; FP64 CPU arrival_f macro diagnostics, not candidate CUDA execution or performance", + "tstar": 0.48, + "cases": [ + { + "n": 32, + "historical_gpu_report": "CelerisLab_9805383_tg_matched_re100_n32_d2.json", + "historical_gpu_sha256": "6e84837781f2c4eb001721d6f828ea0d927b255fd4c46efdef443a9ae59af730", + "constant_max_abs_cpu_gpu_metric_difference": { + "velocity_l2": 6.928351422436085e-08, + "pressure_l2": 3.5717157327030957e-06, + "energy_relative": 2.1743219684455806e-07 + }, + "runs": [ + { + "explosion": "constant", + "final": { + "velocity_l2": 0.03249421441793122, + "pressure_l2": 0.08986961517973119, + "energy_relative": 0.044970292392050686, + "amplitude_bias": -0.02307193176461997, + "orthogonal_residual_l2": 0.022881432107435998, + "mass_signed_relative": -7.616129948928574e-14, + "momentum_change": [ + -3.26405569239796e-14, + -2.7000103541841014e-14 + ] + }, + "max_abs_mass_relative": 7.616129948928574e-14, + "max_abs_momentum_change": 3.608224830031759e-14 + }, + { + "explosion": "projected_linear", + "final": { + "velocity_l2": 0.015226371368103809, + "pressure_l2": 0.048560123195715434, + "energy_relative": 0.02285059695090361, + "amplitude_bias": -0.011601431258836059, + "orthogonal_residual_l2": 0.009861499773660909, + "mass_signed_relative": -7.616129948928574e-14, + "momentum_change": [ + -4.030109579389318e-14, + -2.3444787777826548e-15 + ] + }, + "max_abs_mass_relative": 7.616129948928574e-14, + "max_abs_momentum_change": 4.030109579389318e-14 + } + ] + }, + { + "n": 64, + "historical_gpu_report": "CelerisLab_9805383_tg_re100_n64_d2.json", + "historical_gpu_sha256": "251f743d6adc33dff1d829db91fcc5f982c9ac299c177eb2c94e313249b25b9f", + "constant_max_abs_cpu_gpu_metric_difference": { + "velocity_l2": 1.7330029265794078e-08, + "pressure_l2": 8.634212284933107e-07, + "energy_relative": 8.74280938097105e-08 + }, + "runs": [ + { + "explosion": "constant", + "final": { + "velocity_l2": 0.01131134574353476, + "pressure_l2": 0.03550495463085949, + "energy_relative": 0.012514764591008386, + "amplitude_bias": -0.006336459299048802, + "orthogonal_residual_l2": 0.009369942693596348, + "mass_signed_relative": -1.412203687323199e-13, + "momentum_change": [ + -2.1049828546892968e-13, + -1.632548263241773e-14 + ] + }, + "max_abs_mass_relative": 1.412203687323199e-13, + "max_abs_momentum_change": 2.1049828546892968e-13 + }, + { + "explosion": "projected_linear", + "final": { + "velocity_l2": 0.004930396085095945, + "pressure_l2": 0.016580125701009935, + "energy_relative": 0.0056291365148441195, + "amplitude_bias": -0.0028419317091766816, + "orthogonal_residual_l2": 0.004028924138812507, + "mass_signed_relative": -1.4099832412739488e-13, + "momentum_change": [ + -1.709743457922741e-13, + 7.568945470382005e-14 + ] + }, + "max_abs_mass_relative": 1.4099832412739488e-13, + "max_abs_momentum_change": 1.709743457922741e-13 + } + ] + } + ] +} diff --git a/tests/validation/taylor_green/run_cpu_explosion.py b/tests/validation/taylor_green/run_cpu_explosion.py new file mode 100644 index 0000000..e8969b7 --- /dev/null +++ b/tests/validation/taylor_green/run_cpu_explosion.py @@ -0,0 +1,85 @@ +"""Diagnostic-only rectangular adapter to CPUREF-0.3. No production changes.""" +import sys,json,time,argparse,hashlib +from pathlib import Path +import numpy as np +parser=argparse.ArgumentParser(description=__doc__) +parser.add_argument('--repo',type=Path,required=True,help='Candidate checkout at 43c527b') +parser.add_argument('--n',type=int,choices=(32,64),default=32) +parser.add_argument('--output-dir',type=Path,required=True) +args=parser.parse_args() +ROOT=args.repo.resolve() +sys.path[:0]=[str(ROOT/'src'),str(ROOT/'tests/validation/multigrid'),str(ROOT/'tests/validation/taylor_green'),str(ROOT/'tests/unit')] +import cpu_reference as cpu +import oracle +from run_case import provenance +from test_projected_explosion import reconstruct +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.refinement import FixedPatchSet2D +from CelerisLab.lbm.refinement_tables import projected_explosion_tables +from CelerisLab.lbm.blocks import BlockLayout2D +class RectSpec(cpu.Spec): + @property + def n_fine(self): return 2*self.n_coarse + +def setup(n): + nu=oracle.U*n/100 + cfg=LBMConfig(nx=n,ny=n,viscosity=nu) + bounds=oracle.rectangles(n,2)[0] + p=FixedPatchSet2D.from_config(cfg,bounds) + spec=RectSpec(n_coarse=n,patch_lo=n//8,patch_hi=7*n//8,halo=n//4,tau_coarse=p.tau_coarse,tau_fine=p.tau_fine,collision='srt') + q,k=np.divmod(p.coalesce_coarse.astype(int),n*n); y,x=np.divmod(k,n) + rows=np.column_stack((q,y,x)); links=[] + for qi,yi,xi in rows: + sx,sy=xi-cpu.CX[qi],yi-cpu.CY[qi] + if not p.covered[sy,sx]: links.append(cpu.CornerLink(f'{xi}_{yi}_q{qi}',int(qi),(int(sx),int(sy)),(int(xi),int(yi)))) + assert len(links)==8 + maps=cpu.Maps(~p.covered,p.covered,p.fine_kind==0,p.fine_kind==1,p.fine_kind==2,p.fine_kind,np.arange(2*n)//2,(np.arange(2*n)+.5)/2,rows,tuple(links)) + # Verify the CPU coalescence child locations equal actual production tables. + children=np.array([[qi*(2*n)**2+(2*yi+dy)*2*n+2*xi+dx for dy in (0,1) for dx in (0,1)] for qi,yi,xi in rows]) + np.testing.assert_array_equal(children,p.coalesce_fine) + layout=BlockLayout2D.build(n,n,8,8,periodic=(True,True)) + neighbors,offsets=projected_explosion_tables(p,layout) + inverse=np.full(layout.nblocks*np.prod(layout.local_shape),-1,dtype=int) + inverse[layout.owner.ravel()]=np.arange(n*n) + neighbors=inverse[neighbors] + def explode(gc,cv,gf,fv,m,writes): + # Only actual transfer footprint, including inner and outer ghost slabs. + gy,gx=np.divmod(p.ghost_fine.astype(int),2*n) + if mode[0]=='constant': values=gc.reshape(9,-1)[:,p.ghost_coarse] + else: values=reconstruct(gc,p,neighbors,offsets)[0] + assert np.isfinite(values).all() + gf[:,gy,gx]=values; fv[:,gy,gx]=True; writes[gy,gx]+=1 + mode=['constant']; cpu._explode=explode + return spec,maps,p,mode,nu + +def run(n,which): + spec,m,p,mode,nu=setup(n); mode[0]=which + masks=[m.coarse_active,m.fine_active] + arrays=[oracle.initial(mask.shape,l,n,nu) for l,mask in enumerate(masks)] + for a,mask in zip(arrays,masks):a[:,~mask]=np.nan + s=cpu.State(*arrays); rows=[]; initmass=None; momentum0=None; start=time.monotonic() + for tstar in oracle.TIMES: + step=round(tstar*n/oracle.U) + while s.tick<2*step: s,_=cpu.advance_coarse_once(spec,m,s) + arrays=[s.coarse_f,s.fine_f] + row=oracle.metrics(arrays,masks,step,n,nu,initmass) + if initmass is None:initmass=row['mass'] + dot=den=res=0.; parts=[]; momentum=np.zeros(2) + for l,(a,mask) in enumerate(zip(arrays,masks)): + v=a[:,mask];rho=v.sum(0);u=cpu.CX@v/rho;w=cpu.CY@v/rho + x,y=oracle.centers(mask.shape,l);ex,ey,_=oracle.exact(x[mask],y[mask],step,n,nu);vol=.25**l + dot+=vol*np.sum(u*ex+w*ey);den+=vol*np.sum(ex*ex+ey*ey) + parts.append((u,w,ex,ey,vol));momentum+=vol*np.array([np.sum(cpu.CX@v),np.sum(cpu.CY@v)]) + amp=dot/den + res=sum(vol*np.sum((u-amp*ex)**2+(w-amp*ey)**2) for u,w,ex,ey,vol in parts) + if momentum0 is None:momentum0=momentum.copy() + row.update(step=step,exact_energy_quadrature_factor=den/(2*n*n*row['exact_energy']),amplitude_ratio=amp,amplitude_bias=amp-1,orthogonal_residual_l2=float(np.sqrt(res/den)),energy_signed_relative=row['energy']/row['exact_energy']-1,mass_signed_relative=row['mass']/initmass-1,momentum=momentum.tolist(),momentum_change=(momentum-momentum0).tolist()) + rows.append(row) + return dict(n=n,explosion=which,seconds=time.monotonic()-start,rows=rows) +if __name__=='__main__': + a=args + a.output_dir.mkdir(parents=True,exist_ok=True) + result=dict(implementation_revision='43c527b',source=provenance(ROOT),diagnostic_input_sha256={name:hashlib.sha256((ROOT/name).read_bytes()).hexdigest() for name in ('tests/validation/multigrid/cpu_reference.py','tests/unit/test_projected_explosion.py')},precision='FP64 physical populations',phase='arrival_f; macro moments invariant under collision, initialization equilibrium',runs=[]) + out=a.output_dir/f'n{a.n}_ab.json' + for mode in ('constant','projected_linear'): + r=run(a.n,mode);result['runs'].append(r);out.write_text(json.dumps(result,indent=2)+'\n');print(json.dumps(r),flush=True)