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

Add receiver-local static cut-link mapping to BlockLayout2D

Browse Source 
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 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-15 02:33:00 +00:00
co-authored by Claude Opus 5.5
parent a64af5bece
commit 2aadffe183
2 changed files with 183 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/CelerisLab/lbm/blocks.py
+74
View File
@@ -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)
tests/unit/test_block_layout.py
+109 -1
View File
@@ -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):