From 2aadffe18328ee2dcf20498729420bb1b58c7bc2 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 15 Oct 2026 02:33:00 +0000 Subject: [PATCH] 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):