Compare commits

..
Author SHA1 Message Date
ClaudeandClaude Opus 5.5 512dcc923b step(n) docstring: raw-exchange ledger holds the final interval only
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-20 03:10:39 +00:00
ClaudeandClaude Opus 5.5 ad1d82cb7f FixedPatchStepper.step(n): wait on the stream once per call
Enqueue all n base intervals in order on the one stream and synchronize
once before returning (step(0): no launch, no wait). step_count, the raw
exchange ledger and an open window are published only after that wait;
any failure leaves step_count at its entry value and invalidates both
ledgers. Async device errors are no longer attributed to an interval;
the docstring says to use step(1) loops to localize one. Launch order and
buffers are unchanged. Adds a mocked-driver host test of the contract.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-20 03:06:45 +00:00
12 changed files with 234 additions and 2813 deletions
-300
View File
@@ -159,303 +159,3 @@ 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 the shear wave;
the default flow and its receipt are the unchanged Taylor–Green path.
### Convention
Integer mode (2,1) on the L=N periodic root domain:
- the wavevector `k = (2*pi/L)*(2,1)` points 26.565° from +x;
- mode (1,2) would point 63.435°, and a wrong-mode test catches a swap.
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
`nu = U*N/100`, the Taylor–Green Re100 viscosity at the same N. This is an
exact incompressible Navier–Stokes solution. The initial condition is FP64
second-order equilibrium at the actual leaf centers, with no
non-equilibrium part.
### Physical caveats
1. *Scope of the cubic-term argument.* The second-order D2Q9 equilibrium
lacks the cubic third moment `rho*u_a*u_b*u_g` of the Maxwellian (Dellar,
https://people.maths.ox.ac.uk/dellar/papers/Galilean.pdf, eq. 3). At
Navier–Stokes order this produces a stress error proportional to
`d_g(rho u_a u_b u_g)`.
For the ideal constant-density, zero-background transverse wave, k.E = 0
makes that contracted derivative vanish. That is the only claim made.
It does not remove other discrete or finite-Mach errors. It is not
guaranteed once AMR introduces density, longitudinal or higher-mode
contamination. The analytic NS field itself is exact.
2. *Startup.* Equilibrium initialization is not the discrete shear
eigenmode (Mei et al., https://doi.org/10.1016/j.compfluid.2005.08.008).
- The decay rate is a free-intercept least-squares fit of ln R over the
samples with t* >= 0.12, where `R = hypot(a_hat, b_hat)`.
- The intercept cancels only a constant amplitude prefactor, and only
once a single exponential dominates.
- Per-interval rates and the maximum log-fit residual are kept, so a
lingering transient stays visible.
- Late sampling alone does not show that transients have gone.
### Metrics
Computed on ACTIVE leaf cells, each weighted by its area `4**-level` in
root-cell units. The partition must cover the domain exactly once, and one
global pressure gauge is used.
- `a_hat`, `b_hat`: weighted least-squares fit of `u.E` onto sin/cos, using
the 2x2 Gram matrix of the actual composite quadrature. A wrong Fourier
mode need not project to exactly zero on this quadrature.
- `amplitude_error = R/A_exact - 1`.
- `phase_drift = atan2(b_hat, a_hat)`, in radians, reported separately.
- Rates are positive (minus the fitted slope).
- Per root step they are compared with `nu*|k|^2`. Uniform 2N uses its
own `nu = 0.0768` at N=256 and doubled step counts.
- In t* they are compared with `5*(2*pi)^2/Re = 1.97392088`.
- `rate_error = |rate/exact - 1|`.
- `velocity_l2`, `residual_l2` (velocity left after removing the fitted
mode) and `longitudinal_l2` (the `u.k_hat` part), all relative to the
exact-field norm. A longitudinal perturbation also raises the
full-vector residual and the velocity L2.
- `spurious_pressure_rms`: the raw gauged rms of `p = (rho-1)/3`.
- `spurious_pressure_scaled`: that rms divided by `A_exact(t)^2/4`, a
continuum reference scale. It is *not* the discrete, globally gauged
Taylor–Green norm on the composite grid. Exact pressure is constant, so
no relative-to-exact claim is made. Max and per-level values are kept.
### Preregistered screens (`shear_wave.assess`)
Run set: uniform N, uniform 2N, and refined depth-2 runs at root N.
- `eps(N) = (|k| h_root)^2 = 5*(2*pi/N)^2`: dimensionless, coefficient
one, 0.01205 at N=128. This is a chosen screen, not a proven error bound.
- Absolute screens:
- rate_error, max |phase_drift| (radians) and max residual_l2 must each
be <= eps;
- velocity_l2 <= 0.05, mass <= 1e-5, spurious_pressure_scaled <= 0.02.
- Budgets keep the `assess.py` form, `refined <= max(2*control_N, floor)`:
- floors are 0.002 for rate and residual, 0.02 for pressure;
- the comparative pressure budget is always redundant:
`max(2*control, 0.02)` is never below the 0.02 absolute ceiling;
- raw ratios to both the root and finest controls, including raw
pressure rms, are always reported.
- Every record is validated before any screen is applied:
- exact scheduled steps and t*, with no duplicates;
- every numeric field finite, including list fields and density;
- error quantities nonnegative and rho_min > 0;
- signed quantities (amplitudes, `amplitude_error`, `phase_drift`,
per-interval rate errors) keep their sign.
- Outcomes are kept distinct:
- INVALID_CONTROL: a control or refined run is missing, a control is
incomplete or invalid, or the records come from mixed sources;
- INVALID_RUN: a refined record fails validation;
- CONTROL_SCREEN_FAIL: complete, valid controls fail their own screens.
Each refined run's `own_screens_pass` is still reported, but
`controlled_pass` is false, so it is never a controlled-comparison pass;
- PASS/FAIL: every refined run's `controlled_pass` (own screens and
budgets against valid, passing controls).
- Finite modal extinction (late R = 0) is a physical failure, not
malformed data. The run gets `rate_fit_valid=false` with null fit
quantities, and its other metrics are kept. Nothing is clipped and no
rate is invented. The result is CONTROL_SCREEN_FAIL for a control, or
FAIL for a refined run.
- Screens are never retuned.
### Limitations
One flow, one phase and one rectangle, SRT only. The rectangle's vertical
faces span 0.625 tangential periods and its horizontal faces 1.5, so this
is not placement robustness. Two resolutions show a trend, not an order.
### Receipt assessment and runner integrity
- Fixed four-case N128 contract:
`python assess_shear_wave.py U128.json U256.json R128_constant.json R128_projected.json --output report.json --expected-revision <sha>`.
It reuses the shared provenance, build, workload and sample-time checks
from `assess_projected_pairs.py`, and adds shear metadata checks (flow,
mode, polarization, angle, phase, eps, limits, disabled early trace, and
pinned `shear_wave.py` hash).
- Missing, extra, duplicate or mixed-source records give `INVALID_SET`,
as does any failed shared check.
- Otherwise the result comes from `shear_wave.assess`, which separates
`INVALID_CONTROL` from `CONTROL_SCREEN_FAIL`.
- `run_case.verify_source(..., flow="shear_wave")` requires the imported
`shear_wave.py` to be the selected harness file with the pinned
`SHEAR_WAVE_SHA256`. The default TG flow does not import it.
- Sample and READY finiteness use `finite_row`, which accepts scalar and
list diagnostics. For TG's all-scalar rows it behaves the same as before.
- Merge note: native 92b8afe adds a `scaled_n256` receipt block. When the
two are merged, keep both blocks and make `--flow shear_wave` and
`scaled_n256` mutually exclusive. Not merged here.
### Scaled shear diagnostic (opt-in `scaled_shear_n256`)
Runner: `run_case.py --flow shear_wave --diagnostic scaled_shear_n256 --re 100`
with one of these four cases:
- `--n 256 --depth 1` (U256)
- `--n 512 --depth 1` (U512)
- `--n 256 --depth 2 --explosion constant` (constant R256)
- `--n 256 --depth 2 --explosion projected_linear` (projected R256)
Everything physical is unchanged from N128:
- `rectangles(256,2)` covers the same physical box;
- same t* schedule (16N root steps);
- `nu = U*N/100`;
- same `eps(N)`, limits and floors.
`scaled_n256` (Taylor–Green, from native 92b8afe) and `scaled_shear_n256`
each require their own flow, so they are mutually exclusive. N512 is
accepted only through one of the two diagnostics. The shear diagnostic
takes no early-pressure trace.
Assessment:
python assess_shear_wave.py U256.json U512.json R256_constant.json R256_projected.json \
--diagnostic scaled_shear_n256 --reference-files <four de5b0cc N128 receipts> \
--output report.json --expected-revision <sha>
- The scaled N256/512 case, build and workload checks come from the
92b8afe assessor. Its TG-only diagnostic identity and early-sampling
checks are not used. Instead every record must carry
`diagnostic=scaled_shear_n256`, and N128 records must carry none.
- The reference N128 receipts are hash-pinned to the audited de5b0cc
evidence (commit cba6bca). They are never relabelled as fresh controls.
- The report always includes, for each case:
- the raw signed rate bias `rate/exact - 1` and the interval errors;
- each refined run's signed excess over *each* uniform control (rate
bias, residual, raw pressure rms);
- the N128→N256 signed changes and ratios. The ratio is null when the
N128 value is zero.
- These are descriptive. Passing or failing another floor-controlled
screen is not the conclusion, and a two-resolution ratio is not a
convergence order.
- Portability: the receipt's recomputed `eps = |k|^2` is compared within
4 ULPs (`close_ulps`), not by exact equality. A cross-runtime replay
differed by one ULP: 0.012047856934923534 on Horsea vs
0.012047856934923532 elsewhere. Source and file hashes, case identity,
the other metadata and every scientific threshold stay exact.
### Nested shear preset (opt-in `nested_shear_n128`)
A fixed three-level, two-region preset, defined in `shear_nested.py`. The
wave, phase, Re100/U=0.03, initialization, five sample times, eps(root N),
limits and floors are the same as before.
- Refinements, given as parent-frame `(x0, y0, x1, y1)` rectangles on
level-wide, zero-origin unions:
- L0→L1: `(16,16,60,112)`, `(68,24,112,120)`
- L1→L2: `(44,44,108,212)`, `(148,60,212,228)`
- Leaves `[7936, 12288, 86016]`, hops `[2048, 4096, 8192]`, work
771227648 per refined case; 3723493376 for all four cases.
- Level ν `[.0384, .0768, .1536]`, τ `[.6152, .7304, .9608]`.
Cases (`run_case.py --flow shear_wave --diagnostic nested_shear_n128 --re 100`):
- `--n 128 --depth 1`
- `--n 512 --depth 1` — the finest control is 4N, not 2N
- `--n 128 --depth 3`
- `--n 128 --depth 3 --explosion projected_linear [--final-snapshot]`
`--final-snapshot` is allowed only for the projected depth-3 case. It writes
`<output>.final_macro.npz`, containing:
- ACTIVE leaves packed as `level`/`i`/`j`;
- `rho`, `ux`, `uy` from the FP64 physical populations at the final root
synchronization;
- embedded metadata: exact time, phase, mode and config.
It contains no populations, restart state or time series. Its SHA-256 and
size go into the receipt. A file over 8 MiB is removed and reported, never
kept. All old paths stay scalar-only.
Assessment: `assess_shear_wave.py U128.json U512.json R128c.json R128p.json
--diagnostic nested_shear_n128 --output report.json --expected-revision <sha>`.
- It reuses the case-agnostic provenance checks from the TG assessor, plus
exact nested scales, builds, patch, workload, metadata and the pinned
`shear_nested.py` hash.
- The physics is `shear_wave.assess` generalized only in its depth-3 / 4N
shape contract; the metrics are unchanged.
- `shear_wave.py` is byte-identical (81d5c8ee…), so historical receipt
identities and verdicts are untouched.
- This preset changes placement and refined volume as well as depth. It
gives no isolated cost of an extra level and no resolution trend.
### Depth-2 twin of the nested preset (opt-in `twin_shear_n128`)
`shear_twin.py` uses the same two L0→L1 boxes as `shear_nested.py`, with no
L2. Leaves are `[7936, 33792]`, hops `[2048, 4096]`, and work is 154664960 per
refined case (611319808 for all four).
The campaign is four fresh runs:
`run_case.py --flow shear_wave --diagnostic twin_shear_n128 --re 100` with
- `--n 128 --depth 1`
- `--n 256 --depth 1` (the finest control is 2N)
- `--n 128 --depth 2`
- `--n 128 --depth 2 --explosion projected_linear`
There is no snapshot option. Older exact controls serve only as consistency
checks; they are not reused as references.
Assessment uses the same preset contract as the nested preset, through
`assess_shear_wave.py ... --diagnostic twin_shear_n128`:
- the physics is the unchanged depth-2 / 2N `shear_wave.assess`;
- receipts must pin `shear_twin.py`, `shear_nested.py` and `shear_wave.py`;
- `shear_nested.py` stays byte-identical (fa6cfe87…), so the 335624c nested
receipts still validate.
Labels:
- *Twin vs nested* is the effect of adding this nested L2 refinement.
Removing L2 also changes the finest spacing and refined volume, so it does
not isolate abstract depth.
- *Twin vs the single-box depth-2 runs* is a placement/geometry contrast at
equal depth.
- Neither contrast is additive by construction.
- The thin-band/acoustic explanation for the nested L1 pressure is a
hypothesis, not established.
### Signed interface-probe traces (opt-in `--probe-trace`, mechanism diagnostic)
`shear_probes.py` (pinned) uses native's validated 24-target layout through
box 1, in root-128 physical coordinates:
- a horizontal line y = 56, crossing the vertical faces x = 16 and 22;
- a vertical line x = 38, crossing the horizontal faces y = 16 and 22;
- normal coordinates 12.5, 14.5, 15.5, 16.25, 17.25, 19.25, 20.75, 21.75,
22.125, 23.125, 24.125, 26.125.
Each target maps to the nearest ACTIVE leaf center. The deterministic tie
rule is: distance, then finer level, then smaller j, then smaller i. This
gives 24 unique cells in every case, and nested has 3 L0 + 5 L1 + 4 L2 per
transect. Maximum displacement from target to center is 0.625 root cells in
U128 (0.177 in U512).
Allowed cases:
- nested U128, nested U512, nested projected depth 3;
- twin projected depth 2.
The option can't be combined with `--final-snapshot`.
Sampling:
- The receipt stores each probe's leaf level, `i`/`j`, actual cell center
and displacement from its target, all in root-physical units (N=128).
- Grids differ between cases, so centers are *not* collocated. Analytic
fields are evaluated at each actual center (`shear_probes.residuals`), and
comparisons across cases are between residuals at nearby but distinct
centers.
- Times are root-physical steps 0–128 and 1920–2048. U512 is sampled every 4
lattice steps to match.
- Each sample stores `rho`, `ux`, `uy` from the FP64 physical populations,
POST-collision G at root synchronization only. No substep content is
observed, and grid-scale alternation can alias.
What's kept and enforced:
- The five qualification samples stay unchanged.
- No field dumps. The trace is validated for exact schedule, finiteness and
a 600 kB cap. If invalid it is dropped and an error is recorded, keeping
the scalar receipt.
- There are no new pass gates; this is a mechanism diagnostic only.
+21 -7
View File
@@ -822,27 +822,41 @@ class FixedPatchStepper:
grid=((pairs + self._block[0] - 1) // self._block[0], 1), **kw)
def step(self, n=1):
"""Advance G→G by n base intervals: level l takes 2**l hops per interval."""
"""Advance G→G by n base intervals: level l takes 2**l hops per interval.
All n intervals are enqueued in order on the one stream, and the host waits
once, before returning; step(0) launches and waits for nothing. On return
the device work is complete, step_count and an open window cover all n
intervals, and the raw-exchange ledger holds the final interval only.
Asynchronous device errors surface at that single wait, so they are not
attributed to the interval that caused them; call step(1) repeatedly to
localize one (same launches, same bits). On any failure the instance is
failed, the exchange ledger and window are invalid, and step_count keeps
its value from entry to the failing call.
"""
if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
raise ValueError("n must be a nonnegative integer")
with self._lock, self.ctx.activate():
self._require("ready")
if n == 0:
return
root = self._levels[0]
kw = dict(block=self._block, stream=self._stream)
try:
self._status = "advancing"
self._exchange_valid = False
for _ in range(n):
if self._observe_raw_exchange:
self._exchange_valid = False
if self._observe_raw_exchange: # stream-ordered: holds the last interval only
for b in self._level_boundaries:
if b is not None:
cuda.memset_d8_async(b.interval_integral, 0,
len(b.tables["cell"]) * 32, self._stream)
self._hop(0, root.current, root.spare, **kw)
self._stream.synchronize()
root.current, root.spare = root.spare, root.current
self._step_count += 1
self._exchange_valid = self._observe_raw_exchange
root.current, root.spare = root.spare, root.current # host handles only
self._stream.synchronize()
self._step_count += n
self._exchange_valid = self._observe_raw_exchange
self._status = "ready"
except BaseException:
self._window_valid = False
+170
View File
@@ -0,0 +1,170 @@
"""Host contract of FixedPatchStepper.step(n): one stream wait per call (mocked driver).
All n base intervals are enqueued in order on the backend stream and the host
waits once before step returns. Launch order and buffer arguments must equal n
step(1) calls; completion metadata, raw-exchange and window ledgers, and failure
invalidation are published only after that wait.
"""
from contextlib import nullcontext
from unittest.mock import Mock
import unittest
import numpy as np
from CelerisLab.config import LBMConfig
from tests.unit.test_raw_exchange import body
from tests.unit.test_refinement_three_level_multipatch import host_module
BODY_LEVEL, HOPS = 3, 8 # body on L3: 2**3 hops per base interval
def backend(window=False):
"""Four-level host backend; every launch, memset and stream wait goes to one log."""
host = host_module()
log, allocations = [], []
def allocate(size):
allocations.append(Mock(size=size))
return allocations[-1]
module = Mock()
module.get_global.return_value = (Mock(), 20)
module.get_function.side_effect = lambda name: Mock(
side_effect=lambda *a, **k: log.append((name, a, k)))
host.compiler.load_module.return_value = module
host.cuda.mem_alloc.side_effect = allocate
host.cuda.memcpy_dtoh.side_effect = lambda array, pointer: array.fill(0)
host.cuda.memset_d8_async.side_effect = lambda *a: log.append(("memset", a, {}))
run = host.FixedPatchStepper(
LBMConfig(nx=48, ny=48, viscosity=.02), context=Mock(activate=nullcontext, sm_arch="sm_70"),
refinements=[[(6, 6, 26, 26)], [(20, 20, 44, 44)], [(56, 56, 72, 72)]],
level_blocks=[(4, 4)] * 4, level_objects=[(), (), (), (body(),)],
observe_raw_exchange=True, observe_window_exchange=window)
run._stream.synchronize.side_effect = lambda: log.append(("sync", (), {}))
run.initialize(*[np.full(level.shape, 1 / 9, np.float32) for level in run._levels])
log.clear()
return run, log, allocations
def signature(log, allocations):
"""Event names with device-pointer arguments replaced by allocation index."""
index = {id(a): i for i, a in enumerate(allocations)}
return [(name, tuple(index.get(id(a), "-") for a in args)) for name, args, _ in log]
def names(log):
return [name for name, _, _ in log]
class StepSync(unittest.TestCase):
def test_step1_waits_once_after_all_launches(self):
run, log, _ = backend()
run.step(1)
self.assertEqual(names(log).count("sync"), 1)
self.assertEqual(names(log)[-1], "sync")
self.assertEqual(names(log)[0], "memset")
self.assertEqual(names(log).count("BlockStep"), 1 + 2 + 4 + 8)
self.assertEqual(run.step_count, 1)
run.close()
def test_step_n_equals_n_step1_launches_with_one_wait(self):
single, single_log, single_alloc = backend()
for _ in range(3):
single.step(1)
batched, batched_log, batched_alloc = backend()
batched.step(3)
expected = [e for e in signature(single_log, single_alloc) if e[0] != "sync"]
got = signature(batched_log, batched_alloc)
self.assertEqual(got[-1][0], "sync")
self.assertEqual(names(batched_log).count("sync"), 1)
self.assertEqual(got[:-1], expected) # same order, same ping-pong buffers
self.assertEqual(names(batched_log).count("memset"), 3)
self.assertEqual(batched.step_count, single.step_count)
index = [{id(a): i for i, a in enumerate(alloc)} for alloc in (single_alloc, batched_alloc)]
self.assertEqual(index[1][id(batched._levels[0].current)],
index[0][id(single._levels[0].current)])
single.close()
batched.close()
def test_step0_launches_and_waits_for_nothing(self):
run, log, _ = backend()
run.step(2)
log.clear()
run.step(0)
self.assertEqual(log, [])
self.assertEqual(run.step_count, 2)
value = run.read_raw_exchange(0, level=BODY_LEVEL) # step(0) keeps the ledger valid
self.assertEqual((value["time_start"], value["time_end"]), (1, 2))
run.close()
def test_completion_metadata_covers_all_intervals(self):
run, _, _ = backend()
with self.assertRaisesRegex(RuntimeError, "completed"):
run.read_raw_exchange(0, level=BODY_LEVEL)
run.step(3)
value = run.read_raw_exchange(0, level=BODY_LEVEL)
self.assertEqual((value["time_start"], value["time_end"], value["sample_count"]),
(2, 3, HOPS)) # step(n) retains only the final base interval
run.step(2)
self.assertEqual(run.read_raw_exchange(0, level=BODY_LEVEL)["time_end"], 5)
run.close()
def test_window_accounting_across_batched_calls(self):
windows = []
for calls in ((1, 3), (4,), (1, 1, 1, 1)):
run, log, _ = backend(window=True)
run.step(1)
run.start_exchange_window()
log.clear()
for n in calls:
run.step(n)
self.assertEqual(names(log).count("BlockAccumulateExchange"), 2 * 4 * HOPS)
self.assertEqual(names(log).count("sync"), len(calls))
out = run.read_exchange_window(0, level=BODY_LEVEL)
windows.append((out["time_start"], out["time_end"], out["sample_count"]))
run.close()
self.assertEqual(windows, [(1., 5, 4 * HOPS)] * 3)
def test_failed_wait_publishes_nothing(self):
run, _, _ = backend(window=True)
run.step(2)
run.start_exchange_window()
run._stream.synchronize.side_effect = RuntimeError("device failure")
with self.assertRaisesRegex(RuntimeError, "device failure"):
run.step(3)
self.assertEqual(run.step_count, 2) # value at entry to the failing call
self.assertFalse(run._exchange_valid)
self.assertFalse(run._window_valid)
for read in (lambda: run.read_raw_exchange(0, level=BODY_LEVEL),
lambda: run.read_exchange_window(0, level=BODY_LEVEL), run.download,
lambda: run.step(1)):
with self.assertRaisesRegex(RuntimeError, "failed"):
read()
run._stream.synchronize.side_effect = None
run.close()
def test_launch_error_mid_batch_skips_wait_and_fails(self):
run, log, _ = backend()
run.step(1)
log.clear()
hop, calls = run._hop, []
def failing(index, *args, **kw):
if index == 0:
calls.append(1)
if len(calls) == 2:
raise RuntimeError("launch failure")
return hop(index, *args, **kw)
run._hop = failing
with self.assertRaisesRegex(RuntimeError, "launch failure"):
run.step(3)
self.assertNotIn("sync", names(log))
self.assertEqual(run.step_count, 1)
self.assertFalse(run._exchange_valid)
with self.assertRaisesRegex(RuntimeError, "failed"):
run.read_raw_exchange(0, level=BODY_LEVEL)
run.close()
if __name__ == "__main__":
unittest.main()
@@ -1,7 +1,6 @@
"""New-source four-case repeat of 11901a7; frozen oracle and gates are unchanged."""
import argparse
import hashlib
import json
import math
from pathlib import Path
@@ -9,25 +8,16 @@ from pathlib import Path
from assess import BASELINE_FACTOR, FLOORS
from oracle import LIMITS, TIMES, U, rectangles
from run_case import (INITIALIZATION, sample_steps, NUMERICAL_REVISION, SOURCE_SHA256,
BUILD_CAPTURE_SHA256, REFERENCE_HARNESS, SCALED_DIAGNOSTIC)
BUILD_CAPTURE_SHA256, REFERENCE_HARNESS)
# Host-generated snapshots of a135252 at the frozen local configurations.
# Host-generated snapshots of a135252 at the frozen two local configurations.
# These describe expected inputs; actual sidecars are captured before cleanup.
BUILD_SOURCES = {
128: "803177416b7599ab89cc5afd35cf4a28a673b2f4bb8e8a6a1abc0b6d77610873",
256: "17ec5e60f87b840d3d13b8778bc8a3aecf485469c5a061af46c5cc19b73396cf",
512: "4e85e38098858ec95c137d1b1e96846a58d9cf3cc242b349a3ad7c8b8e3ce0aa",
}
CASES = ((100, 128, 1, "constant"), (100, 256, 1, "constant"),
(100, 128, 2, "constant"), (100, 128, 2, "projected_linear"))
REFERENCE_REVISION = "2cf9853952e52f59bea63cf5c1d3a1cd76e53eed"
# Audited 2cf9853 scalar receipts, archived at evidence commit 90849a55cad4a0da9b72747da7bfe6eb75962a44.
REFERENCE_RECEIPTS = {
"uniform128": "724efafc267eaa59dcc59e8d4c16177b67f6c6a08fd2732a481833237bb19cee",
"uniform256": "3c312a99ebf747a135674dd481b77455a3d2c9a406ca388f30524b37dad9bbe9",
"constant128": "2ee28efd1e81b2e28b7902b6c4f95bb2c1ea4e5f0ec484638fa5ae83dca34bba",
"projected128": "ee67b1d4d90a048c0d3d1a16447c6eab3d89c94d151f791306e373f96dd8b963",
}
NUMERICS = dict(lattice_model="D2Q9", collision="SRT", store_precision="FP32",
ddf_shifting=True, streaming="double_buffer", les_enabled=False,
boundary="periodic", block_shape=[8, 8],
@@ -83,11 +73,7 @@ def valid_builds(record, n, depth):
return True
def assess(records, expected_revision=None, *, diagnostic=None):
if diagnostic not in (None, SCALED_DIAGNOSTIC):
raise ValueError("Unknown diagnostic")
base_n = 256 if diagnostic else 128
cases = tuple((re, n * base_n // 128, depth, scheme) for re, n, depth, scheme in CASES)
def assess(records, expected_revision=None):
indexed, duplicates = {}, []
for record in records:
key = tuple(record.get(k) for k in ("re_L", "n", "depth", "explosion"))
@@ -95,11 +81,9 @@ def assess(records, expected_revision=None, *, diagnostic=None):
duplicates.append(str(key))
else:
indexed[key] = record
missing = [str(k) for k in cases if k not in indexed]
checks = dict(complete_matrix=not missing, no_extra_cases=set(indexed) <= set(cases),
missing = [str(k) for k in CASES if k not in indexed]
checks = dict(complete_matrix=not missing, no_extra_cases=set(indexed) <= set(CASES),
no_duplicate_cases=not duplicates, finite_receipts=finite_tree(records))
if diagnostic or any("diagnostic" in r for r in records):
checks["diagnostic_identity"] = all(r.get("diagnostic") == diagnostic for r in records)
sources = [r.get("source", {}) for r in records]
checks["same_source_and_harness"] = bool(sources) and len({
json.dumps(s, sort_keys=True) for s in sources}) == 1
@@ -122,11 +106,9 @@ def assess(records, expected_revision=None, *, diagnostic=None):
) and len(set(devices)) == 1
enabled = [r.get("early_pressure", {}).get("enabled") for r in records]
checks["same_early_sampling"] = bool(enabled) and all(type(v) is bool for v in enabled) and len(set(enabled)) == 1
if diagnostic:
checks["diagnostic_early_sampling"] = all(v is True for v in enabled)
shared_valid = all(checks.values())
complete, case_valid, early = {}, {}, []
for key in cases:
for key in CASES:
if key not in indexed:
continue
r = indexed[key]
@@ -153,7 +135,7 @@ def assess(records, expected_revision=None, *, diagnostic=None):
(("module", "__init__.py"), ("compiler", "cuda/compiler_v2.py"), ("kernel_root", "lbm/kernels")))
expected_patch = [] if depth == 1 else rectangles(n, depth)
checks[label + " patch"] = json.dumps(r.get("refinements")) == json.dumps(expected_patch)
cells = [n*n] if depth == 1 else [17*n*n//32, 15*n*n//8]
cells = [n*n] if depth == 1 else [8704, 30720]
hops = [16*n * 2**l for l in range(depth)]
work = r.get("workload", {})
checks[label + " workload"] = (
@@ -186,22 +168,19 @@ def assess(records, expected_revision=None, *, diagnostic=None):
comparisons, scheme_results = [], {}
for scheme in ("constant", "projected_linear"):
key = (100, base_n, 2, scheme)
key = (100, 128, 2, scheme)
for metric, floor in FLOORS.items():
amr, coarse, fine = (worst(k, metric) for k in (key, cases[0], cases[1]))
amr, coarse, fine = (worst(k, metric) for k in (key, CASES[0], CASES[1]))
budget = None if coarse is None else max(BASELINE_FACTOR * coarse, floor)
budget = budget if finite(budget) else None
ratio = None if amr is None or fine is None else amr / max(fine, 1e-30)
passed = amr is not None and budget is not None and amr <= budget
checks[f"baseline {scheme}/{metric}"] = passed
comparisons.append(dict(explosion=scheme, metric=metric, amr=amr,
baseline_budget=budget, passed=passed,
**{f"uniform{base_n}": coarse, f"uniform{2*base_n}": fine,
f"amr_to_uniform{2*base_n}": ratio if finite(ratio) else None}))
if diagnostic:
comparisons[-1]["floor_dominates"] = None if coarse is None else floor >= BASELINE_FACTOR * coarse
uniform128=coarse, baseline_budget=budget, passed=passed,
uniform256=fine, amr_to_uniform256=ratio if finite(ratio) else None))
baseline_passed = all(checks[f"baseline {scheme}/{metric}"] for metric in FLOORS)
valid = shared_valid and all(case_valid.get(k, False) for k in (cases[0], cases[1], key))
valid = shared_valid and all(case_valid.get(k, False) for k in (CASES[0], CASES[1], key))
scheme_results[scheme] = dict(passed=valid and baseline_passed,
valid_case_and_controls=valid, baseline_passed=baseline_passed)
return dict(passed=all(checks.values()), missing=missing, duplicates=duplicates, checks=checks,
@@ -211,89 +190,19 @@ def assess(records, expected_revision=None, *, diagnostic=None):
floors=FLOORS, baseline_factor=BASELINE_FACTOR,
main_statistic="maximum of four post-initial samples at t*=.12,.24,.36,.48",
early_statistic="sampled pressure maximum only, t*=k*.001875, k=0..64; not an all-step or continuous maximum",
finest_control=f"uniform{2*base_n} and AMR/ uniform{2*base_n} ratios are descriptive, not acceptance gates",
finest_control="uniform256 and AMR/ uniform256 ratios are descriptive, not acceptance gates",
scope="2D periodic Re100 smooth-decay diagnostic; no new algorithm, convergence-order or general AMR qualification")
def load_reference(paths):
"""Use the four audited historical files verbatim, never relabel them fresh controls."""
blobs = [path.read_bytes() for path in paths]
hashes = [hashlib.sha256(blob).hexdigest() for blob in blobs]
if len(hashes) != 4 or set(hashes) != set(REFERENCE_RECEIPTS.values()):
raise ValueError("Historical N128 receipt hashes differ from the audited 2cf9853 evidence")
return [json.loads(blob) for blob in blobs]
def resolution_trend(records, reference, expected_revision):
"""Signed changes are descriptive; they never replace the frozen acceptance budgets."""
old = assess(reference, REFERENCE_REVISION)
new = assess(records, expected_revision, diagnostic=SCALED_DIAGNOSTIC)
# Numerical screen failures still have useful raw errors; only malformed inputs block trends.
valid = all(v for result in (old, new) for k, v in result["checks"].items()
if not k.startswith("baseline ") and not k.endswith(" absolute"))
valid = valid and all(maximum(r["rows"], m) is not None for r in reference + records for m in LIMITS)
oracle = reference[0]["source"]["harness_sha256"].get("oracle.py") if valid else None
valid = valid and bool(oracle) and records[0]["source"]["harness_sha256"].get("oracle.py") == oracle
if not valid:
return dict(valid=False, per_time=[], worst_post_initial=[], uniform256_bridge=[])
key = lambda r: (r["re_L"], r["n"], r["depth"], r["explosion"])
before, after = {key(r): r for r in reference}, {key(r): r for r in records}
def change(role, metric, low, high, **extra):
ratio = high / low if low else None
percent = 100*(ratio-1) if ratio is not None else None
return dict(role=role, metric=metric, n128_campaign_error=low, n256_campaign_error=high,
signed_change=high-low, new_to_old=ratio if finite(ratio) else None,
ratio_unavailable_reason=None if finite(ratio) else "zero historical error" if not low else "overflow",
change_percent=percent if finite(percent) else None, **extra)
per_time, worst, bridge = [], [], []
for old_key, role in zip(CASES, ("uniform_root", "uniform_fine", "constant", "projected_linear")):
re, n, depth, scheme = old_key
low, high = before[old_key], after[(re, 2*n, depth, scheme)]
for metric in LIMITS:
per_time.extend(change(role, metric, a[metric], b[metric], tstar=a["tstar"])
for a, b in zip(low["rows"], high["rows"]))
worst.append(change(role, metric, maximum(low["rows"][1:], metric),
maximum(high["rows"][1:], metric)))
for a, b in zip(before[CASES[1]]["rows"], after[CASES[1]]["rows"]):
bridge.extend(dict(metric=m, tstar=a["tstar"], historical=a[m], fresh=b[m], signed_change=b[m]-a[m])
for m in LIMITS)
low, high = before[CASES[1]], after[CASES[1]]
bridge_identity = dict(main_rows_identical=low["rows"] == high["rows"],
early_pressure_rows_identical=low["early_pressure"]["rows"] == high["early_pressure"]["rows"])
same_runtime = reference[0]["runtime_identity"] == records[0]["runtime_identity"]
return dict(valid=True, reference_revision=REFERENCE_REVISION, reference_receipt_sha256=REFERENCE_RECEIPTS,
reference_scheme_results=old["scheme_results"], per_time=per_time, worst_post_initial=worst,
uniform256_bridge=bridge, uniform256_bridge_identity=bridge_identity,
same_recorded_runtime=same_runtime,
trend_interpretation_ready=all(bridge_identity.values()) and same_runtime,
scope="Historical same-production/oracle reference versus fresh controls; same-placement two-resolution trend only, not convergence order or phase robustness")
if __name__ == "__main__":
p = argparse.ArgumentParser(description=__doc__)
p.add_argument("files", type=Path, nargs="+")
p.add_argument("--output", type=Path, required=True)
p.add_argument("--expected-revision", required=True)
p.add_argument("--diagnostic", choices=[SCALED_DIAGNOSTIC])
p.add_argument("--reference-files", type=Path, nargs=4)
args = p.parse_args()
if args.output.exists():
p.error("Output already exists; choose a new path to retain prior evidence")
if bool(args.diagnostic) != bool(args.reference_files):
p.error("--diagnostic scaled_n256 requires --reference-files; references are diagnostic-only")
records = [json.loads(path.read_text()) for path in args.files]
result = assess(records, args.expected_revision, diagnostic=args.diagnostic)
if args.diagnostic:
try:
reference = load_reference(args.reference_files)
except ValueError as exc:
p.error(str(exc))
result["diagnostic"] = args.diagnostic
result["resolution_trend"] = resolution_trend(records, reference, args.expected_revision)
result["checks"]["valid_historical_reference_and_trend"] = result["resolution_trend"]["valid"]
result["passed"] = all(result["checks"].values())
result = assess([json.loads(path.read_text()) for path in args.files], args.expected_revision)
args.output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n")
print(json.dumps(result, indent=2, allow_nan=False))
raise SystemExit(0 if result["passed"] else 1)
@@ -1,336 +0,0 @@
"""Fixed four-case shear-wave contracts.
Default: N128 (U128, U256, constant R128, projected R128). Opt-in
--diagnostic scaled_shear_n256: N256 (U256, U512, constant R256, projected
R256), same physical rectangle, t* schedule, nu = U N / 100, eps(N), limits
and floors, plus a descriptive comparison with the hash-pinned N128 receipts.
Shared provenance, build, workload and sample-time checks are reused from
assess_projected_pairs (same cases, source pins and geometry); only the
Taylor-Green-specific checks are replaced by shear-wave metadata and physics.
"""
import argparse
import hashlib
import json
import math
from pathlib import Path
import assess_projected_pairs as tg
import shear_wave as sw
import shear_nested as nested
import shear_twin as twin
from run_case import (NESTED_SHEAR, SCALED_DIAGNOSTIC, SCALED_SHEAR, SHEAR_NESTED_SHA256,
SHEAR_TWIN_SHA256, SHEAR_WAVE_SHA256, TWIN_SHEAR, sample_steps)
# Audited de5b0cc N128 receipts (evidence commit cba6bca, shear-wave-de5b0cc zip).
REFERENCE_REVISION = "de5b0cc080794d0ce44d37f845aa8e8dbf8c2016"
REFERENCE_RECEIPTS = {
"uniform128": "d4a60fb806d04bfa8098a862cceba51e00b49979a6a250fef17a5fc8e3ea3255",
"uniform256": "e0b58e734724be5f82ee5d6df84395e3181b4e96e4246dd5a574ea23a8d12386",
"constant128": "c70d7cd14b1257a01cc49f8edfcc7685d430e5301e2f2d60ed426936f8b96c51",
"projected128": "980161458b9463af09a8e3178063fea0f4bbeb465329e1d3becdb0f224537b54",
}
ROLES = ("uniform_root", "uniform_fine", "constant", "projected_linear")
SHARED = ("complete_matrix", "no_extra_cases", "no_duplicate_cases", "finite_receipts",
"same_source_and_harness", "clean_a135252_numerical_source",
"same_recorded_runtime", "same_explicit_device", "same_early_sampling")
PER_CASE = ("complete_times", "scales", "actual_builds", "reference", "runtime_source",
"patch", "workload")
METADATA_ULPS = 4
def close_ulps(recorded, expected, ulps=METADATA_ULPS):
"""Portable check for recomputed float metadata only (e.g. eps = |k|^2).
The recorded value comes from another runtime/libm, so exact equality can
fail by one ULP (observed: 0.012047856934923534 Horsea vs ...532 here).
A few ULPs is far below any meaningful change; hashes, case identity and
scientific thresholds stay exact.
"""
return (isinstance(recorded, float) and math.isfinite(recorded)
and abs(recorded - expected) <= ulps * math.ulp(max(abs(recorded), abs(expected))))
def metadata_ok(record, n):
return (record.get("flow") == "shear_wave" and record.get("mode") == list(sw.MODE)
and record.get("polarization") == sw.E.tolist()
and record.get("angle_deg") == sw.ANGLE_DEG and record.get("phase") == sw.PHI
and close_ulps(record.get("eps"), sw.eps(n))
and record.get("numerics") == tg.NUMERICS
and record.get("initialization") == sw.INITIALIZATION
and record.get("limits") == sw.LIMITS
and record.get("early_pressure") == dict(enabled=False, rows=[],
sampled_max_pressure_l2=None)
and (record.get("source") or {}).get("shear_wave_sha256") == SHEAR_WAVE_SHA256)
def cases_for(scaled):
factor = 2 if scaled else 1
return tuple((re, n * factor, depth, scheme) for re, n, depth, scheme in tg.CASES)
def assess(records, expected_revision=None, *, scaled=False):
if not isinstance(records, list) or not all(isinstance(r, dict) for r in records):
raise ValueError("records must be a list of receipt objects")
# The scaled TG assessor supplies N256/512 cases, builds and workload; its
# TG-only diagnostic identity and early-sampling checks are not used here.
base_checks = tg.assess(records, expected_revision,
diagnostic=SCALED_DIAGNOSTIC if scaled else None)["checks"]
checks = {name: base_checks.get(name, False) for name in SHARED}
checks["shear_diagnostic_identity"] = all(
r.get("diagnostic") == (SCALED_SHEAR if scaled else None) for r in records)
cases = cases_for(scaled)
indexed = {}
for r in records:
indexed.setdefault(tuple(r.get(k) for k in ("re_L", "n", "depth", "explosion")), r)
for key in cases:
_, n, depth, scheme = key
label = f"{n}/{depth}/{scheme}"
for name in PER_CASE:
checks[f"{label} {name}"] = base_checks.get(f"{label} {name}", False)
checks[f"{label} shear_metadata"] = key in indexed and metadata_ok(indexed[key], n)
contract_valid = all(checks.values())
result = dict(checks=checks, contract_valid=contract_valid,
cases=[list(k) for k in cases], expected_revision=expected_revision,
diagnostic=SCALED_SHEAR if scaled else None,
scope="2D periodic Re100 oblique (2,1) shear wave; one phase and one "
"rectangle; no placement robustness or convergence order")
if not contract_valid:
result.update(status="INVALID_SET", passed=False)
return result
try:
physics = sw.assess(indexed[cases[0]], indexed[cases[1]],
[indexed[cases[2]], indexed[cases[3]]])
except (KeyError, ValueError) as error:
result.update(status="INVALID_SET", passed=False, reason=str(error))
return result
result.update(physics=physics, status=physics["status"],
passed=physics["status"] == "PASS")
try:
result["signed"] = signed_summary([indexed[k] for k in cases])
except (KeyError, TypeError, ValueError) as error:
result["signed_unavailable"] = str(error)
return result
CASE_AGNOSTIC = ("finite_receipts", "same_source_and_harness", "clean_a135252_numerical_source",
"same_recorded_runtime", "same_explicit_device", "same_early_sampling")
# Fixed multi-box presets sharing one contract; only geometry, controls, pins differ.
PRESETS = {
NESTED_SHEAR: dict(
module=nested, physics=nested.assess, snapshot=True,
pins={"shear_nested_sha256": SHEAR_NESTED_SHA256},
finest_control="uniform 4N (U512)",
scope="2D periodic Re100 oblique (2,1) shear wave; fixed three-level/two-region "
"placement; changes placement and refined volume as well as depth, so no "
"isolated-level cost or resolution trend"),
TWIN_SHEAR: dict(
module=twin, physics=sw.assess, snapshot=False,
pins={"shear_nested_sha256": SHEAR_NESTED_SHA256,
"shear_twin_sha256": SHEAR_TWIN_SHA256},
finest_control="uniform 2N (U256)",
scope="2D periodic Re100 oblique (2,1) shear wave; depth-2 twin of the nested "
"preset (same two L0->L1 boxes, no L2). Twin vs nested is the effect of adding "
"this nested refinement; it also changes finest spacing and refined volume, so it "
"does not isolate abstract depth"),
}
def _nested_case_ok(r, n, depth, diagnostic=NESTED_SHEAR):
"""Per-case identity of a fixed multi-box preset (scales, builds, patch, work, pins)."""
preset = PRESETS[diagnostic]
geometry = preset["module"]
runtime = r.get("runtime_source", {})
package = Path(runtime.get("checkout", "")) / "src/CelerisLab"
cells = [n * n] if depth == 1 else geometry.LEAVES
hops = [16 * n * 2**level for level in range(depth)]
work = r.get("workload", {})
snapshot = r.get("final_snapshot")
expected_patch = [] if depth == 1 else geometry.REFINEMENTS
return {
"complete_times": r.get("completed") is True and tg.at_steps(
r.get("rows", []), sample_steps(n)[0], n),
"scales": (r.get("U") == sw.U and tg.finite(r.get("viscosity"))
and abs(r["viscosity"] - sw.U * n / 100) < 1e-15
and r.get("tau_by_level") == [0.5 + 3 * sw.U * n / 100 * 2**l
for l in range(depth)]
and r.get("mach") == sw.U * math.sqrt(3)),
"actual_builds": tg.valid_builds(r, n, depth),
"reference": (r.get("numerical_revision") == tg.NUMERICAL_REVISION
and r.get("reference_harness") == tg.REFERENCE_HARNESS),
"runtime_source": bool(runtime.get("checkout")) and all(
runtime.get(k) == str(package / tail) for k, tail in
(("module", "__init__.py"), ("compiler", "cuda/compiler_v2.py"),
("kernel_root", "lbm/kernels"))),
"patch": json.dumps(r.get("refinements")) == json.dumps(expected_patch),
"workload": (work.get("active_cells_by_level") == cells and work.get("level_hops") == hops
and work.get("root_intervals") == 16 * n
and work.get("active_collision_updates_by_level") == [
c * h for c, h in zip(cells, hops)]
and work.get("active_collision_updates") == sum(
c * h for c, h in zip(cells, hops))),
"shear_metadata": metadata_ok(r, n),
"nested_module_pin": all((r.get("source") or {}).get(key) == pin
for key, pin in preset["pins"].items()),
# Optional, nested projected depth-3 only; bounded and hashed when present.
"snapshot": snapshot is None or (
preset["snapshot"] and depth == nested.DEPTH and r.get("explosion") == "projected_linear"
and snapshot.get("kept") is True and isinstance(snapshot.get("sha256"), str)
and len(snapshot["sha256"]) == 64
and 0 < snapshot.get("bytes", 0) <= nested.SNAPSHOT_MAX_BYTES
and snapshot.get("leaf_counts") == nested.LEAVES
and snapshot.get("step") == sample_steps(n)[0][-1]),
}
def assess_nested(records, expected_revision=None, diagnostic=NESTED_SHEAR):
"""Fixed four-case preset contract: nested_shear_n128 (U128, U512, constant/projected
R128 depth 3) or twin_shear_n128 (U128, U256, constant/projected R128 depth 2)."""
preset = PRESETS[diagnostic]
geometry = preset["module"]
if not isinstance(records, list) or not all(isinstance(r, dict) for r in records):
raise ValueError("records must be a list of receipt objects")
base_checks = tg.assess(records, expected_revision)["checks"]
checks = {name: base_checks.get(name, False) for name in CASE_AGNOSTIC}
indexed, duplicates = {}, []
for r in records:
key = tuple(r.get(k) for k in ("re_L", "n", "depth", "explosion"))
(duplicates.append(key) if key in indexed else indexed.__setitem__(key, r))
checks.update(complete_matrix=all(k in indexed for k in geometry.CASES),
no_extra_cases=set(indexed) <= set(geometry.CASES),
no_duplicate_cases=not duplicates,
shear_diagnostic_identity=all(r.get("diagnostic") == diagnostic
for r in records))
for key in geometry.CASES:
_, n, depth, scheme = key
label = f"{n}/{depth}/{scheme}"
case = _nested_case_ok(indexed[key], n, depth, diagnostic) if key in indexed else {}
for name in ("complete_times", "scales", "actual_builds", "reference", "runtime_source",
"patch", "workload", "shear_metadata", "nested_module_pin", "snapshot"):
checks[f"{label} {name}"] = bool(case.get(name, False))
contract_valid = all(checks.values())
result = dict(checks=checks, contract_valid=contract_valid, diagnostic=diagnostic,
cases=[list(k) for k in geometry.CASES], expected_revision=expected_revision,
finest_control=preset["finest_control"], scope=preset["scope"])
if not contract_valid:
result.update(status="INVALID_SET", passed=False)
return result
try:
cases = geometry.CASES
physics = preset["physics"](indexed[cases[0]], indexed[cases[1]],
[indexed[cases[2]], indexed[cases[3]]])
except (KeyError, ValueError) as error:
result.update(status="INVALID_SET", passed=False, reason=str(error))
return result
result.update(physics=physics, status=physics["status"], passed=physics["status"] == "PASS")
try:
# Role "uniform_fine" is the preset's finest control (4N nested, 2N twin).
result["signed"] = signed_summary([indexed[k] for k in geometry.CASES])
except (KeyError, TypeError, ValueError) as error:
result["signed_unavailable"] = str(error)
return result
def _signed(record):
"""Raw signed quantities of one valid record (None when the fit is unavailable)."""
s = sw.summarize(record)
bias = None if not s["rate_fit_valid"] else s["rate"] / s["exact_rate"] - 1
return dict(n=record["n"], depth=record["depth"], explosion=record["explosion"],
rate_fit_valid=s["rate_fit_valid"], signed_rate_bias=bias,
rate_tstar=s["rate_tstar"], exact_rate_tstar=s["exact_rate_tstar"],
interval_rate_errors=s["interval_rate_errors"],
residual_l2=s["residual_l2"], spurious_pressure_scaled=s["spurious_pressure_scaled"],
spurious_pressure_rms=s["spurious_pressure_rms"],
final_amplitude_error=record["rows"][-1]["amplitude_error"])
def signed_summary(records):
"""Signed bias per case and refined excess over each uniform control."""
rows = dict(zip(ROLES, (_signed(r) for r in records)))
def minus(a, b):
return None if a is None or b is None else a - b
for role in ("constant", "projected_linear"):
for control in ("uniform_root", "uniform_fine"):
rows[role][f"excess_over_{control}"] = {
key: minus(rows[role][key], rows[control][key])
for key in ("signed_rate_bias", "residual_l2", "spurious_pressure_rms")}
return rows
def load_reference(paths):
"""The four audited N128 receipts verbatim; never relabelled as fresh controls."""
blobs = [Path(path).read_bytes() for path in paths]
hashes = {hashlib.sha256(blob).hexdigest() for blob in blobs}
if len(blobs) != 4 or hashes != set(REFERENCE_RECEIPTS.values()):
raise ValueError("N128 receipt hashes differ from the audited de5b0cc evidence")
return [json.loads(blob) for blob in blobs]
def resolution_trend(new_result, reference, reference_revision=REFERENCE_REVISION):
"""Descriptive N128 -> N256 signed changes; never an acceptance gate or an order."""
old_result = assess(reference, reference_revision)
if "signed" not in old_result or "signed" not in new_result:
return dict(valid=False, reason="N128 reference or N256 set lacks a valid signed summary")
out = {}
for role in ROLES:
old, new = old_result["signed"][role], new_result["signed"][role]
out[role] = {}
for key in ("signed_rate_bias", "residual_l2", "spurious_pressure_rms",
*(f"excess_over_{c}" for c in ("uniform_root", "uniform_fine"))):
if key not in old:
continue
pairs = old[key].items() if isinstance(old[key], dict) else [(None, old[key])]
for sub, low in pairs:
high = new[key][sub] if sub else new[key]
ok = low is not None and high is not None
out[role][key if sub is None else f"{key}.{sub}"] = dict(
n128=low, n256=high, signed_change=high - low if ok else None,
n256_to_n128=high / low if ok and low else None)
return dict(valid=True, reference_revision=reference_revision,
reference_receipt_sha256=REFERENCE_RECEIPTS,
reference_status=old_result["status"], changes=out,
scope="same-placement two-resolution trend; not a convergence order, "
"placement robustness or redesign decision")
if __name__ == "__main__":
p = argparse.ArgumentParser(description=__doc__)
p.add_argument("files", type=Path, nargs=4)
p.add_argument("--output", type=Path, required=True)
p.add_argument("--expected-revision", required=True)
p.add_argument("--diagnostic", choices=[SCALED_SHEAR, NESTED_SHEAR, TWIN_SHEAR])
p.add_argument("--reference-files", type=Path, nargs=4)
args = p.parse_args()
if args.output.exists():
p.error("Output already exists; choose a new path to retain prior evidence")
if args.diagnostic in PRESETS:
if args.reference_files:
p.error(f"{args.diagnostic} takes no reference files (no resolution trend)")
result = assess_nested([json.loads(path.read_text()) for path in args.files],
args.expected_revision, args.diagnostic)
text = json.dumps(result, indent=2, allow_nan=False) + "\n"
args.output.write_text(text)
print(text, end="")
raise SystemExit(0 if result["passed"] else 1)
if bool(args.diagnostic) != bool(args.reference_files):
p.error("--diagnostic scaled_shear_n256 requires --reference-files (N128 receipts)")
result = assess([json.loads(path.read_text()) for path in args.files], args.expected_revision,
scaled=bool(args.diagnostic))
if args.diagnostic:
try:
reference = load_reference(args.reference_files)
except ValueError as exc:
p.error(str(exc))
result["resolution_trend"] = resolution_trend(result, reference)
text = json.dumps(result, indent=2, allow_nan=False) + "\n"
args.output.write_text(text)
print(text, end="")
raise SystemExit(0 if result["passed"] else 1)
+18 -165
View File
@@ -18,37 +18,6 @@ SOURCE_SHA256 = "a5dbfd804befce868067d79a23e5fd47dce7f505b3f9b1bb2145fb57cb5f217
BUILD_CAPTURE_SHA256 = "3dbaa6547cb0a4cb1541054f3b64ad8b6e18cc70527ac6166e56b944e9bd464a"
REFERENCE_HARNESS = "11901a7f77557a901e017a171eeae1758d584fa5"
INITIALIZATION = "physical FP64 equilibrium G at actual cell centers; rho=1+3p"
SHEAR_WAVE_SHA256 = "81d5c8ee60864372887cddf9239580c82fb6750d956038f1e3ee9de61aa4e877"
SCALED_DIAGNOSTIC = "scaled_n256"
SCALED_SHEAR = "scaled_shear_n256" # mutually exclusive with scaled_n256 (one flow each)
NESTED_SHEAR = "nested_shear_n128" # fixed three-level/two-region shear preset
SHEAR_NESTED_SHA256 = "fa6cfe87fb4df8c6d732c637c885276833c81973b950c7724d734f8fc350bcb0"
TWIN_SHEAR = "twin_shear_n128" # depth-2 twin: nested L0->L1 boxes without L2
SHEAR_TWIN_SHA256 = "1cbf37f18e89aea577581631c4565f51e68d22a3739b698816ef64acbf2b2bb4"
SHEAR_PROBES_SHA256 = "53ba8af8ac9c8078c40f84488c7ade4c11ccc91dd1fc275ff806042e9af30166" # signed interface-probe mechanism diagnostic
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 finite_row(row):
"""Every scalar or list-valued diagnostic is finite (TG rows are all scalars)."""
return all(np.isfinite(np.asarray(value, dtype=float)).all() for value in row.values())
def sample(run, masks, step, n, viscosity, initial_mass, metrics, depth, label="sample"):
arrays = (run.download(),) if depth == 1 else run.download()
row = metrics(arrays, masks, step, n, viscosity, initial_mass)
if not finite_row(row):
raise FloatingPointError(f"nonfinite {label} diagnostics")
return row
def sample_steps(n, early_pressure=False):
@@ -95,7 +64,7 @@ def provenance(repo):
)
def verify_source(repo, revision, flow="taylor_green", probes=False):
def verify_source(repo, revision):
"""Exact clean harness commit, fixed production tree, and reviewed capture helper."""
import oracle
from tests.validation.multigrid import check_projected_pairs as gate
@@ -114,34 +83,6 @@ def verify_source(repo, revision, flow="taylor_green", probes=False):
if result["source_tree_sha256"] != SOURCE_SHA256:
raise RuntimeError("Wrong numerical source fingerprint")
result["build_capture_sha256"] = BUILD_CAPTURE_SHA256
if probes:
import shear_probes
module = harness / "shear_probes.py"
if (Path(shear_probes.__file__).resolve() != module
or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_PROBES_SHA256):
raise RuntimeError("Wrong shear_probes module checkout or content")
result["shear_probes_sha256"] = SHEAR_PROBES_SHA256
if flow == "shear_twin":
import shear_twin
module = harness / "shear_twin.py"
if (Path(shear_twin.__file__).resolve() != module
or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_TWIN_SHA256):
raise RuntimeError("Wrong shear_twin module checkout or content")
result["shear_twin_sha256"] = SHEAR_TWIN_SHA256
if flow in ("shear_nested", "shear_twin"): # the twin reuses the nested geometry
import shear_nested
module = harness / "shear_nested.py"
if (Path(shear_nested.__file__).resolve() != module
or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_NESTED_SHA256):
raise RuntimeError("Wrong shear_nested module checkout or content")
result["shear_nested_sha256"] = SHEAR_NESTED_SHA256
if flow in ("shear_wave", "shear_nested", "shear_twin"):
import shear_wave
module = harness / "shear_wave.py"
if (Path(shear_wave.__file__).resolve() != module
or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_WAVE_SHA256):
raise RuntimeError("Wrong shear_wave module checkout or content")
result["shear_wave_sha256"] = SHEAR_WAVE_SHA256
return result
@@ -204,79 +145,26 @@ def main():
p.add_argument("--run-gpu", action="store_true")
p.add_argument("--dry-run", action="store_true")
p.add_argument("--re", type=int, choices=[100, 1000], required=True)
p.add_argument("--n", type=int, choices=[32, 64, 128, 256, 512], required=True)
p.add_argument("--n", type=int, choices=[32, 64, 128, 256], required=True)
p.add_argument("--depth", type=int, choices=[1, 2, 3], required=True)
p.add_argument("--output", type=Path, required=True)
p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant")
p.add_argument("--diagnostic",
choices=[SCALED_DIAGNOSTIC, SCALED_SHEAR, NESTED_SHEAR, TWIN_SHEAR],
help="explicit four-case N256/fine512 same-placement diagnostic "
"(scaled_n256: Taylor-Green; scaled_shear_n256: shear wave) or the "
"fixed three-level/two-region shear preset (nested_shear_n128) "
"or its depth-2 twin (twin_shear_n128)")
p.add_argument("--probe-trace", action="store_true",
help="signed interface-probe trace (nested U128/U512/projected d3, "
"twin projected d2); mechanism diagnostic, no gates")
p.add_argument("--final-snapshot", action="store_true",
help="nested_shear_n128 projected depth-3 only: final macro snapshot")
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()
nested = args.diagnostic == NESTED_SHEAR
if args.flow == "shear_wave" and (args.re != 100 or (args.depth == 3 and not nested)
or args.early_pressure):
p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure")
probe_cases = {(NESTED_SHEAR, 128, 1, "constant"), (NESTED_SHEAR, 512, 1, "constant"),
(NESTED_SHEAR, 128, 3, "projected_linear"),
(TWIN_SHEAR, 128, 2, "projected_linear")}
if args.probe_trace and (args.final_snapshot or (
args.diagnostic, args.n, args.depth, args.explosion) not in probe_cases):
p.error("--probe-trace requires nested U128/U512/projected depth 3 or twin projected "
"depth 2, without --final-snapshot")
if args.final_snapshot and not (nested and args.depth == 3
and args.explosion == "projected_linear"):
p.error("--final-snapshot requires nested_shear_n128 projected depth 3")
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 and not nested:
if args.explosion != "constant" and args.depth != 2:
p.error("projected_linear requires depth 2")
twin = args.diagnostic == TWIN_SHEAR
if nested or twin:
import shear_nested
import shear_twin
preset = shear_nested if nested else shear_twin
if args.flow != "shear_wave":
p.error(f"--diagnostic {args.diagnostic} requires --flow shear_wave")
if (args.re, args.n, args.depth, args.explosion) not in preset.CASES:
p.error(f"{args.diagnostic} requires one of its four Re100 cases")
elif args.diagnostic:
cases = {(256, 1, "constant"), (512, 1, "constant"),
(256, 2, "constant"), (256, 2, "projected_linear")}
flow = "shear_wave" if args.diagnostic == SCALED_SHEAR else "taylor_green"
if args.flow != flow:
p.error(f"--diagnostic {args.diagnostic} requires --flow {flow}")
if args.diagnostic == SCALED_DIAGNOSTIC and (
args.re != 100 or (args.n, args.depth, args.explosion) not in cases
or not args.early_pressure):
p.error("scaled_n256 requires one of its four Re100 cases and --early-pressure")
if args.diagnostic == SCALED_SHEAR and (args.n, args.depth, args.explosion) not in cases:
p.error("scaled_shear_n256 requires one of its four Re100 cases")
elif args.n == 512:
p.error("N512 requires --diagnostic scaled_n256, scaled_shear_n256 or nested_shear_n128")
selected = os.environ.get("CUDA_VISIBLE_DEVICES", "")
if not args.dry_run and (not args.run_gpu or not selected or "," in selected or selected == "-1"):
p.error("--run-gpu and exactly one explicit CUDA_VISIBLE_DEVICES are required")
repo = args.repo.resolve()
sys.path[:0] = [str(repo / "src"), str(repo)]
source = verify_source(repo, args.expected_revision,
"shear_nested" if nested else "shear_twin" if twin else args.flow,
probes=args.probe_trace)
source = verify_source(repo, args.expected_revision)
runtime = runtime_source(repo)
if args.dry_run:
print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True,
**({"diagnostic": args.diagnostic} if args.diagnostic else {}))))
print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True)))
return 0
from CelerisLab.cuda.context import CudaContext
from CelerisLab.lbm.blocks import BlockLayout2D
@@ -286,12 +174,7 @@ def main():
from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition
cfg = config(args.n, args.re)
if args.depth == 1:
refinements = []
elif nested or twin:
refinements = [list(level) for level in preset.REFINEMENTS]
else:
refinements = rectangles(args.n, args.depth)
refinements = [] if args.depth == 1 else rectangles(args.n, args.depth)
h = build_hierarchy(cfg, refinements)
np.testing.assert_array_equal(leaf_partition(h), np.ones((args.n, args.n)))
masks = [k == ACTIVE for k in h.kinds]
@@ -314,7 +197,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,
@@ -322,20 +205,13 @@ 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))
if args.diagnostic:
report["diagnostic"] = args.diagnostic
args.output.parent.mkdir(parents=True, exist_ok=True)
def save():
temporary = args.output.with_suffix(args.output.suffix + ".tmp")
@@ -370,37 +246,17 @@ def main():
]
)
initial_mass = None
probe_steps = set()
if args.probe_trace:
import shear_probes
probe_root = dict(zip(shear_probes.case_steps(args.n), shear_probes.STEPS))
probe_steps = set(probe_root)
probes = shear_probes.select(masks, args.n)
report["probe_trace"] = shear_probes.new_trace(
probes, cfg.viscosity, [args.diagnostic, args.n, args.depth, args.explosion])
for step in sorted(set(main_steps + early_steps) | probe_steps):
for step in sorted(set(main_steps + early_steps)):
run.step(step - run.step_count)
if step in probe_steps:
trace = report["probe_trace"]
trace["steps"].append(probe_root[step]) # root-physical step
trace["values"].append(shear_probes.sample(
(run.download(),) if args.depth == 1 else run.download(), probes))
if step not in main_steps and step not in early_steps:
continue # probe-only step: no qualification metrics
row = sample(run, masks, step, args.n, cfg.viscosity, initial_mass,
metrics, args.depth)
arrays = (run.download(),) if args.depth == 1 else run.download()
row = metrics(arrays, masks, step, args.n, cfg.viscosity, initial_mass)
if not all(np.isfinite(value) for value in row.values()):
raise FloatingPointError("nonfinite sample diagnostics")
row["step"] = step
if initial_mass is None:
initial_mass = row["mass"]
if step in main_steps:
report["rows"].append(row)
if args.final_snapshot and step == main_steps[-1]:
# Optional macro-only snapshot; scalar receipt stays authoritative.
report["final_snapshot"] = shear_nested.write_snapshot(
args.output.with_name(args.output.stem + ".final_macro.npz"),
run.download(), masks,
shear_nested.snapshot_meta(step, cfg.viscosity, args.explosion,
source.get("revision")))
if step in early_steps:
trace = report["early_pressure"]
trace["rows"].append({key: row[key] for key in ("step", "tstar", "pressure_l2")})
@@ -422,8 +278,10 @@ def main():
# Exercise real kernels/halos/transfers on every level before READY.
# This is the first measured trajectory step, not discarded warmup.
run.step(1)
sample(run, masks, 1, args.n, cfg.viscosity, initial_mass,
metrics, args.depth, "READY")
probe = (run.download(),) if args.depth == 1 else run.download()
probe_metrics = metrics(probe, masks, 1, args.n, cfg.viscosity, initial_mass)
if not all(np.isfinite(value) for value in probe_metrics.values()):
raise FloatingPointError("nonfinite READY diagnostics")
report["ready_after_step"] = 1
save()
print(
@@ -439,16 +297,11 @@ def main():
),
flush=True,
)
if args.probe_trace:
try:
report["probe_trace"]["bytes"] = shear_probes.validate(report["probe_trace"])
except ValueError as error: # keep the scalar receipt; drop an invalid trace
report["probe_trace"] = dict(error=str(error))
report.update(
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()
@@ -1,174 +0,0 @@
"""Fixed three-level/two-region shear preset (opt-in nested_shear_n128).
Four fresh cases: U128, U512 (finest control = 4N), constant R128 depth 3,
projected R128 depth 3. Same wave, phase, Re100/U=0.03, initialization, five
sample times, eps(root N), limits and floors as shear_wave. This changes
placement and refined volume as well as depth: no isolated-level cost or
resolution-trend label. shear_wave.py is reused unchanged (historical receipt
identities stay valid); only the depth/finest-control contract is generalized
here.
"""
import hashlib
import json
from pathlib import Path
import numpy as np
import shear_wave as sw
from oracle import CX, CY, U, centers
ROOT_N = 128
DEPTH = 3
FINEST_FACTOR = 2 ** (DEPTH - 1) # finest-spacing uniform control is 4N
# Parent-frame (x0, y0, x1, y1) rectangles; level-wide, zero-origin unions.
REFINEMENTS = [
[(16, 16, 60, 112), (68, 24, 112, 120)], # L0 -> L1, root coordinates
[(44, 44, 108, 212), (148, 60, 212, 228)], # L1 -> L2, L1 coordinates
]
CASES = ((100, ROOT_N, 1, "constant"), (100, FINEST_FACTOR * ROOT_N, 1, "constant"),
(100, ROOT_N, DEPTH, "constant"), (100, ROOT_N, DEPTH, "projected_linear"))
ROOT_STEPS = 16 * ROOT_N
LEAVES = [7936, 12288, 86016]
HOPS = [ROOT_STEPS * 2**level for level in range(DEPTH)] # [2048, 4096, 8192]
WORK = sum(c * h for c, h in zip(LEAVES, HOPS)) # 771227648 per refined case
NU = [U * ROOT_N / sw.RE * 2**level for level in range(DEPTH)] # [.0384, .0768, .1536]
TAU = [0.5 + 3 * nu for nu in NU] # [.6152, .7304, .9608]
SNAPSHOT_MAX_BYTES = 8 * 1024 * 1024
SNAPSHOT_KEYS = ("level", "i", "j", "rho", "ux", "uy", "meta")
def leaf_masks():
"""Independent numpy leaf masks of the fixed preset (root, L1, L2 level-wide grids)."""
masks, covered = [], None
for level in range(DEPTH):
size = ROOT_N * 2**level
inside = np.ones((size, size), bool) if level == 0 else covered
if level < DEPTH - 1:
child = np.zeros((size, size), bool)
for x0, y0, x1, y1 in REFINEMENTS[level]:
child[y0:y1, x0:x1] = True
masks.append(inside & ~child)
covered = np.kron(child, np.ones((2, 2), bool))
else:
masks.append(inside)
return masks
def assess(base_control, finest_control, refined):
"""shear_wave.assess generalized to depth 3 with a 4N finest control.
Same summaries, screens, budgets and outcome semantics; only the
(depth, finest-control) shape contract differs.
"""
if base_control is None or finest_control is None or not refined:
return dict(status="INVALID_CONTROL", reason="missing control or refined run")
n = base_control["n"]
if (base_control["depth"], finest_control["depth"], finest_control["n"]) != (
1, 1, FINEST_FACTOR * n):
raise ValueError("controls must be uniform N and uniform 4N")
if any(r["depth"] != DEPTH or r["n"] != n for r in refined):
raise ValueError("refined runs must be depth 3 at root N")
if len({r.get("explosion") for r in refined}) != len(refined):
raise ValueError("duplicate refined explosion scheme")
sources = {repr(r.get("source")) for r in [base_control, finest_control, *refined]}
base, fine = sw._usable(base_control), sw._usable(finest_control)
if base is None or fine is None or len(sources) != 1:
return dict(status="INVALID_CONTROL", reason="incomplete/invalid control or mixed source")
controls_pass = base["own_screens_pass"] and fine["own_screens_pass"]
report = dict(base_control=base, finest_control=fine, refined={},
controls_pass=controls_pass)
statuses = []
for r in refined:
try:
s = sw.summarize(r)
except (KeyError, TypeError, ValueError) as error:
report["refined"][r["explosion"]] = dict(valid=False, reason=str(error))
statuses.append("INVALID_RUN")
continue
s["valid"] = True
s["budget"] = {
key: (None if base[key] is None else max(sw.BASELINE_FACTOR * base[key], floor))
for key, floor in sw.FLOORS.items()
}
s["budget_pass"] = all(
s[key] is not None and s["budget"][key] is not None and s[key] <= s["budget"][key]
for key in sw.FLOORS)
s["control_ratios"] = {
name: {key: (s[key] / c[key] if s[key] is not None and c[key] else None)
for key in (*sw.FLOORS, "spurious_pressure_rms")}
for name, c in (("root", base), ("finest", fine))
}
s["controlled_pass"] = controls_pass and s["own_screens_pass"] and s["budget_pass"]
report["refined"][r["explosion"]] = s
statuses.append("PASS" if s["controlled_pass"] else "FAIL")
if "INVALID_RUN" in statuses:
report["status"] = "INVALID_RUN"
elif not controls_pass:
report["status"] = "CONTROL_SCREEN_FAIL"
else:
report["status"] = "PASS" if all(v == "PASS" for v in statuses) else "FAIL"
return report
def write_snapshot(path, arrays, masks, meta):
"""Final macro snapshot: packed ACTIVE leaves (level, i, j), rho, ux, uy only.
No populations, no restart state, no time series. Returns the receipt block;
a file over SNAPSHOT_MAX_BYTES is removed and reported, never kept.
"""
path = Path(path)
if path.exists():
raise FileExistsError(path)
cols = {k: [] for k in ("level", "i", "j", "rho", "ux", "uy")}
counts = []
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)
rho = values.sum(axis=0)
j, i = np.nonzero(mask) # row = y index, column = x index
cols["level"].append(np.full(rho.size, level, np.uint8))
cols["i"].append(i.astype(np.uint16))
cols["j"].append(j.astype(np.uint16))
cols["rho"].append(rho)
cols["ux"].append(CX @ values / rho)
cols["uy"].append(CY @ values / rho)
counts.append(int(rho.size))
packed = {k: np.concatenate(v) for k, v in cols.items()}
if not all(np.isfinite(packed[k]).all() for k in ("rho", "ux", "uy")):
raise FloatingPointError("nonfinite snapshot macros")
meta = dict(meta, leaf_counts=counts, cell_center="h*(index+0.5), h=2**-level",
fields="rho, ux, uy from FP64 physical populations at ACTIVE leaves")
with path.open("xb") as handle:
np.savez_compressed(handle, meta=np.array(json.dumps(meta, sort_keys=True)), **packed)
size = path.stat().st_size
block = dict(file=path.name, bytes=size, leaf_counts=counts, max_bytes=SNAPSHOT_MAX_BYTES,
step=meta.get("step"), tstar=meta.get("tstar"))
if size > SNAPSHOT_MAX_BYTES:
path.unlink()
return dict(block, kept=False, sha256=None, reason="exceeds snapshot byte cap")
block.update(kept=True, sha256=hashlib.sha256(path.read_bytes()).hexdigest())
return block
def read_snapshot(path):
"""Validate a snapshot: exact keys, no populations, consistent packed arrays."""
with np.load(path, allow_pickle=False) as data:
if tuple(sorted(data.files)) != tuple(sorted(SNAPSHOT_KEYS)):
raise ValueError("unexpected snapshot contents")
out = {k: data[k] for k in data.files}
meta = json.loads(str(out.pop("meta")))
n = len(out["rho"])
if any(len(v) != n for v in out.values()) or sum(meta["leaf_counts"]) != n:
raise ValueError("inconsistent snapshot arrays")
return meta, out
def snapshot_meta(step, viscosity, explosion, source_revision):
"""Exact time/phase/config identity stored inside the snapshot."""
return dict(step=step, tstar=step * U / ROOT_N, phase=sw.PHI, mode=list(sw.MODE),
polarization=sw.E.tolist(), amplitude=U, re_L=sw.RE, n=ROOT_N, depth=DEPTH,
explosion=explosion, viscosity=viscosity, tau_by_level=TAU,
refinements=REFINEMENTS, source_revision=source_revision,
origin="level-wide zero-origin grids", time="root synchronization, postcollision")
@@ -1,123 +0,0 @@
"""Signed interface-probe traces for the nested/twin shear presets (mechanism diagnostic).
No pass gates. Native-validated 24-target layout through box 1, in root-128
physical coordinates, crossing L0->L1 and (nested only) L1->L2 away from corners:
horizontal line y = 56 across the vertical faces x = 16 and x = 22;
vertical line x = 38 across the horizontal faces y = 16 and y = 22;
normal coordinates NORMAL_COORDS on both (nested: 3 L0 + 5 L1 + 4 L2 each).
Each target maps to the nearest ACTIVE leaf center (tie rule: distance, then
finer level, then smaller j, then smaller i). Grids differ between cases
(U128, U512, nested, twin), so actual centers, levels and displacements are
stored and analytic fields are evaluated at the actual centers.
Samples are POST-collision G at root synchronization only (root-step units);
no substep content is observed, and grid-scale alternation can alias.
"""
import math
import numpy as np
import shear_wave as sw
from oracle import CX, CY, U
ROOT_N = 128
NORMAL_COORDS = (12.5, 14.5, 15.5, 16.25, 17.25, 19.25, 20.75, 21.75,
22.125, 23.125, 24.125, 26.125)
TRANSECTS = {"horizontal_y56": dict(axis=0, other=56.0), # x varies, y fixed
"vertical_x38": dict(axis=1, other=38.0)} # y varies, x fixed
STEPS = tuple(range(0, 129)) + tuple(range(1920, 2049)) # root steps, both inclusive
TRACE_MAX_BYTES = 600_000
def points():
out = []
for name, t in TRANSECTS.items():
for along in NORMAL_COORDS:
xy = (along, t["other"]) if t["axis"] == 0 else (t["other"], along)
out.append((name, along, xy))
return out
def case_steps(n):
"""Case lattice steps at the fixed root-physical sample times (U512: 4 per root step)."""
return [s * (n // ROOT_N) for s in STEPS]
def select(masks, n=ROOT_N):
"""Nearest ACTIVE leaf center for each target (deterministic tie rule).
Targets and stored centers are in root-physical units (root cell = 1); a case
with root N = n has cell size h = ROOT_N / (n * 2**level) there (U512: 0.25).
"""
if n % ROOT_N:
raise ValueError("case root N must be a multiple of 128")
probes = []
for name, along, (x, y) in points():
best = None
for level, mask in enumerate(masks):
h = ROOT_N / (n * 2 ** level)
i0, j0 = int(math.floor(x / h)), int(math.floor(y / h))
for j in range(j0 - 2, j0 + 3):
for i in range(i0 - 2, i0 + 3):
if 0 <= j < mask.shape[0] and 0 <= i < mask.shape[1] and mask[j, i]:
cx, cy = h * (i + 0.5), h * (j + 0.5)
key = (math.hypot(cx - x, cy - y), -level, j, i)
if best is None or key < best[0]:
best = (key, level, i, j, cx, cy)
if best is None:
raise ValueError(f"no ACTIVE leaf near target {(x, y)}")
_, level, i, j, cx, cy = best
probes.append(dict(transect=name, normal=along, point=[x, y], level=level, i=i, j=j,
center=[cx, cy], displacement=[cx - x, cy - y]))
if len({(p["level"], p["i"], p["j"]) for p in probes}) != len(probes):
raise ValueError("probe targets map to duplicate cells")
return probes
def sample(arrays, probes):
"""Flat [rho, ux, uy] per probe from FP64 physical populations."""
row = []
for p in probes:
f = np.asarray(arrays[p["level"]][:, p["j"], p["i"]], dtype=float)
rho = float(f.sum())
row += [rho, float(CX @ f) / rho, float(CY @ f) / rho]
return row
def new_trace(probes, case_viscosity, case):
"""Times/coordinates are root-physical (N=128); nu_root = U*128/Re for all cases."""
return dict(schema=1, kind="mechanism diagnostic, no pass gates", case=case,
time_unit="root step (N=128)", sample_phase="postcollision G at root synchronization",
root_n=ROOT_N, viscosity_root=U * ROOT_N / sw.RE, case_viscosity=case_viscosity,
mode=list(sw.MODE), phase=sw.PHI,
amplitude=U, fields=["rho", "ux", "uy"], probes=probes, steps=[], values=[])
def residuals(trace):
"""Signed (drho, dux, duy) per step and probe, analytic field at each actual center."""
out = []
for step, row in zip(trace["steps"], trace["values"]):
res = []
for n, p in enumerate(trace["probes"]):
x, y = p["center"]
ux, uy, _ = sw.exact(x, y, step, trace["root_n"], trace["viscosity_root"])
rho, vx, vy = row[3 * n: 3 * n + 3]
res.append((rho - 1.0, vx - float(ux), vy - float(uy)))
out.append(res)
return out
def validate(trace):
"""Exact schedule, finite values, fixed probe count and byte cap."""
import json
if trace.get("steps") != list(STEPS) or len(trace["values"]) != len(STEPS):
raise ValueError("trace steps differ from the fixed schedule")
if len(trace["probes"]) != len(NORMAL_COORDS) * len(TRANSECTS):
raise ValueError("unexpected probe count")
width = 3 * len(trace["probes"])
if any(len(r) != width or not all(math.isfinite(v) for v in r) for r in trace["values"]):
raise ValueError("malformed or nonfinite trace row")
size = len(json.dumps(trace, allow_nan=False).encode())
if size > TRACE_MAX_BYTES:
raise ValueError(f"trace exceeds {TRACE_MAX_BYTES} bytes ({size})")
return size
@@ -1,36 +0,0 @@
"""Depth-2 twin of the nested shear preset (opt-in twin_shear_n128).
Same two L0->L1 boxes as shear_nested, no L2. Four fresh cases: U128, U256
(finest control = 2N), constant R128 depth 2, projected R128 depth 2; same
wave, phase, Re100/U=0.03, initialization, times, eps(root N), limits, floors.
Twin vs nested measures the effect of adding this nested L2 refinement.
Removing L2 also changes the finest spacing and refined volume, so the
contrast does not isolate abstract depth. Physics is the unchanged depth-2,
2N-control shear_wave.assess; shear_nested.py stays byte-identical.
"""
import numpy as np
import shear_nested as nested
ROOT_N = nested.ROOT_N
DEPTH = 2
FINEST_FACTOR = 2
REFINEMENTS = [nested.REFINEMENTS[0]] # L0 -> L1 only, root coordinates
CASES = ((100, ROOT_N, 1, "constant"), (100, FINEST_FACTOR * ROOT_N, 1, "constant"),
(100, ROOT_N, DEPTH, "constant"), (100, ROOT_N, DEPTH, "projected_linear"))
ROOT_STEPS = 16 * ROOT_N
LEAVES = [7936, 33792]
HOPS = [ROOT_STEPS * 2**level for level in range(DEPTH)] # [2048, 4096]
WORK = sum(c * h for c, h in zip(LEAVES, HOPS)) # 154664960 per refined case
NU = nested.NU[:DEPTH]
TAU = nested.TAU[:DEPTH]
def leaf_masks():
"""Independent numpy leaf masks (root complement, full L1 boxes)."""
child = np.zeros((ROOT_N, ROOT_N), bool)
for x0, y0, x1, y1 in REFINEMENTS[0]:
child[y0:y1, x0:x1] = True
return [~child, np.kron(child, np.ones((2, 2), bool))]
-312
View File
@@ -1,312 +0,0 @@
"""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, spurious_pressure_scaled=0.02, mass_relative=1e-5)
# The comparative pressure budget max(2*control, 0.02) is never below the 0.02
# absolute ceiling, so it is always redundant; raw control ratios are kept.
FLOORS = dict(rate_error=0.002, residual_l2=0.002, spurious_pressure_scaled=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 # continuum reference scale, not the TG norm
modal = math.hypot(a_hat, b_hat)
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),
modal_amplitude=modal,
amplitude_error=modal / amplitude - 1,
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)),
spurious_pressure_rms=float(np.sqrt(w @ (p * p) / volume)),
spurious_pressure_scaled=float(np.sqrt(w @ (p * p) / volume) / scale),
spurious_pressure_max_scaled=float(np.max(np.abs(p)) / scale),
spurious_pressure_scaled_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):
"""Free-intercept least squares of ln R, R = hypot(a_hat, b_hat), late samples.
The intercept absorbs only a constant amplitude prefactor once a single
exponential dominates; per-interval rates and the fit residual are kept so
a lingering transient stays visible. Rates are per root step of the record
(uniform 2N uses its own nu and doubled step counts).
"""
late = [r for r in rows if r["tstar"] >= FIT_FROM_TSTAR - 1e-12]
amp = np.array([math.hypot(r["a_hat"], r["b_hat"]) for r in late])
if len(late) < 3 or not np.isfinite(amp).all():
raise ValueError("need >= 3 late samples with finite R")
exact_rate = decay_rate(length, viscosity)
to_tstar = length / U # root steps per unit t*
if not (amp > 0).all():
# Finite modal extinction: a physical failure, not malformed data. No
# clipping and no invented rate; fit quantities are unavailable.
return dict(rate_fit_valid=False, rate=None, exact_rate=exact_rate,
rate_tstar=None, exact_rate_tstar=exact_rate * to_tstar,
rate_error=None, interval_rate_errors=None, log_fit_residual_max=None)
t = np.array([r["step"] for r in late], dtype=float)
(slope, intercept) = np.polyfit(t, np.log(amp), 1)
residual = np.log(amp) - (slope * t + intercept)
return dict(
rate_fit_valid=True, rate=float(-slope), exact_rate=exact_rate,
rate_tstar=float(-slope * to_tstar),
exact_rate_tstar=exact_rate * to_tstar, # = 5 (2 pi)^2 / Re
rate_error=abs(-slope / exact_rate - 1),
interval_rate_errors=[
float(-np.log(amp[i + 1] / amp[i]) / (t[i + 1] - t[i]) / exact_rate - 1)
for i in range(len(t) - 1)
],
log_fit_residual_max=float(np.max(np.abs(residual))),
)
NONNEGATIVE = ("tstar", "mass", "mass_relative", "amplitude_exact", "modal_amplitude",
"phase_drift_abs", "velocity_l2", "residual_l2", "longitudinal_l2",
"spurious_pressure_rms", "spurious_pressure_scaled",
"spurious_pressure_max_scaled", "rho_min", "rho_max", "max_speed")
SIGNED = ("a_hat", "b_hat", "amplitude_ratio", "amplitude_error", "phase_drift")
def _number(value):
return (isinstance(value, (int, float)) and not isinstance(value, bool)
and math.isfinite(value))
def validate_rows(record):
"""Exact schedule, finite numbers (incl. lists) and nonnegative error fields.
Signed quantities (amplitudes, amplitude/phase errors) keep their sign.
Raises ValueError so callers can report an invalid run before any screen.
"""
n, rows = record["n"], record["rows"]
steps = [round(tstar * n / U) for tstar in TIMES]
if not isinstance(rows, list) or len(rows) != len(steps):
raise ValueError("rows do not match the scheduled samples")
for row, step in zip(rows, steps):
if type(row.get("step")) is not int or row["step"] != step:
raise ValueError(f"unscheduled or duplicate step {row.get('step')!r}")
if not _number(row.get("tstar")) or abs(row["tstar"] - step * U / n) > 3e-17:
raise ValueError(f"wrong sample time at step {step}")
if not all(_number(row.get(key)) for key in (*NONNEGATIVE, *SIGNED)):
raise ValueError(f"missing or nonfinite diagnostic at step {step}")
if any(row[key] < 0 for key in NONNEGATIVE) or row["rho_min"] <= 0:
raise ValueError(f"negative error quantity or density at step {step}")
levels = row.get("spurious_pressure_scaled_by_level")
if (not isinstance(levels, list) or len(levels) != record["depth"]
or not all(_number(v) and v >= 0 for v in levels)):
raise ValueError(f"invalid per-level pressure at step {step}")
def summarize(record):
validate_rows(record)
rows = record["rows"]
n = record["n"]
out = fit_rate(rows, n, record["viscosity"])
out.update(
eps=eps(n),
complete=record["completed"] is True,
residual_l2=max(r["residual_l2"] for r in rows),
phase_drift_abs=max(r["phase_drift_abs"] for r in rows),
spurious_pressure_scaled=max(r["spurious_pressure_scaled"] for r in rows),
spurious_pressure_rms=max(r["spurious_pressure_rms"] for r in rows),
amplitude_error_abs=max(abs(r["amplitude_error"]) for r in rows),
velocity_l2=max(r["velocity_l2"] for r in rows),
mass_relative=max(r["mass_relative"] for r in rows),
)
if out["rate_fit_valid"] and not all(
_number(v) for v in (out["rate"], out["rate_error"], *out["interval_rate_errors"],
out["log_fit_residual_max"])):
raise ValueError("nonfinite rate fit")
out["own_screens_pass"] = out["complete"] and out["rate_fit_valid"] 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 _usable(record):
try:
out = summarize(record)
except (KeyError, TypeError, ValueError):
return None
return out if out["complete"] else None
def assess(base_control, fine_control, refined):
"""Uniform root-spacing and finest-spacing controls versus refined runs.
INVALID_CONTROL: missing, incomplete, invalid-row or mixed-source controls.
INVALID_RUN: a refined record has invalid rows (schedule, NaN, sign).
CONTROL_SCREEN_FAIL: complete valid controls that fail their own screens;
refined own-screen results are kept, but no controlled comparison passes.
PASS/FAIL: every refined run's controlled_pass (own screens and budgets).
"""
if base_control is None or fine_control is None or not refined:
return dict(status="INVALID_CONTROL", reason="missing control or refined run")
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")
if len({r.get("explosion") for r in refined}) != len(refined):
raise ValueError("duplicate refined explosion scheme")
sources = {repr(r.get("source")) for r in [base_control, fine_control, *refined]}
base, fine = _usable(base_control), _usable(fine_control)
if base is None or fine is None or len(sources) != 1:
return dict(status="INVALID_CONTROL",
reason="incomplete/invalid control or mixed source")
controls_pass = base["own_screens_pass"] and fine["own_screens_pass"]
report = dict(base_control=base, fine_control=fine, refined={},
controls_pass=controls_pass)
statuses = []
for r in refined:
try:
s = summarize(r)
except (KeyError, TypeError, ValueError) as error:
report["refined"][r["explosion"]] = dict(valid=False, reason=str(error))
statuses.append("INVALID_RUN")
continue
s["valid"] = True
# rate_error is None when a modal fit is unavailable; never a pass.
s["budget"] = {
key: (None if base[key] is None else max(BASELINE_FACTOR * base[key], floor))
for key, floor in FLOORS.items()
}
s["budget_pass"] = all(
s[key] is not None and s["budget"][key] is not None and s[key] <= s["budget"][key]
for key in FLOORS)
s["control_ratios"] = {
name: {key: (s[key] / c[key] if s[key] is not None and c[key] else None)
for key in (*FLOORS, "spurious_pressure_rms")}
for name, c in (("root", base), ("finest", fine))
}
s["controlled_pass"] = controls_pass and s["own_screens_pass"] and s["budget_pass"]
report["refined"][r["explosion"]] = s
statuses.append("PASS" if s["controlled_pass"] else "FAIL")
if "INVALID_RUN" in statuses:
report["status"] = "INVALID_RUN"
elif not controls_pass:
report["status"] = "CONTROL_SCREEN_FAIL"
else:
report["status"] = "PASS" if all(v == "PASS" for v in statuses) else "FAIL"
return report
@@ -17,7 +17,6 @@ from unittest.mock import patch
import numpy as np
import run_case
import assess_projected_pairs as assessor
from assess_projected_pairs import (BUILD_SOURCES, CASES, NUMERICS, NUMERICAL_REVISION,
SOURCE_SHA256, assess)
from oracle import LIMITS, TIMES, U, rectangles
@@ -32,7 +31,7 @@ RUNTIME = dict(device=dict(name="host mock", pci_bus_id="mock", capability=[7, 0
def source():
return dict(revision=MOCK_REVISION, dirty=False, source_tree_sha256=SOURCE_SHA256,
harness_sha256={"run_case.py": "a" * 64, "oracle.py": "c" * 64},
harness_sha256={"run_case.py": "a" * 64},
build_capture_sha256=run_case.BUILD_CAPTURE_SHA256)
@@ -43,14 +42,12 @@ def build_metadata(size):
arch="sm_70", flags=["--use_fast_math"], ptx_sha256="b" * 64)
def records(early=False, diagnostic=None):
def records(early=False):
from CelerisLab.lbm.hierarchy import build_hierarchy
from CelerisLab.lbm.refinement import ACTIVE
out = []
for re, n, depth, scheme in CASES:
if diagnostic:
n *= 2
cfg = run_case.config(n, re)
refinements = [] if depth == 1 else rectangles(n, depth)
h = build_hierarchy(cfg, refinements)
@@ -69,8 +66,6 @@ def records(early=False, diagnostic=None):
early_pressure=dict(enabled=early, rows=[
dict(step=s, tstar=s*U/n, pressure_l2=.01 if s not in main_steps else 0.0) for s in early_steps],
sampled_max_pressure_l2=.01 if early else None)))
if diagnostic:
out[-1]["diagnostic"] = diagnostic
return out
@@ -256,7 +251,7 @@ class ProjectedPairsTests(unittest.TestCase):
def test_pinned_build_sources_match_host_generated_inputs(self):
from CelerisLab.cuda import compiler_v2
for size in (128, 256, 512):
for size in (128, 256):
cfg = replace(run_case.config(size, 100), nx=10, ny=10)
cfg.validate()
root = compiler_v2.generate_config(cfg)
@@ -268,111 +263,6 @@ class ProjectedPairsTests(unittest.TestCase):
finally:
compiler_v2.cleanup_build_dir(root)
def test_scaled_geometry_work_and_frozen_gates(self):
data = records(True, run_case.SCALED_DIAGNOSTIC)
self.assertEqual(sum(r["workload"]["active_collision_updates"] for r in data), 4714397696)
self.assertEqual([r["viscosity"] for r in data], [.0768, .1536, .0768, .0768])
self.assertEqual(data[-1]["tau_by_level"], [.7303999999999999, .9608])
for r in data:
self.assertEqual([row["tstar"] for row in r["rows"]], TIMES)
self.assertEqual(len(r["early_pressure"]["rows"]), 65)
if r["depth"] == 2:
self.assertEqual(r["refinements"], [[(32, 32, 224, 192)]])
self.assertEqual(r["workload"]["active_cells_by_level"], [34816, 122880])
self.assertEqual(r["workload"]["level_hops"], [4096, 8192])
result = assess(data, MOCK_REVISION, diagnostic=run_case.SCALED_DIAGNOSTIC)
self.assertTrue(result["passed"])
self.assertEqual(result["floors"], dict(velocity_l2=.002, pressure_l2=.02, energy_relative=.004))
self.assertTrue(all(row["floor_dominates"] for row in result["comparisons"]))
data[-1]["rows"][-1]["velocity_l2"] = .00200001
self.assertFalse(assess(data, diagnostic=run_case.SCALED_DIAGNOSTIC)["passed"])
self.assertFalse(assess(data)["passed"])
data[-1]["rows"][-1]["velocity_l2"] = 0.0
for mutate in (lambda r: r.pop("diagnostic"), lambda r: r.update(diagnostic="wrong"),
lambda r: r["source"].update(revision="wrong"),
lambda r: r["early_pressure"].update(enabled=False),
lambda r: r.update(refinements=[[(32, 32, 224, 200)]])):
bad = deepcopy(data)
mutate(bad[-1])
self.assertFalse(assess(bad, MOCK_REVISION, diagnostic=run_case.SCALED_DIAGNOSTIC)["passed"])
def test_raw_resolution_changes_and_reference_identity(self):
old, new = records(True), records(True, run_case.SCALED_DIAGNOSTIC)
for r in old:
r["source"]["revision"] = assessor.REFERENCE_REVISION
old[-2]["rows"][-1]["velocity_l2"] = .004 # historical constant failure remains usable
old[-1]["rows"][-1]["velocity_l2"] = .001
new[-1]["rows"][-1]["velocity_l2"] = .0015 # floor pass, but worse raw error
result = assessor.resolution_trend(new, old, MOCK_REVISION)
self.assertTrue(result["valid"])
self.assertTrue(result["trend_interpretation_ready"])
self.assertFalse(result["reference_scheme_results"]["constant"]["passed"])
self.assertEqual((len(result["per_time"]), len(result["worst_post_initial"]),
len(result["uniform256_bridge"])), (80, 16, 20))
row = next(r for r in result["per_time"] if r["role"] == "projected_linear"
and r["metric"] == "velocity_l2" and r["tstar"] == .48)
self.assertEqual((row["signed_change"], row["new_to_old"], row["change_percent"]), (.0005, 1.5, 50.))
self.assertTrue(assess(new, diagnostic=run_case.SCALED_DIAGNOSTIC)["passed"])
self.assertIsNone(result["per_time"][0]["new_to_old"]) # zero denominator is explicit
self.assertEqual(result["per_time"][0]["ratio_unavailable_reason"], "zero historical error")
new[-1]["rows"][-1]["velocity_l2"] = .06 # retain raw errors even when absolute screen fails
self.assertTrue(assessor.resolution_trend(new, old, MOCK_REVISION)["valid"])
self.assertFalse(assess(new, diagnostic=run_case.SCALED_DIAGNOSTIC)["passed"])
new[-1]["rows"][-1]["velocity_l2"] = .0015
for sample in ("rows", "early_pressure"):
bad = deepcopy(new)
if sample == "rows":
bad[0]["rows"][-1]["velocity_l2"] = .0001
else:
bad[0]["early_pressure"]["rows"][1]["pressure_l2"] = .009
result = assessor.resolution_trend(bad, old, MOCK_REVISION)
self.assertTrue(result["valid"])
self.assertFalse(result["trend_interpretation_ready"])
for mutate in (lambda r: r["source"].update(revision=MOCK_REVISION),
lambda r: r["source"].update(source_tree_sha256="wrong"),
lambda r: r["source"]["harness_sha256"].update(**{"oracle.py": "wrong"}),
lambda r: r.update(diagnostic=run_case.SCALED_DIAGNOSTIC),
lambda r: r.update(refinements=[]), lambda r: r["rows"].pop()):
bad = deepcopy(old)
mutate(bad[-1])
self.assertFalse(assessor.resolution_trend(new, bad, MOCK_REVISION)["valid"])
with tempfile.TemporaryDirectory() as tmp:
paths = [Path(tmp) / f"{i}.json" for i in range(4)]
for path, receipt in zip(paths, old):
path.write_text(json.dumps(receipt))
with self.assertRaisesRegex(ValueError, "receipt hashes"):
assessor.load_reference(paths)
witnesses = {p.name: hashlib.sha256(p.read_bytes()).hexdigest() for p in paths}
with patch.object(assessor, "REFERENCE_RECEIPTS", witnesses):
self.assertEqual(assessor.load_reference(paths), json.loads(json.dumps(old)))
with self.assertRaises(ValueError):
assessor.load_reference(paths[:-1] + paths[:1])
paths[-1].write_text(paths[-1].read_text() + " ")
with self.assertRaises(ValueError):
assessor.load_reference(paths)
def test_scaled_runner_requires_explicit_bounded_diagnostic(self):
with tempfile.TemporaryDirectory() as tmp:
base = ["run_case.py", "--repo", str(REPO), "--expected-revision", MOCK_REVISION,
"--output", str(Path(tmp) / "case.json"), "--dry-run"]
diagnostic = ["--diagnostic", run_case.SCALED_DIAGNOSTIC, "--early-pressure"]
for n, depth, re, extra, valid in (
(512, 1, 100, [], False), (512, 1, 100, diagnostic, True),
(256, 2, 100, diagnostic, True), (128, 2, 100, diagnostic, False),
(256, 3, 100, diagnostic, False), (256, 2, 1000, diagnostic, False),
(256, 2, 100, diagnostic[:2], False), (512, 2, 100, diagnostic, False)):
argv = base + ["--n", str(n), "--depth", str(depth), "--re", str(re)] + extra
with self.subTest(argv=argv), patch.object(sys, "argv", argv), \
patch.object(run_case, "verify_source", return_value=source()) as verify, \
redirect_stdout(io.StringIO()), patch("sys.stderr", io.StringIO()):
if valid:
self.assertEqual(run_case.main(), 0)
else:
with self.assertRaises(SystemExit) as caught:
run_case.main()
self.assertEqual(caught.exception.code, 2)
verify.assert_not_called()
def test_source_preflight_requires_clean_exact_revision_and_tree(self):
import oracle
from tests.validation.multigrid import check_projected_pairs as gate
@@ -429,7 +319,7 @@ class ProjectedPairsTests(unittest.TestCase):
patch.object(run_case, "runtime_identity", side_effect=AssertionError("CUDA reached")), \
redirect_stdout(stdout):
self.assertEqual(run_case.main(), 0)
verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
verify.assert_called_once_with(REPO, MOCK_REVISION)
self.assertTrue(json.loads(stdout.getvalue())["dry_run"])
self.assertFalse(output.exists())
@@ -459,15 +349,9 @@ class ProjectedPairsTests(unittest.TestCase):
def test_mocked_runner_default_early_and_failure(self):
# Replace the GPU context and steppers before import. No CUDA module or flow is used.
cases = [(early, depth, fail, 128, None, "projected_linear" if early else "constant")
for early, depth, fail in
((False, 1, None), (False, 2, None), (True, 2, None),
for early, depth, fail in ((False, 1, None), (False, 2, None), (True, 2, None),
(True, 2, "resource"), (True, 2, "nonfinite"),
(True, 2, "probe_nonfinite"))]
cases += [(True, depth, None, n, run_case.SCALED_DIAGNOSTIC, scheme)
for n, depth, scheme in ((256, 1, "constant"), (512, 1, "constant"),
(256, 2, "constant"), (256, 2, "projected_linear"))]
for early, depth, fail, n, diagnostic, scheme in cases:
(True, 2, "probe_nonfinite")):
calls, selectors = [], []
class FakeRun:
@@ -511,14 +395,11 @@ class ProjectedPairsTests(unittest.TestCase):
with tempfile.TemporaryDirectory() as tmp:
output = Path(tmp) / "case.json"
argv = ["run_case.py", "--repo", str(REPO), "--re", "100", "--n", str(n),
argv = ["run_case.py", "--repo", str(REPO), "--re", "100", "--n", "128",
"--depth", str(depth), "--output", str(output),
"--expected-revision", MOCK_REVISION, "--run-gpu"]
argv += ["--explosion", scheme]
if early:
argv += ["--early-pressure"]
if diagnostic:
argv += ["--diagnostic", diagnostic]
argv += ["--early-pressure", "--explosion", "projected_linear"]
with patch.dict(sys.modules, modules), patch.object(sys, "argv", argv), \
patch.dict("os.environ", CUDA_VISIBLE_DEVICES="mock"), \
patch.object(run_case, "verify_source", return_value=source()) as verify, \
@@ -531,17 +412,16 @@ class ProjectedPairsTests(unittest.TestCase):
else:
self.assertEqual(run_case.main(), 0)
receipt = json.loads(output.read_text())
verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
self.assertEqual(selectors, [scheme] if depth == 2 else [])
verify.assert_called_once_with(REPO, MOCK_REVISION)
self.assertEqual(selectors, ["projected_linear" if early else "constant"] if depth == 2 else [])
self.assertEqual(receipt["source"], source())
self.assertEqual(receipt["runtime_identity"], RUNTIME)
self.assertEqual(receipt["completed"], not fail)
self.assertEqual(receipt.get("diagnostic"), diagnostic)
self.assertEqual(calls[:2], [0, 1]) # READY probe is part of this trajectory
if not fail:
self.assertEqual([r["tstar"] for r in receipt["rows"]], TIMES)
self.assertEqual(len(receipt["early_pressure"]["rows"]), 65 if early else 0)
self.assertEqual(calls[-1], 16*n)
self.assertEqual(calls[-1], 2048)
else:
self.assertEqual(len(receipt["rows"]), 1)
self.assertEqual(len(receipt["early_pressure"]["rows"]),
File diff suppressed because it is too large Load Diff