From a64af5beceece3f2a1415d37671d95bf7a33ae65 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 02:30:40 +0000 Subject: [PATCH 1/5] Add CPU BlockLayout2D: core ownership, halo pairs, blocked packing Fixed-size 2D blocks, block-contiguous Q-major storage, explicit per-axis periodicity; non-periodic exterior halo cells stay unpaired (NaN). Provides owner/halo tables for a later periodic GPU consumer. Cut-link mapping deliberately omitted until the BFL write policy is settled. Tests: ownership, pair invariants, pack round trip, tagged pull vs dense. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/blocks.py | 173 ++++++++++++++++++++++++++++++++ tests/unit/test_block_layout.py | 139 +++++++++++++++++++++++++ 2 files changed, 312 insertions(+) create mode 100644 src/CelerisLab/lbm/blocks.py create mode 100644 tests/unit/test_block_layout.py diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py new file mode 100644 index 0000000..b7885d8 --- /dev/null +++ b/src/CelerisLab/lbm/blocks.py @@ -0,0 +1,173 @@ +"""CPU layout for fixed-size 2D blocks: core ownership, halo pairs and packing. + +Dense arrays are Q-major (q, y, x) with cell y*nx+x, as in core/layout.cuh. +Blocked storage is block-contiguous and Q-major inside each block: +(nblocks, q, by+2*halo, bx+2*halo). Block b = jy*nbx + jx has core origin +(jy*by, jx*bx). Every dense cell, fluid or solid, has exactly one core owner. + +A "block cell" id is b*L + c with L=(by+2h)*(bx+2h) and c the local row-major +cell; its population q lives at flat slot b*NQ*L + q*L + c. Halo pairs copy +all NQ populations from an owner core cell into another block's halo cell. +Periodicity is explicit per axis. Across a non-periodic exterior, halo cells +have no pair and stay invalid (NaN after packing): no wrap, no implied BC. + +Preparation only: no CUDA, no boundary conditions, no refinement levels. +""" +from dataclasses import dataclass +from numbers import Integral + +import numpy as np + +NQ = 9 +FACE, CORNER = 0, 1 + + +def _readonly(values, dtype): + array = np.ascontiguousarray(values, dtype=dtype) + array.setflags(write=False) + return array + + +def _positive_int(value, name): + if isinstance(value, (bool, np.bool_)) or not isinstance(value, Integral) or value <= 0: + raise ValueError(f"{name} must be a positive integer") + return int(value) + + +@dataclass(frozen=True, eq=False) +class BlockLayout2D: + """Fixed-size block partition of an (ny, nx) grid; construct with build(). + + ``origins`` rows are (y0, x0) core origins. ``owner`` maps each dense cell + to its owning block cell id. ``halo_dst``/``halo_src`` are block cell ids + (destination halo cell, source owner core cell) sorted by destination; + ``halo_kind`` marks FACE or CORNER pairs. Sources are never halo cells, so + one gather in any order fills every valid halo cell. + """ + shape: tuple[int, int] + block: tuple[int, int] + halo: int + periodic: tuple[bool, bool] + origins: np.ndarray + owner: np.ndarray + halo_dst: np.ndarray + halo_src: np.ndarray + halo_kind: np.ndarray + + @classmethod + def build(cls, ny, nx, by, bx, *, halo=1, periodic=(True, True)) -> "BlockLayout2D": + """``periodic`` is (periodic_x, periodic_y); block sizes must divide the grid.""" + ny, nx, by, bx, halo = (_positive_int(v, n) for v, n in + ((ny, "ny"), (nx, "nx"), (by, "by"), (bx, "bx"), (halo, "halo"))) + if ny % by or nx % bx: + raise ValueError("block sizes must divide the grid") + if halo > min(by, bx): + raise ValueError("halo may not exceed the block size") + if len(periodic) != 2 or not all(isinstance(p, (bool, np.bool_)) for p in periodic): + raise ValueError("periodic must be two booleans (x, y)") + px, py = map(bool, periodic) + nby, nbx = ny // by, nx // bx + ly, lx = by + 2 * halo, bx + 2 * halo + cells = ly * lx + if nby * nbx * NQ * cells > np.iinfo(np.int64).max // 2: + raise ValueError("blocked storage exceeds 64-bit indexing") + + jy, jx = np.divmod(np.arange(nby * nbx), nbx) + origins = np.stack([jy * by, jx * bx], axis=1) + + gy, gx = np.mgrid[0:ny, 0:nx] + owner = ((gy // by) * nbx + gx // bx) * cells + (gy % by + halo) * lx + gx % bx + halo + + # Halo cells of one block in local row-major order, then all blocks. + hy, hx = np.mgrid[0:ly, 0:lx] + out_y = (hy < halo) | (hy >= by + halo) + out_x = (hx < halo) | (hx >= bx + halo) + ring = out_y | out_x + hy, hx, kind = hy[ring], hx[ring], np.where(out_y[ring] & out_x[ring], CORNER, FACE) + dst, src, kinds = [], [], [] + for b, (y0, x0) in enumerate(origins): + sy, sx = y0 + hy - halo, x0 + hx - halo + if py: + sy %= ny + if px: + sx %= nx + keep = (sy >= 0) & (sy < ny) & (sx >= 0) & (sx < nx) + dst.append(b * cells + hy[keep] * lx + hx[keep]) + src.append(owner[sy[keep], sx[keep]]) + kinds.append(kind[keep]) + return cls((ny, nx), (by, bx), halo, (px, py), + _readonly(origins, np.int64), _readonly(owner, np.uint64), + _readonly(np.concatenate(dst), np.uint64), + _readonly(np.concatenate(src), np.uint64), + _readonly(np.concatenate(kinds), np.uint8)) + + @property + def nblocks(self) -> int: + return len(self.origins) + + @property + def local_shape(self) -> tuple[int, int]: + return self.block[0] + 2 * self.halo, self.block[1] + 2 * self.halo + + @property + def blocked_shape(self) -> tuple[int, int, int, int]: + return (self.nblocks, NQ, *self.local_shape) + + def slot(self, block_cell, q): + """Flat blocked population slot of block cell id(s) for direction q.""" + cells = int(np.prod(self.local_shape)) + b, c = np.divmod(np.asarray(block_cell, dtype=np.uint64), np.uint64(cells)) + return (b * np.uint64(NQ) + np.asarray(q, dtype=np.uint64)) * np.uint64(cells) + c + + def block_cell(self, block, y, x): + """Address of in-grid cell (y, x) in ``block``'s own storage, core or halo. + + Periodic axes resolve to the image inside the block's halo window. + Raises if the cell is not represented there or lies across an exterior. + """ + by, bx = self.block + ly, lx = self.local_shape + ny, nx = self.shape + y0, x0 = self.origins[block] + ry, rx = int(y) - int(y0) + self.halo, int(x) - int(x0) + self.halo + if self.periodic[1] and not 0 <= ry < ly: + ry = next((r for r in (ry - ny, ry + ny) if 0 <= r < ly), ry) + if self.periodic[0] and not 0 <= rx < lx: + rx = next((r for r in (rx - nx, rx + nx) if 0 <= r < lx), rx) + if not (0 <= ry < ly and 0 <= rx < lx and 0 <= y < ny and 0 <= x < nx): + raise ValueError("cell is not represented in this block's storage") + return int(block) * ly * lx + ry * lx + rx + + def exchange(self, blocked) -> None: + """Fill every valid halo cell in place from its owner's core.""" + blocked = self._check_blocked(blocked) + flat = blocked.reshape(self.nblocks, NQ, -1) + cells = flat.shape[2] + db, dc = np.divmod(self.halo_dst.astype(np.intp), cells) + sb, sc = np.divmod(self.halo_src.astype(np.intp), cells) + flat[db, :, dc] = flat[sb, :, sc] + + def to_blocked(self, dense): + """Pack (NQ, ny, nx) floats; halo cells without a pair stay NaN.""" + dense = np.asarray(dense) + if dense.shape != (NQ, *self.shape) or not np.issubdtype(dense.dtype, np.floating): + raise ValueError(f"expected floating (NQ, ny, nx) = {(NQ, *self.shape)}") + blocked = np.full(self.blocked_shape, np.nan, dtype=dense.dtype) + flat = blocked.reshape(self.nblocks, NQ, -1) + ob, oc = np.divmod(self.owner.ravel().astype(np.intp), flat.shape[2]) + flat[ob, :, oc] = dense.reshape(NQ, -1).T + self.exchange(blocked) + return blocked + + def to_dense(self, blocked): + """Unpack owner cores only; halo contents are ignored.""" + blocked = self._check_blocked(blocked) + flat = blocked.reshape(self.nblocks, NQ, -1) + ob, oc = np.divmod(self.owner.ravel().astype(np.intp), flat.shape[2]) + return np.ascontiguousarray(flat[ob, :, oc].T).reshape(NQ, *self.shape) + + def _check_blocked(self, blocked): + if not isinstance(blocked, np.ndarray) or blocked.shape != self.blocked_shape \ + or not blocked.flags.c_contiguous: + raise ValueError(f"expected C-contiguous blocked array {self.blocked_shape}") + return blocked diff --git a/tests/unit/test_block_layout.py b/tests/unit/test_block_layout.py new file mode 100644 index 0000000..39d8623 --- /dev/null +++ b/tests/unit/test_block_layout.py @@ -0,0 +1,139 @@ +"""CPU block layout contracts: ownership, halo pairs, packing and tagged pull.""" +import sys +import unittest + +import numpy as np + +from CelerisLab.lbm.blocks import CORNER, FACE, NQ, BlockLayout2D +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY + +CASES = ( # (ny, nx, by, bx, halo) + (6, 8, 3, 4, 1), # 2x2 split + (6, 12, 3, 4, 1), # 2x3 split + (8, 8, 4, 4, 2), # wider halo + (4, 8, 4, 4, 1), # one block along y +) +PERIODIC = ((True, True), (False, False), (True, False), (False, True)) + + +def tags(ny, nx): + """Unique, exactly representable value per dense population slot.""" + return np.arange(NQ * ny * nx, dtype=np.float64).reshape(NQ, ny, nx) + + +def dense_pull(a, periodic): + """f[q](x) = a[q](x - c_q); NaN where a non-periodic exterior is crossed.""" + px, py = periodic + _, ny, nx = a.shape + gy, gx = np.mgrid[0:ny, 0:nx] + out = np.empty_like(a) + for q in range(NQ): + sy, sx = gy - D2Q9_EY[q], gx - D2Q9_EX[q] + valid = ((py | ((sy >= 0) & (sy < ny))) & (px | ((sx >= 0) & (sx < nx)))) + out[q] = np.where(valid, a[q, sy % ny, sx % nx], np.nan) + return out + + +def blocked_pull(layout, blocked): + """Exchange, then pull every core cell from its own block storage only.""" + layout.exchange(blocked) + out = np.full_like(blocked, np.nan) + (by, bx), h = layout.block, layout.halo + for q in range(NQ): + ey, ex = D2Q9_EY[q], D2Q9_EX[q] + out[:, q, h:h + by, h:h + bx] = blocked[:, q, h - ey:h - ey + by, h - ex:h - ex + bx] + return out + + +class BlockLayoutTests(unittest.TestCase): + def test_cpu_import_boundary(self): + self.assertFalse(any(name == "pycuda" or name.startswith("pycuda.") for name in sys.modules)) + + def test_origins_and_unique_core_ownership(self): + layout = BlockLayout2D.build(6, 8, 3, 4) + np.testing.assert_array_equal(layout.origins, [[0, 0], [0, 4], [3, 0], [3, 4]]) + self.assertEqual(layout.blocked_shape, (4, NQ, 5, 6)) + for ny, nx, by, bx, h in CASES: + layout = BlockLayout2D.build(ny, nx, by, bx, halo=h) + ly, lx = layout.local_shape + core = [b * ly * lx + (y + h) * lx + x + h + for b in range(layout.nblocks) for y in range(by) for x in range(bx)] + self.assertEqual(sorted(layout.owner.ravel().tolist()), sorted(core)) + for y in range(ny): + for x in range(nx): + b = (y // by) * (nx // bx) + x // bx + self.assertEqual(layout.block_cell(b, y, x), layout.owner[y, x]) + + def test_halo_pairs_are_owner_core_to_unique_halo(self): + for ny, nx, by, bx, h in CASES: + for periodic in PERIODIC: + layout = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic) + core = set(layout.owner.ravel().tolist()) + dst = layout.halo_dst.tolist() + self.assertEqual(dst, sorted(dst)) + self.assertEqual(len(dst), len(set(dst))) + self.assertTrue(core.isdisjoint(dst)) + self.assertTrue(set(layout.halo_src.tolist()) <= core) + self.assertTrue(set(layout.halo_kind.tolist()) <= {FACE, CORNER}) + full = BlockLayout2D.build(6, 8, 3, 4) + per_block = np.bincount(full.halo_kind, minlength=2) // full.nblocks + np.testing.assert_array_equal(per_block, [2 * (3 + 4), 4]) + closed = BlockLayout2D.build(6, 8, 3, 4, periodic=(False, False)) + # Each 2x2 block keeps one inner face column, one face row and one corner. + self.assertEqual(len(closed.halo_dst), 4 * (3 + 4 + 1)) + + def test_pack_round_trip_and_exterior_stays_invalid(self): + for ny, nx, by, bx, h in CASES: + for periodic in PERIODIC: + layout = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic) + a = tags(ny, nx) + blocked = layout.to_blocked(a) + np.testing.assert_array_equal(layout.to_dense(blocked), a) + filled = np.isfinite(blocked[:, 0]).sum() + self.assertEqual(filled, ny * nx + len(layout.halo_dst)) + flat = blocked.ravel() + for q in (0, 5): + np.testing.assert_array_equal(flat[layout.slot(layout.halo_dst, q).astype(np.intp)], + flat[layout.slot(layout.halo_src, q).astype(np.intp)]) + + def test_tagged_pull_matches_dense_pull(self): + for ny, nx, by, bx, h in CASES: + for periodic in PERIODIC: + with self.subTest(case=(ny, nx, by, bx, h), periodic=periodic): + layout = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic) + a = tags(ny, nx) + got = layout.to_dense(blocked_pull(layout, layout.to_blocked(a))) + np.testing.assert_array_equal(got, dense_pull(a, periodic)) + + def test_block_cell_resolves_periodic_halo_images(self): + layout = BlockLayout2D.build(6, 8, 3, 4) + ly, lx = layout.local_shape + # Block 0 sees x=7 (wrapped) in its west halo and y=5 in its south halo. + self.assertEqual(layout.block_cell(0, 0, 7), 1 * lx + 0) + self.assertEqual(layout.block_cell(0, 5, 0), 0 * lx + 1) + closed = BlockLayout2D.build(6, 8, 3, 4, periodic=(False, False)) + with self.assertRaises(ValueError): + closed.block_cell(0, 0, 7) + with self.assertRaises(ValueError): + layout.block_cell(0, 4, 0) # two rows past the one-cell halo + + def test_invalid_inputs_rejected(self): + for args, kw in (((6, 8, 4, 4), {}), ((6, 8, 3, 4), {"halo": 0}), + ((6, 8, 3, 4), {"halo": 4}), ((6, 8, 3, 4), {"periodic": (1, 1)}), + ((6, 8, 3, True), {}), ((6, 8.0, 3, 4), {})): + with self.assertRaises(ValueError): + BlockLayout2D.build(*args, **kw) + layout = BlockLayout2D.build(6, 8, 3, 4) + with self.assertRaises(ValueError): + layout.to_blocked(np.zeros((NQ, 6, 8), dtype=np.int64)) + with self.assertRaises(ValueError): + layout.to_dense(np.zeros((4, NQ, 5, 5))) + + def test_tables_are_readonly(self): + layout = BlockLayout2D.build(6, 8, 3, 4) + for array in (layout.origins, layout.owner, layout.halo_dst, layout.halo_src, layout.halo_kind): + self.assertFalse(array.flags.writeable) + + +if __name__ == "__main__": + unittest.main() From 2aadffe18328ee2dcf20498729420bb1b58c7bc2 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 02:33:00 +0000 Subject: [PATCH 2/5] Add receiver-local static cut-link mapping to BlockLayout2D Each physical link belongs once to the block owning its fluid receiver core; the reflected source population is written into that block's own core/halo copy (no halo copy between BFL and pull). Rejects repeated links and repeated local writers; one row per link so force counts once. Documents stale remote solid-core slots vs unique core ownership, and that the 4-cell refinement halo is separate. Static geometry only. Tests: address equivalence incl. cross-block q<.5 donors, and static BFL+pull bit-equal to the dense operator on fluid cells. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/blocks.py | 74 +++++++++++++++++++++ tests/unit/test_block_layout.py | 110 +++++++++++++++++++++++++++++++- 2 files changed, 183 insertions(+), 1 deletion(-) diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py index b7885d8..3b9759d 100644 --- a/src/CelerisLab/lbm/blocks.py +++ b/src/CelerisLab/lbm/blocks.py @@ -12,12 +12,25 @@ Periodicity is explicit per axis. Across a non-periodic exterior, halo cells have no pair and stay invalid (NaN after packing): no wrap, no implied BC. Preparation only: no CUDA, no boundary conditions, no refinement levels. +The one-cell same-level halo here is unrelated to the four-cell fine halo of +the fixed refinement patch (refinement.py), whose semantics are unchanged. + +Static cut links (map_cut_links) use RECEIVER-LOCAL staging: each physical +link belongs once to the block whose core holds its fluid receiver. Order per +step: gather halos -> BFL writes the reflected source population into that +block's own core/halo copy of the solid cell -> pull/collide core cells, with +NO halo copy between BFL and pull. Reflected populations are temporary, +receiver-local boundary values: the solid's owner-core copy of that slot may +be stale, which is distinct from the unique core ownership of stored state. +No moving geometry is implied. """ from dataclasses import dataclass from numbers import Integral import numpy as np +from .descriptors import D2Q9_EX, D2Q9_EY + NQ = 9 FACE, CORNER = 0, 1 @@ -34,6 +47,23 @@ def _positive_int(value, name): return int(value) +@dataclass(frozen=True, eq=False) +class BlockCutLinks: + """One row per input link, in input order (so force is counted once). + + ``block`` owns the receiver core cell; ``receiver`` and ``donor`` are block + cell ids in that block (donor = NO_DONOR unless q<.5 with a legal donor); + ``write`` is the flat blocked slot of the reflected population. + """ + block: np.ndarray + receiver: np.ndarray + donor: np.ndarray + write: np.ndarray + + +NO_DONOR = np.iinfo(np.uint64).max + + @dataclass(frozen=True, eq=False) class BlockLayout2D: """Fixed-size block partition of an (ny, nx) grid; construct with build(). @@ -138,6 +168,50 @@ class BlockLayout2D: raise ValueError("cell is not represented in this block's storage") return int(block) * ly * lx + ry * lx + rx + def _neighbor(self, y, x, dy, dx): + ny, nx = self.shape + y, x = y + dy, x + dx + if self.periodic[1]: + y %= ny + if self.periodic[0]: + x %= nx + if not (0 <= y < ny and 0 <= x < nx): + raise ValueError("cut-link stencil crosses a non-periodic exterior") + return y, x + + def map_cut_links(self, links) -> BlockCutLinks: + """Map dense static cut links (fluid_idx, dir, q, fallback) to blocks. + + Rejects repeated physical links and repeated local writers. Addresses + are geometry-free: q, fallback and BFL math stay as built on the dense + grid. Halo must be at least one cell (always true here). + """ + ny, nx = self.shape + cells = int(np.prod(self.local_shape)) + rows, seen, writers = [], set(), set() + for link in links: + r, q = int(link.fluid_idx), int(link.dir) + if not (0 <= r < ny * nx and 1 <= q < NQ) or (r, q) in seen: + raise ValueError("cut link is invalid or repeated") + seen.add((r, q)) + y, x = divmod(r, nx) + receiver = int(self.owner[y, x]) + block = receiver // cells + sy, sx = self._neighbor(y, x, D2Q9_EY[q], D2Q9_EX[q]) + opp = q + 1 if q % 2 else q - 1 + write = int(self.slot(self.block_cell(block, sy, sx), opp)) + if write in writers: + raise ValueError("duplicate local boundary writer") + writers.add(write) + donor = NO_DONOR + if float(link.q) < .5 and int(link.fallback) == 0: + dy, dx = self._neighbor(y, x, -D2Q9_EY[q], -D2Q9_EX[q]) + donor = self.block_cell(block, dy, dx) + rows.append((block, receiver, donor, write)) + table = np.array(rows, dtype=np.uint64).reshape(-1, 4) + return BlockCutLinks(_readonly(table[:, 0], np.uint32), + *(_readonly(table[:, i], np.uint64) for i in (1, 2, 3))) + def exchange(self, blocked) -> None: """Fill every valid halo cell in place from its owner's core.""" blocked = self._check_blocked(blocked) diff --git a/tests/unit/test_block_layout.py b/tests/unit/test_block_layout.py index 39d8623..eb3be45 100644 --- a/tests/unit/test_block_layout.py +++ b/tests/unit/test_block_layout.py @@ -4,7 +4,9 @@ import unittest import numpy as np -from CelerisLab.lbm.blocks import CORNER, FACE, NQ, BlockLayout2D +from CelerisLab.body.geometry.base import CutLink +from CelerisLab.body.geometry.circle import CircleGeometry +from CelerisLab.lbm.blocks import CORNER, FACE, NO_DONOR, NQ, BlockLayout2D from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY CASES = ( # (ny, nx, by, bx, halo) @@ -45,6 +47,58 @@ def blocked_pull(layout, blocked): return out +def opp(q): + return q + 1 if q % 2 else q - 1 + + +def reflect(q, toward, other): + """Static linear Bouzidi stand-in; equivalence holds for any read->write map.""" + if q < .5: + return 2 * q * toward + (1 - 2 * q) * other + return toward / (2 * q) + (2 * q - 1) / (2 * q) * other + + +def dense_bfl(a, links): + out = a.copy() + _, ny, nx = a.shape + for link in links: + y, x = divmod(link.fluid_idx, nx) + q, d = link.q, link.dir + if q < .5 and link.fallback == 0: + other = a[d, (y - D2Q9_EY[d]) % ny, (x - D2Q9_EX[d]) % nx] + elif q < .5: + other = 0.0 + else: + other = a[opp(d), y, x] + out[opp(d), (y + D2Q9_EY[d]) % ny, (x + D2Q9_EX[d]) % nx] = reflect(q, a[d, y, x], other) + return out + + +def blocked_bfl(layout, blocked, links, mapped): + """Reads and writes only the receiver block's storage; no exchange inside.""" + flat = blocked.ravel() + source = flat.copy() + for link, receiver, donor, write in zip(links, mapped.receiver, mapped.donor, mapped.write): + q, d = link.q, link.dir + toward = source[int(layout.slot(receiver, d))] + if donor != NO_DONOR: + other = source[int(layout.slot(donor, d))] + elif q < .5: + other = 0.0 + else: + other = source[int(layout.slot(receiver, opp(d)))] + flat[int(write)] = reflect(q, toward, other) + + +def seam_circle(n=16): + """Circle on the shared corner of a 2x2 split (integer-node frame).""" + geometry = CircleGeometry(n / 2 - .5, n / 2 - .5, 3.3) + links = geometry.build_cut_links(n, n) + gy, gx = np.mgrid[0:n, 0:n] + solid = (gx - (n / 2 - .5)) ** 2 + (gy - (n / 2 - .5)) ** 2 < 3.3 ** 2 + return links, solid + + class BlockLayoutTests(unittest.TestCase): def test_cpu_import_boundary(self): self.assertFalse(any(name == "pycuda" or name.startswith("pycuda.") for name in sys.modules)) @@ -129,6 +183,60 @@ class BlockLayoutTests(unittest.TestCase): with self.assertRaises(ValueError): layout.to_dense(np.zeros((4, NQ, 5, 5))) + def test_cut_links_are_receiver_local_and_unique(self): + links, _ = seam_circle() + for periodic in ((True, True), (False, False)): + layout = BlockLayout2D.build(16, 16, 8, 8, periodic=periodic) + mapped = layout.map_cut_links(links) + self.assertEqual(len(mapped.write), len(links)) # one row, one force per link + self.assertEqual(len(set(mapped.write.tolist())), len(links)) + self.assertEqual(len(set(mapped.block.tolist())), 4) # seams are crossed + cells = int(np.prod(layout.local_shape)) + core = set(layout.owner.ravel().tolist()) + a = tags(16, 16) + flat = layout.to_blocked(a).ravel() + halo_reads = 0 + for link, b, r, d, w in zip(links, mapped.block, mapped.receiver, + mapped.donor, mapped.write): + y, x = divmod(link.fluid_idx, 16) + self.assertEqual(r, layout.owner[y, x]) + self.assertEqual(int(r) // cells, b) + self.assertEqual(int(w) // (NQ * cells), b) + sy, sx = y + D2Q9_EY[link.dir], x + D2Q9_EX[link.dir] + self.assertEqual(flat[int(w)], a[opp(link.dir), sy, sx]) + if d != NO_DONOR: + self.assertEqual(int(d) // cells, b) + self.assertEqual(flat[int(layout.slot(d, link.dir))], + a[link.dir, y - D2Q9_EY[link.dir], x - D2Q9_EX[link.dir]]) + halo_reads += int(d) not in core + self.assertGreater(halo_reads, 0) # some q<.5 donors come from a neighbor's core + + def test_static_bfl_pull_matches_dense_on_fluid_cells(self): + links, solid = seam_circle() + a = tags(16, 16) / 7.0 + expected = dense_pull(dense_bfl(a, links), (True, True)) + layout = BlockLayout2D.build(16, 16, 8, 8) + blocked = layout.to_blocked(a) # gather halos + blocked_bfl(layout, blocked, links, layout.map_cut_links(links)) + # Pull core cells without another halo copy (it would undo halo writes). + out = np.full_like(blocked, np.nan) + (by, bx), h = layout.block, layout.halo + for q in range(NQ): + ey, ex = D2Q9_EY[q], D2Q9_EX[q] + out[:, q, h:h + by, h:h + bx] = blocked[:, q, h - ey:h - ey + by, h - ex:h - ex + bx] + got = layout.to_dense(out) + np.testing.assert_array_equal(got[:, ~solid], expected[:, ~solid]) + + def test_cut_link_rejections(self): + layout = BlockLayout2D.build(16, 16, 8, 8, periodic=(False, False)) + link = CutLink(fluid_idx=5 * 16 + 5, dir=1, q=.6, body_id=0) + with self.assertRaises(ValueError): + layout.map_cut_links([link, link]) + with self.assertRaises(ValueError): + layout.map_cut_links([CutLink(fluid_idx=15, dir=1, q=.6, body_id=0)]) + with self.assertRaises(ValueError): + layout.map_cut_links([CutLink(fluid_idx=5 * 16 + 5, dir=0, q=.6, body_id=0)]) + def test_tables_are_readonly(self): layout = BlockLayout2D.build(6, 8, 3, 4) for array in (layout.origins, layout.owner, layout.halo_dst, layout.halo_src, layout.halo_kind): From 2a68d7e7c73708eb1340666a70bf05b5e005bc96 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 02:38:13 +0000 Subject: [PATCH 3/5] Fix periodic-seam cut-link addressing in map_cut_links Source (write) and reverse-donor addresses are now the receiver's local core address +/- (ey*local_nx + ex), i.e. the window image the next local pull reads. block_cell preferred the in-window core image on a one-block periodic axis, so BFL wrote a copy the pull never reads. Physical neighbors are still formed only to reject non-periodic exteriors. block_cell semantics unchanged; docstring notes the image ambiguity. Adds an end-to-end synthetic seam test (x=0/x=nx-1/corner wraps) on one-block-x, one-block-y and 1x1 periodic layouts. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/blocks.py | 17 ++++++++++++----- tests/unit/test_block_layout.py | 34 ++++++++++++++++++++++++++------- 2 files changed, 39 insertions(+), 12 deletions(-) diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py index 3b9759d..fabdaa1 100644 --- a/src/CelerisLab/lbm/blocks.py +++ b/src/CelerisLab/lbm/blocks.py @@ -152,7 +152,9 @@ class BlockLayout2D: def block_cell(self, block, y, x): """Address of in-grid cell (y, x) in ``block``'s own storage, core or halo. - Periodic axes resolve to the image inside the block's halo window. + Periodic axes resolve to the image inside the block's halo window; + if two images fit, the core-side one is returned (use local offsets, + not this lookup, when a stencil must follow a specific image). Raises if the cell is not represented there or lies across an exterior. """ by, bx = self.block @@ -187,6 +189,7 @@ class BlockLayout2D: grid. Halo must be at least one cell (always true here). """ ny, nx = self.shape + lx = self.local_shape[1] cells = int(np.prod(self.local_shape)) rows, seen, writers = [], set(), set() for link in links: @@ -197,16 +200,20 @@ class BlockLayout2D: y, x = divmod(r, nx) receiver = int(self.owner[y, x]) block = receiver // cells - sy, sx = self._neighbor(y, x, D2Q9_EY[q], D2Q9_EX[q]) + # Local stencil offsets from the core receiver: the window image the + # next local pull reads, even where a periodic seam wraps physically. + # Physical neighbors are formed only to reject non-periodic exteriors. + step = D2Q9_EY[q] * lx + D2Q9_EX[q] + self._neighbor(y, x, D2Q9_EY[q], D2Q9_EX[q]) opp = q + 1 if q % 2 else q - 1 - write = int(self.slot(self.block_cell(block, sy, sx), opp)) + write = int(self.slot(receiver + step, opp)) if write in writers: raise ValueError("duplicate local boundary writer") writers.add(write) donor = NO_DONOR if float(link.q) < .5 and int(link.fallback) == 0: - dy, dx = self._neighbor(y, x, -D2Q9_EY[q], -D2Q9_EX[q]) - donor = self.block_cell(block, dy, dx) + self._neighbor(y, x, -D2Q9_EY[q], -D2Q9_EX[q]) + donor = receiver - step rows.append((block, receiver, donor, write)) table = np.array(rows, dtype=np.uint64).reshape(-1, 4) return BlockCutLinks(_readonly(table[:, 0], np.uint32), diff --git a/tests/unit/test_block_layout.py b/tests/unit/test_block_layout.py index eb3be45..f38e6cb 100644 --- a/tests/unit/test_block_layout.py +++ b/tests/unit/test_block_layout.py @@ -39,6 +39,11 @@ def dense_pull(a, periodic): def blocked_pull(layout, blocked): """Exchange, then pull every core cell from its own block storage only.""" layout.exchange(blocked) + return blocked_pull_core(layout, blocked) + + +def blocked_pull_core(layout, blocked): + """Pull core cells with no halo copy (a copy would undo BFL halo writes).""" out = np.full_like(blocked, np.nan) (by, bx), h = layout.block, layout.halo for q in range(NQ): @@ -218,15 +223,30 @@ class BlockLayoutTests(unittest.TestCase): layout = BlockLayout2D.build(16, 16, 8, 8) blocked = layout.to_blocked(a) # gather halos blocked_bfl(layout, blocked, links, layout.map_cut_links(links)) - # Pull core cells without another halo copy (it would undo halo writes). - out = np.full_like(blocked, np.nan) - (by, bx), h = layout.block, layout.halo - for q in range(NQ): - ey, ex = D2Q9_EY[q], D2Q9_EX[q] - out[:, q, h:h + by, h:h + bx] = blocked[:, q, h - ey:h - ey + by, h - ex:h - ex + bx] - got = layout.to_dense(out) + got = layout.to_dense(blocked_pull_core(layout, blocked)) np.testing.assert_array_equal(got[:, ~solid], expected[:, ~solid]) + def test_periodic_seam_links_follow_the_pulled_halo_image(self): + # Receivers on x=0 / x=nx-1 / a corner whose sources wrap physically. + links = [CutLink(fluid_idx=1 * 6 + 0, dir=2, q=.6, body_id=0), + CutLink(fluid_idx=2 * 6 + 0, dir=2, q=.3, body_id=0), + CutLink(fluid_idx=3 * 6 + 5, dir=1, q=.3, body_id=0), + CutLink(fluid_idx=0 * 6 + 0, dir=6, q=.4, body_id=0)] + sources = np.zeros((8, 6), dtype=bool) + for y, x in ((1, 5), (2, 5), (3, 0), (7, 5)): + sources[y, x] = True + a = tags(8, 6) / 7.0 + expected = dense_pull(dense_bfl(a, links), (True, True)) + # One block along x, along y, along both (periodic self-images). + for by, bx in ((4, 6), (8, 3), (8, 6)): + with self.subTest(block=(by, bx)): + layout = BlockLayout2D.build(8, 6, by, bx) + mapped = layout.map_cut_links(links) + blocked = layout.to_blocked(a) + blocked_bfl(layout, blocked, links, mapped) + got = layout.to_dense(blocked_pull_core(layout, blocked)) + np.testing.assert_array_equal(got[:, ~sources], expected[:, ~sources]) + def test_cut_link_rejections(self): layout = BlockLayout2D.build(16, 16, 8, 8, periodic=(False, False)) link = CutLink(fluid_idx=5 * 16 + 5, dir=1, q=.6, body_id=0) From 118635aebd442572624c393ae256783b47e77962 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 02:39:15 +0000 Subject: [PATCH 4/5] Reject non-integral, bool and out-of-range block-layout indices One _index helper guards block_cell (block, y, x), slot (block_cell, q) and map_cut_links (fluid_idx, dir) before addressing; previously e.g. block_cell(-1, 0, 0) returned -49 and fractional values were truncated. Adds focused rejection tests and Astra's one-block seam reproducer (build(4,8,4,8), receiver (2,0), dir 2 -> west-halo slot 90). Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/blocks.py | 29 +++++++++++++++++++++-------- tests/unit/test_block_layout.py | 22 ++++++++++++++++++++++ 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py index fabdaa1..3b7ac30 100644 --- a/src/CelerisLab/lbm/blocks.py +++ b/src/CelerisLab/lbm/blocks.py @@ -41,6 +41,14 @@ def _readonly(values, dtype): return array +def _index(value, upper, name, lower=0): + """Integer index or index array in [lower, upper); rejects bools and floats.""" + array = np.asarray(value) + if array.dtype.kind not in "iu" or (array.size and (array.min() < lower or array.max() >= upper)): + raise ValueError(f"{name} must be integer index(es) in [{lower}, {upper})") + return array + + def _positive_int(value, name): if isinstance(value, (bool, np.bool_)) or not isinstance(value, Integral) or value <= 0: raise ValueError(f"{name} must be a positive integer") @@ -146,8 +154,10 @@ class BlockLayout2D: def slot(self, block_cell, q): """Flat blocked population slot of block cell id(s) for direction q.""" cells = int(np.prod(self.local_shape)) - b, c = np.divmod(np.asarray(block_cell, dtype=np.uint64), np.uint64(cells)) - return (b * np.uint64(NQ) + np.asarray(q, dtype=np.uint64)) * np.uint64(cells) + c + block_cell = _index(block_cell, self.nblocks * cells, "block_cell").astype(np.uint64) + q = _index(q, NQ, "q").astype(np.uint64) + b, c = np.divmod(block_cell, np.uint64(cells)) + return (b * np.uint64(NQ) + q) * np.uint64(cells) + c def block_cell(self, block, y, x): """Address of in-grid cell (y, x) in ``block``'s own storage, core or halo. @@ -160,15 +170,17 @@ class BlockLayout2D: by, bx = self.block ly, lx = self.local_shape ny, nx = self.shape + block = int(_index(block, self.nblocks, "block")) + y, x = int(_index(y, ny, "y")), int(_index(x, nx, "x")) y0, x0 = self.origins[block] - ry, rx = int(y) - int(y0) + self.halo, int(x) - int(x0) + self.halo + ry, rx = y - int(y0) + self.halo, x - int(x0) + self.halo if self.periodic[1] and not 0 <= ry < ly: ry = next((r for r in (ry - ny, ry + ny) if 0 <= r < ly), ry) if self.periodic[0] and not 0 <= rx < lx: rx = next((r for r in (rx - nx, rx + nx) if 0 <= r < lx), rx) - if not (0 <= ry < ly and 0 <= rx < lx and 0 <= y < ny and 0 <= x < nx): + if not (0 <= ry < ly and 0 <= rx < lx): raise ValueError("cell is not represented in this block's storage") - return int(block) * ly * lx + ry * lx + rx + return block * ly * lx + ry * lx + rx def _neighbor(self, y, x, dy, dx): ny, nx = self.shape @@ -193,9 +205,10 @@ class BlockLayout2D: cells = int(np.prod(self.local_shape)) rows, seen, writers = [], set(), set() for link in links: - r, q = int(link.fluid_idx), int(link.dir) - if not (0 <= r < ny * nx and 1 <= q < NQ) or (r, q) in seen: - raise ValueError("cut link is invalid or repeated") + r = int(_index(link.fluid_idx, ny * nx, "fluid_idx")) + q = int(_index(link.dir, NQ, "dir", lower=1)) + if (r, q) in seen: + raise ValueError("cut link is repeated") seen.add((r, q)) y, x = divmod(r, nx) receiver = int(self.owner[y, x]) diff --git a/tests/unit/test_block_layout.py b/tests/unit/test_block_layout.py index f38e6cb..1c4c45c 100644 --- a/tests/unit/test_block_layout.py +++ b/tests/unit/test_block_layout.py @@ -257,6 +257,28 @@ class BlockLayoutTests(unittest.TestCase): with self.assertRaises(ValueError): layout.map_cut_links([CutLink(fluid_idx=5 * 16 + 5, dir=0, q=.6, body_id=0)]) + def test_public_indices_are_rejected_not_truncated(self): + layout = BlockLayout2D.build(4, 8, 4, 8) # one block, nx=8, ny=4 + for args in ((-1, 0, 0), (1, 0, 0), (True, 0, 0), (0.0, 0, 0), (0, 0.5, 0), + (0, 1, 2.0), (0, -1, 0), (0, 4, 0), (0, 0, 8), (0, False, 0)): + with self.subTest(block_cell=args), self.assertRaises(ValueError): + layout.block_cell(*args) + cells = int(np.prod(layout.local_shape)) + for args in ((-1, 1), (cells, 1), (True, 1), (1.0, 1), (0, NQ), (0, -1), (0, 1.0), + (0, True), (np.array([0, cells]), 1)): + with self.subTest(slot=args), self.assertRaises(ValueError): + layout.slot(*args) + for idx, d in ((5.7, 1), (True, 1), (-1, 1), (32, 1), (5, 1.0), (5, 9), (5, True), (5, 0)): + with self.subTest(fluid_idx=idx, dir=d), self.assertRaises(ValueError): + layout.map_cut_links([CutLink(fluid_idx=idx, dir=d, q=.6, body_id=0)]) + + def test_one_block_seam_reproducer(self): + # Receiver (2,0) -x wraps to x=7; the write must land in the west halo. + layout = BlockLayout2D.build(4, 8, 4, 8) + mapped = layout.map_cut_links([CutLink(fluid_idx=2 * 8, dir=2, q=.6, body_id=0)]) + lx = layout.local_shape[1] + self.assertEqual(int(mapped.write[0]), int(layout.slot(3 * lx + 0, 1))) + def test_tables_are_readonly(self): layout = BlockLayout2D.build(6, 8, 3, 4) for array in (layout.origins, layout.owner, layout.halo_dst, layout.halo_src, layout.halo_kind): From d27d344f43bb61de1c80c297e8013a1183dcbf11 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 02:39:58 +0000 Subject: [PATCH 5/5] Add source-path header comment to blocks.py Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/blocks.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/CelerisLab/lbm/blocks.py b/src/CelerisLab/lbm/blocks.py index 3b7ac30..4fc9975 100644 --- a/src/CelerisLab/lbm/blocks.py +++ b/src/CelerisLab/lbm/blocks.py @@ -1,3 +1,4 @@ +# CelerisLab/lbm/blocks.py """CPU layout for fixed-size 2D blocks: core ownership, halo pairs and packing. Dense arrays are Q-major (q, y, x) with cell y*nx+x, as in core/layout.cuh.