ClaudeandClaude Opus 5.5
19db7d362f
Restores the PeriodicBlockStepper safeguard dropped in 8e1feff, for all
FixedPatchStepper levels: both G buffers, halo pairs and boundary columns
of every level are compared with cuda.mem_get_info() before the stream,
any module or any buffer exists, raising MemoryError('level buffers exceed
free device memory'). Transfer tables are not counted (lower bound).
Tests: the single-level MemoryError case is back beside the param-layout
check (no stream, no allocation); host test_single_level covers depth 1
and depth 2 (no stream, build or allocation) and the wait failure after a
good step (count stays 1). Driver-mocking fixtures report ample free memory.
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
384 lines
21 KiB
Python
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"), patch.object(cuda, "mem_get_info", return_value=(1 << 40, 1 << 40)), 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(cuda, "mem_get_info", return_value=(1 << 40, 1 << 40)), 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)
|