Files
CelerisLab/tests/validation/multigrid/test_refinement_followup.py
ClaudeandClaude Opus 5.5 ee7074a53d tests: restore the negative checks masked by stale fixtures (WP-B groups B-E)
Test-side only; no production, manifest or pin change. Group A identity
checks (test_source_pin and the jets.py pins) are unchanged and stay
failing as recorded.

B test_refined_dfg: receipt/plan fixtures load the frozen manifest JSON
  that verify_source returns for a matching tree (frozen_manifest), and
  the two manifest-selection assertions call verify_source with only its
  working-tree hash comparison neutralised (selected_manifest), so the
  malformed/collision/observation-mode/checkout receipt checks run again.
  test_source_pin remains the single identity check.
C test_channel_root: int(owner)-1 (uint64-int became float on NumPy 1);
  test_shifted_storage: FP32 NaN fill so the arrays stay FP32. Both reach
  the intended immutable/representation/FP32/density checks.
D test_block_stepper fake gets _ddf_shifting; test_fine_surface_gpu sets
  _boundary_hosts/_level_jets instead of the read-only fine_boundary view;
  test_refinement_followup fake_backend sets _layouts/_ddf_shifting/
  _level_cfgs instead of the read-only fine_layout view (as reviewed in
  32da2d73).
E test_kernel_contracts_reject_before_upload: an inconsistent pair tau is
  now first rejected by the hierarchy chain check; the test asserts that
  rejection, then bypasses _check_chain to prove the per-level omega/tau
  guard still rejects, and asserts no mem_alloc/memcpy_htod in either case.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-21 07:04:47 +00:00

384 lines
21 KiB
Python

"""Bounded fixed-patch checks: fault guards, periodic geometry and coarse slot BC.
Use the existing bounded runner/environment; no audit-build rerun. The scalar
oracle covers the fixed tau pair and two boundary modes used by these cases.
"""
from contextlib import nullcontext
from dataclasses import replace
import json
from pathlib import Path
import resource
from threading import RLock
import unittest
from unittest.mock import MagicMock, Mock, patch
import numpy as np
import pycuda.driver as cuda
import CelerisLab
from CelerisLab.lbm import refinement_stepper as implementation
from test_refinement_gpu import (ROOT, POP_ATOL, MOMENT_ATOL, ref, postcollision,
totals, LBMConfig, CudaContext, compiler,
PeriodicPatchStepper)
def scalar_interval(coarse, fine, bounds, slot=None):
"""Independent G→G interval; optional (amplitude, width) bottom slot/top rho=1."""
_, ny, nx = coarse.shape
_, fy, fx = fine.shape
x0, y0, x1, y1 = bounds
cx, cy = (0, 1, -1, 0, 0, 1, -1, 1, -1), (0, 0, 0, 1, -1, 1, -1, -1, 1)
weights = (4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36)
def collide(values, tau):
rho = sum(values)
ux = sum(c * v for c, v in zip(cx, values)) / rho
uy = sum(c * v for c, v in zip(cy, values)) / rho
equilibrium = [w * rho * (1 + 3*(a*ux+b*uy) + 4.5*(a*ux+b*uy)**2
- 1.5*(ux*ux+uy*uy)) for w, a, b in zip(weights, cx, cy)]
return [v - (v-e)/tau for v, e in zip(values, equilibrium)]
current = fine.astype(np.float64).copy()
for y in range(fy):
for x in range(fx):
if not (4 <= x < fx-4 and 4 <= y < fy-4):
current[:, y, x] = coarse[:, y0+(y-4)//2, x0+(x-4)//2]
for _ in range(2):
out = np.full_like(current, np.nan)
for y in range(2, fy-2):
for x in range(2, fx-2):
values = [current[q, y-cy[q], x-cx[q]] for q in range(9)]
out[:, y, x] = (collide(values, 1.1)
if 4 <= x < fx-4 and 4 <= y < fy-4 else values)
current = out
result = np.full(coarse.shape, np.nan)
for y in range(ny):
for x in range(nx):
if (x0 <= x < x1 and y0 <= y < y1) or (slot is not None and y in (0, ny-1)):
continue
values = []
for q in range(9):
mx, my = x+.5-cx[q]/2, y+.5-cy[q]/2
if x0 <= mx <= x1 and y0 <= my <= y1:
a, b = 4+2*(x-x0), 4+2*(y-y0)
value = sum(current[q, b+dy, a+dx] for dy in (0, 1) for dx in (0, 1))/4
else:
value = coarse[q, (y-cy[q]) % ny, (x-cx[q]) % nx]
values.append(value)
if slot is not None and y == 1:
amplitude, width = slot
s = 2*((x+.5)-nx/2)/width
v = amplitude * max(0.0, 1-s*s)**2
rho = (values[0]+values[1]+values[2]
+ 2*(values[4]+values[6]+values[7]))/(1-v)
values[3] = values[4] + 2*rho*v/3
values[5] = values[6] + (values[2]-values[1])/2 + rho*v/6
values[8] = values[7] + (values[1]-values[2])/2 + rho*v/6
elif slot is not None and y == ny-2:
jy = values[0]+values[1]+values[2]+2*(values[3]+values[5]+values[8])-1
values[4] = values[3] - 2*jy/3
values[6] = values[5] + (values[1]-values[2])/2 - jy/6
values[7] = values[8] + (values[2]-values[1])/2 - jy/6
result[:, y, x] = collide(values, .8)
return result, current
class FaultGuardsCPU(unittest.TestCase):
def test_boundary_mode_and_margin_contract(self):
cfg = LBMConfig(nx=32, ny=24, viscosity=.1, velocity=0, wall_normal_velocity=.001,
wall_slot_width=8, y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
context = Mock()
context.activate.side_effect = RuntimeError("validated before CUDA")
for mode, bounds in (("unknown", (12, 4, 20, 10)),
("slot_reservoir", (12, 3, 20, 10)),
("slot_reservoir", (12, 4, 20, 21))):
with self.subTest(mode=mode, bounds=bounds):
with self.assertRaises(ValueError):
implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[bounds])
context.activate.assert_not_called()
# Both exact channel limits and the original periodic two-cell margin
# must reach context activation, with no allocation or device required.
for mode, bounds in (("slot_reservoir", (12, 4, 20, 20)),
("periodic", (2, 2, 6, 5))):
with self.subTest(mode=mode, bounds=bounds):
with self.assertRaisesRegex(RuntimeError, "validated before CUDA"):
implementation.FixedPatchStepper(cfg, context=context, coarse_boundary=mode,
refinements=[bounds])
def fake_backend(self):
backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper)
backend._layouts = [None, None] # dense storage on both levels
backend._ddf_shifting = False # unshifted on every level, as _check_storage expects
backend._level_cfgs = [LBMConfig(nx=16, ny=16, viscosity=.1)] * 2 # ddf_shifting off
backend.ctx = Mock(activate=nullcontext)
backend._lock = RLock()
backend._status = "ready"
backend._stream = Mock()
backend._levels = [Mock(shape=(9, 1, 1), layout=None), Mock(shape=(9, 1, 1), layout=None)]
backend._allocations = [Mock(), Mock()]
return backend
def test_failed_download_blocks_read_and_step(self):
# Fail on each level separately, including after coarse data was read.
for failed_read in (0, 1):
backend = self.fake_backend()
reads = [None] * failed_read + [RuntimeError("device read failed")]
with patch.object(cuda, "memcpy_dtoh", side_effect=reads) as read:
with self.assertRaisesRegex(RuntimeError, "device read failed"):
backend.download()
calls = read.call_count
for operation in (backend.download, backend.step):
with self.assertRaisesRegex(RuntimeError, "failed"):
operation()
self.assertEqual(read.call_count, calls)
backend.close()
def test_activation_failure_retains_resources_for_retry(self):
backend = self.fake_backend()
levels, allocations, stream = backend._levels[:], backend._allocations[:], backend._stream
unavailable = MagicMock()
unavailable.__enter__.side_effect = RuntimeError("activation failed")
backend.ctx.activate = Mock(side_effect=[unavailable, nullcontext()])
with self.assertRaisesRegex(RuntimeError, "activation failed"):
backend.close()
self.assertEqual(backend._status, "failed")
self.assertEqual(backend._levels, levels)
self.assertEqual(backend._allocations, allocations)
self.assertIs(backend._stream, stream)
for allocation in allocations:
allocation.free.assert_not_called()
backend.close()
self.assertEqual(backend._status, "closed")
for allocation in allocations:
allocation.free.assert_called_once()
def test_kernel_contracts_reject_before_upload(self):
cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
context = Mock(activate=nullcontext, sm_arch="sm_70")
with patch.object(implementation, "HALO", 2):
with self.assertRaisesRegex(RuntimeError, "four-cell"):
PeriodicPatchStepper(cfg, context=context, refinements=[(4, 4, 12, 12)])
p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12))
from CelerisLab.lbm import hierarchy
bad = patch.object(implementation.FixedPatch2D, "from_config",
return_value=replace(p, tau_coarse=.9))
# A pair tau inconsistent with its level config is rejected twice over, both before
# any device allocation or upload: first by the hierarchy chain check, and, if that
# is bypassed, by the per-level omega/tau check in _build_level.
with patch.object(cuda, "Stream"), bad, patch.object(cuda, "mem_alloc") as alloc, \
patch.object(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "config chain is inconsistent"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
with patch.object(hierarchy, "_check_chain"), \
self.assertRaisesRegex(RuntimeError, "omega/tau"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
alloc.assert_not_called()
upload.assert_not_called()
module = Mock()
module.get_global.return_value = (Mock(), 16)
with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \
patch.object(compiler, "compile_kernel"), patch.object(compiler, "cleanup_build_dir"), \
patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "layout mismatch"):
PeriodicPatchStepper(cfg, context=context, refinements=[p.bounds])
upload.assert_not_called()
class NormalBuildGPU(unittest.TestCase):
def setUp(self):
self.assertEqual(Path(CelerisLab.__file__).resolve(), ROOT / "src/CelerisLab/__init__.py")
self.assertEqual(Path(compiler._KERNEL_DIR).resolve(), ROOT / "src/CelerisLab/lbm/kernels")
self.assertEqual(Path(implementation.__file__).resolve(),
ROOT / "src/CelerisLab/lbm/refinement_stepper.py")
def tearDown(self):
rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
print(f"host_peak_rss={rss:.3f} MiB", flush=True)
def test_slot_reservoir_eight_steps(self):
nx, ny, bounds, amplitude, width = 32, 24, (12, 4, 20, 10), .001, 8
x0, y0, x1, y1 = bounds
fx, fy = 2*(x1-x0)+8, 2*(y1-y0)+8
cfg = LBMConfig(nx=nx, ny=ny, viscosity=.1, velocity=0, rho=1,
wall_normal_velocity=amplitude, wall_slot_width=width,
y_wall_bc="on_site_velocity", top_wall_bc="reservoir")
cx = np.array((0, 1, -1, 0, 0, 1, -1, 1, -1))
cy = np.array((0, 0, 0, 1, -1, 1, -1, -1, 1))
weights = np.array((4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36))
def slot_velocity(x):
return amplitude * np.maximum(0, 1-(2*(x-nx/2)/width)**2)**2
def equilibrium(shape, origin, dx):
ys = origin[1] + dx*(np.arange(shape[0])+.5)
xs = origin[0] + dx*(np.arange(shape[1])+.5)
v = slot_velocity(xs)[None, :] * (1-(ys[:, None]-1.5)/(ny-3))
cu = cy[:, None, None] * v
return np.array(weights[:, None, None]*(1+3*cu+4.5*cu*cu-1.5*v*v), np.float32)
coarse_leaf = np.ones((ny, nx), bool)
coarse_leaf[[0, -1]] = False
coarse_leaf[y0:y1, x0:x1] = False
fine_leaf = np.zeros((fy, fx), bool)
fine_leaf[4:-4, 4:-4] = True
initial = [equilibrium((ny, nx), (0, 0), 1),
equilibrium((fy, fx), (x0-2, y0-2), .5)]
initial[0][:, ~coarse_leaf] = np.nan
initial[1][:, ~fine_leaf] = np.nan
# A postcollision zero-mass/zero-momentum perturbation at a grazing
# source distinguishes closed-corner coupling from ordinary coarse pull.
initial[0][[5, 6], y0-1, x0] += 1e-5
initial[0][[7, 8], y0-1, x0] -= 1e-5
expected = tuple(a.astype(np.float64) for a in initial)
def discrete_mass(arrays):
return (arrays[0][:, coarse_leaf].sum(dtype=np.float64)
+ .25*arrays[1][:, fine_leaf].sum(dtype=np.float64))
# Observe the actual generated headers, without changing generation or
# compilation. In particular the fine module must have slot BC disabled.
generated = []
generate = compiler.generate_config
def record_config(*args, **kwargs):
build = generate(*args, **kwargs)
metadata = json.loads((Path(build)/".celeris-build").read_text())
definitions = {}
for line in (Path(build)/"config/config_method.h").read_text().splitlines():
parts = line.split()
if len(parts) == 3 and parts[0] == "#define":
definitions[parts[1]] = parts[2]
generated.append((metadata["config"], definitions))
return build
max_error = max_bc_error = 0.0
with CudaContext(0) as context:
with patch.object(compiler, "generate_config", side_effect=record_config):
backend = implementation.FixedPatchStepper(cfg, context=context, refinements=[bounds],
coarse_boundary="slot_reservoir")
with backend:
self.assertEqual(len(generated), 2)
fine_config, fine_defines = generated[1]
self.assertEqual((fine_config["nx"], fine_config["ny"]), (fx, fy))
for name in ("Y_WALL_BC", "TOP_WALL_BC", "WALL_SLOT_WIDTH",
"WALL_NORMAL_VELOCITY", "TOP_WALL_NORMAL_VELOCITY"):
self.assertEqual(float(fine_defines[name].removesuffix("f")), 0, name)
self.assertEqual(int(generated[0][1]["Y_WALL_BC"]), 2)
self.assertEqual(int(generated[0][1]["TOP_WALL_BC"]), 1)
# Independently derive all target links, including the eight
# whose coarse sources are uncovered, from physical midpoints.
targets, corners = set(), set()
for y, x in np.argwhere(coarse_leaf):
for q in range(1, 9):
if x0 <= x+.5-cx[q]/2 <= x1 and y0 <= y+.5-cy[q]/2 <= y1:
target = q*nx*ny+y*nx+x
targets.add(target)
sx, sy = x-cx[q], y-cy[q]
if not (x0 <= sx < x1 and y0 <= sy < y1):
corners.add(target)
self.assertEqual(len(corners), 8)
self.assertEqual(set(backend.patch.coalesce_coarse.tolist()), targets)
backend.initialize(*initial)
before = discrete_mass(initial)
for step in range(1, 9):
expected = scalar_interval(*expected, bounds, slot=(amplitude, width))
backend.step()
self.assertEqual(backend.step_count, step)
raw = [np.empty(a.shape, np.float32) for a in initial]
for level, array in zip(backend._levels, raw):
cuda.memcpy_dtoh(array, level.current)
# Read raw storage before masking: active NaNs cannot be
# hidden by the public download's non-leaf cleanup.
self.assertTrue(np.isfinite(raw[0][:, coarse_leaf]).all())
self.assertTrue(np.isfinite(raw[1][:, fine_leaf]).all())
self.assertTrue(np.isnan(raw[0][:, ~coarse_leaf]).all())
for got, want in zip(raw, expected):
np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
max_error = max(max_error, float(np.nanmax(np.abs(got-want))))
self.assertLessEqual(max_error, POP_ATOL)
public = backend.download()
np.testing.assert_array_equal(public[0], raw[0])
np.testing.assert_array_equal(public[1][:, fine_leaf], raw[1][:, fine_leaf])
self.assertTrue(np.isnan(public[1][:, ~fine_leaf]).all())
walls = raw[0][:, [1, ny-2], :].astype(np.float64)
rho = walls.sum(axis=0)
ux = np.einsum("q,qyx->yx", cx, walls)/rho
vy = np.einsum("q,qyx->yx", cy, walls)/rho
max_bc_error = max(max_bc_error, float(np.abs(ux).max()),
float(np.abs(vy[0]-slot_velocity(np.arange(nx)+.5)).max()),
float(np.abs(rho[1]-1).max()))
self.assertLessEqual(max_bc_error, 2e-6)
change = discrete_mass(raw)-before
bottom_flux, top_flux = (rho*vy).sum(axis=1)
print(f"slot/reservoir steps=8 max_population={max_error:.9g} "
f"max_BC_error={max_bc_error:.9g} discrete_leaf_mass_change={change:+.9g} "
f"final_sum_rho_v_bottom={bottom_flux:.9g} top={top_flux:.9g}; "
"raw halos/corners/fine BC isolation PASS", flush=True)
def test_nonequilibrium_seed1_eight_steps(self):
spec = ref.make_spec(collision="srt")
maps = ref.build_maps(spec)
initial = ref.make_state(spec, velocity=(.03, -.02))
rng = np.random.default_rng(1)
state = ref.State(*(a * (1 + 1e-3*rng.standard_normal(a.shape))
for a in (initial.coarse_f, initial.fine_f)))
max_pop = max_drift = max_moment = 0.0
with CudaContext(0) as context, PeriodicPatchStepper(
LBMConfig(nx=16, ny=16, viscosity=.1), context=context,
refinements=[(4, 4, 12, 12)]) as backend:
backend.initialize(*postcollision(state, spec, maps))
before = totals(backend.download(), maps)
for step in range(1, 9):
state, _ = ref.advance_coarse_once(spec, maps, state)
backend.step()
actual, expected = backend.download(), postcollision(state, spec, maps)
self.assertEqual(backend.step_count, step)
for got, want in zip(actual, expected):
np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
max_pop = max(max_pop, float(np.nanmax(np.abs(got-want))))
max_drift = max(max_drift, float(np.max(np.abs(totals(actual, maps)-before))))
max_moment = max(max_moment, float(np.max(np.abs(totals(actual, maps)-totals(expected, maps)))))
self.assertLessEqual(max_pop, POP_ATOL)
self.assertLessEqual(max_drift, MOMENT_ATOL)
self.assertLessEqual(max_moment, MOMENT_ATOL)
print(f"normal seed1 steps=8: max_population={max_pop:.9g} "
f"max_leaf_drift={max_drift:.9g} max_leaf_error={max_moment:.9g}", flush=True)
def test_offcenter_rectangle_margin_two_one_step(self):
cfg, bounds = LBMConfig(nx=17, ny=13, viscosity=.1), (2, 2, 6, 5)
rng = np.random.default_rng(7)
# Explicit asymmetric, nonequilibrium G0, not an arrival array to collide again.
initial = [ref.W[:, None, None] * (1+1e-3*rng.standard_normal((9, ny, nx)))
for ny, nx in ((13, 17), (14, 16))]
initial = [np.array(a, dtype=np.float32) for a in initial]
initial[0][:, 2:5, 2:6] = np.nan
fine_mask = np.zeros((14, 16), bool)
fine_mask[4:-4, 4:-4] = True
initial[1][:, ~fine_mask] = np.nan
expected = scalar_interval(*[a.astype(np.float64) for a in initial], bounds)
with CudaContext(0) as context, PeriodicPatchStepper(cfg, context=context, refinements=[bounds]) as backend:
self.assertEqual(backend.patch.fine_shape, (14, 16))
p = backend.patch
allocated = 2 * 9 * (17*13+16*14) * 4 + sum(a.nbytes for a in
(p.ghost_fine, p.ghost_coarse, p.coalesce_coarse, p.coalesce_fine))
self.assertLess(allocated, 128*1024)
backend.initialize(*initial)
backend.step()
coarse, fine = backend.download()
raw_fine = np.empty((9, 14, 16), np.float32)
cuda.memcpy_dtoh(raw_fine, backend._levels[1].current)
max_error = 0.0
for got, want in zip((coarse, raw_fine), expected):
np.testing.assert_array_equal(np.isnan(got), np.isnan(want))
max_error = max(max_error, float(np.nanmax(np.abs(got-want))))
self.assertLessEqual(max_error, POP_ATOL)
self.assertTrue(np.isnan(fine[:, ~fine_mask]).all())
print(f"normal rectangle 17x13 bounds={bounds} fine=16x14 steps=1: "
f"max_population={max_error:.9g} device_arrays={allocated}; halo/stride PASS", flush=True)
if __name__ == "__main__":
unittest.main(verbosity=2)