217 lines
7.8 KiB
Python
217 lines
7.8 KiB
Python
"""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))
|