"""Independent tensor/oracle tests against host-compiled exact copied headers.""" import argparse, ctypes, hashlib, itertools, json, math, pathlib, subprocess, sys, tempfile HERE = pathlib.Path(__file__).resolve().parent from audit import directions, wrapper def tensor(dim, kind, p, a=0, b=1): t = [[0.0] * dim for _ in range(dim)] if kind == "shear": t[a][b] = t[b][a] = p elif kind == "extension": t[a][a] = p t[b][b] = -p elif kind == "isotropic": for i in range(dim): t[i][i] = p elif kind != "equilibrium": raise ValueError(kind) return t def state(dim, rho, t): # Second-order Hermite projection: sum c_i c_i fneq_i = prescribed tensor. # D2Q9/D3Q19 fourth-order isotropy makes this a distinct construction from audit.py. c = directions(dim) w = [ ( (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 ] eq = [rho * wi for wi in w] neq = [ 4.5 * wi * sum( (v[a] * v[b] - (1 / 3 if a == b else 0)) * t[a][b] for a in range(dim) for b in range(dim) ) for wi, v in zip(w, c) ] assert abs(sum(neq)) < 1e-12 for a in range(dim): assert abs(sum(v[a] * f for v, f in zip(c, neq))) < 1e-12 for b in range(dim): assert abs(sum(v[a] * v[b] * f for v, f in zip(c, neq)) - t[a][b]) < 1e-12 return [x + y for x, y in zip(eq, neq)], eq def norm(t): return math.sqrt(sum(x * x for row in t for x in row)) def standard_tau(tau0, cs, p, rho): # Solve tau*(tau-tau0)=9 Cs^2 P/(sqrt(2) rho), independently via bisection. k = 9 * cs * cs * p / (math.sqrt(2) * rho) lo = tau0 hi = tau0 + max(1.0, k / tau0) for _ in range(90): mid = (lo + hi) / 2 if mid * (mid - tau0) > k: hi = mid else: lo = mid return (lo + hi) / 2 def build(root, mode, dim, cs, closed, iters, tmp): name = f"{mode}_d{dim}_cs{cs}_c{closed}_i{iters}" src = tmp / (name + ".cpp") lib = tmp / (name + ".so") src.write_text(wrapper(dim, closed, cs, iters)) subprocess.run( [ "g++", "-std=c++11", "-O2", "-shared", "-fPIC", "-I", str(root / mode / "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, rho, tau0): ar = ctypes.c_float * len(f) return float(fn(ar(*f), ar(*eq), rho, 1 / tau0)) return call def run(root): n = 0 maxerr = 0.0 maxres = 0.0 examples = [] with tempfile.TemporaryDirectory(prefix="les-candidate-") as tmp: for dim, cs in itertools.product((2, 3), (0.0, 0.08, 0.16, 0.24)): kernels = { (mode, c, i): build(root, mode, dim, cs, c, i, pathlib.Path(tmp)) for mode in ("baseline", "patched") for c, i in ((1, 3), (0, 3), (0, 200)) } for rho, tau0 in itertools.product((0.7, 1.0, 1.3), (0.6, 0.8, 1.2)): # Pure shear/extension all coordinate planes and both signs; Pi/rho fixed under density scaling. cases = [("equilibrium", 0.0, 0, 1), ("isotropic", 0.01 * rho, 0, 1)] cases += [ (kind, sign * 0.02 * rho, a, b) for kind in ("shear", "extension") for sign in (-1, 1) for a, b in itertools.combinations(range(dim), 2) ] for kind, p, a, b in cases: t = tensor(dim, kind, p, a, b) f, eq = state(dim, rho, t) P = norm(t) out = {key: fn(f, eq, rho, tau0) for key, fn in kernels.items()} want = 1 / standard_tau(tau0, cs, P, rho) got = out["patched", 1, 3] err = abs(got - want) maxerr = max(maxerr, err) assert err < 3e-6, (dim, cs, rho, tau0, kind, p, a, b, got, want) tau = 1 / got res = abs( tau * (tau - tau0) - 9 * cs * cs * P / (math.sqrt(2) * rho) ) maxres = max(maxres, res) assert res < 2e-6, (res, dim, kind) # No edits to iterative path: exact float results must remain identical. for it in (3, 200): assert out["baseline", 0, it] == out["patched", 0, it] if kind in ("shear", "extension", "equilibrium") or cs == 0: assert abs(out["patched", 0, 200] - want) < 3e-6 # Three iterations need not equal closed form; they approach the same root. assert ( abs(out["patched", 0, 200] - want) <= abs(out["patched", 0, 3] - want) + 3e-7 ) if kind == "equilibrium" or cs == 0: assert abs(out["baseline", 1, 3] - got) < 3e-6 elif kind in ("shear", "extension"): assert out["baseline", 1, 3] > got if ( rho == 1.0 and tau0 == 0.8 and cs == 0.16 and p >= 0 and a == 0 and b == 1 ): examples.append( dict( dim=dim, kind=kind, baseline_closed=out["baseline", 1, 3], patched_closed=got, iterative3=out["patched", 0, 3], iterative200=out["patched", 0, 200], standard=want, ) ) n += 1 # Explicit fixed-P/nonunit-rho test, distinct from scaling Pi and rho together. for rho in (0.7, 1.3): t = tensor(dim, "shear", 0.02) f, eq = state(dim, rho, t) got = kernels["patched", 1, 3](f, eq, rho, 0.8) assert abs(got - 1 / standard_tau(0.8, cs, norm(t), rho)) < 3e-6 n += 1 # Cap path and trace behavior unchanged; arbitrary stress is a helper-only diagnostic. t = tensor(dim, "shear", 0.2) f, eq = state(dim, 1.0, t) for it in (3, 200): assert kernels["baseline", 0, it](f, eq, 1.0, 0.503) == kernels[ "patched", 0, it ](f, eq, 1.0, 0.503) return dict( parameter_cases=n, max_closed_omega_error=maxerr, max_quadratic_residual=maxres, examples=examples, gpu_run=False, ) if __name__ == "__main__": from prepare import verify_candidate ap = argparse.ArgumentParser() ap.add_argument("--candidate-dir", type=pathlib.Path, required=True) ap.add_argument("--output", type=pathlib.Path, required=True) args = ap.parse_args() if args.output.exists(): ap.error("Refusing existing evidence file") verify_candidate(args.candidate_dir) result = run(args.candidate_dir.resolve()) args.output.write_text(json.dumps(result, indent=2) + "\n") print(json.dumps(result, indent=2))