"""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, bounds, context=context, coarse_boundary=mode) 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, bounds, context=context, coarse_boundary=mode) def fake_backend(self): backend = PeriodicPatchStepper.__new__(PeriodicPatchStepper) backend.fine_layout = None # dense fine storage backend.ctx = Mock(activate=nullcontext) backend._lock = RLock() backend._status = "ready" backend._stream = Mock() backend._levels = [Mock(shape=(9, 1, 1)), Mock(shape=(9, 1, 1))] 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, (4, 4, 12, 12), context=context) p = implementation.FixedPatch2D.from_config(cfg, (4, 4, 12, 12)) with patch.object(cuda, "Stream"), patch.object(implementation.FixedPatch2D, "from_config", return_value=replace(p, tau_coarse=.9)): with self.assertRaisesRegex(RuntimeError, "omega/tau"): PeriodicPatchStepper(cfg, p.bounds, context=context) 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, p.bounds, context=context) 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, bounds, context=context, 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), (4, 4, 12, 12), context=context) 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, bounds, context=context) 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)