"""Bounded host-only audit of exact production Smagorinsky helper; stdlib only.""" import argparse, ctypes, hashlib, json, math, pathlib, subprocess, tempfile HERE = pathlib.Path(__file__).resolve().parent DEFAULT_ROOT = HERE.parents[2] HEADERS = ["core/descriptors.cuh", "operators/turbulence_smag.cuh"] def directions(dim): result = [(0,) * dim] for axis in range(dim): v = [0] * dim v[axis] = 1 result += [tuple(v), tuple(-x for x in v)] for sign in (1, -1): for a in range(dim): for b in range(a + 1, dim): v = [0] * dim v[a] = 1 v[b] = sign result += [tuple(v), tuple(-x for x in v)] return result def populations(dim, kind, p=0.02): c = directions(dim) eq = [ ( (4 / 9 if dim == 2 else 1 / 3) if not any(v) else ( (1 / 9 if dim == 2 else 1 / 18) if sum(x * x for x in v) == 1 else 1 / 36 ) ) for v in c ] neq = [0.0] * len(c) if kind == "shear": for i, v in enumerate(c): if v[0] and v[1]: neq[i] = p * v[0] * v[1] / 4 elif kind == "isotropic": for i, v in enumerate(c): if sum(x * x for x in v) == 1: neq[i] = p / 2 neq[0] = -dim * p elif kind != "equilibrium": raise ValueError(kind) return [a + b for a, b in zip(eq, neq)], eq def stress(f, eq, dim): c = directions(dim) return [ [sum((a - b) * v[i] * v[j] for a, b, v in zip(f, eq, c)) for j in range(dim)] for i in range(dim) ] def norm(matrix, dev=False): n = len(matrix) tr = sum(matrix[i][i] for i in range(n)) / n if dev else 0 return math.sqrt( sum( (matrix[i][j] - (tr if i == j else 0)) ** 2 for i in range(n) for j in range(n) ) ) def tau_closed(tau0, cs, pnorm, rho=1.0, standard=True): factor = math.sqrt(2) if standard else 1.0 return 0.5 * (tau0 + math.sqrt(tau0 * tau0 + 18 * factor * cs * cs * pnorm / rho)) def iterative(tau0, cs, pnorm_dev, rho=1.0, iters=3, cap=20.0, relax=0.7): tau = tau0 nu0 = (tau0 - 0.5) / 3 for _ in range(iters): s = 3 * pnorm_dev / (math.sqrt(2) * rho * tau) nut = min(cs * cs * s, cap * nu0) new = min(100.0, max(0.500001, 0.5 + 3 * (nu0 + nut))) tau = relax * new + (1 - relax) * tau return min(1.99, max(0.01, 1 / tau)) def wrapper(dim, closed, cs, iters=3, gpu=False): prefix = ( "" if gpu else "#define __device__\n#define __constant__\n#define __forceinline__ inline\n" ) declarations = f"""#include {prefix} #define LATTICE_D2Q9 0 #define LATTICE_D3Q19 1 #define LATTICE_MODEL {0 if dim==2 else 1} #define DIM {dim} #define NQ {9 if dim==2 else 19} #define OMEGA_COLLISION_MIN .01f #define OMEGA_COLLISION_MAX 1.99f #define LES_CS {cs:.9f}f #define LES_CLOSED_FORM {int(closed)} #define LES_FP_ITERS {iters} #include "core/descriptors.cuh" #include "operators/turbulence_smag.cuh" """ if gpu: return ( declarations + """extern "C" __global__ void probe(const float* f,const float* eq,float rho,float omega,float* out) { if(threadIdx.x==0 && blockIdx.x==0) out[0]=compute_omega_smag(f,eq,rho,omega); }\n""" ) return ( declarations + """extern "C" float probe(const float* f,const float* eq,float rho,float omega) {return compute_omega_smag(f,eq,rho,omega);}\n""" ) def compile_host(root, directory, dim, closed, cs, iters): stem = f"d{dim}_closed{closed}_cs{cs}_i{iters}" src = directory / (stem + ".cpp") lib = directory / (stem + ".so") src.write_text(wrapper(dim, closed, cs, iters)) subprocess.run( [ "g++", "-std=c++11", "-O2", "-shared", "-fPIC", "-I", str(root / "src/CelerisLab/lbm/kernels"), str(src), "-o", str(lib), ], check=True, capture_output=True, ) dll = ctypes.CDLL(str(lib)) fn = dll.probe fn.argtypes = [ctypes.POINTER(ctypes.c_float)] * 2 + [ctypes.c_float] * 2 fn.restype = ctypes.c_float def call(f, eq, tau0): arr = ctypes.c_float * len(f) return float(fn(arr(*f), arr(*eq), 1.0, 1 / tau0)) return call def run(root): rows = [] checks = 0 with tempfile.TemporaryDirectory(prefix="les-host-") as tmp: for dim in (2, 3): for cs in (0.16, 0.0): kernels = { (closed, iters): compile_host( root, pathlib.Path(tmp), dim, closed, cs, iters ) for closed, iters in ((1, 3), (0, 3), (0, 200)) } for case, kind, p, tau0 in [ ("equilibrium", "equilibrium", 0.0, 0.8), ("pure_shear", "shear", 0.02, 0.8), ("isotropic", "isotropic", 0.02, 0.8), ("cap_active", "shear", 0.2, 0.503), ]: f, eq = populations(dim, kind, p) tensor = stress(f, eq, dim) pf, pd = norm(tensor), norm(tensor, True) actual = { f"closed{c}_iterations{i}": fn(f, eq, tau0) for (c, i), fn in kernels.items() } expected = { "closed1_iterations3": min( 1.99, max(0.01, 1 / tau_closed(tau0, cs, pf, standard=False)), ), "closed0_iterations3": iterative(tau0, cs, pd), "closed0_iterations200": iterative(tau0, cs, pd, iters=200), } for key in expected: assert abs(actual[key] - expected[key]) < 3e-6, ( dim, cs, case, key, actual, expected, ) checks += 1 standard = 1 / tau_closed(tau0, cs, pd) if case == "pure_shear" and cs: assert abs(actual["closed0_iterations200"] - standard) < 3e-6 assert abs(actual["closed1_iterations3"] - standard) > 1e-3 checks += 2 if case == "isotropic" and cs: assert abs(actual["closed0_iterations3"] - 1 / tau0) < 3e-6 assert actual["closed1_iterations3"] < 1 / tau0 - 1e-3 checks += 2 rows.append( dict( dim=dim, cs=cs, case=case, tau0=tau0, stress=tensor, actual_omega=actual, expected_omega=expected, standard_uncapped_deviatoric_omega=standard, ) ) hashes = { name: hashlib.sha256( (root / "src/CelerisLab/lbm/kernels" / name).read_bytes() ).hexdigest() for name in HEADERS } return dict( source_root=str(root), source_headers=hashes, compiler=subprocess.check_output(["g++", "--version"], text=True).splitlines()[ 0 ], host_checks=checks, gpu_executed=False, results=rows, ) if __name__ == "__main__": ap = argparse.ArgumentParser() ap.add_argument("--root", type=pathlib.Path, default=DEFAULT_ROOT) ap.add_argument("--output", type=pathlib.Path, required=True) args = ap.parse_args() report = run(args.root.resolve()) args.output.write_text(json.dumps(report, indent=2) + "\n") print( f"PASS: {report['host_checks']} host checks; exact production header, no GPU; {args.output}" )