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