Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

lbm: delete PeriodicBlockStepper; a single blocked level is FixedPatchStepper(refinements=[]) (S2)

Browse Source 
The second owner of one blocked level goes. Its capability is the
zero-refinement FixedPatchStepper added in f764d68/985db14 (periodic SRT
level, stationary passive circles, shifted storage, NaN spare, single wait),
which matched it array for array on GPU at 1318fc2 (3 splits x +-circle x
shifted/unshifted) and passes the host lifecycle and fault contracts.

API change for callers:
  PeriodicBlockStepper(cfg, BlockLayout2D.build(ny, nx, by, bx), context=c, objects=o)
  -> FixedPatchStepper(cfg, context=c, refinements=[], level_blocks=[(by, bx)],
                       level_objects=[tuple(o)])
  download() returns a 1-tuple; static-circle tables are
  _level_boundaries[0].tables (was .boundary); no free-memory preflight.

- lbm/block_stepper.py deleted; MAX_GRID_DIM lives in refinement_stepper.
- tests/validation/multigrid/test_block_stepper.py: same frozen GPU gates
  (split identity, FP64 oracle, circle across seams) and pre-CUDA rejection
  cases on the single level; the old-vs-new bitwise test is retired with the
  old class (its result is the 1318fc2 record).
- tests/unit/test_block_step_sync.py folded into test_single_level.py
  (step(n) vs n x step(1) roles and one wait, step(0), mid-batch launch
  failure), which also gains bad-state rejection before upload, download
  failure and activation-failure retry on close.
- test_shifted_storage and test_raw_exchange_gpu use the single level; the
  raw-exchange reference feeds the AMR Jet column and raw ledger through the
  root level boundary instead of patching a private kernel field.
- Pinned historical runners (taylor_green/run_case.py, run_jet_circle_cv.py)
  are unchanged: they only run on trees that still ship PeriodicBlockStepper.
  No alias.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-27 01:02:08 +00:00
co-authored by Claude Opus 5.5
parent d58f06f0fa
commit 8e1feffcad
12 changed files with 153 additions and 536 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/lifecycle.md
+3 -3
View File
@@ -1,4 +1,4 @@
# Device state lifecycle (uniform Simulation, PeriodicBlockStepper, FixedPatchStepper)
# Device state lifecycle (uniform Simulation, FixedPatchStepper)
State is grouped by what a change invalidates. This file describes the code as
it is; it adds no framework. Line references are to candidate 0f8fdbc3.
@@ -13,8 +13,8 @@ it is; it adds no framework. Line references are to candidate 0f8fdbc3.
Storage routes are separate concerns:
- Original uniform grid (`uniform.cu`, `Simulation`): whole-domain `NX`/`NY`.
- Blocked levels (`blocked.cu`; `PeriodicBlockStepper` and every blocked
`FixedPatchStepper` level): `NX`/`NY` are one block window including halo.
- Blocked levels (`blocked.cu`; every `FixedPatchStepper` level, including the
single level of `refinements=[]`): `NX`/`NY` are one block window including halo.
- Dense refinement levels (`refinement.cu`, `level_blocks[l] is None`).
A single-rectangle refinement (`FixedPatch2D`) is stored in its own fine
docs/supported_settings.md
+1 -1
View File
@@ -15,7 +15,7 @@ measured; they are not general bounds.
| Collision on refinement levels | SRT or MRT; TRT and LES rejected | `lbm/refinement.py` `FixedPatch2D.from_config` (re-checked for fine boundaries in `prepare_fine_boundary`) |
| Shifted (ddf_shifting) storage | only when every level is blocked; dense or mixed rejected | `lbm/refinement_stepper.py` `FixedPatchStepper.__init__`, `FixedPatchStepper._check_storage`; `#error` guard in `lbm/kernels/refinement.cu` (dense patch kernel unshifted only) |
| Fine relaxation | 1/tau_fine must lie inside the configured omega guards | `lbm/refinement.py` `FixedPatch2D.from_config` |
| Uniform blocked backend `PeriodicBlockStepper` | doubly periodic, SRT only, stationary passive circles | `lbm/block_stepper.py` `PeriodicBlockStepper.__init__`, `_static_circles` |
| Single blocked level `FixedPatchStepper(refinements=[], level_blocks=[(by, bx)], level_objects=[circles])` (replaces `PeriodicBlockStepper` after 888c8038) | doubly periodic, SRT only, no LES/Jets/projection/exchange observation, stationary passive circles; `download()` returns a 1-tuple | `lbm/refinement_stepper.py` `_check_levels`, `FixedPatchStepper.__init__`; `lbm/blocks.py` `static_circle_tables` |
Tests: `tests/unit/test_refinement_mrt.py`, `tests/validation/multigrid/test_refinement_mrt_gpu.py`,
`tests/unit/test_shifted_storage.py`, `tests/unit/test_shifted_block_operators.py`.
src/CelerisLab/cuda/compiler.py
+3 -3
View File
@@ -281,9 +281,9 @@ def build_module(cfg: LBMConfig, source: str, *, arch: str, n_objects: int = 0,
block_root_bc: int = 0):
"""Generate, compile and load one private module; the build directory is removed.
The single build path of every owner (Simulation, PeriodicBlockStepper,
FixedPatchStepper levels). Each call returns a NEW module: callers must not
share it, because its globals (d_params, ...) are mutable per module.
The single build path of every owner (Simulation, FixedPatchStepper levels).
Each call returns a NEW module: callers must not share it, because its
globals (d_params, ...) are mutable per module.
"""
build_dir = generate_config(cfg, n_objects=n_objects, block_root_bc=block_root_bc)
try:
src/CelerisLab/lbm/block_stepper.py
-230
View File
@@ -1,230 +0,0 @@
# CelerisLab/lbm/block_stepper.py
"""Opt-in periodic backend for one BlockLayout2D level, optional static circles.
D2Q9/SRT/FP32 unforced double buffers. Each step refreshes every halo from
its owner core, then pulls/collides core cells into the spare buffer, then
swaps. Only halo pair tables are uploaded; dense ownership stays on the host.
The caller owns the CudaContext and must close this backend before it.
Optional ``objects`` are stationary passive circles (static linear BFL via
receiver-local staging: halo refresh -> cut links -> pull, no copy between).
No force/actuation readout, motion, walls, refinement, checkpoint or runtime
parameter mutation; wall/slot/inlet fields in cfg are ignored.
"""
from dataclasses import replace
from numbers import Integral
from threading import RLock
import numpy as np
import pycuda.driver as cuda
from ..cuda import compiler
from .blocks import (NQ, BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config,
static_circle_tables)
from .population_storage import _StoredPopulations, _encode_physical, _decode_stored
MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks)
class PeriodicBlockStepper:
"""G->G stepping of Q-major (9, ny, nx) postcollision populations.
On a device/launch failure the instance becomes unusable: close and
recreate it rather than assuming a partially advanced buffer rolled back.
"""
def __init__(self, cfg, layout, *, context, objects=()):
if not isinstance(layout, BlockLayout2D) or layout.periodic != (True, True):
raise ValueError("PeriodicBlockStepper requires a doubly periodic BlockLayout2D")
checked = replace(cfg)
checked.validate()
if (not checked.is_d2q9 or checked.collision != "SRT"
or checked.streaming != "double_buffer" or checked.store_precision != "FP32"
or checked.les_enabled
or (checked.ny, checked.nx) != layout.shape):
raise ValueError("PeriodicBlockStepper requires matching D2Q9/SRT/FP32/"
"double_buffer config without LES")
ly, lx = layout.local_shape
pairs = len(layout.halo_dst)
tpb = checked.threads_per_block
if layout.nblocks > MAX_GRID_DIM or ly > MAX_GRID_DIM or not 0 < pairs < 2**31:
raise ValueError("block count, block rows or halo pairs exceed launch limits")
objects = tuple(objects) # an empty iterable takes the body-free path
self.boundary = static_circle_tables(layout, objects) if objects else None
links = 0 if self.boundary is None else len(self.boundary["cell"])
if links >= 2**31:
raise ValueError("cut links exceed launch limits")
self._ddf_shifting = bool(checked.ddf_shifting)
self.layout = layout
self.ctx = context
self._links = links
self._local_cfg = level_config(checked, layout, 3 * checked.viscosity + .5)
self._block = (tpb, 1, 1)
self._grid = ((lx + tpb - 1) // tpb, ly, layout.nblocks)
self._bytes = int(np.prod(layout.blocked_shape)) * 4
self._lock = RLock()
self._status = "new"
self._step_count = 0
self._allocations = []
self._stream = self._module = self._copy = self._step = self._curved = None
self._tables = {}
with self.ctx.activate():
try:
free, _ = cuda.mem_get_info()
tables = sum(a.nbytes for k, a in (self.boundary or {}).items() if k != "solid")
if 2 * self._bytes + layout.halo_dst.nbytes + layout.halo_src.nbytes + tables > free:
raise MemoryError("blocked buffers exceed free device memory")
self._stream = cuda.Stream()
config = self._local_cfg
module = compiler.build_module(config, "blocked.cu", arch=self.ctx.sm_arch)
compiler.set_params(module, config.omega)
self._module = module
self._current = self._allocate(self._bytes)
self._spare = self._allocate(self._bytes)
self._dst = self._upload(layout.halo_dst)
self._src = self._upload(layout.halo_src)
self._copy = module.get_function("BlockHaloCopy")
self._step = module.get_function("BlockStep")
if self.boundary is not None:
self._tables = {k: self._upload(a) for k, a in self.boundary.items() if k != "solid"}
self._curved = module.get_function("BlockCurvedBoundary")
except BaseException:
try:
self.close()
except Exception:
pass # Preserve the original construction failure.
raise
def _allocate(self, size):
pointer = cuda.mem_alloc(size)
self._allocations.append(pointer)
return pointer
def _upload(self, array):
array = np.ascontiguousarray(array)
pointer = self._allocate(array.nbytes)
cuda.memcpy_htod(pointer, array)
return pointer
def _require(self, *states):
if self._status not in states:
raise RuntimeError(f"Block stepper is {self._status}; expected {' or '.join(states)}")
@property
def step_count(self):
return self._step_count
def initialize(self, g):
"""Upload dense G_0; the spare buffer starts NaN so stale halos cannot hide."""
with self._lock, self.ctx.activate():
self._require("new", "ready")
dense = np.array(g, dtype=np.float64 if self._ddf_shifting else np.float32,
order="C", copy=True)
if dense.shape != (NQ, *self.layout.shape):
raise ValueError("Expected a Q-major (9, ny, nx) array")
fluid = np.ones(self.layout.shape, bool) if self.boundary is None else ~self.boundary["solid"]
if not np.isfinite(dense[:, fluid]).all() or (dense[:, fluid].sum(axis=0) <= 0).any():
raise ValueError("Fluid populations must be finite with positive density")
dense[:, ~fluid] = np.nan # Solid cores are deliberately invalid.
encoded = _encode_physical(dense, self._ddf_shifting)
values = encoded[:, fluid]
if (not np.isfinite(values).all()
or (values.sum(axis=0, dtype=np.float64) + int(self._ddf_shifting) <= 0).any()):
raise ValueError("Stored fluid populations must be finite with positive density")
blocked = self.layout.to_blocked(encoded)
spare = np.full_like(blocked, np.nan)
try:
cuda.memcpy_htod(self._current, blocked)
cuda.memcpy_htod(self._spare, spare)
self._step_count = 0
self._status = "ready"
except BaseException:
self._status = "failed"
raise
def step(self, n=1):
"""Advance n steps: halo refresh, then core pull/collide, per step.
All n steps are enqueued in order on the one stream and the host waits
once before returning (same launches and bits as n step(1) calls).
Asynchronous device errors surface at that single wait, so they are not
attributed to the step that caused them; call step(1) repeatedly to
localize one. On any failure the instance is failed and step_count keeps
its entry value.
"""
if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
raise ValueError("n must be a nonnegative integer")
with self._lock, self.ctx.activate():
self._require("ready")
if n == 0:
return
kw = dict(block=self._block, stream=self._stream)
try:
self._status = "advancing"
for _ in range(n):
launch_halo_copy(self._copy, self._current, self._dst, self._src,
len(self.layout.halo_dst), **kw)
if self._curved is not None: # stationary: no actuation or observation
launch_curved_boundary(self._curved, self._current, self._tables,
self._links, **kw)
self._step(self._current, self._spare, np.uint32(self.layout.halo),
self._tables.get("flags", np.uintp(0)), np.uintp(0), # no kind
grid=self._grid, **kw)
self._current, self._spare = self._spare, self._current # host handles only
self._stream.synchronize()
self._step_count += n
self._status = "ready"
except BaseException:
self._status = "failed"
raise
def download(self):
"""Return physical core populations (FP64 when shifted); no host cache."""
snapshot = self._download_raw()
return _decode_stored(snapshot.values[0], snapshot.shifted)
def _download_raw(self):
"""Return one tagged dense FP32 population array without conversion."""
with self._lock, self.ctx.activate():
self._require("ready")
blocked = np.empty(self.layout.blocked_shape, np.float32)
try:
cuda.memcpy_dtoh(blocked, self._current)
except BaseException:
self._status = "failed"
raise
return _StoredPopulations((self.layout.to_dense(blocked),), self._ddf_shifting)
def close(self):
"""Release owned resources; retain them for retry if context activation fails."""
with self._lock:
if self._status == "closed":
return
self._status = "failed"
with self.ctx.activate():
self._status = "closed"
error = None
try:
if self._stream is not None:
self._stream.synchronize()
except Exception as exc:
error = exc
self._copy = self._step = self._curved = self._module = None
self._tables = {}
for pointer in self._allocations:
try:
pointer.free()
except Exception as exc:
error = error or exc
self._allocations.clear()
self._stream = None
if error is not None:
raise error
def __enter__(self):
return self
def __exit__(self, *exc):
try:
self.close()
except Exception:
if exc[0] is None:
raise
src/CelerisLab/lbm/blocks.py
+1 -2
View File
@@ -375,8 +375,7 @@ def blocked_link_tables(layout, flags, fluid_idx, dir, q, fallback):
return tables, mapped
# One blocked level as both owners (PeriodicBlockStepper, FixedPatchStepper) build and
# launch it. Raw argument orders of the blocked.cu kernels live only here.
# One blocked level as FixedPatchStepper builds and launches it (any depth). Raw argument orders of the blocked.cu kernels live only here.
def level_config(cfg, layout, tau):
"""Kernel config of one blocked level: NX/NY = block window incl. halo.
src/CelerisLab/lbm/refinement_stepper.py
+3 -2
View File
@@ -19,7 +19,6 @@ from .blocks import (BlockLayout2D, launch_curved_boundary, launch_halo_copy, le
static_circle_tables)
from .channel import channel_masks, validate_channel_config, validate_channel_transfers
from ..body.jets import jet_link_velocity
from .block_stepper import MAX_GRID_DIM
from .hierarchy import build_hierarchy, is_rectangle
from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary
from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noqa: F401 (re-export)
@@ -67,6 +66,8 @@ def _level_name(level):
return ("base", "fine", "child")[level] if level < 3 else f"level {level}"
MAX_GRID_DIM = 65535 # CUDA gridDim.y/z limit (rows, blocks)
def _check_levels(refinements, blocks, objects, jets):
"""Candidate limits of the level lists, checked before any device work."""
depth = len(refinements) + 1
@@ -104,7 +105,7 @@ class FixedPatchStepper:
hold one entry per level, L0 included (L0 bodies/Jets are not supported).
``refinements=[]`` is one blocked level: periodic SRT root without LES,
Jets or exchange observation; its ``level_objects[0]`` are stationary
passive circles (static BFL, as PeriodicBlockStepper).
passive circles (static BFL).
``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
a bottom on-site slot and top fixed-density reservoir from ``cfg``.
tests/unit/test_block_step_sync.py
-104
View File
@@ -1,104 +0,0 @@
"""Host contract of PeriodicBlockStepper.step(n): one stream wait per call (mocked driver)."""
import ast
from contextlib import nullcontext
import importlib.util
from pathlib import Path
import types
from unittest.mock import Mock
import unittest
import numpy as np
from CelerisLab.config import LBMConfig
from CelerisLab.lbm.blocks import BlockLayout2D
def host_module():
"""block_stepper with pycuda/compiler replaced by mocks; no GPU context."""
# Locate the source without importing it: importing block_stepper loads pycuda.driver,
# which CPU-only import-boundary tests (test_build_contracts) must never see.
path = Path(importlib.util.find_spec("CelerisLab.lbm.block_stepper").origin)
tree = ast.parse(path.read_text())
tree.body = [node for node in tree.body if not (
isinstance(node, ast.Import) and any(a.name == "pycuda.driver" for a in node.names)
or isinstance(node, ast.ImportFrom) and node.module == "cuda")]
module = types.ModuleType("CelerisLab.lbm._test_host_block_stepper")
module.__package__ = "CelerisLab.lbm"
module.__dict__.update(cuda=Mock(), compiler=Mock())
exec(compile(tree, str(path), "exec"), module.__dict__)
return module
def backend():
host, log = host_module(), []
host.cuda.mem_get_info.return_value = (1 << 40, 1 << 40)
host.cuda.mem_alloc.side_effect = lambda size: Mock(size=size) # distinct buffers
module = Mock()
module.get_global.return_value = (Mock(), 20)
module.get_function.side_effect = lambda name: Mock(
side_effect=lambda *a, **k: log.append((name, a[:2])))
host.compiler.build_module.return_value = module
cfg = LBMConfig(nx=16, ny=16, viscosity=.1)
run = host.PeriodicBlockStepper(cfg, BlockLayout2D.build(16, 16, 8, 8),
context=Mock(activate=nullcontext, sm_arch="sm_70"))
run._stream.synchronize.side_effect = lambda: log.append(("sync", ()))
run.initialize(np.full((9, 16, 16), 1 / 9))
log.clear()
return run, log
class BlockStepSync(unittest.TestCase):
def test_step_n_waits_once_and_matches_n_single_steps(self):
many, log_many = backend()
many.step(3)
single, log_single = backend()
for _ in range(3):
single.step(1)
names = [name for name, _ in log_many]
self.assertEqual(names.count("sync"), 1)
self.assertEqual(names[-1], "sync")
self.assertEqual(names.count("BlockStep"), 3)
# Same launch order and buffer roles as three step(1) calls.
strip = lambda log: [(n, tuple(id(x) for x in a)) for n, a in log if n != "sync"]
roles = lambda run, log: [(n, tuple({id(run._current): "c", id(run._spare): "s"}.get(x, "-")
for x in ids)) for n, ids in strip(log)]
self.assertEqual([n for n, _ in strip(log_many)], [n for n, _ in strip(log_single)])
self.assertEqual(roles(many, log_many), roles(single, log_single))
self.assertIsNot(many._current, many._spare)
self.assertEqual((many.step_count, single.step_count), (3, 3))
def test_step_zero_launches_and_waits_for_nothing(self):
run, log = backend()
run.step(0)
self.assertEqual((log, run.step_count), ([], 0))
def test_failure_keeps_entry_step_count(self):
run, log = backend()
run.step(1)
run._stream.synchronize.side_effect = RuntimeError("device fault")
with self.assertRaises(RuntimeError):
run.step(2)
self.assertEqual((run.step_count, run._status), (1, "failed"))
def test_launch_failure_mid_batch_fails_closed_without_waiting(self):
run, log = backend()
calls = []
step = run._step
def failing(*args, **kwargs):
calls.append(1)
if len(calls) == 2:
raise RuntimeError("launch failure")
return step(*args, **kwargs)
run._step = failing
with self.assertRaises(RuntimeError):
run.step(3)
self.assertNotIn("sync", [name for name, _ in log])
self.assertEqual((run.step_count, run._status), (0, "failed"))
with self.assertRaises(RuntimeError):
run.download()
if __name__ == "__main__":
unittest.main()
tests/unit/test_refinement_three_level_multipatch.py
+2 -2
View File
@@ -22,11 +22,11 @@ def host_module():
tree.body=[node for node in tree.body if not (
isinstance(node,ast.Import) and any(a.name=='pycuda.driver' for a in node.names)
or isinstance(node,ast.ImportFrom) and node.module in (
'cuda','block_stepper','curved_links'))]
'cuda','curved_links'))]
module=types.ModuleType('CelerisLab.lbm._test_host_refinement')
module.__package__='CelerisLab.lbm'
from CelerisLab.body.jets import surface_link_velocity
module.__dict__.update(cuda=Mock(),compiler=Mock(),MAX_GRID_DIM=65535,
module.__dict__.update(cuda=Mock(),compiler=Mock(),
surface_link_velocity=surface_link_velocity)
exec(compile(tree,str(path),'exec'),module.__dict__)
return module
tests/unit/test_shifted_storage.py
+5 -17
View File
@@ -1,7 +1,4 @@
"""CPU host contracts; mocked CUDA executes real lifecycle, not numerical kernels."""
import ast
import types
from pathlib import Path
from contextlib import nullcontext
from dataclasses import replace
import unittest
@@ -138,24 +135,14 @@ class ShiftedStorage(unittest.TestCase):
module.compiler.build_module.assert_not_called()
def test_uniform_block_codec_and_raw_snapshot(self):
path = Path(backend.r.__file__).with_name('block_stepper.py')
tree = ast.parse(path.read_text())
tree.body = [n for n in tree.body if not (
isinstance(n, ast.Import) and any(a.name == 'pycuda.driver' for a in n.names)
or isinstance(n, ast.ImportFrom) and n.module == 'cuda')]
module = types.ModuleType('CelerisLab.lbm._test_block_host')
module.__package__ = 'CelerisLab.lbm'
cuda, compiler = Mock(), Mock()
module.__dict__.update(cuda=cuda,compiler=compiler)
exec(compile(tree,str(path),'exec'),module.__dict__)
cuda.mem_get_info.return_value=(10**9,10**9)
module = backend.host_module() # single blocked level: FixedPatchStepper(refinements=[])
cuda, compiler = module.cuda, module.compiler
cuda.mem_alloc.side_effect=lambda size: Mock()
gpu=Mock();gpu.get_global.return_value=(Mock(),20)
compiler.build_module.return_value=gpu
context=Mock(activate=nullcontext,sm_arch='sm_70')
layout=module.BlockLayout2D.build(16,16,4,4)
cfg=replace(backend.CFG,nx=16,ny=16,ddf_shifting=True)
run=module.PeriodicBlockStepper(cfg,layout,context=context)
run=module.FixedPatchStepper(cfg,context=context,refinements=[],level_blocks=[(4,4)])
storage={}
cuda.memcpy_htod.side_effect=lambda p,a: storage.__setitem__(id(p),np.array(a,copy=True))
cuda.memcpy_dtoh.side_effect=lambda a,p: np.copyto(a,storage[id(p)])
@@ -166,7 +153,8 @@ class ShiftedStorage(unittest.TestCase):
self.assertTrue(raw.shifted)
self.assertEqual(len(raw.values),1)
self.assertTrue((raw.values[0][1]==np.float32(2**-35)).all())
np.testing.assert_array_equal(run.download(),physical)
(down,)=run.download()
np.testing.assert_array_equal(down,physical)
finally: run.close()
def test_passive_shifted_body_raw_remap_preserves_support(self):
tests/unit/test_single_level.py
+84 -2
View File
@@ -2,10 +2,11 @@
Per step: halo refresh -> static BFL (circles only) -> BlockStep with a null
kind, one stream wait per step(n) call, NaN spare after initialize. An upload,
launch or wait failure fails the instance and keeps its step count.
launch, wait or download failure fails the instance and keeps its step count;
close retains resources for retry when context activation fails.
"""
from contextlib import nullcontext
from unittest.mock import Mock
from unittest.mock import MagicMock, Mock
import unittest
import numpy as np
@@ -108,6 +109,87 @@ class SingleLevel(unittest.TestCase):
with self.assertRaisesRegex(RuntimeError, "failed"):
run.step(1)
def test_step_n_waits_once_and_matches_n_single_steps(self):
many, log_many, _ = backend()
many.step(3)
single, log_single, _ = backend()
for _ in range(3):
single.step(1)
names = [name for name, _ in log_many]
self.assertEqual(names.count("sync"), 1)
self.assertEqual(names[-1], "sync")
self.assertEqual(names.count("BlockStep"), 3)
def roles(run, log): # same launch order and buffer roles as three step(1) calls
level = run._levels[0]
tags = {id(level.current): "c", id(level.spare): "s"}
return [(n, tuple(tags.get(id(x), "-") for x in a)) for n, a in log if n != "sync"]
self.assertEqual(roles(many, log_many), roles(single, log_single))
self.assertIsNot(many._levels[0].current, many._levels[0].spare)
self.assertEqual((many.step_count, single.step_count), (3, 3))
def test_step_zero_launches_and_waits_for_nothing(self):
run, log, _ = backend()
log.clear()
run.step(0)
self.assertEqual((log, run.step_count), ([], 0))
def test_launch_failure_mid_batch_fails_without_waiting(self):
run, log, host = backend()
step = run._levels[0].kernel
calls = []
def failing(*args, **kwargs):
calls.append(1)
if len(calls) == 2:
raise RuntimeError("launch failure")
return step(*args, **kwargs)
run._levels[0].kernel = failing
log.clear()
with self.assertRaises(RuntimeError):
run.step(3)
self.assertNotIn("sync", [name for name, _ in log])
self.assertEqual((run.step_count, run._status), (0, "failed"))
with self.assertRaises(RuntimeError):
run.download()
def test_initialize_rejects_bad_state_before_upload(self):
run, log, _ = backend(initialize=False)
good = state()
for bad in (good[:, :-1], np.where(np.arange(32) == 3, np.nan, good), -good):
with self.assertRaises(ValueError):
run.initialize(bad)
self.assertNotIn("upload", [name for name, _ in log])
run.initialize(good)
self.assertEqual([name for name, _ in log].count("upload"), 2)
def test_download_failure_fails_instance(self):
run, _, host = backend()
host.cuda.memcpy_dtoh.side_effect = RuntimeError("read failed")
with self.assertRaisesRegex(RuntimeError, "read failed"):
run.download()
for operation in (run.download, run.step):
with self.assertRaisesRegex(RuntimeError, "failed"):
operation()
def test_activation_failure_retains_resources_for_retry(self):
run, _, _ = backend()
allocations, stream = run._allocations[:], run._stream
unavailable = MagicMock()
unavailable.__enter__.side_effect = RuntimeError("activation failed")
run.ctx.activate = Mock(side_effect=[unavailable, nullcontext()])
with self.assertRaisesRegex(RuntimeError, "activation failed"):
run.close()
self.assertEqual(run._status, "failed")
self.assertEqual(run._allocations, allocations)
self.assertIs(run._stream, stream)
run.close()
self.assertEqual(run._status, "closed")
self.assertEqual(run._levels, [])
for allocation in allocations:
allocation.free.assert_called_once()
if __name__ == "__main__":
unittest.main()
tests/validation/multigrid/test_block_stepper.py
+40 -159
View File
@@ -1,4 +1,4 @@
"""Periodic body-free block backend: CPU contract doubles plus one short GPU case.
"""Single blocked level, FixedPatchStepper(refinements=[]): CPU contracts plus short GPU cases.
Frozen GPU gates (32x24, 8 steps, positive low-Mach non-equilibrium state):
finite cores; bit-identical 1x1, 2x2 and 4x2 splits; max |G - FP64 oracle|
@@ -11,9 +11,8 @@ Mass/momentum drift is reported, not gated.
from contextlib import nullcontext
from pathlib import Path
import resource
from threading import RLock
import unittest
from unittest.mock import MagicMock, Mock, patch
from unittest.mock import Mock, patch
import numpy as np
import pycuda.driver as cuda
@@ -22,9 +21,8 @@ from CelerisLab.body.geometry.circle import CircleGeometry
from CelerisLab.body.objects import SimObject
from CelerisLab.config import LBMConfig
from CelerisLab.cuda import compiler
from CelerisLab.lbm import block_stepper as implementation
from CelerisLab.lbm.blocks import NO_DONOR, BlockLayout2D
from CelerisLab.lbm.block_stepper import PeriodicBlockStepper
from CelerisLab.lbm import refinement_stepper as implementation
from CelerisLab.lbm.blocks import NO_DONOR, BlockLayout2D, static_circle_tables
from CelerisLab.lbm.refinement_stepper import FixedPatchStepper
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
@@ -113,92 +111,72 @@ def totals(g):
return g.sum(), (CX[:, None, None] * g).sum(), (CY[:, None, None] * g).sum()
def single(cfg=None, *, context, block=(12, 16), objects=(), **kw):
"""One periodic blocked level (the former PeriodicBlockStepper capability)."""
return FixedPatchStepper(cfg or config(), context=context, refinements=[], level_blocks=[block],
level_objects=[tuple(objects)], **kw)
class CpuContracts(unittest.TestCase):
"""No device: driver doubles prove rejection before allocation and lifecycle."""
"""No device: driver doubles prove rejection before allocation."""
def forbidden_context(self):
return Mock(activate=Mock(side_effect=AssertionError("CUDA reached")), sm_arch="sm_70")
def test_single_level_fixed_patch_contract_before_cuda(self):
def single(cfg=None, **kw):
options = dict(context=self.forbidden_context(), refinements=[], level_blocks=[(12, 16)])
options.update(kw)
return FixedPatchStepper(cfg or config(), **options)
for cfg, kw in ((None, dict(level_blocks=[None])),
def test_unsupported_inputs_rejected_before_cuda(self):
for cfg, kw in ((None, dict(block=None)),
(None, dict(coarse_boundary="channel")),
(None, dict(coarse_boundary="slot_reservoir")),
(config(collision="MRT"), {}),
(config(les_enabled=True), {}),
(None, dict(explosion="projected_linear")),
(None, dict(observe_raw_exchange=True)),
(None, dict(observe_raw_exchange=True, observe_window_exchange=True)),
(config(les_enabled=True), {}),
(None, dict(level_jets=[(object(),)])),
(None, dict(explosion="projected_linear")),
(None, dict(level_objects=[(circle((15.5, 4.0), 4.3),)]))): # touches edge
(None, dict(level_jets=[(object(),)]))):
with self.subTest(kw=kw), self.assertRaises(ValueError):
single(cfg, **kw)
with self.assertRaisesRegex(AssertionError, "CUDA reached"):
single(level_objects=[(circle(),)])
def test_unsupported_inputs_rejected_before_cuda(self):
layout = BlockLayout2D.build(NY, NX, 12, 16)
cases = ((config(), BlockLayout2D.build(NY, NX, 12, 16, periodic=(True, False))),
(config(), BlockLayout2D.build(NY, NX * 2, 12, 16)),
(config(collision="MRT"), layout),
(config(), "layout"))
for cfg, lay in cases:
with self.subTest(cfg=cfg.collision, layout=getattr(lay, "periodic", lay)):
with self.assertRaises(ValueError):
PeriodicBlockStepper(cfg, lay, context=self.forbidden_context())
single(cfg, context=self.forbidden_context(), **kw)
moving = circle()
moving.state.vx = .01
build_links = CircleGeometry.build_cut_links
for objects in ([moving], [circle((1.5, 11.5), 2.0)], [circle(), circle()]):
with self.assertRaises(ValueError):
PeriodicBlockStepper(config(), layout, context=self.forbidden_context(), objects=objects)
single(context=self.forbidden_context(), objects=objects)
for center, radius in (((16.0, 12.0), -.1), ((16.0, 12.0), -.5), ((15.5, 11.5), np.nan),
((np.inf, 11.5), 4.3), ((15.5, np.nan), 4.3)):
with self.subTest(center=center, radius=radius), self.assertRaises(ValueError):
PeriodicBlockStepper(config(), layout, context=self.forbidden_context(),
objects=[circle(center, radius)])
single(context=self.forbidden_context(), objects=[circle(center, radius)])
with patch.object(implementation, "static_circle_tables") as prepare:
with self.assertRaisesRegex(AssertionError, "CUDA reached"): # body-free path
PeriodicBlockStepper(config(), layout, context=self.forbidden_context(), objects=iter(()))
single(context=self.forbidden_context(), objects=iter(()))
prepare.assert_not_called()
with self.assertRaisesRegex(AssertionError, "CUDA reached"):
single(context=self.forbidden_context(), objects=[circle()])
with patch.object(CircleGeometry, "build_cut_links",
lambda geometry, *a, **k: build_links(geometry, *a, **k)[1:]):
with self.assertRaisesRegex(ValueError, "cover every fluid pull"):
PeriodicBlockStepper(config(), layout, context=self.forbidden_context(),
objects=[circle()])
single(context=self.forbidden_context(), objects=[circle()])
with patch.object(implementation, "MAX_GRID_DIM", 3):
with self.assertRaisesRegex(ValueError, "launch limits"):
PeriodicBlockStepper(config(), layout, context=self.forbidden_context())
single(context=self.forbidden_context())
def test_memory_and_param_layout_rejected_before_allocation(self):
layout = BlockLayout2D.build(NY, NX, 12, 16)
def test_param_layout_rejected_before_allocation(self):
context = Mock(activate=nullcontext, sm_arch="sm_70")
with patch.object(cuda, "mem_get_info", return_value=(1024, 1 << 30)), \
patch.object(cuda, "Stream") as stream, patch.object(cuda, "mem_alloc") as alloc:
with self.assertRaises(MemoryError):
PeriodicBlockStepper(config(), layout, context=context)
stream.assert_not_called()
alloc.assert_not_called()
module = Mock()
module.get_global.return_value = (Mock(), 16)
with patch.object(cuda, "mem_get_info", return_value=(1 << 30, 1 << 30)), \
patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \
with patch.object(cuda, "Stream"), patch.object(compiler, "generate_config"), \
patch.object(compiler, "compile_kernel") as compile_kernel, \
patch.object(compiler, "cleanup_build_dir"), \
patch.object(compiler, "load_module", return_value=module), \
patch.object(cuda, "mem_alloc") as alloc, patch.object(cuda, "memcpy_htod") as upload:
with self.assertRaisesRegex(RuntimeError, "layout mismatch"):
PeriodicBlockStepper(config(), layout, context=context)
single(context=context)
self.assertEqual(compile_kernel.call_args.kwargs["source"], "blocked.cu")
alloc.assert_not_called()
upload.assert_not_called()
def test_static_circle_tables(self):
layout = BlockLayout2D.build(NY, NX, 12, 16)
tables = implementation.static_circle_tables(layout, (circle(),))
tables = static_circle_tables(layout, (circle(),))
links = CircleGeometry(*CENTER, RADIUS).build_cut_links(NX, NY, domain_flags=dense_flags())
self.assertEqual(len(tables["cell"]), len(links))
np.testing.assert_array_equal(tables["q"], np.array([l.q for l in links], np.float32))
@@ -206,85 +184,11 @@ class CpuContracts(unittest.TestCase):
self.assertTrue((tables["flags"][solid_cells] & 2).all())
self.assertEqual(int(np.count_nonzero(tables["flags"] & 2)), int(tables["solid"].sum()))
def fake_backend(self):
backend = PeriodicBlockStepper.__new__(PeriodicBlockStepper)
backend.layout = BlockLayout2D.build(NY, NX, 12, 16)
backend.ctx = Mock(activate=nullcontext)
backend._lock = RLock()
backend._status = "ready"
backend._step_count = 0
backend._ddf_shifting = False # unshifted storage, as for the default config
backend._stream = Mock()
backend._module = Mock()
backend._copy = backend._step = Mock()
backend.boundary, backend._curved, backend._tables, backend._links = None, None, {}, 0
backend._dst, backend._src = Mock(), Mock()
backend._block, backend._grid = (32, 1, 1), (1, 1, 1)
backend._allocations = [Mock(), Mock()]
backend._current, backend._spare = backend._allocations
return backend
def test_initialize_rejects_bad_state_before_upload(self):
backend = self.fake_backend()
good = initial_state()
bad = (good[:, :-1], np.where(np.arange(NX) == 3, np.nan, good), -good)
with patch.object(cuda, "memcpy_htod") as upload:
for g in bad:
with self.assertRaises(ValueError):
backend.initialize(g)
upload.assert_not_called()
backend.initialize(good)
self.assertEqual(upload.call_count, 2)
self.assertTrue(np.isnan(upload.call_args_list[1].args[1]).all())
backend.close()
def test_failed_launch_and_download_make_backend_unusable(self):
for failure in ("step", "download"):
backend = self.fake_backend()
if failure == "step":
backend._step = Mock(side_effect=RuntimeError("launch failed"))
with self.assertRaisesRegex(RuntimeError, "launch failed"):
backend.step()
else:
with patch.object(cuda, "memcpy_dtoh", side_effect=RuntimeError("read failed")):
with self.assertRaisesRegex(RuntimeError, "read failed"):
backend.download()
for operation in (backend.download, backend.step):
with self.assertRaisesRegex(RuntimeError, "failed"):
operation()
backend.close()
for allocation in backend_allocations(backend):
allocation.free.assert_called_once()
def test_activation_failure_retains_resources_for_retry(self):
backend = self.fake_backend()
allocations, stream, module = backend._allocations[:], backend._stream, backend._module
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._allocations, allocations)
self.assertIs(backend._stream, stream)
self.assertIs(backend._module, module)
backend.close()
self.assertEqual(backend._status, "closed")
self.assertIsNone(backend._module)
self.assertIsNone(backend._copy)
self.assertIsNone(backend._step)
for allocation in allocations:
allocation.free.assert_called_once()
def backend_allocations(backend):
return (backend._current, backend._spare)
class PeriodicBlockGPU(unittest.TestCase):
class SingleLevelGPU(unittest.TestCase):
def setUp(self):
self.assertEqual(Path(implementation.__file__).resolve(),
ROOT / "src/CelerisLab/lbm/block_stepper.py")
ROOT / "src/CelerisLab/lbm/refinement_stepper.py")
def tearDown(self):
rss = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / 1024
@@ -295,13 +199,13 @@ class PeriodicBlockGPU(unittest.TestCase):
g0 = initial_state()
results = []
with CudaContext(0) as context:
for by, bx in SPLITS:
layout = BlockLayout2D.build(NY, NX, by, bx)
with PeriodicBlockStepper(config(), layout, context=context) as backend:
for split in SPLITS:
with single(context=context, block=split) as backend:
backend.initialize(g0)
backend.step(STEPS)
self.assertEqual(backend.step_count, STEPS)
results.append(backend.download())
(g,) = backend.download()
results.append(g)
for (by, bx), g in zip(SPLITS, results):
self.assertTrue(np.isfinite(g).all(), (by, bx))
np.testing.assert_array_equal(g, results[0], err_msg=f"split {(by, bx)}")
@@ -328,16 +232,17 @@ class PeriodicBlockGPU(unittest.TestCase):
g0 = initial_state()
results = []
with CudaContext(0) as context:
for by, bx in SPLITS:
layout = BlockLayout2D.build(NY, NX, by, bx)
with PeriodicBlockStepper(config(), layout, context=context, objects=[circle()]) as backend:
for split in SPLITS:
with single(context=context, block=split, objects=[circle()]) as backend:
tables = backend._level_boundaries[0].tables
for key, values, dtype in (("q", [l.q for l in links], np.float32),
("dir", [l.dir for l in links], np.uint8),
("fallback", [l.fallback for l in links], np.uint8)):
np.testing.assert_array_equal(backend.boundary[key], np.array(values, dtype))
np.testing.assert_array_equal(tables[key], np.array(values, dtype))
backend.initialize(g0)
backend.step(STEPS)
results.append(backend.download())
(g,) = backend.download()
results.append(g)
for (by, bx), g in zip(SPLITS, results):
self.assertTrue(np.isfinite(g[:, ~solid]).all(), (by, bx))
self.assertTrue(np.isnan(g[:, solid]).all(), (by, bx))
@@ -352,30 +257,6 @@ class PeriodicBlockGPU(unittest.TestCase):
*(a - b for a, b in zip(after, before))), flush=True)
self.assertLessEqual(error, ORACLE_ATOL)
def test_single_level_fixed_patch_is_bit_identical(self):
"""FixedPatchStepper(refinements=[]) reproduces PeriodicBlockStepper bit for bit."""
from CelerisLab.cuda.context import CudaContext
g0 = initial_state()
with CudaContext(0) as context:
for (by, bx), objects, shifted in [(split, objects, shifted) for split in SPLITS
for objects in ((), (circle(),))
for shifted in (False, True)]:
with self.subTest(split=(by, bx), bodies=len(objects), shifted=shifted):
layout = BlockLayout2D.build(NY, NX, by, bx)
with PeriodicBlockStepper(config(ddf_shifting=shifted), layout, context=context,
objects=objects) as backend:
backend.initialize(g0)
backend.step(STEPS)
expected = backend.download()
with FixedPatchStepper(config(ddf_shifting=shifted), context=context, refinements=[],
level_blocks=[(by, bx)],
level_objects=[objects]) as backend:
backend.initialize(g0)
backend.step(STEPS)
self.assertEqual(backend.step_count, STEPS)
(actual,) = backend.download()
np.testing.assert_array_equal(actual, expected)
if __name__ == "__main__":
unittest.main()
tests/validation/multigrid/test_raw_exchange_gpu.py
+11 -11
View File
@@ -224,14 +224,12 @@ class RawExchangeGPU(unittest.TestCase):
Eight local hops are a short-time equivalence check, not convergence.
The test adapter only enables the standalone kernel's optional inputs.
"""
from CelerisLab.lbm.blocks import BlockLayout2D
from CelerisLab.lbm.block_stepper import PeriodicBlockStepper
fine_cfg = replace(CFG, nx=64, ny=48, viscosity=.2)
obj = SimObject(0, geometry=CircleGeometry(31.5, 23.5, 4.3),
center=(31.5, 23.5), radius=4.3)
layout = BlockLayout2D.build(48, 64, 8, 8)
with CudaContext(0) as ctx:
with PeriodicBlockStepper(fine_cfg, layout, context=ctx, objects=[obj]) as uniform, \
with FixedPatchStepper(fine_cfg, context=ctx, refinements=[], level_blocks=[(8, 8)],
level_objects=[(obj,)]) as uniform, \
FixedPatchStepper(CFG, context=ctx, **fine((16, 20)), observe_raw_exchange=True) as amr:
amr.initialize(*(state(l.shape) for l in amr._levels))
uniform.initialize(state((9, 48, 64)))
@@ -240,18 +238,20 @@ class RawExchangeGPU(unittest.TestCase):
links = obj.geometry.build_cut_links(64, 48)
np.testing.assert_array_equal([v.dir for v in links], b.tables['dir'])
np.testing.assert_allclose([v.q for v in links], b.tables['q'], atol=2e-7, rtol=0)
# Same fine-grid force arms after translation to the same center.
rx = uniform._upload(b.tables['rx'])
ry = uniform._upload(b.tables['ry'])
# Same fine-grid force arms after translation to the same center: the
# single level's static BFL gets the AMR Jet column and raw ledger inputs.
root = uniform._level_boundaries[0]
velocity = uniform._upload(b.tables['surface_velocity'])
raw = uniform._allocate(count*16)
original = uniform._curved
root.device.update(surface_velocity=velocity, rx=uniform._upload(b.tables['rx']),
ry=uniform._upload(b.tables['ry']))
root.last_raw = raw = uniform._allocate(count*16)
original = root.kernel
samples = []
def measured_boundary(*args, **kw):
original(*args[:7], velocity, rx, ry, raw, **kw)
original(*args, **kw)
uniform._stream.synchronize()
samples.append(download(raw, (count, 4)))
uniform._curved = measured_boundary
root.kernel = measured_boundary
for speed in (0., .002, -.001, 0.):
amr.set_jet_amplitude(0, speed)
cuda.memcpy_htod(velocity, b.tables['surface_velocity'])