Files
CelerisLab/tests/unit/test_block_layout.py
T
ClaudeandClaude Opus 5.5 aecee19c19 BlockLayout2D: copy caller-supplied origins
An int64 C-contiguous origins array was kept by reference (and frozen in
the caller's hands), so mutating it later changed layout.origins and
block_cell without changing owner/halos. Copy before marking read-only;
regression checks origins, owner and block_cell invariance.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
2026-10-16 03:15:37 +00:00

439 lines
22 KiB
Python

"""CPU block layout contracts: ownership, halo pairs, packing and tagged pull."""
import sys
import unittest
from types import SimpleNamespace
import numpy as np
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, blocked_link_tables
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)
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):
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
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))
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_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))
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)
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_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):
self.assertFalse(array.flags.writeable)
def all_origins(ny, nx, by, bx):
return [(y, x) for y in range(0, ny, by) for x in range(0, nx, bx)]
def link(y, x, d, q, nx=12, fallback=0):
return SimpleNamespace(fluid_idx=y * nx + x, dir=d, q=q, fallback=fallback)
# 6x12 grid of 3x4 blocks (2x3 block grid); stored L: (0,0), (0,4), (3,4).
L_ORIGINS = [(0, 0), (0, 4), (3, 4)]
SPARSE = ( # (ny, nx, by, bx, halo, origins)
(6, 12, 3, 4, 1, L_ORIGINS),
(6, 12, 3, 4, 1, [(3, 8)]), # N=1 corner block
(8, 8, 4, 4, 2, [(0, 4), (4, 0)]), # diagonal pair, wider halo
(6, 8, 3, 4, 1, [(0, 0), (0, 4), (3, 0)]), # missing one of 2x2
)
class SparseBlockLayoutTests(unittest.TestCase):
"""Explicit selected block origins: absent blocks are holes, not wrapped -1."""
def test_all_origins_reproduce_full_layout(self):
for ny, nx, by, bx, h in CASES + ((4, 8, 4, 8, 1),): # last: N=1 full grid
for periodic in PERIODIC:
with self.subTest(case=(ny, nx, by, bx, h), periodic=periodic):
full = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic)
explicit = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic,
origins=all_origins(ny, nx, by, bx))
for name in ("origins", "owner", "halo_dst", "halo_src", "halo_kind"):
np.testing.assert_array_equal(getattr(explicit, name), getattr(full, name))
self.assertEqual(getattr(explicit, name).dtype, getattr(full, name).dtype)
self.assertTrue(full.owned.all())
self.assertEqual(explicit.hole, full.hole)
def test_holes_and_ownership(self):
layout = BlockLayout2D.build(6, 12, 3, 4, periodic=(False, False), origins=L_ORIGINS)
np.testing.assert_array_equal(layout.origins, L_ORIGINS)
self.assertEqual(layout.nblocks, 3)
cells = int(np.prod(layout.local_shape))
self.assertEqual(layout.hole, 3 * cells)
expected = np.zeros((6, 12), dtype=bool)
expected[0:3, 0:8] = expected[3:6, 4:8] = True
np.testing.assert_array_equal(layout.owned, expected)
self.assertTrue((layout.owner[~expected] == layout.hole).all())
for b, (y0, x0) in enumerate(L_ORIGINS):
for y in range(y0, y0 + 3):
for x in range(x0, x0 + 4):
self.assertEqual(layout.block_cell(b, y, x), layout.owner[y, x])
owned = layout.owner[expected].tolist()
self.assertEqual(len(owned), len(set(owned)))
with self.assertRaises(ValueError): # the hole value is never a valid slot
layout.slot(layout.hole, 0)
def test_halo_pairs_skip_absent_blocks(self):
layout = BlockLayout2D.build(6, 12, 3, 4, periodic=(False, False), origins=L_ORIGINS)
# Block 0: 3 east faces + SE corner; block 1: 3 west + 4 south faces;
# block 2: 4 north faces + NW corner. Nothing reads (3..5, 0..3) or x>=8.
self.assertEqual(len(layout.halo_dst), 16)
np.testing.assert_array_equal(np.bincount(layout.halo_kind, minlength=2), [14, 2])
for ny, nx, by, bx, h, origins in SPARSE:
for periodic in PERIODIC:
with self.subTest(origins=origins, periodic=periodic):
layout = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic,
origins=origins)
core = set(layout.owner[layout.owned].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.assertLess(max(dst + [-1]), layout.hole)
def test_pack_unpack_and_pull_with_holes(self):
for ny, nx, by, bx, h, origins in SPARSE:
for periodic in PERIODIC:
with self.subTest(origins=origins, periodic=periodic):
layout = BlockLayout2D.build(ny, nx, by, bx, halo=h, periodic=periodic,
origins=origins)
owned = layout.owned
a = tags(ny, nx)
holed = a.copy()
holed[:, ~owned] = np.nan
blocked = layout.to_blocked(a)
np.testing.assert_array_equal(layout.to_dense(blocked), holed)
poisoned = a.copy()
poisoned[:, ~owned] = np.inf # hole values are ignored by packing
np.testing.assert_array_equal(layout.to_blocked(poisoned), blocked)
self.assertEqual(np.isfinite(blocked[:, 0]).sum(),
owned.sum() + len(layout.halo_dst))
expected = dense_pull(holed, periodic)
expected[:, ~owned] = np.nan
got = layout.to_dense(blocked_pull(layout, blocked))
np.testing.assert_array_equal(got, expected)
def test_cut_links_in_or_next_to_absent_blocks_are_rejected(self):
layout = BlockLayout2D.build(6, 12, 3, 4, periodic=(False, False), origins=L_ORIGINS)
ok = [link(1, 1, 3, .6), link(1, 7, 2, .6), link(1, 4, 2, .3)]
mapped = layout.map_cut_links(ok) # fully owned stencils, incl. across a seam
cells = int(np.prod(layout.local_shape))
np.testing.assert_array_equal(mapped.block, [0, 1, 1])
np.testing.assert_array_equal(mapped.receiver, [layout.owner[1, 1], layout.owner[1, 7],
layout.owner[1, 4]])
self.assertEqual(int(mapped.donor[2]), int(layout.owner[1, 5]))
# Same physical link with q<.5 is mapped separately (repeats are rejected).
donor = layout.map_cut_links([link(1, 1, 3, .3)]).donor
self.assertEqual(int(donor[0]), int(layout.owner[0, 1]))
self.assertTrue((mapped.write < np.uint64(layout.nblocks * NQ * cells)).all())
for bad in (link(4, 1, 1, .6), # receiver in an absent block
link(2, 1, 3, .6), # solid neighbor (3,1) absent
link(1, 7, 2, .3), # q<.5 donor (1,8) absent
link(2, 7, 1, .6)): # solid neighbor (2,8) absent
with self.subTest(fluid=bad.fluid_idx, dir=bad.dir, q=bad.q):
with self.assertRaises(ValueError):
layout.map_cut_links([bad])
def test_link_tables_store_owned_flags_only(self):
layout = BlockLayout2D.build(6, 12, 3, 4, periodic=(False, False), origins=L_ORIGINS)
flags = np.arange(1, 6 * 12 + 1, dtype=np.uint16)
tables, _ = blocked_link_tables(layout, flags, [], [], [], [])
owned = layout.owned.ravel()
stored = tables["flags"]
self.assertEqual(stored.shape, (layout.hole,))
np.testing.assert_array_equal(stored[layout.owner.ravel()[owned].astype(np.intp)], flags[owned])
self.assertEqual(np.count_nonzero(stored), owned.sum())
def test_caller_origins_are_copied(self):
origins = np.array([[0, 0], [0, 4]], dtype=np.int64)
layout = BlockLayout2D.build(6, 8, 3, 4, origins=origins)
owner, address = layout.owner.copy(), layout.block_cell(0, 0, 0)
self.assertTrue(origins.flags.writeable) # caller array is not frozen
origins[0] = [3, 0]
np.testing.assert_array_equal(layout.origins, [[0, 0], [0, 4]])
np.testing.assert_array_equal(layout.owner, owner)
self.assertEqual(layout.block_cell(0, 0, 0), address)
self.assertEqual(address, layout.owner[0, 0])
def test_invalid_origins_rejected(self):
for origins in ([(0, 4), (0, 0)], # not row-major
[(3, 0), (0, 4)], # not row-major across rows
[(0, 0), (0, 0)], # repeated
[(1, 0)], [(0, 2)], # misaligned
[(6, 0)], [(0, 12)], [(-3, 0)], [(0, -4)], # outside
np.zeros((0, 2), dtype=np.int64), # empty
[(0.0, 0.0)], [(True, False)], [0, 0], [(0, 0, 0)]):
with self.subTest(origins=origins), self.assertRaises(ValueError):
BlockLayout2D.build(6, 12, 3, 4, origins=origins)
if __name__ == "__main__":
unittest.main()