diff --git a/tests/validation/experimental_full9/README.md b/tests/validation/experimental_full9/README.md new file mode 100644 index 0000000..f119c1a --- /dev/null +++ b/tests/validation/experimental_full9/README.md @@ -0,0 +1,38 @@ +# Isolated experimental full-nine inlet comparison + +This is preparation for one separately admitted, at-most320-step bodyless robustness comparison. It is not public support, a default change, proof of root cause, or boundary-accuracy acceptance. + +## Source and isolation +Numerical commit: `926ca27407c26ca11283282a6e53249925725376`, based on history diagnostic `e92de3742476ea5afd6b84ff9538992df48548f7`. The new manifest pins numerical files and separately hashes the two runtime harness files. Actual checkout HEAD/dirty state is recorded on execution; --source must resolve to the executing harness checkout. The historical manifest is unchanged and intentionally rejects the new numerical tree. The original Zou–He default and IDs0..3 remain unchanged. + +Explicit scheme `experimental_regularized_full9` is ID4. Configuration guards restrict it to D2Q9, SRT, double-buffer, FP32 storage, bounce-back y walls, and no LES; shifted and unshifted storage are tested. The bodyless runner uses exactly the old config factory, changes only inlet scheme, and removes only the original cylinder. Original initialization, profile, boundary planes, outlet/all-nine damping, viscosity, timestep and health gates remain unchanged. + +The new helper computes rho from the existing west known-population closure. At the prescribed velocity, known Kq=Fq−Feq gives Pi_xx=2(K2+K6+K8), Pi_yy=K3+K4+2(K6+K8), Pi_xy=2(K6−K8). It reconstructs all nine pre-collision populations from equilibrium plus the second-order stress projection. The existing step kernel then applies SRT exactly once. Old-buffer timing is unchanged. This preserves prescribed physical velocity moments while intentionally changing the boundary's kinetic reconstruction; it is not numerically identical to Zou–He. + +Primary mapping: [Palabos stress construction](https://unigespc.gitlab.io/palabos/docs/html/d0/dfa/boundary_templates_8h_source.html), [velocity completion before base collision](https://unigespc.gitlab.io/palabos/docs/html/d2/d68/regularized_boundary_dynamics_8hh_source.html), [Latt–Chopard equations10–12](https://arxiv.org/pdf/physics/0506157). This independently implemented mapping does not use Celeris's differently defined existing incoming-only donor-damped `regularized` option. + +Near-wall inlet known inputs are unchanged. This comparison includes their existing interaction; it does not silently move walls or repair initialization. An EsoPull dispatch arm is present because that function is preprocessed even in double-buffer builds; the configuration guard rejects running that streaming variant. + +## Tests and bounded analytical check + +```sh +python tests/validation/experimental_full9/test_full9.py +python tests/validation/experimental_full9/test_runner.py +OPENBLAS_NUM_THREADS=1 OMP_NUM_THREADS=1 python tests/validation/experimental_full9/check_linear_strip.py +python tests/validation/experimental_full9/run_bodyless.py --record-inlet-history --output /tmp/new-full9-preflight +``` + +The host tests compile the actual helper and canonical macro/SRT headers with g++, in shifted/unshifted modes. They check equilibrium, mass/momentum, prescribed stress, ignored wrapped q1/5/7, exactly-one collision, source ordering, legacy IDs/defaults and rejected configurations. These tests do not certify CUDA compilation. + +The fixed homogeneous linear comparison uses the original five k samples, width32, rho1, U=.0028664576821029186, tau=.524 and zero far incoming perturbation. The prior3π/4 Zou–He radius1.0186489 becomes .9802484 with full9; all five sampled full9 radii are below1. This is a finite sampled surrogate result, not a continuous stability proof or a match to the actual dominant≈2π/3 history spectrum. No parameter tuning or new flow is part of this check. + +## Exact admitted command +From the final clean reviewed checkout, after explicit GPU admission: + +```sh +CUDA_VISIBLE_DEVICES= timeout --signal=TERM --kill-after=30s 240s python tests/validation/experimental_full9/run_bodyless.py --record-inlet-history --execute --gpu-admitted --wall-seconds 180 --output output/new-experimental-full9-bodyless +``` + +The shared diagnostic still stops on its unchanged first-health/320-step/wall ceiling. Optional raw history records the fixed0..196 window only, with the existing12MiB bound and completeness receipt; evolution can continue beyond196. Use a new output directory and preserve old evidence. + +A healthy320-step result would support only bounded robustness to this exact startup. Failure remains evidence, not a reason to tune the scheme after seeing results. Compare step-count/health, finite fields, density/Mach, inlet target moments and retained modal history; conservation and long-time/cylinder accuracy require later separate tests. No default/dev integration is authorized by this experiment. diff --git a/tests/validation/experimental_full9/check_linear_strip.py b/tests/validation/experimental_full9/check_linear_strip.py new file mode 100644 index 0000000..3c95f6d --- /dev/null +++ b/tests/validation/experimental_full9/check_linear_strip.py @@ -0,0 +1,31 @@ +"""Bounded fixed-parameter linear comparison; no nonlinear flow or GPU run.""" +import json +from pathlib import Path +import numpy as np +from scipy.linalg import eigvals +cx=np.array([0,1,-1,0,0,1,-1,1,-1]);cy=np.array([0,0,0,1,-1,1,-1,-1,1]);w=np.array([4/9]+[1/9]*4+[1/36]*4) +U=.0028664576821029186;omega=1/.524;cu=cx*U +base=w*(1+3*cu+4.5*cu**2-1.5*U**2) +E=w[:,None]*((1+3*cu+4.5*cu**2-1.5*U*U)[:,None]+(3*cx+9*cu*cx-3*U)[:,None]*(cx-U)[None,:]+(3*cy+9*cu*cy)[:,None]*cy[None,:]) +C=(1-omega)*np.eye(9)+omega*E +r=np.array([1,0,2,1,1,0,2,0,2])/(1-U);I=np.eye(9) +B=I.copy();B[1]=I[2]+2/3*U*r;B[5]=I[6]+.5*(I[4]-I[3])+U/6*r;B[7]=I[8]+.5*(I[3]-I[4])+U/6*r +T=np.array([[0,0,2,0,0,0,2,0,2],[0,0,0,1,1,0,2,0,2],[0,0,0,0,0,0,2,0,-2]]) +R=w[:,None]*np.array([4.5*(cx*cx-1/3),4.5*(cy*cy-1/3),9*cx*cy]).T +Bfull=base[:,None]*r+R@T@(I-base[:,None]*r) +assert np.max(abs(Bfull[:,[1,5,7]]))<1e-15 +assert np.max(abs(np.array([np.ones(9),cx,cy])@Bfull-np.array([r,U*r,np.zeros(9)])))<1e-14 + +def radius(n,k,b): + A=np.zeros((9*n,9*n),complex) + for x in range(n): + local=C@b if x==0 else C + for q in range(9): + donor=x-int(cx[q]) + if 0<=donor +#define __device__ +#define __forceinline__ inline +#define __fdiv_rn(a,b) ((a)/(b)) +#define LATTICE_D2Q9 0 +#define LATTICE_MODEL 0 +#define NQ 9 +#define CS2 (1.0f/3.0f) +#define NY 330 +#define U0 .01f +#define W0_VAL (4.0f/9.0f) +#define WS_VAL (1.0f/9.0f) +#define WE_VAL (1.0f/36.0f) +#define USE_DDF_SHIFTING {shifted} +#define INLET_SCHEME 4 +#include "operators/macro.cuh" +#include "operators/collision_srt.cuh" +#include "boundary/inlet/common.cuh" +#include "boundary/inlet/experimental_full9.cuh" +extern "C" int apply(float* f,float u,float v,int collision) {{ + apply_experimental_full9_inlet_d2q9(f,u,v); + int count=0; + if (collision && inlet_scheme_uses_post_collision_ghost()) {{ + float rho,ux,uy,eq[9]; compute_rho_u(f,rho,ux,uy); compute_feq(rho,ux,uy,eq); + collide_srt_no_force(f,eq,1.0f/.524f); ++count; + }} + return count; +}} +''' + src=Path(TEMP.name)/f'test{shifted}.cpp';src.write_text(code);lib=src.with_suffix('.so') + subprocess.run(['g++','-std=c++11','-shared','-fPIC','-O2','-ffp-contract=off','-I',str(ROOT/'src/CelerisLab/lbm/kernels'),str(src),'-o',str(lib)],check=True) + dll=ctypes.CDLL(str(lib));dll.apply.argtypes=[ctypes.POINTER(ctypes.c_float),ctypes.c_float,ctypes.c_float,ctypes.c_int];dll.apply.restype=ctypes.c_int;libs.append(dll) + +def eq(rho,u,v): + cu=cx*u+cy*v;return w*rho*(1+3*cu+4.5*cu*cu-1.5*(u*u+v*v)) +def reference(f,u,v): + rho=(f[0]+f[3]+f[4]+2*(f[2]+f[6]+f[8]))/(1-u);e=eq(rho,u,v);k=f-e + pi=np.array([2*(k[2]+k[6]+k[8]),k[3]+k[4]+2*(k[6]+k[8]),2*(k[6]-k[8])]) + stress=w*(4.5*((cx*cx-1/3)*pi[0]+(cy*cy-1/3)*pi[1])+9*cx*cy*pi[2]) + return e,stress,pi + +def apply(f,u,v,shifted,collision=0): + raw=np.asarray(f-shifted*w,dtype=np.float32);count=libs[shifted].apply(raw.ctypes.data_as(ctypes.POINTER(ctypes.c_float)),u,v,collision) + return raw.astype(float)+shifted*w,count + +def test_moments_stress_and_shift(): + rng=np.random.default_rng(92) + for _ in range(20): + f=eq(1.02,.003,0)+rng.normal(0,.001,9);e,s,pi=reference(f,.003,.001) + outputs=[] + for shifted in [0,1]: + out,count=apply(f,.003,.001,shifted);outputs.append(out) + assert count==0 + np.testing.assert_allclose(out,e+s,atol=2e-7,rtol=0) + np.testing.assert_allclose(M@out,M@e,atol=3e-7,rtol=0) + np.testing.assert_allclose(Q@(out-e),pi,atol=2e-7,rtol=0) + np.testing.assert_allclose(*outputs,atol=2e-7,rtol=0) + +def test_unknown_inputs_ignored_and_equilibrium(): + for shifted in [0,1]: + f=eq(1.01,.004,0);out,_=apply(f,.004,0,shifted) + np.testing.assert_allclose(out,f,atol=2e-7,rtol=0) + f[[1,5,7]]=np.nan;ignored,_=apply(f,.004,0,shifted) + assert np.isfinite(ignored).all() + np.testing.assert_array_equal(out,ignored) + +def test_exactly_once_srt(): + f=eq(1,.003,0)+np.array([.001,0,.002,.003,-.001,0,.002,0,-.002]);e,s,_=reference(f,.003,0) + for shifted in [0,1]: + out,count=apply(f,.003,0,shifted,1) + assert count==1 + np.testing.assert_allclose(out,e+(1-1/.524)*s,atol=3e-7,rtol=0) + assert np.max(abs(out-(e+(1-1/.524)**2*s)))>1e-4 + source=(ROOT/'src/CelerisLab/lbm/kernels/step/one_step_double.cu').read_text() + part=source.split('#if DIM == 2',2)[-1].split('#elif DIM == 3')[0] + assert part.index('apply_boundary_pull(f,')0 else '-Inf')) -def source_config(source): - manifest=json.loads((HERE/'source_manifest.json').read_text()) - for name,digest in manifest['sha256'].items(): +def source_config(source, manifest_path=None): + manifest=json.loads((HERE/'source_manifest.json' if manifest_path is None else manifest_path).read_text()) + for name,digest in {**manifest['sha256'],**manifest.get('runtime_sha256',{})}.items(): if hashlib.sha256((source/name).read_bytes()).hexdigest()!=digest: raise ValueError('Pinned source mismatch: '+name) sys.path.insert(0,str(source/'src')); sys.path.insert(0,str(source/'tests/validation'))