ClaudeandClaude Opus 5.5
10069feadb
Modal magnitude R=hypot(a,b) with amplitude error and separate phase; free-intercept ln R rate fit with per-interval rates, log-fit residual, t* units and finite-positive check; spurious-pressure naming with raw rms; INVALID_CONTROL vs CONTROL_SCREEN_FAIL (refined results kept), mixed-source check, raw ratios to both controls; docs narrow the cubic-term scope, startup and placement limitations. oracle.py still byte-identical. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
280 lines
13 KiB
Python
280 lines
13 KiB
Python
"""Standalone CPU manufactured checks for the oblique shear wave.
|
|
|
|
Run from this directory with PYTHONPATH=<checkout>/src (same as test_oracle.py).
|
|
"""
|
|
|
|
import hashlib
|
|
import math
|
|
import unittest
|
|
from pathlib import Path
|
|
|
|
import numpy as np
|
|
import oracle
|
|
import run_case
|
|
import shear_wave as sw
|
|
from oracle import CX, CY, U, TIMES, centers, rectangles
|
|
|
|
N = 32
|
|
NU = U * N / sw.RE
|
|
HERE = Path(__file__).resolve().parent
|
|
ORACLE_SHA256_EEC4677 = "83e4d7425d45538278e9fb181af450a36358bd9a5ed6be57ff611dd9ae5eb35b"
|
|
|
|
|
|
def composite_masks(n):
|
|
"""Level-wide union masks for rectangles(n, 2): root complement + fine box."""
|
|
(x0, y0, x1, y1), = rectangles(n, 2)[0]
|
|
root = np.ones((n, n), dtype=bool)
|
|
root[y0:y1, x0:x1] = False
|
|
fine = np.zeros((2 * n, 2 * n), dtype=bool)
|
|
fine[2 * y0:2 * y1, 2 * x0:2 * x1] = True
|
|
return [root, fine]
|
|
|
|
|
|
def populations(masks, field, t=0.0):
|
|
"""Equilibrium G of an arbitrary (ux, uy, rho) field on every level."""
|
|
out = []
|
|
for level, m in enumerate(masks):
|
|
x, y = centers(m.shape, level)
|
|
ux, uy, rho = field(x, y, t)
|
|
out.append(sw.equilibrium(rho, ux, uy))
|
|
return out
|
|
|
|
|
|
def exact_field(phase=0.0, rate_factor=1.0, mode=sw.MODE, longitudinal=0.0, ripple=0.0,
|
|
offset=0.0):
|
|
k = 2 * np.pi / N * np.array(mode, dtype=float)
|
|
e = np.array([-mode[1], mode[0]], dtype=float) / math.hypot(*mode)
|
|
rate = rate_factor * sw.decay_rate(N, NU)
|
|
kh = k / np.linalg.norm(k)
|
|
|
|
def field(x, y, t):
|
|
theta = k[0] * x + k[1] * y + phase
|
|
a = U * np.exp(-rate * t)
|
|
ux = a * e[0] * np.sin(theta) + longitudinal * kh[0] * np.sin(theta)
|
|
uy = a * e[1] * np.sin(theta) + longitudinal * kh[1] * np.sin(theta)
|
|
rho = 1 + offset + ripple * np.cos(theta)
|
|
return ux, uy, rho
|
|
|
|
return field
|
|
|
|
|
|
class ConventionTests(unittest.TestCase):
|
|
def test_mode_angle_and_polarization(self):
|
|
k = sw.wavevector(N)
|
|
self.assertAlmostEqual(sw.ANGLE_DEG, 26.565051177077994)
|
|
self.assertAlmostEqual(math.degrees(math.atan2(k[1], k[0])), sw.ANGLE_DEG)
|
|
self.assertAlmostEqual(float(sw.E @ k), 0.0)
|
|
self.assertAlmostEqual(float(sw.E @ sw.E), 1.0)
|
|
self.assertAlmostEqual(sw.eps(128), 5 * (2 * math.pi / 128) ** 2)
|
|
|
|
def test_navier_stokes_and_vanishing_advection(self):
|
|
# Independent central differences: du/dt = nu lap u, (u.grad)u = 0, div u = 0.
|
|
x = np.array([2.7, 8.1, 13.3])
|
|
y = np.array([3.2, 7.4, 21.6])
|
|
t, h = 10.0, 1e-3
|
|
f = lambda x, y, t: np.array(sw.exact(x, y, t, N, NU)[:2])
|
|
z = f(x, y, t)
|
|
dt = (f(x, y, t + h) - f(x, y, t - h)) / (2 * h)
|
|
dx = (f(x + h, y, t) - f(x - h, y, t)) / (2 * h)
|
|
dy = (f(x, y + h, t) - f(x, y - h, t)) / (2 * h)
|
|
lap = (f(x + h, y, t) + f(x - h, y, t) + f(x, y + h, t) + f(x, y - h, t) - 4 * z) / h**2
|
|
np.testing.assert_allclose(dt, NU * lap, atol=1e-9)
|
|
np.testing.assert_allclose(z[0] * dx + z[1] * dy, 0.0, atol=1e-12)
|
|
np.testing.assert_allclose(dx[0] + dy[1], 0.0, atol=1e-12)
|
|
self.assertTrue(np.all(sw.exact(x, y, t, N, NU)[2] == 0))
|
|
|
|
|
|
class InitialTests(unittest.TestCase):
|
|
def test_moments_on_each_level(self):
|
|
for level, m in enumerate(composite_masks(N)):
|
|
g = sw.initial(m.shape, level, N, NU)
|
|
x, y = centers(m.shape, level)
|
|
ux, uy, _ = sw.exact(x, y, 0.0, N, NU)
|
|
rho = g.sum(axis=0)
|
|
np.testing.assert_allclose(rho, 1.0, atol=1e-14)
|
|
np.testing.assert_allclose(np.tensordot(CX, g, 1) / rho, ux, atol=1e-15)
|
|
np.testing.assert_allclose(np.tensordot(CY, g, 1) / rho, uy, atol=1e-15)
|
|
|
|
def test_equilibrium_matches_frozen_oracle_formula(self):
|
|
# Same equilibrium as oracle.initial when fed the TG field.
|
|
shape = (8, 8)
|
|
x, y = centers(shape, 1)
|
|
ux, uy, p = oracle.exact(x, y, 0.0, 4.0, 0.01)
|
|
np.testing.assert_array_equal(
|
|
sw.equilibrium(1 + 3 * p, ux, uy), oracle.initial(shape, 1, 4.0, 0.01)
|
|
)
|
|
|
|
|
|
class MetricTests(unittest.TestCase):
|
|
def setUp(self):
|
|
self.masks = composite_masks(N)
|
|
|
|
def sample(self, field, t):
|
|
return sw.metrics(populations(self.masks, field, t), self.masks, t, N, NU, None)
|
|
|
|
def test_exact_field_on_composite_partition(self):
|
|
for t in (0.0, 100.0, 512.0):
|
|
r = self.sample(exact_field(), t)
|
|
self.assertLess(r["velocity_l2"], 1e-12)
|
|
self.assertLess(r["residual_l2"], 1e-12)
|
|
self.assertLess(r["longitudinal_l2"], 1e-12)
|
|
self.assertLess(r["phase_drift_abs"], 1e-12)
|
|
self.assertAlmostEqual(r["amplitude_ratio"], 1.0, places=12)
|
|
self.assertLess(abs(r["amplitude_error"]), 1e-12)
|
|
self.assertLess(r["spurious_pressure_scaled"], 1e-10)
|
|
self.assertEqual(len(r["spurious_pressure_scaled_by_level"]), 2)
|
|
self.assertEqual(r["leaf_cells"], int(sum(m.sum() for m in self.masks)))
|
|
|
|
def test_phase_control(self):
|
|
r = self.sample(exact_field(phase=0.01), 0.0)
|
|
self.assertAlmostEqual(r["phase_drift"], 0.01, delta=1e-9)
|
|
self.assertLess(abs(r["amplitude_error"]), 1e-12) # R = hypot(a, b) is phase-free
|
|
self.assertLess(r["residual_l2"], 1e-12) # a pure phase shift stays modal
|
|
|
|
def test_rate_control_and_fit(self):
|
|
for factor, expected in ((1.0, 0.0), (1.01, 0.01)):
|
|
field = exact_field(rate_factor=factor)
|
|
rows = []
|
|
for tstar in TIMES:
|
|
step = round(tstar * N / U)
|
|
rows.append(dict(self.sample(field, step), step=step))
|
|
fit = sw.fit_rate(rows, N, NU)
|
|
self.assertAlmostEqual(fit["rate_error"], expected, delta=1e-9)
|
|
np.testing.assert_allclose(fit["interval_rate_errors"], expected, atol=1e-9)
|
|
self.assertLess(fit["log_fit_residual_max"], 1e-9)
|
|
self.assertAlmostEqual(fit["exact_rate_tstar"], 5 * (2 * math.pi) ** 2 / sw.RE)
|
|
self.assertAlmostEqual(fit["exact_rate_tstar"], 1.97392088, places=8)
|
|
|
|
def test_fit_ignores_startup_sample(self):
|
|
rows = [dict(tstar=0.0, step=0, a_hat=2 * U, b_hat=0.0)]
|
|
rows += [dict(tstar=s, step=round(s * N / U), b_hat=0.0,
|
|
a_hat=U * math.exp(-sw.decay_rate(N, NU) * round(s * N / U)))
|
|
for s in TIMES[1:]]
|
|
self.assertLess(sw.fit_rate(rows, N, NU)["rate_error"], 1e-12)
|
|
rows[-1]["a_hat"] = float("nan")
|
|
with self.assertRaises(ValueError):
|
|
sw.fit_rate(rows, N, NU)
|
|
|
|
def test_wrong_mode_control(self):
|
|
# Swapping (2,1) for (1,2), i.e. 26.6 vs 63.4 deg, is not silently modal.
|
|
# Composite quadrature need not give an exactly zero projection.
|
|
r = self.sample(exact_field(mode=(1, 2)), 0.0)
|
|
self.assertLess(abs(r["amplitude_ratio"]), 0.05)
|
|
self.assertGreater(r["residual_l2"], 0.95)
|
|
|
|
def test_longitudinal_control(self):
|
|
r = self.sample(exact_field(longitudinal=1e-3), 0.0)
|
|
self.assertAlmostEqual(r["longitudinal_l2"], 1e-3 / U, delta=1e-9)
|
|
# The full-vector residual and velocity error see it too.
|
|
self.assertAlmostEqual(r["residual_l2"], 1e-3 / U, delta=1e-9)
|
|
self.assertAlmostEqual(r["velocity_l2"], 1e-3 / U, delta=1e-9)
|
|
self.assertLess(r["phase_drift_abs"], 1e-12)
|
|
self.assertAlmostEqual(r["amplitude_ratio"], 1.0, places=12)
|
|
|
|
def test_pressure_gauge_and_ripple(self):
|
|
r = self.sample(exact_field(offset=1e-4), 0.0)
|
|
self.assertLess(r["spurious_pressure_scaled"], 1e-10) # one global gauge
|
|
ripple = 3e-5
|
|
r = self.sample(exact_field(ripple=ripple), 0.0)
|
|
# p = ripple/3 cos(theta); rms over a whole number of periods = ripple/(3 sqrt 2).
|
|
# Composite midpoint quadrature is not exactly periodic: 3% tolerance.
|
|
expected = ripple / (3 * math.sqrt(2)) / (U * U / 4)
|
|
self.assertAlmostEqual(r["spurious_pressure_scaled"], expected, delta=0.03 * expected)
|
|
self.assertAlmostEqual(r["spurious_pressure_rms"], expected * U * U / 4,
|
|
delta=0.03 * expected * U * U / 4)
|
|
|
|
def test_volume_weighting_and_double_count(self):
|
|
bad = [np.ones_like(m) for m in self.masks]
|
|
with self.assertRaises(ValueError):
|
|
sw.metrics(populations(bad, exact_field()), bad, 0.0, N, NU)
|
|
# Unweighted (count-based) averaging would change the fit; weights matter.
|
|
r = self.sample(exact_field(), 0.0)
|
|
self.assertEqual(r["leaf_cells"], N * N - self.masks[1].sum() // 4 + self.masks[1].sum())
|
|
|
|
|
|
def record(n, depth, rows, explosion="constant", completed=True, source="s"):
|
|
return dict(n=n, depth=depth, viscosity=U * n / sw.RE, rows=rows, completed=completed,
|
|
explosion=explosion, source=source)
|
|
|
|
|
|
def rows_for(n, rate_error=0.0, residual=1e-4, pressure=1e-3, phase=0.0):
|
|
rate = (1 + rate_error) * sw.decay_rate(n, U * n / sw.RE)
|
|
out = []
|
|
for tstar in TIMES:
|
|
step = round(tstar * n / U)
|
|
out.append(dict(tstar=tstar, step=step, a_hat=U * math.exp(-rate * step), b_hat=0.0,
|
|
amplitude_error=0.0, residual_l2=residual, phase_drift_abs=phase,
|
|
spurious_pressure_scaled=pressure,
|
|
spurious_pressure_rms=pressure * U * U / 4,
|
|
velocity_l2=residual, mass_relative=0.0))
|
|
return out
|
|
|
|
|
|
class AssessTests(unittest.TestCase):
|
|
def test_budget_form_is_max_not_sum(self):
|
|
base = record(128, 1, rows_for(128, rate_error=1e-3, pressure=0.015))
|
|
fine = record(256, 1, rows_for(256, rate_error=2e-4))
|
|
# rate budget = max(2*1e-3, 0.002) = 0.002; 0.0025 would pass a sum form.
|
|
refined = record(128, 2, rows_for(128, rate_error=0.0025))
|
|
report = sw.assess(base, fine, [refined])
|
|
self.assertEqual(report["status"], "FAIL")
|
|
self.assertAlmostEqual(report["refined"]["constant"]["budget"]["rate_error"], 0.002)
|
|
# pressure budget = max(2*0.015, 0.02) = 0.03 (control-dominated).
|
|
self.assertAlmostEqual(
|
|
report["refined"]["constant"]["budget"]["spurious_pressure_scaled"], 0.03)
|
|
self.assertIn("finest", report["refined"]["constant"]["control_ratios"])
|
|
ok = record(128, 2, rows_for(128, rate_error=0.0015), explosion="projected_linear")
|
|
self.assertEqual(sw.assess(base, fine, [ok])["status"], "PASS")
|
|
|
|
def test_control_outcomes_are_distinct(self):
|
|
fine = record(256, 1, rows_for(256))
|
|
refined = record(128, 2, rows_for(128))
|
|
failing = record(128, 1, rows_for(128, residual=0.5)) # complete, > eps(128)
|
|
report = sw.assess(failing, fine, [refined])
|
|
self.assertEqual(report["status"], "CONTROL_SCREEN_FAIL")
|
|
self.assertIn("constant", report["refined"]) # refined results preserved
|
|
good = record(128, 1, rows_for(128))
|
|
self.assertEqual(sw.assess(None, fine, [refined])["status"], "INVALID_CONTROL")
|
|
incomplete = record(128, 1, rows_for(128), completed=False)
|
|
self.assertEqual(sw.assess(incomplete, fine, [refined])["status"], "INVALID_CONTROL")
|
|
mixed = record(256, 1, rows_for(256), source="other")
|
|
self.assertEqual(sw.assess(good, mixed, [refined])["status"], "INVALID_CONTROL")
|
|
|
|
def test_wrong_shapes(self):
|
|
base = record(128, 1, rows_for(128))
|
|
fine = record(256, 1, rows_for(256))
|
|
refined = record(128, 2, rows_for(128))
|
|
with self.assertRaises(ValueError):
|
|
sw.assess(base, record(128, 1, rows_for(128)), [refined])
|
|
with self.assertRaises(ValueError):
|
|
sw.assess(record(128, 1, rows_for(128)), fine, [record(256, 2, rows_for(256))])
|
|
|
|
def test_incomplete_record_fails(self):
|
|
base = record(128, 1, rows_for(128))
|
|
fine = record(256, 1, rows_for(256))
|
|
refined = record(128, 2, rows_for(128), completed=False)
|
|
self.assertEqual(sw.assess(base, fine, [refined])["status"], "FAIL")
|
|
|
|
|
|
class DefaultPathTests(unittest.TestCase):
|
|
def test_frozen_oracle_is_byte_identical(self):
|
|
digest = hashlib.sha256((HERE / "oracle.py").read_bytes()).hexdigest()
|
|
self.assertEqual(digest, ORACLE_SHA256_EEC4677)
|
|
|
|
def test_default_flow_is_frozen_taylor_green(self):
|
|
init, met, limits, text = run_case.flow_functions()
|
|
self.assertIs(init, oracle.initial)
|
|
self.assertIs(met, oracle.metrics)
|
|
self.assertIs(limits, oracle.LIMITS)
|
|
self.assertEqual(text, run_case.INITIALIZATION)
|
|
self.assertEqual(run_case.flow_functions("taylor_green"), (init, met, limits, text))
|
|
shear = run_case.flow_functions("shear_wave")
|
|
self.assertIs(shear[0], sw.initial)
|
|
self.assertIs(shear[1], sw.metrics)
|
|
with self.assertRaises(ValueError):
|
|
run_case.flow_functions("tg")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
unittest.main()
|