refinement: FixedPatchStepper(refinements=[]) owns a single blocked level (S1)

A zero-refinement FixedPatchStepper is the PeriodicBlockStepper capability:
one doubly periodic blocked SRT level (no LES, Jets or exchange
observation) whose level_objects[0] are stationary passive circles. The
circle tables come from the same checks, now blocks.static_circle_tables
(moved from block_stepper, pycuda-free); BFL, halo and BlockStep launches
are the shared helpers, with an all-ACTIVE (null) kind. Base-level bodies
stay rejected for depth >= 2. PeriodicBlockStepper and its callers are
unchanged (deleted only after S1's GPU result).

Tests: host contract rejections before CUDA for the single level; GPU
bitwise comparison against PeriodicBlockStepper for every split, with and
without the circle (the S2 gate).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/CelerisLab/lbm/block_stepper.py

@@ -18,61 +18,13 @@ import numpy as np
 import pycuda.driver as cuda
 
 from ..cuda import compiler
-from .blocks import (NQ, BlockLayout2D, blocked_link_tables, launch_curved_boundary,
+from .blocks import (NQ, BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config,
-                     launch_halo_copy, level_config)
+                     static_circle_tables)
-from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID
 from .population_storage import _StoredPopulations, _encode_physical, _decode_stored
 
 MAX_GRID_DIM = 65535  # CUDA gridDim.y/z limit (rows, blocks)
 
 
-def _static_circles(layout, objects):
-    """Host tables for stationary circles; rejects before any CUDA work.
-
-    Every periodic neighbor of every fluid core is inspected, including domain
-    edges that CircleGeometry does not enumerate, so clipped or incomplete
-    geometry cannot pass. Each physical link has one receiver-owning block.
-    """
-    from ..body.geometry.circle import CircleGeometry
-
-    ny, nx = layout.shape
-    flags = np.full(ny * nx, FLUID, np.uint16)
-    for obj in objects:
-        if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region
-                or obj.surface_actuation is not None
-                or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega"))
-                or any(vars(obj.control).values())):
-            raise ValueError("Block bodies must be stationary passive circles")
-        cx, cy = obj.geometry.center
-        radius = obj.geometry._radius
-        if not np.isfinite([cx, cy, radius]).all() or radius <= 0:
-            raise ValueError("Block circle center/radius must be finite with radius > 0")
-        mask = obj.geometry.build_flag_mask(nx, ny)
-        solid = (mask & SOLID) != 0
-        if not solid.any() or np.any(flags[solid] & SOLID):
-            raise ValueError("Block circles must be resolved and disjoint")
-        flags[solid] = mask[solid]
-    solid = ((flags & SOLID) != 0).reshape(ny, nx)
-    if solid[[0, -1]].any() or solid[:, [0, -1]].any():
-        raise ValueError("Block circles must not touch the domain edge")
-    links = [link for obj in objects
-             for link in obj.geometry.build_cut_links(nx, ny, domain_flags=flags)]
-    def source(r, q):
-        y, x = divmod(r, nx)
-        return (y + D2Q9_EY[q]) % ny, (x + D2Q9_EX[q]) % nx
-    receivers = {(int(link.fluid_idx), int(link.dir)) for link in links}
-    required = {(int(r), q) for r in np.flatnonzero(~solid.ravel())
-                for q in range(1, NQ) if solid[source(r, q)]}
-    if len(receivers) != len(links) or receivers != required:
-        raise ValueError("Cut links must cover every fluid pull from a solid exactly once")
-    tables, _ = blocked_link_tables(
-        layout, flags, np.array([link.fluid_idx for link in links], np.uint32),
-        np.array([link.dir for link in links], np.uint8),
-        np.array([link.q for link in links], np.float32),
-        np.array([link.fallback for link in links], np.uint8))
-    return dict(solid=solid, **tables)
-
-
 class PeriodicBlockStepper:
     """G->G stepping of Q-major (9, ny, nx) postcollision populations.
 
@@ -96,7 +48,7 @@ class PeriodicBlockStepper:
         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_circles(layout, objects) if objects else None
+        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")

src/CelerisLab/lbm/blocks.py

@@ -407,3 +407,51 @@ def launch_curved_boundary(kernel, buffer, tables, count, *, block, stream,
            np.uint32(count), none if velocity is None else velocity,
            none if rx is None else rx, none if ry is None else ry, none if raw is None else raw,
            grid=((count + block[0] - 1) // block[0], 1), block=block, stream=stream)
+
+
+def static_circle_tables(layout, objects):
+    """Host tables for stationary circles; rejects before any CUDA work.
+
+    Every periodic neighbor of every fluid core is inspected, including domain
+    edges that CircleGeometry does not enumerate, so clipped or incomplete
+    geometry cannot pass. Each physical link has one receiver-owning block.
+    """
+    from ..body.geometry.circle import CircleGeometry
+    from .descriptors import FLUID, SOLID
+
+    ny, nx = layout.shape
+    flags = np.full(ny * nx, FLUID, np.uint16)
+    for obj in objects:
+        if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region
+                or obj.surface_actuation is not None
+                or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega"))
+                or any(vars(obj.control).values())):
+            raise ValueError("Block bodies must be stationary passive circles")
+        cx, cy = obj.geometry.center
+        radius = obj.geometry._radius
+        if not np.isfinite([cx, cy, radius]).all() or radius <= 0:
+            raise ValueError("Block circle center/radius must be finite with radius > 0")
+        mask = obj.geometry.build_flag_mask(nx, ny)
+        solid = (mask & SOLID) != 0
+        if not solid.any() or np.any(flags[solid] & SOLID):
+            raise ValueError("Block circles must be resolved and disjoint")
+        flags[solid] = mask[solid]
+    solid = ((flags & SOLID) != 0).reshape(ny, nx)
+    if solid[[0, -1]].any() or solid[:, [0, -1]].any():
+        raise ValueError("Block circles must not touch the domain edge")
+    links = [link for obj in objects
+             for link in obj.geometry.build_cut_links(nx, ny, domain_flags=flags)]
+    def source(r, q):
+        y, x = divmod(r, nx)
+        return (y + D2Q9_EY[q]) % ny, (x + D2Q9_EX[q]) % nx
+    receivers = {(int(link.fluid_idx), int(link.dir)) for link in links}
+    required = {(int(r), q) for r in np.flatnonzero(~solid.ravel())
+                for q in range(1, NQ) if solid[source(r, q)]}
+    if len(receivers) != len(links) or receivers != required:
+        raise ValueError("Cut links must cover every fluid pull from a solid exactly once")
+    tables, _ = blocked_link_tables(
+        layout, flags, np.array([link.fluid_idx for link in links], np.uint32),
+        np.array([link.dir for link in links], np.uint8),
+        np.array([link.q for link in links], np.float32),
+        np.array([link.fallback for link in links], np.uint8))
+    return dict(solid=solid, **tables)

src/CelerisLab/lbm/refinement_stepper.py

@@ -15,7 +15,8 @@ import numpy as np
 import pycuda.driver as cuda
 
 from ..cuda import compiler
-from .blocks import BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config
+from .blocks import (BlockLayout2D, launch_curved_boundary, launch_halo_copy, level_config,
+                     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
@@ -23,6 +24,7 @@ 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)
                                 blocked_kind, projected_explosion_tables, _prepare_blocked_boundary,
+                                _BlockedBoundary,
                                 BLOCKED_BOUNDARY_DEVICE_COLUMNS)
 from .descriptors import SOLID
 from .raw_exchange import summarize_exchange
@@ -68,14 +70,16 @@ def _level_name(level):
 def _check_levels(refinements, blocks, objects, jets):
     """Candidate limits of the level lists, checked before any device work."""
     depth = len(refinements) + 1
-    if depth < 2:
-        raise ValueError("at least one refinement level is required")
     for name, values in (("level_blocks", blocks), ("level_objects", objects), ("level_jets", jets)):
         if len(values) != depth:
             raise ValueError(f"{name} needs one entry per level ({depth})")
     for level, block in enumerate(blocks):
         if block is not None and len(block) != 2:
             raise ValueError(f"level_blocks[{level}] must be (by, bx)")
+    if depth == 1:  # one blocked level: the root carries the static circles
+        if blocks[0] is None or jets[0]:
+            raise ValueError("a single level must be blocked and has no Jets")
+        return
     if objects[0] or jets[0]:
         raise ValueError("base-level bodies/Jets are not supported")
     if any(blocks[level - 1] is not None and blocks[level] is None for level in range(1, depth)):
@@ -98,6 +102,9 @@ class FixedPatchStepper:
     level-wide grid, public fine arrays (2*ny, 2*nx), NaN outside ACTIVE cells);
     see lbm.hierarchy. ``level_blocks``, ``level_objects`` and ``level_jets``
     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).
 
     ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
     a bottom on-site slot and top fixed-density reservoir from ``cfg``.
@@ -196,6 +203,10 @@ class FixedPatchStepper:
         self._ddf_shifting = bool(cfg.ddf_shifting)
         if self._ddf_shifting and any(b is None for b in level_blocks):
             raise ValueError("shifted populations require all-blocked levels; dense/mixed is unsupported")
+        if depth == 1 and (coarse_boundary != "periodic" or cfg.collision != "SRT" or cfg.les_enabled
+                           or observe_raw_exchange or observe_window_exchange):
+            raise ValueError("a single level requires a periodic SRT root without LES or exchange "
+                             "observation")
         self._channel = coarse_boundary == "channel"
         if self._channel:
             validate_channel_config(cfg, level_blocks)
@@ -229,12 +240,14 @@ class FixedPatchStepper:
             base = blocked_base_tables(pairs[0], self._layouts[0], np.argsort(pairs[0].coalesce_coarse))
             validate_channel_transfers(pairs[0], self._layouts[0], base[1])
         # Kernel configurations per level (geometry configs live in self.hierarchy).
-        fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].fine_shape[0],
+        self._level_cfgs = [replace(cfg)]
-                           viscosity=2 * cfg.viscosity)
+        if pairs:
-        if self._wall_aware:
+            fine_cfg = replace(cfg, nx=pairs[0].fine_shape[1], ny=pairs[0].fine_shape[0],
-            fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                               viscosity=2 * cfg.viscosity)
-                               wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+            if self._wall_aware:
-        self._level_cfgs = [replace(cfg), fine_cfg]
+                fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                                   wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+            self._level_cfgs.append(fine_cfg)
         for pair in pairs[1:]:
             parent = self._level_cfgs[-1]
             self._level_cfgs.append(replace(parent, nx=pair.fine_shape[1], ny=pair.fine_shape[0],
@@ -275,12 +288,18 @@ class FixedPatchStepper:
         self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects)
                                        for objects in level_objects)
         self._level_jets = level_jets
-        self._level_boundaries = [None] + [
+        root = None
+        if depth == 1 and level_objects[0]:
+            tables = static_circle_tables(self._layouts[0], level_objects[0])
+            root = _BlockedBoundary(None, tables, (), False)
+        self._level_boundaries = [root] + [
             _prepare_blocked_boundary(self._boundary_hosts[level], self._layouts[level],
                                       self._level_cfgs[level], level_jets[level],
                                       _shifted_blocked=self._ddf_shifting)
             for level in range(1, depth)]
         self._leaf = [kind == ACTIVE for kind in self.hierarchy.kinds]
+        if root is not None:
+            self._leaf[0] &= ~root.tables["solid"]
         for level in range(1, depth):
             host = self._boundary_hosts[level]
             if host is not None:
@@ -324,7 +343,9 @@ class FixedPatchStepper:
                                 b.window_integral = self._allocate(count * 4 * 8)
                             b.accumulate = level.module.get_function("BlockAccumulateExchange")
                         b.kernel = level.module.get_function("BlockCurvedBoundary")
-                if self._level_boundaries[1] is not None:
+                if depth == 1:
+                    pass
+                elif self._level_boundaries[1] is not None:
                     b = self._level_boundaries[1]
                     self._fine_flags, self._link_tables, self._curved = b.flags, b.device, b.kernel
                 elif self._boundary_hosts[1] is not None:
@@ -350,7 +371,8 @@ class FixedPatchStepper:
         """Compile one level's step module and allocate its two G buffers."""
         layout, config = self._layouts[index], self._level_cfgs[index]
         pairs = self.hierarchy.pairs
-        tau = pairs[0].tau_coarse if index == 0 else pairs[index - 1].tau_fine
+        tau = (pairs[0].tau_coarse if pairs else 3 * config.viscosity + .5) if index == 0 \
+            else pairs[index - 1].tau_fine
         name, source = ("PatchCoarseStep" if index == 0 else "PatchFineStep"), "refinement.cu"
         if layout is not None:
             config = level_config(config, layout, tau)
@@ -732,7 +754,8 @@ class FixedPatchStepper:
                              np.float64(h), np.float64(h),
                              grid=((count + self._block[0] - 1) // self._block[0], 1), **kw)
         flags = np.uintp(0) if level.boundary is None else level.boundary.flags
-        level.kernel(src, dst, np.uint32(level.layout.halo), flags, level.kind, grid=level.grid, **kw)
+        kind = np.uintp(0) if level.kind is None else level.kind  # single level: all ACTIVE
+        level.kernel(src, dst, np.uint32(level.layout.halo), flags, kind, grid=level.grid, **kw)
 
     def _gather_halo(self, level, buffer, **kw):
         pairs = len(level.layout.halo_dst)

tests/validation/multigrid/test_block_stepper.py

@@ -25,6 +25,7 @@ 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.refinement_stepper import FixedPatchStepper
 from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
 
 ROOT = Path(__file__).resolve().parents[3]
@@ -118,6 +119,23 @@ class CpuContracts(unittest.TestCase):
     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])),
+                        (None, dict(coarse_boundary="channel")),
+                        (None, dict(coarse_boundary="slot_reservoir")),
+                        (config(collision="MRT"), {}),
+                        (None, dict(observe_raw_exchange=True)),
+                        (None, dict(level_jets=[(object(),)])),
+                        (None, dict(level_objects=[(circle((15.5, 4.0), 4.3),)]))):  # touches edge
+            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))),
@@ -139,7 +157,7 @@ class CpuContracts(unittest.TestCase):
             with self.subTest(center=center, radius=radius), self.assertRaises(ValueError):
                 PeriodicBlockStepper(config(), layout, context=self.forbidden_context(),
                                      objects=[circle(center, radius)])
-        with patch.object(implementation, "_static_circles") as prepare:
+        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(()))
             prepare.assert_not_called()
@@ -177,7 +195,7 @@ class CpuContracts(unittest.TestCase):
 
     def test_static_circle_tables(self):
         layout = BlockLayout2D.build(NY, NX, 12, 16)
-        tables = implementation._static_circles(layout, (circle(),))
+        tables = implementation.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))
@@ -331,6 +349,29 @@ 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 in [(split, objects) for split in SPLITS
+                                      for objects in ((), (circle(),))]:
+                with self.subTest(split=(by, bx), bodies=len(objects)):
+                    layout = BlockLayout2D.build(NY, NX, by, bx)
+                    with PeriodicBlockStepper(config(), layout, context=context,
+                                              objects=objects) as backend:
+                        backend.initialize(g0)
+                        backend.step(STEPS)
+                        expected = backend.download()
+                    with FixedPatchStepper(config(), 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()