Integrate claude/block-layout for bounded validation
This commit is contained in:
@@ -0,0 +1,268 @@
|
||||
# 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.
|
||||
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.
|
||||
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
|
||||
|
||||
|
||||
def _readonly(values, dtype):
|
||||
array = np.ascontiguousarray(values, dtype=dtype)
|
||||
array.setflags(write=False)
|
||||
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")
|
||||
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().
|
||||
|
||||
``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))
|
||||
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.
|
||||
|
||||
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
|
||||
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 = 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):
|
||||
raise ValueError("cell is not represented in this block's storage")
|
||||
return 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
|
||||
lx = self.local_shape[1]
|
||||
cells = int(np.prod(self.local_shape))
|
||||
rows, seen, writers = [], set(), set()
|
||||
for link in links:
|
||||
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])
|
||||
block = receiver // cells
|
||||
# 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(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:
|
||||
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),
|
||||
*(_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)
|
||||
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
|
||||
@@ -0,0 +1,289 @@
|
||||
"""CPU block layout contracts: ownership, halo pairs, packing and tagged pull."""
|
||||
import sys
|
||||
import unittest
|
||||
|
||||
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
|
||||
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)
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user