Compare commits
8
Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
|
5fb95f5f4d | ||
|
|
46100e1d53 | ||
|
|
de5b0cc080 | ||
|
|
f337652683 | ||
|
|
3554cb2fde | ||
|
|
10069feadb | ||
|
|
ffff742b69 | ||
|
|
92b8afeeed |
No files matched your search
@@ -159,3 +159,181 @@ 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.
|
||||
@@ -1,6 +1,7 @@
|
||||
"""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
|
||||
@@ -8,16 +9,25 @@ 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)
|
||||
BUILD_CAPTURE_SHA256, REFERENCE_HARNESS, SCALED_DIAGNOSTIC)
|
||||
|
||||
# Host-generated snapshots of a135252 at the frozen two local configurations.
|
||||
# Host-generated snapshots of a135252 at the frozen 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],
|
||||
@@ -73,7 +83,11 @@ def valid_builds(record, n, depth):
|
||||
return True
|
||||
|
||||
|
||||
def assess(records, expected_revision=None):
|
||||
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)
|
||||
indexed, duplicates = {}, []
|
||||
for record in records:
|
||||
key = tuple(record.get(k) for k in ("re_L", "n", "depth", "explosion"))
|
||||
@@ -81,9 +95,11 @@ def assess(records, expected_revision=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
|
||||
@@ -106,9 +122,11 @@ def assess(records, expected_revision=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]
|
||||
@@ -135,7 +153,7 @@ def assess(records, expected_revision=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 [8704, 30720]
|
||||
cells = [n*n] if depth == 1 else [17*n*n//32, 15*n*n//8]
|
||||
hops = [16*n * 2**l for l in range(depth)]
|
||||
work = r.get("workload", {})
|
||||
checks[label + " workload"] = (
|
||||
@@ -168,19 +186,22 @@ def assess(records, expected_revision=None):
|
||||
|
||||
comparisons, scheme_results = [], {}
|
||||
for scheme in ("constant", "projected_linear"):
|
||||
key = (100, 128, 2, scheme)
|
||||
key = (100, base_n, 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,
|
||||
uniform128=coarse, baseline_budget=budget, passed=passed,
|
||||
uniform256=fine, amr_to_uniform256=ratio if finite(ratio) else None))
|
||||
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
|
||||
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,
|
||||
@@ -190,19 +211,89 @@ def assess(records, expected_revision=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="uniform256 and AMR/ uniform256 ratios are descriptive, not acceptance gates",
|
||||
finest_control=f"uniform{2*base_n} and AMR/ uniform{2*base_n} 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")
|
||||
result = assess([json.loads(path.read_text()) for path in args.files], args.expected_revision)
|
||||
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())
|
||||
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)
|
||||
@@ -0,0 +1,188 @@
|
||||
"""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
|
||||
from pathlib import Path
|
||||
|
||||
import assess_projected_pairs as tg
|
||||
import shear_wave as sw
|
||||
from run_case import SCALED_DIAGNOSTIC, SCALED_SHEAR, SHEAR_WAVE_SHA256
|
||||
|
||||
# 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")
|
||||
|
||||
|
||||
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 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
|
||||
|
||||
|
||||
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])
|
||||
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_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,6 +18,32 @@ 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)
|
||||
|
||||
|
||||
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):
|
||||
@@ -64,7 +90,7 @@ def provenance(repo):
|
||||
)
|
||||
|
||||
|
||||
def verify_source(repo, revision):
|
||||
def verify_source(repo, revision, flow="taylor_green"):
|
||||
"""Exact clean harness commit, fixed production tree, and reviewed capture helper."""
|
||||
import oracle
|
||||
from tests.validation.multigrid import check_projected_pairs as gate
|
||||
@@ -83,6 +109,13 @@ def verify_source(repo, revision):
|
||||
if result["source_tree_sha256"] != SOURCE_SHA256:
|
||||
raise RuntimeError("Wrong numerical source fingerprint")
|
||||
result["build_capture_sha256"] = BUILD_CAPTURE_SHA256
|
||||
if flow == "shear_wave":
|
||||
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
|
||||
|
||||
|
||||
@@ -145,26 +178,48 @@ 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], required=True)
|
||||
p.add_argument("--n", type=int, choices=[32, 64, 128, 256, 512], 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],
|
||||
help="explicit four-case N256/fine512 same-placement diagnostic "
|
||||
"(scaled_n256: Taylor-Green; scaled_shear_n256: shear wave)")
|
||||
p.add_argument("--early-pressure", action="store_true",
|
||||
help="also save 65 scalar pressure samples at t*=k*.001875, k=0..64")
|
||||
p.add_argument("--flow", choices=["taylor_green", "shear_wave"], default="taylor_green")
|
||||
args = p.parse_args()
|
||||
if args.flow == "shear_wave" and (args.re != 100 or args.depth == 3 or args.early_pressure):
|
||||
p.error("shear_wave requires --re 100, depth 1 or 2 and no --early-pressure")
|
||||
initial, metrics, limits, initialization = flow_functions(args.flow)
|
||||
if args.output.exists():
|
||||
p.error("Output already exists; choose a new path to retain prior evidence")
|
||||
if args.explosion != "constant" and args.depth != 2:
|
||||
p.error("projected_linear requires depth 2")
|
||||
if 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 or scaled_shear_n256")
|
||||
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)
|
||||
source = verify_source(repo, args.expected_revision, args.flow)
|
||||
runtime = runtime_source(repo)
|
||||
if args.dry_run:
|
||||
print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True)))
|
||||
print(json.dumps(dict(source=source, runtime_source=runtime, dry_run=True,
|
||||
**({"diagnostic": args.diagnostic} if args.diagnostic else {}))))
|
||||
return 0
|
||||
from CelerisLab.cuda.context import CudaContext
|
||||
from CelerisLab.lbm.blocks import BlockLayout2D
|
||||
@@ -197,7 +252,7 @@ def main():
|
||||
mach=U * np.sqrt(3),
|
||||
viscosity=cfg.viscosity,
|
||||
tau_by_level=[0.5 + 3 * c.viscosity for c in h.configs],
|
||||
initialization=INITIALIZATION,
|
||||
initialization=initialization,
|
||||
numerics=dict(lattice_model=cfg.lattice_model, collision=cfg.collision,
|
||||
store_precision=cfg.store_precision, ddf_shifting=cfg.ddf_shifting,
|
||||
streaming=cfg.streaming, les_enabled=cfg.les_enabled,
|
||||
@@ -205,13 +260,20 @@ 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")
|
||||
@@ -248,10 +310,8 @@ def main():
|
||||
initial_mass = None
|
||||
for step in sorted(set(main_steps + early_steps)):
|
||||
run.step(step - run.step_count)
|
||||
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 = sample(run, masks, step, args.n, cfg.viscosity, initial_mass,
|
||||
metrics, args.depth)
|
||||
row["step"] = step
|
||||
if initial_mass is None:
|
||||
initial_mass = row["mass"]
|
||||
@@ -278,10 +338,8 @@ 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)
|
||||
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")
|
||||
sample(run, masks, 1, args.n, cfg.viscosity, initial_mass,
|
||||
metrics, args.depth, "READY")
|
||||
report["ready_after_step"] = 1
|
||||
save()
|
||||
print(
|
||||
@@ -301,7 +359,7 @@ def main():
|
||||
completed=True,
|
||||
seconds=time.monotonic() - started,
|
||||
screening_pass=all(
|
||||
row[key] <= limit for row in report["rows"] for key, limit in LIMITS.items()
|
||||
row[key] <= limit for row in report["rows"] for key, limit in limits.items()
|
||||
),
|
||||
)
|
||||
save()
|
||||
|
||||
@@ -0,0 +1,312 @@
|
||||
"""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,6 +17,7 @@ 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
|
||||
@@ -31,7 +32,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},
|
||||
harness_sha256={"run_case.py": "a" * 64, "oracle.py": "c" * 64},
|
||||
build_capture_sha256=run_case.BUILD_CAPTURE_SHA256)
|
||||
|
||||
|
||||
@@ -42,12 +43,14 @@ def build_metadata(size):
|
||||
arch="sm_70", flags=["--use_fast_math"], ptx_sha256="b" * 64)
|
||||
|
||||
|
||||
def records(early=False):
|
||||
def records(early=False, diagnostic=None):
|
||||
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)
|
||||
@@ -66,6 +69,8 @@ def records(early=False):
|
||||
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
|
||||
|
||||
|
||||
@@ -251,7 +256,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):
|
||||
for size in (128, 256, 512):
|
||||
cfg = replace(run_case.config(size, 100), nx=10, ny=10)
|
||||
cfg.validate()
|
||||
root = compiler_v2.generate_config(cfg)
|
||||
@@ -263,6 +268,111 @@ 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
|
||||
@@ -319,7 +429,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)
|
||||
verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
|
||||
self.assertTrue(json.loads(stdout.getvalue())["dry_run"])
|
||||
self.assertFalse(output.exists())
|
||||
|
||||
@@ -349,9 +459,15 @@ 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.
|
||||
for early, depth, fail in ((False, 1, None), (False, 2, None), (True, 2, None),
|
||||
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),
|
||||
(True, 2, "resource"), (True, 2, "nonfinite"),
|
||||
(True, 2, "probe_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:
|
||||
calls, selectors = [], []
|
||||
|
||||
class FakeRun:
|
||||
@@ -395,11 +511,14 @@ class ProjectedPairsTests(unittest.TestCase):
|
||||
|
||||
with tempfile.TemporaryDirectory() as tmp:
|
||||
output = Path(tmp) / "case.json"
|
||||
argv = ["run_case.py", "--repo", str(REPO), "--re", "100", "--n", "128",
|
||||
argv = ["run_case.py", "--repo", str(REPO), "--re", "100", "--n", str(n),
|
||||
"--depth", str(depth), "--output", str(output),
|
||||
"--expected-revision", MOCK_REVISION, "--run-gpu"]
|
||||
argv += ["--explosion", scheme]
|
||||
if early:
|
||||
argv += ["--early-pressure", "--explosion", "projected_linear"]
|
||||
argv += ["--early-pressure"]
|
||||
if diagnostic:
|
||||
argv += ["--diagnostic", diagnostic]
|
||||
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, \
|
||||
@@ -412,16 +531,17 @@ 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)
|
||||
self.assertEqual(selectors, ["projected_linear" if early else "constant"] if depth == 2 else [])
|
||||
verify.assert_called_once_with(REPO, MOCK_REVISION, "taylor_green")
|
||||
self.assertEqual(selectors, [scheme] 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], 2048)
|
||||
self.assertEqual(calls[-1], 16*n)
|
||||
else:
|
||||
self.assertEqual(len(receipt["rows"]), 1)
|
||||
self.assertEqual(len(receipt["early_pressure"]["rows"]),
|
||||
|
||||
@@ -0,0 +1,560 @@
|
||||
"""Standalone CPU manufactured checks for the oblique shear wave.
|
||||
|
||||
Run from this directory with PYTHONPATH=<checkout>/src (same as test_oracle.py).
|
||||
"""
|
||||
|
||||
import hashlib
|
||||
import io
|
||||
import json
|
||||
import math
|
||||
import sys
|
||||
import tempfile
|
||||
from contextlib import redirect_stdout
|
||||
from unittest.mock import patch
|
||||
import unittest
|
||||
from pathlib import Path
|
||||
|
||||
import numpy as np
|
||||
import oracle
|
||||
import run_case
|
||||
import shear_wave as sw
|
||||
from oracle import CX, CY, U, TIMES, centers, rectangles
|
||||
|
||||
N = 32
|
||||
NU = U * N / sw.RE
|
||||
HERE = Path(__file__).resolve().parent
|
||||
ORACLE_SHA256_EEC4677 = "83e4d7425d45538278e9fb181af450a36358bd9a5ed6be57ff611dd9ae5eb35b"
|
||||
|
||||
|
||||
def composite_masks(n):
|
||||
"""Level-wide union masks for rectangles(n, 2): root complement + fine box."""
|
||||
(x0, y0, x1, y1), = rectangles(n, 2)[0]
|
||||
root = np.ones((n, n), dtype=bool)
|
||||
root[y0:y1, x0:x1] = False
|
||||
fine = np.zeros((2 * n, 2 * n), dtype=bool)
|
||||
fine[2 * y0:2 * y1, 2 * x0:2 * x1] = True
|
||||
return [root, fine]
|
||||
|
||||
|
||||
def populations(masks, field, t=0.0):
|
||||
"""Equilibrium G of an arbitrary (ux, uy, rho) field on every level."""
|
||||
out = []
|
||||
for level, m in enumerate(masks):
|
||||
x, y = centers(m.shape, level)
|
||||
ux, uy, rho = field(x, y, t)
|
||||
out.append(sw.equilibrium(rho, ux, uy))
|
||||
return out
|
||||
|
||||
|
||||
def exact_field(phase=0.0, rate_factor=1.0, mode=sw.MODE, longitudinal=0.0, ripple=0.0,
|
||||
offset=0.0):
|
||||
k = 2 * np.pi / N * np.array(mode, dtype=float)
|
||||
e = np.array([-mode[1], mode[0]], dtype=float) / math.hypot(*mode)
|
||||
rate = rate_factor * sw.decay_rate(N, NU)
|
||||
kh = k / np.linalg.norm(k)
|
||||
|
||||
def field(x, y, t):
|
||||
theta = k[0] * x + k[1] * y + phase
|
||||
a = U * np.exp(-rate * t)
|
||||
ux = a * e[0] * np.sin(theta) + longitudinal * kh[0] * np.sin(theta)
|
||||
uy = a * e[1] * np.sin(theta) + longitudinal * kh[1] * np.sin(theta)
|
||||
rho = 1 + offset + ripple * np.cos(theta)
|
||||
return ux, uy, rho
|
||||
|
||||
return field
|
||||
|
||||
|
||||
class ConventionTests(unittest.TestCase):
|
||||
def test_mode_angle_and_polarization(self):
|
||||
k = sw.wavevector(N)
|
||||
self.assertAlmostEqual(sw.ANGLE_DEG, 26.565051177077994)
|
||||
self.assertAlmostEqual(math.degrees(math.atan2(k[1], k[0])), sw.ANGLE_DEG)
|
||||
self.assertAlmostEqual(float(sw.E @ k), 0.0)
|
||||
self.assertAlmostEqual(float(sw.E @ sw.E), 1.0)
|
||||
self.assertAlmostEqual(sw.eps(128), 5 * (2 * math.pi / 128) ** 2)
|
||||
|
||||
def test_navier_stokes_and_vanishing_advection(self):
|
||||
# Independent central differences: du/dt = nu lap u, (u.grad)u = 0, div u = 0.
|
||||
x = np.array([2.7, 8.1, 13.3])
|
||||
y = np.array([3.2, 7.4, 21.6])
|
||||
t, h = 10.0, 1e-3
|
||||
f = lambda x, y, t: np.array(sw.exact(x, y, t, N, NU)[:2])
|
||||
z = f(x, y, t)
|
||||
dt = (f(x, y, t + h) - f(x, y, t - h)) / (2 * h)
|
||||
dx = (f(x + h, y, t) - f(x - h, y, t)) / (2 * h)
|
||||
dy = (f(x, y + h, t) - f(x, y - h, t)) / (2 * h)
|
||||
lap = (f(x + h, y, t) + f(x - h, y, t) + f(x, y + h, t) + f(x, y - h, t) - 4 * z) / h**2
|
||||
np.testing.assert_allclose(dt, NU * lap, atol=1e-9)
|
||||
# Analytic identities: div u = A cos(theta) (E.k) and (u.grad)u = A sin(theta) (E.k) u/A.
|
||||
k = sw.wavevector(N)
|
||||
self.assertLess(abs(float(sw.E @ k)), 1e-17)
|
||||
# Central differences do not cancel exactly: each derivative carries
|
||||
# -(k_j h)^2/6 relative truncation, so div_h = -A cos(theta) h^2/6 sum_j E_j k_j^3
|
||||
# + O(h^4) (nonzero for an oblique k). Bound it with that leading term plus roundoff.
|
||||
a = U * math.exp(-sw.decay_rate(N, NU) * t)
|
||||
trunc = a * h * h / 6 * float(np.abs(sw.E * k**3).sum()) * 1.01
|
||||
# Argument rounding: eps*max|x|/h*U*|k| ~ 2.2e-16*21.6/1e-3*0.03*0.44 ~ 6e-14.
|
||||
roundoff = 1e-13
|
||||
self.assertLess(np.max(np.abs(dx[0] + dy[1])), trunc + roundoff)
|
||||
# Discrete advection: same truncation scaled by |u| <= a.
|
||||
self.assertLess(np.max(np.abs(z[0] * dx + z[1] * dy)), a * trunc + roundoff)
|
||||
# Halving h reduces the divergence truncation about fourfold.
|
||||
hh = h / 2
|
||||
div2 = ((f(x + hh, y, t) - f(x - hh, y, t))[0] + (f(x, y + hh, t) - f(x, y - hh, t))[1]) / (2 * hh)
|
||||
np.testing.assert_allclose(div2 * 4, dx[0] + dy[1], rtol=1e-2, atol=6 * roundoff)
|
||||
self.assertTrue(np.all(sw.exact(x, y, t, N, NU)[2] == 0))
|
||||
|
||||
|
||||
class InitialTests(unittest.TestCase):
|
||||
def test_moments_on_each_level(self):
|
||||
for level, m in enumerate(composite_masks(N)):
|
||||
g = sw.initial(m.shape, level, N, NU)
|
||||
x, y = centers(m.shape, level)
|
||||
ux, uy, _ = sw.exact(x, y, 0.0, N, NU)
|
||||
rho = g.sum(axis=0)
|
||||
np.testing.assert_allclose(rho, 1.0, atol=1e-14)
|
||||
np.testing.assert_allclose(np.tensordot(CX, g, 1) / rho, ux, atol=1e-15)
|
||||
np.testing.assert_allclose(np.tensordot(CY, g, 1) / rho, uy, atol=1e-15)
|
||||
|
||||
def test_equilibrium_matches_frozen_oracle_formula(self):
|
||||
# Same equilibrium as oracle.initial when fed the TG field.
|
||||
shape = (8, 8)
|
||||
x, y = centers(shape, 1)
|
||||
ux, uy, p = oracle.exact(x, y, 0.0, 4.0, 0.01)
|
||||
np.testing.assert_array_equal(
|
||||
sw.equilibrium(1 + 3 * p, ux, uy), oracle.initial(shape, 1, 4.0, 0.01)
|
||||
)
|
||||
|
||||
|
||||
class MetricTests(unittest.TestCase):
|
||||
def setUp(self):
|
||||
self.masks = composite_masks(N)
|
||||
|
||||
def sample(self, field, t):
|
||||
return sw.metrics(populations(self.masks, field, t), self.masks, t, N, NU, None)
|
||||
|
||||
def test_exact_field_on_composite_partition(self):
|
||||
for t in (0.0, 100.0, 512.0):
|
||||
r = self.sample(exact_field(), t)
|
||||
self.assertLess(r["velocity_l2"], 1e-12)
|
||||
self.assertLess(r["residual_l2"], 1e-12)
|
||||
self.assertLess(r["longitudinal_l2"], 1e-12)
|
||||
self.assertLess(r["phase_drift_abs"], 1e-12)
|
||||
self.assertAlmostEqual(r["amplitude_ratio"], 1.0, places=12)
|
||||
self.assertLess(abs(r["amplitude_error"]), 1e-12)
|
||||
self.assertLess(r["spurious_pressure_scaled"], 1e-10)
|
||||
self.assertEqual(len(r["spurious_pressure_scaled_by_level"]), 2)
|
||||
self.assertEqual(r["leaf_cells"], int(sum(m.sum() for m in self.masks)))
|
||||
|
||||
def test_phase_control(self):
|
||||
r = self.sample(exact_field(phase=0.01), 0.0)
|
||||
self.assertAlmostEqual(r["phase_drift"], 0.01, delta=1e-9)
|
||||
self.assertLess(abs(r["amplitude_error"]), 1e-12) # R = hypot(a, b) is phase-free
|
||||
self.assertLess(r["residual_l2"], 1e-12) # a pure phase shift stays modal
|
||||
|
||||
def test_rate_control_and_fit(self):
|
||||
for factor, expected in ((1.0, 0.0), (1.01, 0.01)):
|
||||
field = exact_field(rate_factor=factor)
|
||||
rows = []
|
||||
for tstar in TIMES:
|
||||
step = round(tstar * N / U)
|
||||
rows.append(dict(self.sample(field, step), step=step))
|
||||
fit = sw.fit_rate(rows, N, NU)
|
||||
self.assertAlmostEqual(fit["rate_error"], expected, delta=1e-9)
|
||||
np.testing.assert_allclose(fit["interval_rate_errors"], expected, atol=1e-9)
|
||||
self.assertLess(fit["log_fit_residual_max"], 1e-9)
|
||||
self.assertAlmostEqual(fit["exact_rate_tstar"], 5 * (2 * math.pi) ** 2 / sw.RE)
|
||||
self.assertAlmostEqual(fit["exact_rate_tstar"], 1.97392088, places=8)
|
||||
|
||||
def test_fit_ignores_startup_sample(self):
|
||||
rows = [dict(tstar=0.0, step=0, a_hat=2 * U, b_hat=0.0)]
|
||||
rows += [dict(tstar=s, step=round(s * N / U), b_hat=0.0,
|
||||
a_hat=U * math.exp(-sw.decay_rate(N, NU) * round(s * N / U)))
|
||||
for s in TIMES[1:]]
|
||||
self.assertLess(sw.fit_rate(rows, N, NU)["rate_error"], 1e-12)
|
||||
rows[-1]["a_hat"] = float("nan")
|
||||
with self.assertRaises(ValueError):
|
||||
sw.fit_rate(rows, N, NU)
|
||||
|
||||
def test_wrong_mode_control(self):
|
||||
# Swapping (2,1) for (1,2), i.e. 26.6 vs 63.4 deg, is not silently modal.
|
||||
# Composite quadrature need not give an exactly zero projection.
|
||||
r = self.sample(exact_field(mode=(1, 2)), 0.0)
|
||||
self.assertLess(abs(r["amplitude_ratio"]), 0.05)
|
||||
self.assertGreater(r["residual_l2"], 0.95)
|
||||
|
||||
def test_longitudinal_control(self):
|
||||
r = self.sample(exact_field(longitudinal=1e-3), 0.0)
|
||||
self.assertAlmostEqual(r["longitudinal_l2"], 1e-3 / U, delta=1e-9)
|
||||
# The full-vector residual and velocity error see it too.
|
||||
self.assertAlmostEqual(r["residual_l2"], 1e-3 / U, delta=1e-9)
|
||||
self.assertAlmostEqual(r["velocity_l2"], 1e-3 / U, delta=1e-9)
|
||||
self.assertLess(r["phase_drift_abs"], 1e-12)
|
||||
self.assertAlmostEqual(r["amplitude_ratio"], 1.0, places=12)
|
||||
|
||||
def test_pressure_gauge_and_ripple(self):
|
||||
r = self.sample(exact_field(offset=1e-4), 0.0)
|
||||
self.assertLess(r["spurious_pressure_scaled"], 1e-10) # one global gauge
|
||||
ripple = 3e-5
|
||||
r = self.sample(exact_field(ripple=ripple), 0.0)
|
||||
# p = ripple/3 cos(theta); rms over a whole number of periods = ripple/(3 sqrt 2).
|
||||
# Composite midpoint quadrature is not exactly periodic: 3% tolerance.
|
||||
expected = ripple / (3 * math.sqrt(2)) / (U * U / 4)
|
||||
self.assertAlmostEqual(r["spurious_pressure_scaled"], expected, delta=0.03 * expected)
|
||||
self.assertAlmostEqual(r["spurious_pressure_rms"], expected * U * U / 4,
|
||||
delta=0.03 * expected * U * U / 4)
|
||||
|
||||
def test_volume_weighting_and_double_count(self):
|
||||
bad = [np.ones_like(m) for m in self.masks]
|
||||
with self.assertRaises(ValueError):
|
||||
sw.metrics(populations(bad, exact_field()), bad, 0.0, N, NU)
|
||||
# Unweighted (count-based) averaging would change the fit; weights matter.
|
||||
r = self.sample(exact_field(), 0.0)
|
||||
self.assertEqual(r["leaf_cells"], N * N - self.masks[1].sum() // 4 + self.masks[1].sum())
|
||||
|
||||
|
||||
def record(n, depth, rows, explosion="constant", completed=True, source="s"):
|
||||
for row in rows:
|
||||
row["spurious_pressure_scaled_by_level"] = [row["spurious_pressure_scaled"]] * depth
|
||||
return dict(n=n, depth=depth, viscosity=U * n / sw.RE, rows=rows, completed=completed,
|
||||
explosion=explosion, source=source)
|
||||
|
||||
|
||||
def rows_for(n, rate_error=0.0, residual=1e-4, pressure=1e-3, phase=0.0):
|
||||
rate = (1 + rate_error) * sw.decay_rate(n, U * n / sw.RE)
|
||||
out = []
|
||||
for tstar in TIMES:
|
||||
step = round(tstar * n / U)
|
||||
amp = U * math.exp(-rate * step)
|
||||
out.append(dict(tstar=step * U / n, step=step, a_hat=amp, b_hat=0.0,
|
||||
amplitude_ratio=1.0, modal_amplitude=amp, amplitude_exact=amp,
|
||||
phase_drift=phase, mass=float(n * n), rho_min=0.999, rho_max=1.001,
|
||||
max_speed=amp, longitudinal_l2=0.0, spurious_pressure_max_scaled=pressure,
|
||||
amplitude_error=0.0, residual_l2=residual, phase_drift_abs=phase,
|
||||
spurious_pressure_scaled=pressure,
|
||||
spurious_pressure_rms=pressure * U * U / 4,
|
||||
velocity_l2=residual, mass_relative=0.0))
|
||||
return out
|
||||
|
||||
|
||||
class AssessTests(unittest.TestCase):
|
||||
def test_budget_form_is_max_not_sum(self):
|
||||
base = record(128, 1, rows_for(128, rate_error=1e-3, pressure=0.015))
|
||||
fine = record(256, 1, rows_for(256, rate_error=2e-4))
|
||||
# rate budget = max(2*1e-3, 0.002) = 0.002; 0.0025 would pass a sum form.
|
||||
refined = record(128, 2, rows_for(128, rate_error=0.0025))
|
||||
report = sw.assess(base, fine, [refined])
|
||||
self.assertEqual(report["status"], "FAIL")
|
||||
self.assertAlmostEqual(report["refined"]["constant"]["budget"]["rate_error"], 0.002)
|
||||
# pressure budget = max(2*0.015, 0.02) = 0.03 (control-dominated).
|
||||
self.assertAlmostEqual(
|
||||
report["refined"]["constant"]["budget"]["spurious_pressure_scaled"], 0.03)
|
||||
self.assertIn("finest", report["refined"]["constant"]["control_ratios"])
|
||||
ok = record(128, 2, rows_for(128, rate_error=0.0015), explosion="projected_linear")
|
||||
self.assertEqual(sw.assess(base, fine, [ok])["status"], "PASS")
|
||||
|
||||
def test_control_outcomes_are_distinct(self):
|
||||
fine = record(256, 1, rows_for(256))
|
||||
refined = record(128, 2, rows_for(128))
|
||||
failing = record(128, 1, rows_for(128, residual=0.5)) # complete, > eps(128)
|
||||
report = sw.assess(failing, fine, [refined])
|
||||
self.assertEqual(report["status"], "CONTROL_SCREEN_FAIL")
|
||||
own = report["refined"]["constant"]
|
||||
self.assertTrue(own["own_screens_pass"]) # refined own screens kept
|
||||
self.assertFalse(own["controlled_pass"]) # but no controlled comparison passes
|
||||
good = record(128, 1, rows_for(128))
|
||||
self.assertEqual(sw.assess(None, fine, [refined])["status"], "INVALID_CONTROL")
|
||||
incomplete = record(128, 1, rows_for(128), completed=False)
|
||||
self.assertEqual(sw.assess(incomplete, fine, [refined])["status"], "INVALID_CONTROL")
|
||||
mixed = record(256, 1, rows_for(256), source="other")
|
||||
self.assertEqual(sw.assess(good, mixed, [refined])["status"], "INVALID_CONTROL")
|
||||
|
||||
def test_invalid_rows_are_rejected_before_screens(self):
|
||||
base = record(128, 1, rows_for(128))
|
||||
fine = record(256, 1, rows_for(256))
|
||||
for mutate in (
|
||||
lambda rows: rows[2].update(rho_min=float("nan")),
|
||||
lambda rows: rows[2].update(spurious_pressure_scaled_by_level=[0.0, float("nan")]),
|
||||
lambda rows: rows[2].update(step=rows[1]["step"], tstar=rows[1]["tstar"]),
|
||||
lambda rows: rows[2].update(residual_l2=-1e-6),
|
||||
lambda rows: rows[2].pop("longitudinal_l2"),
|
||||
):
|
||||
refined = record(128, 2, rows_for(128))
|
||||
mutate(refined["rows"])
|
||||
self.assertEqual(sw.assess(base, fine, [refined])["status"], "INVALID_RUN")
|
||||
# Legitimate signs stay allowed: negative amplitude and phase errors.
|
||||
refined = record(128, 2, rows_for(128))
|
||||
for row in refined["rows"]:
|
||||
row.update(amplitude_error=-1e-3, phase_drift=-1e-4, phase_drift_abs=1e-4)
|
||||
report = sw.assess(base, fine, [refined])
|
||||
self.assertEqual(report["status"], "PASS")
|
||||
self.assertEqual(sw.assess(base, fine, [])["status"], "INVALID_CONTROL")
|
||||
with self.assertRaises(ValueError):
|
||||
sw.assess(base, fine, [record(128, 2, rows_for(128)), record(128, 2, rows_for(128))])
|
||||
|
||||
def test_modal_extinction_is_a_failure_not_malformed(self):
|
||||
def extinct(rows):
|
||||
for row in rows[1:]:
|
||||
row.update(a_hat=0.0, b_hat=0.0, modal_amplitude=0.0, amplitude_ratio=0.0,
|
||||
amplitude_error=-1.0)
|
||||
return rows
|
||||
base = record(128, 1, rows_for(128))
|
||||
fine = record(256, 1, rows_for(256))
|
||||
refined = record(128, 2, extinct(rows_for(128)))
|
||||
report = sw.assess(base, fine, [refined])
|
||||
self.assertEqual(report["status"], "FAIL")
|
||||
own = report["refined"]["constant"]
|
||||
self.assertTrue(own["valid"])
|
||||
self.assertFalse(own["rate_fit_valid"])
|
||||
self.assertIsNone(own["rate_error"])
|
||||
self.assertIsNone(own["control_ratios"]["root"]["rate_error"])
|
||||
self.assertFalse(own["own_screens_pass"] or own["controlled_pass"] or own["budget_pass"])
|
||||
self.assertEqual(own["residual_l2"], 1e-4) # other finite metrics retained
|
||||
json.loads(json.dumps(report, allow_nan=False))
|
||||
control = sw.assess(record(128, 1, extinct(rows_for(128))), fine,
|
||||
[record(128, 2, rows_for(128))])
|
||||
self.assertEqual(control["status"], "CONTROL_SCREEN_FAIL")
|
||||
self.assertIsNone(control["refined"]["constant"]["budget"]["rate_error"])
|
||||
self.assertFalse(control["refined"]["constant"]["controlled_pass"])
|
||||
json.loads(json.dumps(control, allow_nan=False))
|
||||
|
||||
def test_wrong_shapes(self):
|
||||
base = record(128, 1, rows_for(128))
|
||||
fine = record(256, 1, rows_for(256))
|
||||
refined = record(128, 2, rows_for(128))
|
||||
with self.assertRaises(ValueError):
|
||||
sw.assess(base, record(128, 1, rows_for(128)), [refined])
|
||||
with self.assertRaises(ValueError):
|
||||
sw.assess(record(128, 1, rows_for(128)), fine, [record(256, 2, rows_for(256))])
|
||||
|
||||
def test_incomplete_record_fails(self):
|
||||
base = record(128, 1, rows_for(128))
|
||||
fine = record(256, 1, rows_for(256))
|
||||
refined = record(128, 2, rows_for(128), completed=False)
|
||||
self.assertEqual(sw.assess(base, fine, [refined])["status"], "FAIL")
|
||||
|
||||
|
||||
class FakeRun:
|
||||
def __init__(self, arrays):
|
||||
self.arrays = arrays
|
||||
|
||||
def download(self):
|
||||
return self.arrays if len(self.arrays) > 1 else self.arrays[0]
|
||||
|
||||
|
||||
class RunnerSampleTests(unittest.TestCase):
|
||||
def test_refined_shear_sample_with_list_diagnostics(self):
|
||||
masks = composite_masks(N)
|
||||
arrays = [sw.initial(m.shape, level, N, NU) for level, m in enumerate(masks)]
|
||||
row = run_case.sample(FakeRun(arrays), masks, 0, N, NU, None, sw.metrics, 2)
|
||||
self.assertIsInstance(row["spurious_pressure_scaled_by_level"], list)
|
||||
self.assertTrue(run_case.finite_row(row))
|
||||
bad = [a.copy() for a in arrays]
|
||||
bad[1][0, masks[1]] = np.nan
|
||||
with self.assertRaises(FloatingPointError):
|
||||
run_case.sample(FakeRun(bad), masks, 0, N, NU, None, sw.metrics, 2, "READY")
|
||||
self.assertFalse(run_case.finite_row(dict(a=1.0, b=[0.0, float("inf")])))
|
||||
|
||||
def test_taylor_green_sample_unchanged(self):
|
||||
masks = composite_masks(N)
|
||||
arrays = [oracle.initial(m.shape, level, N, NU) for level, m in enumerate(masks)]
|
||||
row = run_case.sample(FakeRun(arrays), masks, 0, N, NU, None, oracle.metrics, 2)
|
||||
self.assertEqual(row, oracle.metrics(arrays, masks, 0, N, NU, None))
|
||||
uniform = [oracle.initial((N, N), 0, N, NU)]
|
||||
row = run_case.sample(FakeRun(uniform), [np.ones((N, N), bool)], 0, N, NU, None,
|
||||
oracle.metrics, 1)
|
||||
self.assertEqual(row["leaf_cells"], N * N)
|
||||
|
||||
|
||||
def receipt(n, depth, explosion="constant", diagnostic=None):
|
||||
"""Synthetic receipt satisfying every contract field of the fixed four-case set."""
|
||||
import assess_projected_pairs as tg
|
||||
steps = run_case.sample_steps(n)[0]
|
||||
rows = rows_for(n)
|
||||
for row, step in zip(rows, steps):
|
||||
row.update(step=step, tstar=step * U / n,
|
||||
spurious_pressure_scaled_by_level=[row["spurious_pressure_scaled"]] * depth)
|
||||
cells = [n * n] if depth == 1 else [17 * n * n // 32, 15 * n * n // 8]
|
||||
hops = [16 * n * 2**level for level in range(depth)]
|
||||
builds = []
|
||||
for level in range(depth):
|
||||
size = n * 2**level
|
||||
builds.append(dict(source_sha256=tg.BUILD_SOURCES[size], arch="sm_70",
|
||||
flags=["--use_fast_math"], n_objects=0, ptx_sha256="a" * 64,
|
||||
config=dict(nx=10, ny=10, nz=1, lattice_model="D2Q9", collision="SRT",
|
||||
store_precision="FP32", streaming="double_buffer",
|
||||
ddf_shifting=True, les_enabled=False, velocity=0.0,
|
||||
viscosity=U * size / 100, rho=1.0,
|
||||
omega=1 / (.5 + 3 * U * size / 100),
|
||||
threads_per_block=256)))
|
||||
checkout = "/x/CelerisLab"
|
||||
return dict(
|
||||
n=n, depth=depth, re_L=100, explosion=explosion, U=U, mach=U * math.sqrt(3),
|
||||
viscosity=U * n / 100, tau_by_level=[0.5 + 3 * U * n / 100 * 2**l for l in range(depth)],
|
||||
source=dict(revision="r", dirty=False, source_tree_sha256=run_case.SOURCE_SHA256,
|
||||
harness_sha256={"x.py": "h"}, build_capture_sha256=run_case.BUILD_CAPTURE_SHA256,
|
||||
shear_wave_sha256=run_case.SHEAR_WAVE_SHA256),
|
||||
runtime_identity=dict(device=dict(name="V100", pci_bus_id="0", capability=[7, 0],
|
||||
driver=1, runtime=[12, 0]),
|
||||
nvcc_path="/nvcc", nvcc_version="v", python="3", numpy="1",
|
||||
pycuda="2"),
|
||||
runtime_source=dict(checkout=checkout, module=checkout + "/src/CelerisLab/__init__.py",
|
||||
compiler=checkout + "/src/CelerisLab/cuda/compiler_v2.py",
|
||||
kernel_root=checkout + "/src/CelerisLab/lbm/kernels"),
|
||||
cuda_visible_devices="0", numerical_revision=run_case.NUMERICAL_REVISION,
|
||||
reference_harness=run_case.REFERENCE_HARNESS, numerics=dict(tg.NUMERICS),
|
||||
initialization=sw.INITIALIZATION, limits=dict(sw.LIMITS),
|
||||
early_pressure=dict(enabled=False, rows=[], sampled_max_pressure_l2=None),
|
||||
refinements=[] if depth == 1 else rectangles(n, depth),
|
||||
workload=dict(active_cells_by_level=cells, level_hops=hops, root_intervals=16 * n,
|
||||
active_collision_updates_by_level=[c * h for c, h in zip(cells, hops)],
|
||||
active_collision_updates=sum(c * h for c, h in zip(cells, hops))),
|
||||
actual_build_sidecars=builds, rows=rows, completed=True,
|
||||
flow="shear_wave", mode=list(sw.MODE), polarization=sw.E.tolist(),
|
||||
angle_deg=sw.ANGLE_DEG, phase=sw.PHI, eps=sw.eps(n),
|
||||
**({"diagnostic": diagnostic} if diagnostic else {}))
|
||||
|
||||
|
||||
class ContractTests(unittest.TestCase):
|
||||
def setUp(self):
|
||||
import assess_shear_wave
|
||||
self.assess = assess_shear_wave.assess
|
||||
self.records = [receipt(128, 1), receipt(256, 1), receipt(128, 2),
|
||||
receipt(128, 2, "projected_linear")]
|
||||
|
||||
def test_complete_contract_passes(self):
|
||||
result = self.assess(self.records, "r")
|
||||
self.assertTrue(result["contract_valid"], [k for k, v in result["checks"].items() if not v])
|
||||
self.assertEqual(result["status"], "PASS")
|
||||
|
||||
def test_reproduced_false_passes_are_rejected(self):
|
||||
a = self.assess
|
||||
self.assertEqual(a(self.records[:2], "r")["status"], "INVALID_SET") # empty refined
|
||||
dup = [*self.records[:3], receipt(128, 2)]
|
||||
self.assertEqual(a(dup, "r")["status"], "INVALID_SET") # duplicate constant-only
|
||||
for mutate in (
|
||||
lambda r: r.pop("source"),
|
||||
lambda r: r.update(flow="taylor_green"),
|
||||
lambda r: r.update(mode=[1, 2]),
|
||||
lambda r: r.update(re_L=1000),
|
||||
lambda r: r.update(refinements=[[(0, 0, 8, 8)]]),
|
||||
lambda r: r["actual_build_sidecars"][0].update(arch="sm_80"),
|
||||
lambda r: r["rows"][0].update(step=1, tstar=U / 128), # wrong initial time
|
||||
lambda r: r["source"].update(shear_wave_sha256="0" * 64),
|
||||
):
|
||||
records = [dict(r, rows=[dict(x) for x in r["rows"]],
|
||||
source=dict(r["source"]),
|
||||
actual_build_sidecars=[dict(b) for b in r["actual_build_sidecars"]])
|
||||
for r in self.records]
|
||||
mutate(records[2])
|
||||
self.assertEqual(a(records, "r")["status"], "INVALID_SET")
|
||||
mixed = [dict(r) for r in self.records]
|
||||
mixed[1] = dict(mixed[1], source=dict(mixed[1]["source"], revision="other"))
|
||||
self.assertEqual(a(mixed, "r")["status"], "INVALID_SET")
|
||||
self.assertEqual(a(self.records, "wrong-revision")["status"], "INVALID_SET")
|
||||
|
||||
|
||||
def scaled_set(rate_error=0.0):
|
||||
d = run_case.SCALED_SHEAR
|
||||
out = [receipt(256, 1, diagnostic=d), receipt(512, 1, diagnostic=d),
|
||||
receipt(256, 2, diagnostic=d), receipt(256, 2, "projected_linear", diagnostic=d)]
|
||||
if rate_error:
|
||||
out[3]["rows"] = rows_for(256, rate_error=rate_error)
|
||||
for row, step in zip(out[3]["rows"], run_case.sample_steps(256)[0]):
|
||||
row.update(step=step, tstar=step * U / 256,
|
||||
spurious_pressure_scaled_by_level=[row["spurious_pressure_scaled"]] * 2)
|
||||
return out
|
||||
|
||||
|
||||
class ScaledContractTests(unittest.TestCase):
|
||||
def setUp(self):
|
||||
import assess_shear_wave
|
||||
self.mod = assess_shear_wave
|
||||
|
||||
def test_scaled_contract_and_signed_excess(self):
|
||||
result = self.mod.assess(scaled_set(1.5e-3), "r", scaled=True)
|
||||
self.assertTrue(result["contract_valid"], [k for k, v in result["checks"].items() if not v])
|
||||
self.assertEqual([c[1] for c in result["cases"]], [256, 512, 256, 256])
|
||||
projected = result["signed"]["projected_linear"]
|
||||
self.assertAlmostEqual(projected["signed_rate_bias"], 1.5e-3, delta=1e-9)
|
||||
self.assertAlmostEqual(projected["excess_over_uniform_root"]["signed_rate_bias"], 1.5e-3,
|
||||
delta=1e-9)
|
||||
self.assertIn("excess_over_uniform_fine", projected)
|
||||
|
||||
def test_scaled_and_default_sets_do_not_mix(self):
|
||||
n128 = [receipt(128, 1), receipt(256, 1), receipt(128, 2), receipt(128, 2, "projected_linear")]
|
||||
self.assertEqual(self.mod.assess(n128, "r", scaled=True)["status"], "INVALID_SET")
|
||||
self.assertEqual(self.mod.assess(scaled_set(), "r")["status"], "INVALID_SET")
|
||||
unlabeled = scaled_set()
|
||||
unlabeled[0].pop("diagnostic")
|
||||
self.assertEqual(self.mod.assess(unlabeled, "r", scaled=True)["status"], "INVALID_SET")
|
||||
tg_label = scaled_set()
|
||||
tg_label[2]["diagnostic"] = run_case.SCALED_DIAGNOSTIC
|
||||
self.assertEqual(self.mod.assess(tg_label, "r", scaled=True)["status"], "INVALID_SET")
|
||||
|
||||
def test_resolution_trend_and_reference_pin(self):
|
||||
reference = [receipt(128, 1), receipt(256, 1), receipt(128, 2),
|
||||
receipt(128, 2, "projected_linear")]
|
||||
new = self.mod.assess(scaled_set(1.5e-3), "r", scaled=True)
|
||||
trend = self.mod.resolution_trend(new, reference, "r")
|
||||
self.assertTrue(trend["valid"])
|
||||
change = trend["changes"]["projected_linear"]["signed_rate_bias"]
|
||||
self.assertAlmostEqual(change["signed_change"], 1.5e-3, delta=1e-9)
|
||||
self.assertIn("excess_over_uniform_fine.residual_l2", trend["changes"]["constant"])
|
||||
with tempfile.TemporaryDirectory() as tmp:
|
||||
paths = []
|
||||
for i, r in enumerate(reference):
|
||||
paths.append(Path(tmp) / f"{i}.json")
|
||||
paths[-1].write_text(json.dumps(r))
|
||||
with self.assertRaises(ValueError):
|
||||
self.mod.load_reference(paths)
|
||||
|
||||
def test_scaled_runner_options(self):
|
||||
repo = Path(__file__).resolve().parents[3]
|
||||
with tempfile.TemporaryDirectory() as tmp:
|
||||
base = ["run_case.py", "--repo", str(repo), "--expected-revision", "f" * 40,
|
||||
"--output", str(Path(tmp) / "case.json"), "--dry-run", "--re", "100"]
|
||||
shear = ["--flow", "shear_wave"]
|
||||
scaled = ["--diagnostic", run_case.SCALED_SHEAR]
|
||||
for extra, valid in (
|
||||
(["--n", "512", "--depth", "1"] + shear + scaled, True),
|
||||
(["--n", "256", "--depth", "2", "--explosion", "projected_linear"] + shear + scaled, True),
|
||||
(["--n", "512", "--depth", "1"] + shear, False),
|
||||
(["--n", "128", "--depth", "2"] + shear + scaled, False),
|
||||
(["--n", "512", "--depth", "1"] + scaled, False), # needs shear flow
|
||||
(["--n", "512", "--depth", "1", "--early-pressure", "--diagnostic",
|
||||
run_case.SCALED_DIAGNOSTIC] + shear, False), # mutually exclusive
|
||||
(["--n", "256", "--depth", "2", "--early-pressure"] + shear + scaled, False)):
|
||||
argv = base + extra
|
||||
with self.subTest(argv=argv), patch.object(sys, "argv", argv), \
|
||||
patch.object(run_case, "verify_source", return_value={}), \
|
||||
patch.object(run_case, "runtime_source", return_value={}), \
|
||||
redirect_stdout(io.StringIO()), patch("sys.stderr", io.StringIO()):
|
||||
if valid:
|
||||
self.assertEqual(run_case.main(), 0)
|
||||
else:
|
||||
with self.assertRaises(SystemExit):
|
||||
run_case.main()
|
||||
|
||||
|
||||
class DefaultPathTests(unittest.TestCase):
|
||||
def test_frozen_oracle_is_byte_identical(self):
|
||||
digest = hashlib.sha256((HERE / "oracle.py").read_bytes()).hexdigest()
|
||||
self.assertEqual(digest, ORACLE_SHA256_EEC4677)
|
||||
|
||||
def test_default_flow_is_frozen_taylor_green(self):
|
||||
init, met, limits, text = run_case.flow_functions()
|
||||
self.assertIs(init, oracle.initial)
|
||||
self.assertIs(met, oracle.metrics)
|
||||
self.assertIs(limits, oracle.LIMITS)
|
||||
self.assertEqual(text, run_case.INITIALIZATION)
|
||||
self.assertEqual(run_case.flow_functions("taylor_green"), (init, met, limits, text))
|
||||
shear = run_case.flow_functions("shear_wave")
|
||||
self.assertIs(shear[0], sw.initial)
|
||||
self.assertIs(shear[1], sw.metrics)
|
||||
with self.assertRaises(ValueError):
|
||||
run_case.flow_functions("tg")
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
Reference in new issue
Block a user