From 3552e6ee92696415491e756824887a3dbbbec492 Mon Sep 17 00:00:00 2001 From: TRON Date: Mon, 19 Oct 2026 05:27:18 -0700 Subject: [PATCH] test: compare fixed centered reconstruction on audited square faces --- docs/validation/taylor_green.md | 56 ++++ tests/validation/taylor_green/run_case.py | 64 +++-- .../validation/taylor_green/shear_centered.py | 229 ++++++++++++++++ .../validation/taylor_green/test_centered.py | 246 ++++++++++++++++++ 4 files changed, 577 insertions(+), 18 deletions(-) create mode 100644 tests/validation/taylor_green/shear_centered.py create mode 100644 tests/validation/taylor_green/test_centered.py diff --git a/docs/validation/taylor_green.md b/docs/validation/taylor_green.md index e9dabaf..4ec3a5c 100644 --- a/docs/validation/taylor_green.md +++ b/docs/validation/taylor_green.md @@ -548,3 +548,59 @@ bottom/top inner fine quartets its mean delta jx / d(uy)/dx is -/+ 1/(576*tau_f), with opposite adjacent-coarse redistribution. Keep raw values, step-zero rounding and fresh controls: sinusoidal curvature and finite-Mach effects prevent treating this coefficient as an exact gate. + +### Test-only centered straight-face comparison + +`centered_face_square_n128` reuses the square axis-wave POST-G initializer, +geometry, native probes and root steps 0–8. It has four fresh cases: U128, +U256, current projected R128, and experimental centered R128. This is a +qualification-ineligible diagnostic; there is no legacy assessor or new physics gate. + +Only the candidate's private generated fine-level `blocks.cu` changes. A +2,112-byte constant mask selects 1,792 rows (448 complete donor quartets). +The diagonal projected correction is doubled on straight faces with exactly +one unsupported normal axis and both legal tangent neighbors. Donor centers +within four root cells of a corner retain the current rule, as do all other +patterns, axial populations and rest. The mask is symmetric under transpose +and reflection. Its hash and the actual modified source/PTX hashes are retained; +production files, restrictions, collision, scheduling and launch arguments stay fixed. +A complete quartet has zero summed offset, so each population's quartet mean +is preserved. The perpendicular gradient also retains equal characteristic +packets; the existing two-hop conservation routing therefore remains applicable. +Host tests additionally check the mixed mask through complete two-hop/collision +cycles, including corner grazing, and execute the actual edited kernel via C++. +The central probes are more than 20 root cells from mask transitions; the +reviewed conservative dependency radius is 2.5 root cells per root interval. + +The separate comparison subtracts each group's step-zero residual and the +fresh uniform control at matching quadrature (U256 inner, U128 outer). It retains +signed density and both momentum components, actual-center gradient normalization, +parities, and sign flips. Step 1 requires bottom-inner cross-momentum magnitude +improvement greater than its deadband, and top-inner nonworsening within it. +Each horizontal face's step-8 RMS combines its inner and outer group means; +both must be nonworse. All four vertical groups must remain unchanged within +the component deadbands at every sampled root step. The norm +`sqrt(drho²/3 + djx² + djy²)` must not worsen beyond its corresponding deadband. +Sampled raw-pressure envelope and step-8 global gauged pressure RMS must not +worsen beyond the pressure deadband. Nonfinite/nonpositive sampled density, +invalid traces, source/build mismatch, or exceeded observed population envelope +reject the comparison. A small top signal or sign reversal remains explicit. + +The predeclared shifted population envelope is 1/32. From per-population oracle +allowance `6e-8 + 2e-7*1/32`, pairwise moment deadbands are +`[1.1925e-6, 7.95e-7, 7.95e-7]`; pressure uses the density deadband divided by 3. +Area averaging cannot enlarge this local allowance. These are practical diagnostic +deadbands, not certified accumulated root-8 rounding bounds. The envelope is +checked on downloaded synchronized ACTIVE populations, not hidden substeps. +Known quadratic fine-cell O(h) terms, corner transitions, repeated coupling and +stability remain unresolved. Passing permits review of a candidate, not promotion. + +Run the existing `run_case.py` with `--diagnostic centered_face_square_n128 +--flow shear_wave --re 100 --probe-trace --test-reconstruction current` for +U128/U256 (depth 1, constant) and current R128 (depth 2, projected_linear). +For candidate R128 use `--test-reconstruction centered_face_v1`. All require +an exact clean `--repo`, `--expected-revision`, fresh `--output`, and either +`--dry-run` or explicit single-device `--run-gpu`. Compare ordered receipts with +`shear_centered.py --repo REPO --expected-revision SHA --output comparison.json +U128.json U256.json current.json candidate.json`. Total work is 1,900,544 ACTIVE +updates, six build sidecars, and 720 cell-time samples; scalar cap stays 1 MiB/case. diff --git a/tests/validation/taylor_green/run_case.py b/tests/validation/taylor_green/run_case.py index 1ef60fa..46ada01 100644 --- a/tests/validation/taylor_green/run_case.py +++ b/tests/validation/taylor_green/run_case.py @@ -1,6 +1,7 @@ """One bounded GPU case, with explicit opt-in experimental explosion.""" import argparse +from contextlib import nullcontext import hashlib import json import os @@ -29,6 +30,8 @@ WIDE_X_SHEAR = "wide_x_shear_n128" # one projected x-gap-width contrast SHEAR_WIDE_X_SHA256 = "13d16ba355117e2d074edccd9dc5d694c2fc13c673a0abdfca6401f6ed917733" AXIS_SQUARE = "axis_square_shear_n128" # short four-face orientation diagnostic SHEAR_AXIS_SQUARE_SHA256 = "0baffa6eb99c268ce998d90300c993d45bb3b7b27c524b636874d75808bd0a8c" +CENTERED_SQUARE = "centered_face_square_n128" # private-build reconstruction comparison +SHEAR_CENTERED_SHA256 = "78bd7360e266e5e28ef4f03496413375551fcdc10ceba61d43bd88c33b4a8f5c" POSTG_SHEAR = "postg_shear_n128" # initialization-only three-case counterfactual SHEAR_POSTG_SHA256 = "9c909fa013c092490ee21532e3125b2a78a4c5102c7a1bdc0e76986d3f0ef903" SHEAR_PROBES_SHA256 = "e2e6490d281cc275f448dc6d6ea7b9b4407cc74a51e9256b86c08c8626b89b76" # signed interface-probe mechanism diagnostic @@ -127,7 +130,14 @@ def verify_source(repo, revision, flow="taylor_green", probes=False): or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_PROBES_SHA256): raise RuntimeError("Wrong shear_probes module checkout or content") result["shear_probes_sha256"] = SHEAR_PROBES_SHA256 - if flow == "shear_axis_square": + if flow == "shear_centered": + import shear_centered + module = harness / "shear_centered.py" + if (Path(shear_centered.__file__).resolve() != module + or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_CENTERED_SHA256): + raise RuntimeError("Wrong shear_centered module checkout or content") + result["shear_centered_sha256"] = SHEAR_CENTERED_SHA256 + if flow in ("shear_axis_square", "shear_centered"): import shear_axis_square module = harness / "shear_axis_square.py" if (Path(shear_axis_square.__file__).resolve() != module @@ -162,7 +172,7 @@ def verify_source(repo, revision, flow="taylor_green", probes=False): or hashlib.sha256(module.read_bytes()).hexdigest() != SHEAR_NESTED_SHA256): raise RuntimeError("Wrong shear_nested module checkout or content") result["shear_nested_sha256"] = SHEAR_NESTED_SHA256 - if flow in ("shear_wave", "shear_nested", "shear_twin", "shear_wide_x", "shear_postg", "shear_axis_square"): + if flow in ("shear_wave", "shear_nested", "shear_twin", "shear_wide_x", "shear_postg", "shear_axis_square", "shear_centered"): import shear_wave module = harness / "shear_wave.py" if (Path(shear_wave.__file__).resolve() != module @@ -237,7 +247,7 @@ def main(): p.add_argument("--explosion", choices=["constant", "projected_linear"], default="constant") p.add_argument("--diagnostic", choices=[SCALED_DIAGNOSTIC, SCALED_SHEAR, NESTED_SHEAR, TWIN_SHEAR, - WIDE_X_SHEAR, POSTG_SHEAR, AXIS_SQUARE], + WIDE_X_SHEAR, POSTG_SHEAR, AXIS_SQUARE, CENTERED_SQUARE], help="explicit four-case N256/fine512 same-placement diagnostic " "(scaled_n256: Taylor-Green; scaled_shear_n256: shear wave) or the " "fixed three-level/two-region shear preset (nested_shear_n128) " @@ -245,6 +255,8 @@ def main(): "x-gap-width contrast (wide_x_shear_n128), or three-case POST-G " "initialization counterfactual (postg_shear_n128), or short square " "orientation diagnostic (axis_square_shear_n128)") + p.add_argument("--test-reconstruction", choices=["current", "centered_face_v1"], + help="required only for centered_face_square_n128; private test build") p.add_argument("--probe-trace", action="store_true", help="signed interface-probe trace (nested U128/U512/projected d3, " "twin projected d2, wide-x projected d3, POST-G three-case set, " @@ -256,7 +268,16 @@ def main(): 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() - axis_square = args.diagnostic == AXIS_SQUARE + centered = args.diagnostic == CENTERED_SQUARE + if centered != (args.test_reconstruction is not None): + p.error("--test-reconstruction is required only for centered_face_square_n128") + if centered: + import shear_centered + if (args.n, args.depth, args.test_reconstruction) not in shear_centered.CASES: + p.error("centered comparison requires one of its four fixed cases") + if args.explosion != ("constant" if args.depth == 1 else "projected_linear"): + p.error("centered comparison uses constant uniforms and projected refined cases") + axis_square = args.diagnostic in (AXIS_SQUARE, CENTERED_SQUARE) nested = args.diagnostic == NESTED_SHEAR wide_x = args.diagnostic == WIDE_X_SHEAR postg = args.diagnostic == POSTG_SHEAR @@ -273,7 +294,7 @@ def main(): (POSTG_SHEAR, 128, 3, "projected_linear")} if axis_square: import shear_axis_square - probe_cases.update((AXIS_SQUARE, n, depth, scheme) + probe_cases.update((args.diagnostic, n, depth, scheme) for _, n, depth, scheme in shear_axis_square.CASES) if args.probe_trace and (args.final_snapshot or ( args.diagnostic, args.n, args.depth, args.explosion) not in probe_cases): @@ -325,7 +346,7 @@ def main(): 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_flow = ("shear_axis_square" if axis_square else "shear_postg" if postg else "shear_wide_x" if wide_x else "shear_nested" if nested + source_flow = ("shear_centered" if centered else "shear_axis_square" if axis_square else "shear_postg" if postg else "shear_wide_x" if wide_x else "shear_nested" if nested else "shear_twin" if twin else args.flow) source = verify_source(repo, args.expected_revision, source_flow, probes=args.probe_trace) runtime = runtime_source(repo) @@ -401,6 +422,10 @@ def main(): report.update(initialization_id=preset.INITIALIZATION_ID, qualification_eligible=False, reference_harness=None, historical_baseline_revision=None, screening_scope="no numerical gates; short orientation diagnostic") + if centered: + report.update(test_reconstruction=args.test_reconstruction, observed_shifted_max=0., + reconstruction=dict(variant="current"), + screening_scope="predeclared centered-witness diagnostic; separate four-case comparison") args.output.parent.mkdir(parents=True, exist_ok=True) def save(): temporary = args.output.with_suffix(args.output.suffix + ".tmp") @@ -415,18 +440,19 @@ def main(): with CudaContext(0) as ctx, build_capture(report, save): report["runtime_identity"] = runtime_identity(ctx) save() - if args.depth == 1: - layout = BlockLayout2D.build(args.n, args.n, 8, 8, periodic=(True, True)) - run = PeriodicBlockStepper(cfg, layout, context=ctx) - else: - run = FixedPatchStepper( - cfg, - context=ctx, - coarse_boundary="periodic", - refinements=refinements, - level_blocks=[(8, 8)] * args.depth, - explosion=args.explosion, - ) + with (shear_centered.private_build(h, report, args.test_reconstruction) if centered else nullcontext()): + if args.depth == 1: + layout = BlockLayout2D.build(args.n, args.n, 8, 8, periodic=(True, True)) + run = PeriodicBlockStepper(cfg, layout, context=ctx) + else: + run = FixedPatchStepper( + cfg, + context=ctx, + coarse_boundary="periodic", + refinements=refinements, + level_blocks=[(8, 8)] * args.depth, + explosion=args.explosion, + ) with run: run.initialize( *[ @@ -452,6 +478,8 @@ def main(): if "error" not in trace: arrays = (run.download(),) if args.depth == 1 else run.download() try: + if centered: + report["observed_shifted_max"] = max(report["observed_shifted_max"], shear_centered.observe(arrays)) values = shear_probes.sample(arrays, probes) shear_probes.validate_row(values, len(probes)) except (ValueError, ZeroDivisionError) as error: diff --git a/tests/validation/taylor_green/shear_centered.py b/tests/validation/taylor_green/shear_centered.py new file mode 100644 index 0000000..698359d --- /dev/null +++ b/tests/validation/taylor_green/shear_centered.py @@ -0,0 +1,229 @@ +"""Test-only straight-face centered witness; private CUDA snapshot, frozen source intact.""" +from contextlib import contextmanager +import hashlib +import json +from pathlib import Path +import numpy as np + +import shear_axis_square as axis +from oracle import CX, CY, W + +DIAGNOSTIC = "centered_face_square_n128" +VARIANT = "centered_face_v1" +BUILD_SHA256 = "81e880017495190209bd867d4b523a2c243be496ed4681892e4cd789b7ba1f21" +CASES = ((128, 1, "current"), (256, 1, "current"), + (128, 2, "current"), (128, 2, VARIANT)) +# Local arithmetic-derived diagnostic deadbands, NOT accumulated trajectory bounds. +POP_ENVELOPE, OP_ATOL, OP_RTOL = 1/32, 6e-8, 2e-7 +MARGIN = 2 * np.array([9, 6, 6]) * (OP_ATOL + OP_RTOL * POP_ENVELOPE) + + +def routing(hierarchy): + """Square-only mask: full donor quartets, one missing normal axis, legal tangent.""" + from CelerisLab.lbm.blocks import BlockLayout2D + from CelerisLab.lbm.refinement_tables import projected_explosion_tables + if hierarchy.depth != 2 or hierarchy.kinds[0].shape != (128, 128): + raise ValueError("centered witness requires the fixed two-level square") + pair = hierarchy.pairs[0] + covered = np.zeros((128, 128), bool); covered[32:96, 32:96] = True + if pair.fine_origin != (0., 0.) or not np.array_equal(pair.covered, covered): + raise ValueError("wrong centered-witness geometry") + layout = BlockLayout2D.build(128, 128, 8, 8, periodic=(True, True)) + neighbors, offsets = projected_explosion_tables(pair, layout, hierarchy.kinds[0]) + donors = layout.owner.ravel()[pair.ghost_coarse] + off = neighbors.reshape(-1, 2, 2) == donors[:, None, None] + disabled = off.all(axis=2) + y, x = np.divmod(pair.ghost_coarse, 128) + vertical = np.isin(x, [30, 31, 96, 97]) & (y >= 36) & (y < 92) + horizontal = np.isin(y, [30, 31, 96, 97]) & (x >= 36) & (x < 92) + mask = ((vertical & disabled[:, 0] & ~off[:, 1].any(axis=1)) | + (horizontal & disabled[:, 1] & ~off[:, 0].any(axis=1))).astype(np.uint8) + for donor in np.unique(pair.ghost_coarse): + rows = pair.ghost_coarse == donor + if rows.sum() != 4 or len(set(mask[rows])) != 1 or sorted(map(tuple, offsets[rows])) != [(-.25,-.25), (-.25,.25), (.25,-.25), (.25,.25)]: + raise ValueError("routing must preserve complete +/- quarter donor quartets") + return mask + + +def metadata(mask): + return dict(variant=VARIANT, mask_sha256=hashlib.sha256(mask.tobytes()).hexdigest(), + ghost_rows=len(mask), selected_rows=int(mask.sum()), + rule="diagonal projected correction x2; audited straight faces only") + + +def patch_source(source, mask): + """Change only the projected function in the private generated build directory.""" + anchor = 'extern "C" __global__ void BlockExplodeProjected(' + position = " to[index_f(t % TOTAL_CELLS, q)] =" + before, kernel = source.split(anchor) + body, after = kernel.split("// Targets are receiver-local SOURCE slots", 1) + if body.count(position) != 1: + raise ValueError("unexpected projected kernel source") + flag = "__device__ __constant__ unsigned char centered_face_mask[] = {" + ",".join(map(str, mask)) + "};\n" + edit = " if (centered_face_mask[row] && ex != 0.f && ey != 0.f) { gx *= 2.f; gy *= 2.f; }\n" + return before + flag + anchor + body.replace(position, edit + position) + "// Targets are receiver-local SOURCE slots" + after + + +def source_digest(build): + digest = hashlib.sha256() + for path in sorted(Path(build).rglob("*")): + if path.is_file() and path.name != ".celeris-build": + digest.update(str(path.relative_to(build)).encode()); digest.update(path.read_bytes()) + return digest.hexdigest() + + +@contextmanager +def private_build(hierarchy, report, variant): + """No kernel upload/launch monkey-patch; same production function and arguments.""" + if variant == "current": + yield + return + from CelerisLab.cuda import compiler_v2 as compiler + mask = routing(hierarchy) + report["reconstruction"] = metadata(mask) + original, calls = compiler.generate_config, [] + def generate(cfg, *args, **kwargs): + build = Path(original(cfg, *args, **kwargs)) + calls.append(cfg.viscosity) + if len(calls) == 2: + if calls != [.0384, .0768]: + raise ValueError("unexpected candidate level configuration") + path = build / "blocks.cu" + path.write_text(patch_source(path.read_text(), mask)) + sidecar = build / ".celeris-build" + data = json.loads(sidecar.read_text()) + data.update(source_sha256=source_digest(build), experimental_reconstruction=metadata(mask)) + sidecar.write_text(json.dumps(data, sort_keys=True)) + elif len(calls) > 2: + raise ValueError("candidate is one adjacent pair only") + return str(build) + compiler.generate_config = generate + try: + yield + if calls != [.0384, .0768]: + raise ValueError("candidate did not build its exact two levels") + finally: + compiler.generate_config = original + + +def observe(arrays): + """Bound observed synchronized ACTIVE values only, not intermediate GPU states.""" + maximum = max(float(np.nanmax(abs(a - W[:, None, None]))) for a in arrays) + if not np.isfinite(maximum) or maximum > POP_ENVELOPE: + raise ValueError("observed shifted-population envelope exceeded") + return maximum + + +def compare(records, expected_source): + """Four fresh receipts; diagnostic screens only, never a physical qualification.""" + import run_case as rc + from assess_projected_pairs import BUILD_SOURCES, NUMERICS, valid_builds, finite_tree + from CelerisLab.lbm.hierarchy import build_hierarchy + from CelerisLab.lbm.refinement import ACTIVE + if not finite_tree(records) or not expected_source.get("revision"): + raise ValueError("nonfinite receipt or missing exact source identity") + keys = [(r['n'], r['depth'], r.get('test_reconstruction')) for r in records] + if keys != list(CASES): + raise ValueError("expected ordered fresh U128/U256/current/candidate cases") + for r in records: + for key in ('source', 'runtime_identity', 'runtime_source', 'cuda_visible_devices'): + if not r.get(key) or r[key] != records[0][key]: + raise ValueError(f"mixed or missing {key}") + if (r['source'] != expected_source or r.get('numerical_revision') != rc.NUMERICAL_REVISION + or r.get('numerics') != NUMERICS or r.get('re_L') != 100 or r.get('U') != .03 + or r.get('polarization') != [0., 1.] or r.get('initialization_id') != axis.INITIALIZATION_ID + or r['source'].get('dirty') is not False or r['source'].get('source_tree_sha256') != rc.SOURCE_SHA256 + or r.get('diagnostic') != DIAGNOSTIC or not r.get('completed') + or r.get('qualification_eligible') is not False or r.get('initialization') != axis.INITIALIZATION + or r.get('mode') != [2, 0] or r.get('phase') != axis.PHI + or r.get('explosion') != ('constant' if r['depth'] == 1 else 'projected_linear') + or not 0 <= r.get('observed_shifted_max', np.inf) <= POP_ENVELOPE): + raise ValueError("wrong candidate receipt identity or observed envelope") + h = build_hierarchy(rc.config(r['n'], 100), [] if r['depth'] == 1 else axis.REFINEMENTS) + if r.get('refinements') != ([] if r['depth'] == 1 else [[list(axis.REFINEMENTS[0][0])]]): + raise ValueError("wrong geometry") + if (r.get('viscosity') != .03*r['n']/100 + or r.get('tau_by_level') != [.5+3*c.viscosity for c in h.configs] + or r.get('workload') != rc.workload([k == ACTIVE for k in h.kinds], 8*(r['n']//128))): + raise ValueError("wrong physical scales or workload") + builds = r.get('actual_build_sidecars', []) + canonical = [] + for level, build in enumerate(builds): + modified = r['test_reconstruction'] == VARIANT and level == 1 + expected = BUILD_SHA256 if modified else BUILD_SOURCES[r['n']*2**level] + if (build.get('source_sha256') != expected + or build.get('experimental_reconstruction') != (metadata(routing(h)) if modified else None)): + raise ValueError("wrong actual build or reconstruction identity") + # Reuse the frozen build-metadata checker after checking the private-source hash above. + canonical.append(dict(build, source_sha256=BUILD_SOURCES[r['n']*2**level])) + if not valid_builds(dict(actual_build_sidecars=canonical), r['n'], r['depth']): + raise ValueError("invalid architecture, flags, config or PTX digest") + if r.get('reconstruction') != (metadata(routing(h)) if r['test_reconstruction'] == VARIANT else dict(variant='current')): + raise ValueError("wrong receipt mask identity") + trace = r['probe_trace'] + declared = axis.new_trace(trace['probes'], r['viscosity'], [DIAGNOSTIC, r['n'], r['depth'], r['explosion']]) + if any(trace.get(k) != v for k,v in declared.items() if k not in ('steps','values')): + raise ValueError("wrong trace phase, clock or physical identity") + if trace['probes'] != axis.select([k == ACTIVE for k in h.kinds], r['n']): + raise ValueError("wrong actual probe cells") + if trace.get("bytes") != axis.validate(trace): + raise ValueError("trace serialized byte count mismatch") + if np.min(np.asarray(trace["values"])[:,0::3]) <= 0: + raise ValueError("nonpositive sampled density") + if [row['step'] for row in r['rows']] != [0, 8*(r['n']//128)] or not all(rc.finite_row(v) for v in r['rows']): + raise ValueError("invalid scalar rows") + if min(row['rho_min'] for row in r['rows']) <= 0: + raise ValueError("nonpositive density") + current, candidate = records[2:] + if current['probe_trace']['values'][0] != candidate['probe_trace']['values'][0]: + raise ValueError("candidate/current initial samples differ; cancellation invalid") + vectors = [] + for r in records: + groups = axis.group_moments(r['probe_trace']) + vectors.append(np.array([[[row['groups'][name][key] for key in ('drho','djx','djy')] + for name in axis.GROUPS] for row in groups])) + # Inner quadrature matches U256; outer quadrature matches U128. + uniform = vectors[1].copy(); uniform[:, 4:] = vectors[0][:, 4:] + residual = [v-v[0]-uniform+uniform[0] for v in vectors[2:]] + old, new = residual + checks = dict(bottom_step1_clear=abs(old[1,2,1])-abs(new[1,2,1]) > MARGIN[1], + top_step1_nonworse=abs(new[1,3,1]) <= abs(old[1,3,1])+MARGIN[1]) + for face, indices in [('bottom',[2,6]), ('top',[3,7])]: + rms = lambda a: float(np.sqrt(np.mean(a[8,indices,1]**2))) + checks[face+'_step8_nonworse'] = rms(new) <= rms(old)+MARGIN[1] + checks['vertical_null'] = bool(np.all(abs(new[:,[0,1,4,5]]-old[:,[0,1,4,5]]) <= MARGIN)) + acoustic = lambda a: np.sqrt(a[...,0]**2/3 + a[...,1]**2 + a[...,2]**2) + acoustic_margin = float(np.sqrt(MARGIN[0]**2/3 + MARGIN[1]**2 + MARGIN[2]**2)) + checks['acoustic_nonworse'] = bool(np.all(acoustic(new) <= acoustic(old)+acoustic_margin)) + raw_pressure = lambda r: float(np.max(abs(np.asarray(r['probe_trace']['values'])[:,0::3]-1))/3) + checks['sampled_pressure_nonworse'] = raw_pressure(candidate) <= raw_pressure(current)+MARGIN[0]/3 + checks['step8_pressure_nonworse'] = candidate['rows'][-1]['spurious_pressure_rms'] <= current['rows'][-1]['spurious_pressure_rms']+MARGIN[0]/3 + return dict(qualified=False, advance_to_review=all(checks.values()), checks=checks, + moment_deadband=MARGIN.tolist(), deadband_scope="local arithmetic-derived diagnostic; no accumulated error guarantee", + groups=list(axis.GROUPS), components=['drho','djx','djy'], + current_signed_residual=old.tolist(), candidate_signed_residual=new.tolist(), + raw_current=vectors[2].tolist(), raw_candidate=vectors[3].tolist(), + group_moments=[axis.group_moments(r['probe_trace']) for r in records], + significant_sign_flips=((old*new < 0) & (abs(old)>MARGIN) & (abs(new)>MARGIN)).tolist(), + top_step1_below_deadband=bool(abs(old[1,3,1]) <= MARGIN[1]), + note="Absolute reduction with sign retained; top below deadband may be inconclusive.") + + +def main(): + import argparse + p = argparse.ArgumentParser(description=__doc__) + p.add_argument("--repo", required=True, type=Path) + p.add_argument("--expected-revision", required=True) + p.add_argument("receipts", nargs=4, type=Path, help="U128 U256 current candidate, in order") + p.add_argument("--output", required=True, type=Path) + args = p.parse_args() + import run_case as rc + source = rc.verify_source(args.repo.resolve(), args.expected_revision, "shear_centered", probes=True) + result = compare([json.loads(path.read_text()) for path in args.receipts], source) + with args.output.open("x") as stream: + json.dump(result, stream, indent=2, allow_nan=False) + return 0 if result["advance_to_review"] else 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/validation/taylor_green/test_centered.py b/tests/validation/taylor_green/test_centered.py new file mode 100644 index 0000000..19764d7 --- /dev/null +++ b/tests/validation/taylor_green/test_centered.py @@ -0,0 +1,246 @@ +"""Independent coordinate oracle, exact kernel host shim and closed pair budgets.""" +import ctypes +from contextlib import nullcontext +from copy import deepcopy +from dataclasses import replace +import hashlib +import json +from pathlib import Path +import subprocess +import tempfile +import unittest +from unittest.mock import patch + +import numpy as np +import run_case as rc +import shear_axis_square as axis +import shear_centered as candidate +import test_shear_wave as legacy +from test_axis_square import masks +from CelerisLab.lbm.hierarchy import build_hierarchy +from CelerisLab.lbm.blocks import BlockLayout2D +from CelerisLab.lbm.refinement_tables import projected_explosion_tables +from CelerisLab.lbm.refinement import ACTIVE, OUTER_GHOST +from tests.unit import test_projected_explosion as old + +REPO = Path(__file__).resolve().parents[3] +SOURCE = REPO/'src/CelerisLab/lbm/kernels/blocks.cu' + + +def oracle(c, pair, centered=True): + """Uses only dense coordinates/legal donors, never production neighbor/mask arrays.""" + out = np.empty((9, len(pair.ghost_fine))) + legal = set(map(int, pair.ghost_coarse)) + for r, (donor, target) in enumerate(zip(pair.ghost_coarse, pair.ghost_fine)): + y, x = divmod(int(donor), 128); fy, fx = divmod(int(target), 256) + delta = np.array([(fx+.5)/2-x-.5, (fy+.5)/2-y-.5]) + usable = [(128*(y-dy)+x-dx in legal and 128*(y+dy)+x+dx in legal) + for dx,dy in [(1,0),(0,1)]] + straight = ((x in (30,31,96,97) and 36 <= y < 92 and usable == [False,True]) + or (y in (30,31,96,97) and 36 <= x < 92 and usable == [True,False])) + for q, e in enumerate(old.E): + grad = np.array([(c[q,y+dy,x+dx]-c[q,y-dy,x-dx])/2 if ok else 0. + for (dx,dy),ok in zip([(1,0),(0,1)],usable)]) + if q: grad -= e*(grad@e)/(e@e) + out[q,r] = c[q,y,x] + (2 if centered and straight and np.all(e) else 1)*(grad@delta) + return out + + +class CenteredTests(unittest.TestCase): + @classmethod + def setUpClass(cls): + cls.h = build_hierarchy(rc.config(128,100), axis.REFINEMENTS) + cls.p = cls.h.pairs[0] + cls.mask = candidate.routing(cls.h) + cls.layout = BlockLayout2D.build(128,128,8,8,periodic=(True,True)) + cls.ids, cls.offsets = projected_explosion_tables(cls.p,cls.layout,cls.h.kinds[0]) + cls.neighbors = old.decode(cls.layout,cls.ids) + + def test_mask_constants_quartets_and_transpose(self): + self.assertEqual(candidate.metadata(self.mask), dict(variant=candidate.VARIANT, + mask_sha256='112de63ce9bafd2e205a9e7029ef6ce5683fe14407b6ac04444579c84960b7ea', + ghost_rows=2112, selected_rows=1792, + rule='diagonal projected correction x2; audited straight faces only')) + c=np.random.default_rng(47).normal(0,.001,(9,128,128)) + actual=oracle(c,self.p); baseline=oracle(c,self.p,False) + base,g=old.reconstruct(c,self.p,self.neighbors,self.offsets) + g[np.ix_(np.flatnonzero(np.all(old.E,axis=1)),np.flatnonzero(self.mask),[0,1])]*=2 + expected=c.reshape(9,-1)[:,self.p.ghost_coarse]+np.einsum('qra,ra->qr',g,self.offsets) + np.testing.assert_allclose(actual,expected,atol=1e-18,rtol=0) + np.testing.assert_array_equal(actual[:,~self.mask.astype(bool)],baseline[:,~self.mask.astype(bool)]) + for d in np.unique(self.p.ghost_coarse): + rows=self.p.ghost_coarse==d + np.testing.assert_allclose(actual[:,rows].mean(axis=1),c.reshape(9,-1)[:,d],atol=1e-18) + const=np.broadcast_to(old.W[:,None,None],c.shape) + np.testing.assert_array_equal(oracle(const,self.p),np.broadcast_to(old.W[:,None],actual.shape)) + qswap=[list(map(tuple,old.E)).index(tuple(e[::-1])) for e in old.E] + index={int(t):r for r,t in enumerate(self.p.ghost_fine)} + swap=[index[int(t%256)*256+int(t//256)] for t in self.p.ghost_fine] + np.testing.assert_array_equal(self.mask,self.mask[swap]) + np.testing.assert_allclose(oracle(c[qswap].transpose(0,2,1),self.p),actual[qswap][:,swap],atol=1e-18) + qreflect=[list(map(tuple,old.E)).index(tuple((-e[0],e[1]))) for e in old.E] + reflect=[index[int(t//256)*256+255-int(t%256)] for t in self.p.ghost_fine] + np.testing.assert_array_equal(self.mask,self.mask[reflect]) + np.testing.assert_allclose(oracle(c[qreflect,:,::-1],self.p),actual[qreflect][:,reflect],atol=1e-18) + # Central dependency corridors remain in a constant mask regime beyond radius8*2.5. + for name,(x,y) in axis.GROUPS.items(): + tangent = y if name.endswith(('left','right')) else x + self.assertGreater(min(tangent+.25-36,92-(tangent+.75)),20) + + def test_actual_generated_kernel_shifted_unshifted(self): + child=self.p.fine_layout((8,4)) # unequal owner strides, unused slots poisoned + patched=candidate.patch_source(SOURCE.read_text(),self.mask) + kernel=patched.split('__device__ __constant__ unsigned char centered_face_mask',1)[1] + kernel='unsigned char centered_face_mask'+kernel.split('// Targets are receiver-local SOURCE slots',1)[0] + kernel=kernel.replace('extern "C" __global__','extern "C"') + code='#include \nusing fpxx=float;\n' + code+=f'constexpr unsigned NQ=9,TOTAL_CELLS={np.prod(child.local_shape)};\n' + code+='struct Id{unsigned x;}; Id blockIdx{0},blockDim{1},threadIdx{0};\n' + code+='unsigned index_f(unsigned k,unsigned q){return q*TOTAL_CELLS+k;}\n' + code+='int d_cx[]={0,1,-1,0,0,1,-1,1,-1},d_cy[]={0,0,0,1,-1,1,-1,-1,1};\n'+kernel + original=SOURCE.read_text().split('extern "C" __global__ void BlockExplodeProjected(',1)[1].split('// Targets are receiver-local SOURCE slots',1)[0] + code+='void Original('+original + code+='extern "C" void invoke(float*c,float*f,uint64_t*t,uint64_t*d,unsigned n,unsigned long l,uint64_t*b,float*o){for(unsigned r=0;r