Files
CelerisLab/tests/unit/test_duetsch_protocol.py
ClaudeandClaude Opus 5.5 48ee0982c9 tests/docs: unique-name driver loading; moving-path lifecycle and settings notes; VIV predeclare renamed
Tests load research drivers via importlib under unique module names
instead of a bare 'import run' (shared-process collision between the
Dütsch and VIV drivers). lifecycle.md qualifies immutable geometry and the
host-only get_flags for translating circles; supported_settings.md points to
the dense-only moving-circle feature; leontini_viv_predeclare.md is renamed
viv_re150_predeclare.md (Re 150 Kolahdouz/Qin case) with references updated.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-30 16:12:05 +00:00

195 lines
8.9 KiB
Python

"""CPU tests of the Dütsch protocol arithmetic (docs/validation/duetsch_re100_kc5_protocol.md §6a)."""
import importlib.util
import math
from pathlib import Path
import numpy as np
import pytest
_spec = importlib.util.spec_from_file_location(
"duetsch_protocol", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/protocol.py")
p = importlib.util.module_from_spec(_spec)
_spec.loader.exec_module(p)
def _load(rel, name):
"""Load a research driver module under a unique name (no bare ``import run``)."""
import sys
path = Path(__file__).resolve().parents[2] / rel
sys.path.insert(0, str(path.parent)) # for the driver's own sibling imports
try:
spec = importlib.util.spec_from_file_location(name, path)
mod = importlib.util.module_from_spec(spec)
spec.loader.exec_module(mod)
finally:
sys.path.remove(str(path.parent))
return mod
T = 2520
def test_phase_and_cycle_bookkeeping():
assert p.phase(0, T) == 270.0 and p.cycle(0, T) == 0 # +A, at rest, phi = -90
assert p.phase(T // 4, T) == 0.0 and p.cycle(T // 4, T) == 1 # mean position, u = -U_max
assert p.position(0, T, 1.0) == 1.0
assert p.cycle(T // 4 - 1, T) == 0 and p.cycle(T // 4 + T - 1, T) == 1
assert p.cycle(T // 4 + T, T) == 2
assert p.phase(T // 6, T) == pytest.approx(330.0) and p.cycle(T // 6, T) == 0 # startup quarter
assert p.phase(7 * T // 6, T) == pytest.approx(330.0) and p.cycle(7 * T // 6, T) == 1
for ph, frac in ((180.0, .75), (210.0, 5 / 6)):
s = round(frac * T) + T
assert p.phase(s, T) == pytest.approx(ph) and p.bin_phase(s, T) == ph
assert p.bin_phase(round(.75 * T) + T - 1, T) is None
assert sum(p.bin_phase(s, T) == 180.0 for s in range(T, 2 * T)) == T // 360
def test_settling_and_common_window():
passes = {k: True for k in range(3, 12)}
assert p.settled_start(passes) == 7 # pairs 3->4, 4->5, 5->6
late = {**passes, 5: False}
assert p.settled_start(late) == 9 # pairs 6->7, 7->8, 8->9
assert p.common_window([7, 9], [passes, passes]) == (9, 10, 11)
assert p.common_window([7, 10], [passes, passes]) is None # K + 2 = 12 beyond the cap
assert p.common_window([7, None], [passes, passes]) is None
assert p.common_window([7, 7], [passes, {**passes, 8: False}]) is None # must stay settled
# a case settled early must keep passing while waiting for the common K
assert p.common_window([7, 9], [{**passes, 8: False}, passes]) is None
long = {k: True for k in range(3, 32)}
assert p.common_window([10], [long]) is None # default cap 11
assert p.common_window([29], [long], last_scored=31) == (29, 30, 31)
assert p.common_window([30], [long], last_scored=31) is None
def test_solid_aware_bilinear():
field = np.arange(16.0).reshape(4, 4)
fluid = np.ones((4, 4), bool)
value, used, total = p.bilinear_fluid(field, fluid, 1.25, 1.5)
assert total == pytest.approx(1.0) and value == pytest.approx(0.75 * 0.5 * (5 + 9) + 0.25 * 0.5 * (6 + 10))
fluid[1, 1] = False # one solid corner
value, used, total = p.bilinear_fluid(field, fluid, 1.25, 1.5)
assert total == pytest.approx(1.0 - 0.75 * 0.5)
assert all(not (i == 1 and j == 1) for i, j, _ in used)
assert sum(w for *_, w in used) == pytest.approx(1.0)
fluid[:] = False
assert p.bilinear_fluid(field, fluid, 1.25, 1.5) == (None, (), 0.0)
def test_interval_bases_and_morison_fit_recover_known_coefficients():
b = np.array([p.interval_bases(n, T) for n in range(T // 4, T // 4 + 3 * T)])
phi = 2 * np.pi * (np.arange(T // 4, T // 4 + 3 * T) + .5) / T - np.pi / 2
assert np.allclose(b[:, 1], np.sin(phi), atol=1e-5) # midpoint ~ interval mean
f = 2.09 * b[:, 0] - (math.pi ** 2 * 1.45 / p.KC) * b[:, 1]
assert p.morison_fit(f, b[:, 0], b[:, 1]) == pytest.approx((2.09, 1.45), abs=1e-9)
def test_geometry_repeats_exactly_across_cycles():
A = p.amplitude(20)
for m in (0, 1, T // 4, T // 2 - 1, 7 * T // 12):
for c in (1, 5, 11):
assert p.position(m + c * T, T, A) == p.position(m, T, A)
assert p.phase(m + c * T, T) == p.phase(m, T)
v0 = p.position(m + 1, T, A) - p.position(m, T, A)
assert p.position(m + 1 + c * T, T, A) - p.position(m + c * T, T, A) == v0
assert p.cycle(T // 4 + 5 * T, T) == 6 # unwrapped index stays separate
assert p.interval_bases(123 + 7 * T, T) == p.interval_bases(123, T)
def test_exact_point_inside_body_is_missing_even_with_fluid_corners():
field = np.zeros((6, 6))
fluid = np.ones((6, 6), bool)
fluid[3, 3] = False # only one corner solid
centre, radius = (3.0, 3.0), 0.3
value, used, total = p.sample(field, fluid, 3.2, 3.1, centre, radius) # inside the circle
assert value is None and used == () and total == 0.0
value, used, total = p.sample(field, fluid, 3.35, 3.0, centre, radius) # just outside
assert value is not None and total > 0.0
def test_generated_config_passes_the_real_loader(tmp_path):
import json
run = _load("research/duetsch_re100_kc5/run.py", "duetsch_run")
from CelerisLab.config import load_lbm_config
cfg_dict, u_max, n = run.config(20, 10.0, 2520)
path = tmp_path / "config_lbm.json"
path.write_text(json.dumps(cfg_dict))
cfg = load_lbm_config(str(path))
assert (cfg.nx, cfg.ny) == (n, n) == (200, 200)
assert cfg.collision == "MRT" and cfg.streaming == "double_buffer" and not cfg.ddf_shifting
assert cfg.velocity == 0.0 and cfg.viscosity == pytest.approx(u_max * 20 / 100)
assert cfg.y_wall_bc == "bounce_back"
def _analyze():
return _load("research/duetsch_re100_kc5/analyze.py", "duetsch_analyze")
def test_hf_energy_is_parseval_weighted():
a = _analyze()
n = 360
t = np.arange(n)
x = 1.0 + np.cos(2 * np.pi * 5 * t / n) + 2.0 * np.cos(2 * np.pi * 30 * t / n) + 0.5 * np.cos(np.pi * t)
e = a.one_sided_energy(np.fft.rfft(x), n) / n
assert e.sum() == pytest.approx(np.sum(x ** 2)) # Parseval
expected = (2.0 ** 2 / 2 + 0.5 ** 2) / np.mean(x ** 2) # bins >= 21 (30 and Nyquist)
assert e[21:].sum() / e.sum() == pytest.approx(expected)
def test_horizon_excludes_future_coverage_and_events(monkeypatch):
a = _analyze()
T, horizon = 360, 4
steps = round((0.25 + 6) * T) # cycles 1..6 recorded
diag = np.zeros((steps, 15))
diag[round((0.25 + horizon) * T):, 13] = 1.0 # synthetic future no_donor
rows = [(s, p.cycle(s, T), 180.0, 0.0, 0.0, 0.1, 0.0, 1.0) for s in range(steps)
if p.bin_phase(s, T) is not None and p.cycle(s, T) > p.DISCARD]
case = dict(path="c", D=20, T=T, box=10, n=200, u_max=.04, nu=.008, tau=.524, A=1.0,
steps=steps, cycles=6, F_ref=1.0, force=np.zeros((steps, 2)), diag=diag,
samples=np.asarray(rows, float))
monkeypatch.setattr(a, "load", lambda path: case)
out = []
monkeypatch.setattr("builtins.print", lambda text: out.append(text))
a.main(["c", "--horizon", str(horizon)])
import json
entry = json.loads(out[0])["cases"][0]
assert sorted(int(k) for k in entry["coverage"]) == [3, 4] # nothing from cycles 5, 6
assert entry["no_donor"] == 0.0
def _capture():
spec = importlib.util.spec_from_file_location(
"duetsch_capture", Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5/capture.py")
mod = importlib.util.module_from_spec(spec)
spec.loader.exec_module(mod)
return mod
def test_capture_crop_downsample_and_trends():
c = _capture()
n, D = 2400, 40
xm = ym = 0.5 * (n - 1)
A, R = p.amplitude(D), 20.0
box = c.crop_box(xm, ym, D, n, (xm - A - R, xm + A + R, ym - R, ym + R))
assert box == (880, 1520, 960, 1440) and box[1] - box[0] + 1 == 641 and box[3] - box[2] + 1 == 481
with pytest.raises(ValueError):
c.crop_box(xm, ym, D, 600, (0, 1, 0, 1)) # leaves the domain
with pytest.raises(ValueError):
c.crop_box(xm, ym, D, n, (xm - 9 * D, xm, ym, ym)) # clips the swept body
field = np.arange(64.0).reshape(8, 8)
fluid = np.ones((8, 8), bool)
fluid[:4, :4] = False
down = c.downsample(field, fluid, 4)
assert np.isnan(down[0, 0]) and down[1, 1] == field[4:, 4:].mean()
rho = np.ones((n, n)); rho[:, 1:81] = 1.01
row = c.trend_row(10, 3, 180.0, rho, np.zeros((n, n)), np.zeros((n, n)),
np.ones((n, n), bool), D, n, box)
assert row[7] == pytest.approx(1.01) and row[9] == pytest.approx(1.0)
def test_dump_steps_are_the_first_step_of_each_bin():
T = 5040
for k in (11, 21, 30, 31):
starts = [s for s in range(T // 4 + (k - 1) * T, T // 4 + k * T)
if p.bin_phase(s, T) is not None and p.phase(s, T) == p.bin_phase(s, T)]
assert len(starts) == 3 and all(p.cycle(s, T) == k for s in starts)