diff --git a/docs/validation/taylor_green.md b/docs/validation/taylor_green.md index d710f4d..f48a817 100644 --- a/docs/validation/taylor_green.md +++ b/docs/validation/taylor_green.md @@ -159,3 +159,87 @@ theoretical stress-consistency certification or production-kernel change. ## Repository integration The four Python harness files remain together under tests/validation/taylor_green, so runtime provenance hashes precisely this harness directory. Formatting changes produce a new harness identity; preparation and integrated receipts must not be mixed. All analytic definitions, numerical floors, assessment gates, sample times, and source-import assertions are unchanged. Production source is unchanged from the cd3 snapshot; these periodic tests do not waive the separate channel GPU gate. + +## Oblique shear wave (test-only diagnostic, not GPU-qualified) + +`shear_wave.py` reuses the frozen oracle's constants, cell centers and +geometry. `oracle.py` is byte-identical, and `test_shear_wave.py` pins its +eec4677 SHA-256. `run_case.py --flow shear_wave` selects it. The default flow +is still the frozen Taylor–Green oracle, and the default receipt is unchanged. + +Convention: integer mode (2,1) on the L=N periodic root domain, so +`k = (2*pi/L)*(2,1)` points 26.565° from +x. Mode (1,2) would point 63.435°. +A wrong-mode test catches a swap between the two. Polarization is +`E = (-1,2)/sqrt(5)`, perpendicular to k: + + u = U*E*sin(k.x + phi)*exp(-nu*|k|^2*t), phi = 0, p = const, rho0 = 1 + +Here `nu = U*N/100`, the Taylor–Green Re100 viscosity at the same N. The +initial condition is FP64 second-order equilibrium at the actual leaf +centers, with no non-equilibrium part. + +Physical caveats: +1. *Finite-Mach error.* The second-order D2Q9 equilibrium has third moment + `sum_i f_eq c_a c_b c_g = rho*cs^2*(u_a d_bg + u_b d_ag + u_g d_ab)`. The + `(c.u)^2` term contributes only even moments. The continuum target also + has `rho*u_a*u_b*u_g`. At Navier–Stokes (first Chapman–Enskog) order, + the resulting stress error is proportional to `d_g(rho u_a u_b u_g)`. + For this field with uniform density, `u.grad = 0` and `div u = 0`, so + that particular term vanishes. + + This does *not* show that all finite-Mach or cubic error vanishes. It + does not cover: + - higher Chapman–Enskog orders; + - lattice truncation for the oblique wavevector; + - FP32 storage; + - any density nonuniformity, which interfaces or anisotropy can create. +2. *Startup.* Equilibrium initialization omits the non-equilibrium stress. + The resulting kinetic transient lasts O(tau) steps and may shift the + in-phase amplitude. The decay rate is fitted from `ln(a_hat)` using + only samples with t* >= 0.12, which is about 1.9N root steps and much + longer than tau. This is intended to reduce startup sensitivity. It + has not been shown to cancel the bias, and that remains unverified + until measured. + +Metrics are computed on ACTIVE leaf cells only, each weighted by its area +`4**-level` in root-cell units. The partition must cover the domain exactly +once. One global pressure gauge is used, never one per level. + +- `a_hat`, `b_hat`: weighted least-squares fit of `u.E` onto the sin and cos + modes, solved with the 2x2 Gram matrix of the actual composite quadrature. +- `phase_drift = atan2(b_hat, a_hat)`, in radians. +- `rate_error = |rate_fit/(nu|k|^2) - 1|`, dimensionless. Rates are per root + step. +- `velocity_l2`: error relative to the exact field. +- `residual_l2`: the velocity left after subtracting the fitted mode, + relative to the exact-field norm. +- `longitudinal_l2`: the `u.k_hat` part of the velocity, relative to the + exact-field norm. +- `pressure_rel = p_rms/(A(t)^2/4)`: the gauged composite rms of + `p = (rho-1)/3`, divided by a *fixed analytic scale*. + - A²/4 is the continuum rms of the Taylor–Green exact pressure at the + same amplitude. That makes the number comparable in magnitude with TG + `pressure_l2`, but it is not that norm. TG normalizes by the discrete, + globally gauged exact-pressure norm on the composite grid. + - Exact pressure here is constant, so no relative-to-exact claim is made. + - Max and per-level values are reported for localization. + +Preregistered screens (`shear_wave.assess`): +- Run set: uniform N, uniform 2N, and refined depth-2 runs at root N, with + constant and/or projected_linear explosion. +- `eps(N) = (|k| h_root)^2 = 5*(2*pi/N)^2`, dimensionless, coefficient one. + It is 0.01205 at N=128. +- Absolute screens: rate_error, max |phase_drift| (radians) and max + residual_l2 must each be <= eps. This is a preregistered diagnostic + choice of the second-order truncation scale, *not* a proven error bound. +- Existing ceilings also apply: velocity_l2 <= 0.05, mass <= 1e-5 and + pressure_rel <= 0.02. +- Budget form, as in `assess.py`: + `refined <= max(2*control_N, floor)`, with floors 0.002 for rate and + residual and 0.02 for pressure. +- The ratio to the 2N control is reported, not gated. +- A control failing an absolute screen makes the set INVALID_CONTROL. The + screens are never retuned. + +Scope: one flow, one rectangle, SRT only. Two resolutions show a trend, +not an order. This does not establish placement robustness. diff --git a/tests/validation/taylor_green/run_case.py b/tests/validation/taylor_green/run_case.py index 2f4954c..df1efdd 100644 --- a/tests/validation/taylor_green/run_case.py +++ b/tests/validation/taylor_green/run_case.py @@ -20,6 +20,16 @@ REFERENCE_HARNESS = "11901a7f77557a901e017a171eeae1758d584fa5" INITIALIZATION = "physical FP64 equilibrium G at actual cell centers; rho=1+3p" +def flow_functions(flow="taylor_green"): + """(initial, metrics, limits, initialization); the default is the frozen TG oracle.""" + if flow == "taylor_green": + return initial, metrics, LIMITS, INITIALIZATION + if flow == "shear_wave": + import shear_wave + return shear_wave.initial, shear_wave.metrics, shear_wave.LIMITS, shear_wave.INITIALIZATION + raise ValueError(f"unknown flow {flow!r}") + + def sample_steps(n, early_pressure=False): """Root synchronization steps; early ladder matches the 707afe6 diagnostic.""" main = [round(tstar * n / U) for tstar in TIMES] @@ -151,7 +161,11 @@ def main(): p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant") p.add_argument("--early-pressure", action="store_true", help="also save 65 scalar pressure samples at t*=k*.001875, k=0..64") + p.add_argument("--flow", choices=["taylor_green", "shear_wave"], default="taylor_green") args = p.parse_args() + if args.flow == "shear_wave" and (args.re != 100 or args.depth == 3 or args.early_pressure): + p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure") + initial, metrics, limits, initialization = flow_functions(args.flow) if args.output.exists(): p.error("Output already exists; choose a new path to retain prior evidence") if args.explosion != "constant" and args.depth != 2: @@ -197,7 +211,7 @@ def main(): mach=U * np.sqrt(3), viscosity=cfg.viscosity, tau_by_level=[0.5 + 3 * c.viscosity for c in h.configs], - initialization=INITIALIZATION, + initialization=initialization, numerics=dict(lattice_model=cfg.lattice_model, collision=cfg.collision, store_precision=cfg.store_precision, ddf_shifting=cfg.ddf_shifting, streaming=cfg.streaming, les_enabled=cfg.les_enabled, @@ -205,13 +219,18 @@ def main(): sample_phase="postcollision G at root synchronization"), cuda_visible_devices=selected, device_id=0, - limits=LIMITS, + limits=limits, workload=workload(masks, round(TIMES[-1] * args.n / U)), rows=[], early_pressure=dict(enabled=args.early_pressure, rows=[], sampled_max_pressure_l2=None), completed=False, ) + if args.flow != "taylor_green": + import shear_wave + report.update(flow=args.flow, mode=list(shear_wave.MODE), + polarization=shear_wave.E.tolist(), angle_deg=shear_wave.ANGLE_DEG, + phase=shear_wave.PHI, eps=shear_wave.eps(args.n)) args.output.parent.mkdir(parents=True, exist_ok=True) def save(): temporary = args.output.with_suffix(args.output.suffix + ".tmp") @@ -301,7 +320,7 @@ def main(): completed=True, seconds=time.monotonic() - started, screening_pass=all( - row[key] <= limit for row in report["rows"] for key, limit in LIMITS.items() + row[key] <= limit for row in report["rows"] for key, limit in limits.items() ), ) save() diff --git a/tests/validation/taylor_green/shear_wave.py b/tests/validation/taylor_green/shear_wave.py new file mode 100644 index 0000000..4864223 --- /dev/null +++ b/tests/validation/taylor_green/shear_wave.py @@ -0,0 +1,197 @@ +"""CPU-only oblique shear-wave reference and composite-leaf diagnostics. + +Convention: integer mode MODE=(2,1) on the L x L periodic root domain, so the +wavevector k = (2*pi/L)*(2,1) points 26.565 deg from +x (mode (1,2) would be +63.435 deg). Velocity is polarized along E = (-1,2)/sqrt(5), perpendicular to k. +The frozen Taylor-Green oracle is reused for constants and geometry only. +""" + +import math + +import numpy as np +from oracle import CX, CY, W, U, TIMES, centers + +MODE = (2, 1) +PHI = 0.0 +E = np.array([-MODE[1], MODE[0]], dtype=float) / math.hypot(*MODE) +ANGLE_DEG = math.degrees(math.atan2(MODE[1], MODE[0])) +RE = 100 # Re = U L / nu, the same nu as Taylor-Green Re100 at equal N +FIT_FROM_TSTAR = 0.12 # late-time samples only; reduces startup sensitivity +INITIALIZATION = ( + "physical FP64 equilibrium G at actual cell centers; rho=1 (exact p const)" +) +# Per-row screens. phase_drift_abs, residual_l2 and the run-level rate screen +# use eps(n) and are preregistered diagnostic choices, not error bounds. +LIMITS = dict(velocity_l2=0.05, pressure_rel=0.02, mass_relative=1e-5) +FLOORS = dict(rate_error=0.002, residual_l2=0.002, pressure_rel=0.02) +BASELINE_FACTOR = 2.0 + + +def wavevector(length): + return 2 * np.pi / length * np.array(MODE, dtype=float) + + +def decay_rate(length, viscosity): + """Exact amplitude decay rate per root step: nu |k|^2.""" + k = wavevector(length) + return viscosity * float(k @ k) + + +def eps(n): + """(|k| h_root)^2 with h_root = 1: dimensionless, coefficient one.""" + k = wavevector(n) + return float(k @ k) + + +def exact(x, y, t, length, viscosity, amplitude=U): + k = wavevector(length) + a = amplitude * np.exp(-decay_rate(length, viscosity) * t) + s = np.sin(k[0] * x + k[1] * y + PHI) + p = np.zeros_like(s) + return a * E[0] * s, a * E[1] * s, p # kinematic pressure, constant + + +def equilibrium(rho, ux, uy): + # Same second-order D2Q9 equilibrium as oracle.initial. + cu = CX[:, None, None] * ux + CY[:, None, None] * uy + return ( + W[:, None, None] + * rho + * (1 + 3 * cu + 4.5 * cu * cu - 1.5 * (ux * ux + uy * uy)) + ) + + +def initial(shape, level, length, viscosity): + x, y = centers(shape, level) + ux, uy, p = exact(x, y, 0.0, length, viscosity) + return equilibrium(1 + 3 * p, ux, uy) + + +def _leaves(arrays, masks): + accum = [] + for level, (g, mask) in enumerate(zip(arrays, masks)): + if g.shape != (9, *mask.shape): + raise ValueError("population/mask shape mismatch") + values = g[:, mask].astype(float) + if not np.isfinite(values).all(): + raise FloatingPointError("nonfinite leaf populations") + rho = values.sum(axis=0) + if (rho <= 0).any(): + raise FloatingPointError("nonpositive leaf density") + x, y = centers(mask.shape, level) + area = np.full(rho.size, 0.25**level) + lev = np.full(rho.size, level) + accum.append((rho, CX @ values / rho, CY @ values / rho, x[mask], y[mask], area, lev)) + return [np.concatenate(v) for v in zip(*accum)] + + +def metrics(arrays, masks, t, length, viscosity, initial_mass=None): + rho, ux, uy, x, y, w, lev = _leaves(arrays, masks) + volume = w.sum() + if not np.isclose(volume, length**2, rtol=0, atol=1e-10): + raise ValueError("leaf partition does not cover domain once") + k = wavevector(length) + theta = k[0] * x + k[1] * y + PHI + s, c = np.sin(theta), np.cos(theta) + v = ux * E[0] + uy * E[1] + # Weighted least squares on the actual composite quadrature (no + # continuous-orthogonality assumption). + gram = np.array([[w @ (s * s), w @ (s * c)], [w @ (s * c), w @ (c * c)]]) + a_hat, b_hat = np.linalg.solve(gram, [w @ (v * s), w @ (v * c)]) + ex, ey, _ = exact(x, y, t, length, viscosity) + exact_norm = w @ (ex * ex + ey * ey) + fit = a_hat * s + b_hat * c + longitudinal = (ux * k[0] + uy * k[1]) / math.sqrt(k @ k) + amplitude = U * math.exp(-decay_rate(length, viscosity) * t) + scale = amplitude * amplitude / 4 # fixed pressure scale, see docs + p = (rho - 1) / 3 + p = p - w @ p / volume # one global gauge + mass = float(w @ rho) + return dict( + tstar=t * U / length, + mass=mass, + mass_relative=( + 0.0 if initial_mass is None else abs(mass - initial_mass) / initial_mass + ), + a_hat=float(a_hat), + b_hat=float(b_hat), + amplitude_exact=amplitude, + amplitude_ratio=float(a_hat / amplitude), + phase_drift=float(math.atan2(b_hat, a_hat)), + phase_drift_abs=abs(float(math.atan2(b_hat, a_hat))), + velocity_l2=float( + np.sqrt(w @ ((ux - ex) ** 2 + (uy - ey) ** 2) / exact_norm) + ), + residual_l2=float( + np.sqrt(w @ ((ux - fit * E[0]) ** 2 + (uy - fit * E[1]) ** 2) / exact_norm) + ), + longitudinal_l2=float(np.sqrt(w @ (longitudinal**2) / exact_norm)), + pressure_rel=float(np.sqrt(w @ (p * p) / volume) / scale), + pressure_max_rel=float(np.max(np.abs(p)) / scale), + pressure_rel_by_level=[ + float(np.sqrt(w[lev == i] @ (p[lev == i] ** 2) / w[lev == i].sum()) / scale) + for i in range(len(masks)) + ], + rho_min=float(rho.min()), + rho_max=float(rho.max()), + max_speed=float(np.max(np.hypot(ux, uy))), + leaf_cells=int(len(rho)), + ) + + +def fit_rate(rows, length, viscosity): + """Least-squares decay rate of ln(a_hat) over late samples.""" + late = [r for r in rows if r["tstar"] >= FIT_FROM_TSTAR - 1e-12] + if len(late) < 2 or any(r["a_hat"] <= 0 for r in late): + raise ValueError("need >= 2 late samples with positive amplitude") + t = np.array([r["step"] for r in late], dtype=float) + slope = np.polyfit(t, np.log([r["a_hat"] for r in late]), 1)[0] + exact_rate = decay_rate(length, viscosity) + return dict(rate=float(-slope), exact_rate=exact_rate, + rate_error=abs(-slope / exact_rate - 1)) + + +def summarize(record): + rows = record["rows"] + n = record["n"] + out = fit_rate(rows, n, record["viscosity"]) + out.update( + eps=eps(n), + complete=record["completed"] and len(rows) == len(TIMES), + residual_l2=max(r["residual_l2"] for r in rows), + phase_drift_abs=max(r["phase_drift_abs"] for r in rows), + pressure_rel=max(r["pressure_rel"] for r in rows), + velocity_l2=max(r["velocity_l2"] for r in rows), + mass_relative=max(r["mass_relative"] for r in rows), + ) + out["absolute_pass"] = out["complete"] and all( + out[key] <= out["eps"] for key in ("rate_error", "phase_drift_abs", "residual_l2") + ) and all(out[key] <= limit for key, limit in LIMITS.items()) + return out + + +def assess(base_control, fine_control, refined): + """Uniform root-spacing and finest-spacing controls versus refined runs.""" + n = base_control["n"] + if (base_control["depth"], fine_control["depth"], fine_control["n"]) != (1, 1, 2 * n): + raise ValueError("controls must be uniform N and uniform 2N") + if any(r["depth"] != 2 or r["n"] != n for r in refined): + raise ValueError("refined runs must be depth 2 at root N") + base, fine = summarize(base_control), summarize(fine_control) + report = dict(base_control=base, fine_control=fine, refined={}, status="PASS") + if not (base["absolute_pass"] and fine["absolute_pass"]): + report["status"] = "INVALID_CONTROL" + for r in refined: + s = summarize(r) + s["budget"] = { + key: max(BASELINE_FACTOR * base[key], floor) for key, floor in FLOORS.items() + } + s["budget_pass"] = all(s[key] <= s["budget"][key] for key in FLOORS) + s["fine_control_ratio"] = { + key: (s[key] / fine[key] if fine[key] > 0 else None) for key in FLOORS + } + s["pass"] = s["absolute_pass"] and s["budget_pass"] + report["refined"][r["explosion"]] = s + if report["status"] == "PASS" and not s["pass"]: + report["status"] = "FAIL" + return report diff --git a/tests/validation/taylor_green/test_shear_wave.py b/tests/validation/taylor_green/test_shear_wave.py new file mode 100644 index 0000000..bf4d3d0 --- /dev/null +++ b/tests/validation/taylor_green/test_shear_wave.py @@ -0,0 +1,247 @@ +"""Standalone CPU manufactured checks for the oblique shear wave. + +Run from this directory with PYTHONPATH=/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(r["pressure_rel"], 1e-10) + self.assertEqual(len(r["pressure_rel_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(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) + + def test_fit_ignores_startup_sample(self): + rows = [dict(tstar=0.0, step=0, a_hat=2 * U)] + rows += [dict(tstar=s, step=round(s * N / U), + 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) + + def test_wrong_mode_control(self): + # Swapping (2,1) for (1,2), i.e. 26.6 vs 63.4 deg, is not silently modal. + 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) + 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["pressure_rel"], 1e-10) # one global gauge removes offsets + 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["pressure_rel"], expected, delta=0.03 * expected) + + 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): + return dict(n=n, depth=depth, viscosity=U * n / sw.RE, rows=rows, completed=completed, + explosion=explosion) + + +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), + residual_l2=residual, phase_drift_abs=phase, pressure_rel=pressure, + 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"]["pressure_rel"], 0.03) + 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_invalid_control_and_wrong_shapes(self): + base = record(128, 1, rows_for(128, residual=0.5)) # > eps(128) + fine = record(256, 1, rows_for(256)) + refined = record(128, 2, rows_for(128)) + self.assertEqual(sw.assess(base, fine, [refined])["status"], "INVALID_CONTROL") + 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()