ClaudeandClaude Opus 5.5
64665319a4
prepare_fine_boundary accepts a FixedPatchSet2D: objects are projected once onto the level-wide fine grid (FixedPatchSet2D.fine_origin = (0, 0), same coarse cell-edge -> fine cell transform), each circle must lie inside one patch, coarse transfer cells are checked against the covered union, and body ids stay object positions, so the existing jet one-body-per-jet, overlap and (receiver, dir) writer checks run once over the whole level (caller obj_id cannot bypass them). The stepper no longer rejects fine_objects/jets with several patches; blocked link tables, BFL and the jet amplitude cache are reused unchanged. No kernel change. README/docstring scope updated. Tests: CPU translation/ids/rejections (outside one patch, one jet over two bodies at zero amplitude); GPU one body via the new path == single-patch path (both modes), two bodies (+jets, amplitude switch) == tiled single (bitwise). Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
315 lines
16 KiB
Python
315 lines
16 KiB
Python
"""CPU geometry/contracts only: no time stepping, CUDA or reference edits."""
|
|
import hashlib
|
|
import importlib.util
|
|
from pathlib import Path
|
|
import sys
|
|
import unittest
|
|
|
|
import numpy as np
|
|
|
|
from CelerisLab.config import LBMConfig
|
|
from CelerisLab.lbm import refinement
|
|
from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY
|
|
from CelerisLab.lbm.refinement import ACTIVE, INNER_GHOST, OUTER_GHOST, FixedPatch2D
|
|
|
|
ROOT = Path(__file__).resolve().parents[2]
|
|
|
|
|
|
def config(**kwargs):
|
|
return LBMConfig(**dict(nx=16, ny=16, viscosity=.1, **kwargs))
|
|
|
|
|
|
class FixedPatchTests(unittest.TestCase):
|
|
def test_source_and_cpu_import_boundary(self):
|
|
self.assertEqual(Path(refinement.__file__).resolve(),
|
|
ROOT / "src/CelerisLab/lbm/refinement.py")
|
|
self.assertFalse(any(name == "pycuda" or name.startswith("pycuda.") for name in sys.modules))
|
|
self.assertNotIn("CelerisLab.simulation", sys.modules)
|
|
self.assertNotIn("CelerisLab.lbm.field", sys.modules)
|
|
|
|
def test_all_tables_match_pinned_reference_without_advancing(self):
|
|
path = ROOT / "tests/validation/multigrid/cpu_reference.py"
|
|
self.assertEqual(hashlib.sha256(path.read_bytes()).hexdigest(),
|
|
"bba2a06a9f6517b32079ffff15e4060d8c3b7ed35f16ca134e1877226c9e4e1e")
|
|
spec = importlib.util.spec_from_file_location("refinement_geometry_reference", path)
|
|
ref = importlib.util.module_from_spec(spec)
|
|
sys.modules[spec.name] = ref
|
|
spec.loader.exec_module(ref)
|
|
maps = ref.build_maps(ref.make_spec(collision="srt"))
|
|
p = FixedPatch2D.from_config(config(), (4, 4, 12, 12))
|
|
np.testing.assert_array_equal(p.fine_kind, maps.fine_class)
|
|
gy, gx = np.nonzero(~maps.fine_active)
|
|
np.testing.assert_array_equal(p.ghost_fine, gy * 24 + gx)
|
|
np.testing.assert_array_equal(p.ghost_coarse,
|
|
maps.fine_parent[gy] * 16 + maps.fine_parent[gx])
|
|
targets = [q * 256 + y * 16 + x for q, y, x in maps.coalesce]
|
|
sources = [[q * 576 + b * 24 + a
|
|
for b in (4 + 2 * (y - 4), 5 + 2 * (y - 4))
|
|
for a in (4 + 2 * (x - 4), 5 + 2 * (x - 4))]
|
|
for q, y, x in maps.coalesce]
|
|
np.testing.assert_array_equal(p.coalesce_coarse, targets)
|
|
np.testing.assert_array_equal(p.coalesce_fine, sources)
|
|
self.assertEqual(len(targets), 100)
|
|
self.assertEqual(len(set(p.coalesce_fine.ravel().tolist())), 400)
|
|
self.assertEqual(np.count_nonzero(p.fine_kind == ACTIVE), 256)
|
|
self.assertEqual(len(p.ghost_fine), 320)
|
|
self.assertEqual(self.corner_pass_count(p), 8)
|
|
|
|
@staticmethod
|
|
def corner_pass_count(p):
|
|
ny, nx = p.coarse_shape
|
|
x0, y0, x1, y1 = p.bounds
|
|
count = 0
|
|
for slot in p.coalesce_coarse.tolist():
|
|
q, cell = divmod(slot, nx * ny)
|
|
y, x = divmod(cell, nx)
|
|
sx, sy = x - D2Q9_EX[q], y - D2Q9_EY[q]
|
|
count += not (x0 <= sx < x1 and y0 <= sy < y1)
|
|
return count
|
|
|
|
def test_rectangular_shifted_patch_against_geometric_definition(self):
|
|
p = FixedPatch2D.from_config(LBMConfig(nx=22, ny=18, viscosity=.1), (7, 6, 13, 10))
|
|
expected = set()
|
|
for y in range(18):
|
|
for x in range(22):
|
|
if 7 <= x < 13 and 6 <= y < 10:
|
|
continue
|
|
for q, (cx, cy) in enumerate(zip(D2Q9_EX, D2Q9_EY)):
|
|
if q and 7 <= x + .5 - cx / 2 <= 13 and 6 <= y + .5 - cy / 2 <= 10:
|
|
expected.add(q * 22 * 18 + y * 22 + x)
|
|
self.assertEqual(set(p.coalesce_coarse.tolist()), expected)
|
|
self.assertEqual(p.fine_shape, (16, 20))
|
|
self.assertEqual(p.fine_origin, (5., 4.))
|
|
fy, fx = p.fine_shape
|
|
for coarse_slot, children in zip(p.coalesce_coarse.tolist(), p.coalesce_fine.tolist()):
|
|
q, cell = divmod(coarse_slot, 22 * 18)
|
|
y, x = divmod(cell, 22)
|
|
centers = []
|
|
for slot in children:
|
|
fq, fine_cell = divmod(slot, fy * fx)
|
|
b, a = divmod(fine_cell, fx)
|
|
self.assertEqual(fq, q)
|
|
self.assertEqual(p.fine_kind[b, a], INNER_GHOST)
|
|
centers.append((p.fine_origin[0] + (a + .5) / 2,
|
|
p.fine_origin[1] + (b + .5) / 2))
|
|
np.testing.assert_array_equal(np.mean(centers, axis=0), (x + .5, y + .5))
|
|
self.assertEqual(self.corner_pass_count(p), 8)
|
|
|
|
def test_smallest_patch_and_extreme_halo_margin(self):
|
|
p = FixedPatch2D.from_config(LBMConfig(nx=5, ny=5, viscosity=.1), (2, 2, 3, 3))
|
|
self.assertEqual(p.fine_shape, (10, 10))
|
|
self.assertEqual(p.fine_origin, (0., 0.))
|
|
self.assertEqual((int(p.ghost_coarse.min()), int(p.ghost_coarse.max())), (0, 24))
|
|
self.assertEqual(np.count_nonzero(p.fine_kind == ACTIVE), 4)
|
|
self.assertEqual(self.corner_pass_count(p), 8)
|
|
self.assertTrue(np.all(p.fine_kind.ravel()[p.coalesce_fine % 100] == INNER_GHOST))
|
|
|
|
def test_config_and_level_scaling_are_isolated(self):
|
|
cfg = config()
|
|
before = vars(cfg).copy()
|
|
p = FixedPatch2D.from_config(cfg, (4, 4, 12, 12))
|
|
self.assertEqual(vars(cfg), before)
|
|
self.assertAlmostEqual(p.tau_coarse, .8)
|
|
self.assertAlmostEqual(p.tau_fine, 1.1)
|
|
self.assertAlmostEqual((p.tau_coarse - .5) / 3,
|
|
(p.tau_fine - .5) / 3 * .5)
|
|
cfg.nx = 999
|
|
self.assertEqual(p.coarse_shape, (16, 16))
|
|
|
|
def test_tables_are_compact_contiguous_and_readonly(self):
|
|
small = FixedPatch2D.from_config(config(), (4, 4, 12, 12))
|
|
large = FixedPatch2D.from_config(LBMConfig(nx=512, ny=256, viscosity=.1), (4, 4, 12, 12))
|
|
self.assertEqual(small.table_bytes, large.table_bytes)
|
|
self.assertEqual(small.table_bytes, 7136)
|
|
for name, dtype in (("fine_kind", np.uint8), ("ghost_fine", np.uint32),
|
|
("ghost_coarse", np.uint32), ("coalesce_coarse", np.uint64),
|
|
("coalesce_fine", np.uint64)):
|
|
array = getattr(small, name)
|
|
self.assertEqual(array.dtype, dtype)
|
|
self.assertTrue(array.flags.c_contiguous)
|
|
with self.assertRaises(ValueError):
|
|
array.flat[0] = 0
|
|
|
|
def test_invalid_geometry_and_index_overflow_rejected_before_allocation(self):
|
|
for bounds in ((4, 4, 4, 12), (12, 4, 4, 12), (1, 4, 12, 12),
|
|
(4, 4, 15, 12), (4., 4, 12, 12), (True, 4, 12, 12),
|
|
(np.bool_(True), 4, 12, 12), (4, 4, 12)):
|
|
with self.subTest(bounds=bounds), self.assertRaises(ValueError):
|
|
FixedPatch2D.from_config(config(), bounds)
|
|
for nx, ny, bounds in ((2**32, 10, (2, 2, 3, 3)),
|
|
(50000, 50000, (2, 2, 49998, 49998))):
|
|
with self.assertRaisesRegex(ValueError, "uint32"):
|
|
FixedPatch2D.from_config(LBMConfig(nx=nx, ny=ny, viscosity=.1), bounds)
|
|
|
|
def test_unsupported_backend_and_fine_relaxation_rejected(self):
|
|
for kwargs in ({"collision": "TRT"}, {"collision": "MRT"}, {"streaming": "esopull"},
|
|
{"ddf_shifting": True}, {"les_enabled": True},
|
|
{"store_precision": "FP16S"}, {"lattice_model": "D3Q19"},
|
|
{"data_type": "FP64"}, {"omega_min": 1.0}):
|
|
cfg = config(**kwargs)
|
|
before = vars(cfg).copy()
|
|
with self.subTest(kwargs=kwargs), self.assertRaises(ValueError):
|
|
FixedPatch2D.from_config(cfg, (4, 4, 12, 12))
|
|
self.assertEqual(vars(cfg), before)
|
|
|
|
def test_handoff_rejects_wrong_phase_and_unsynchronized_counts(self):
|
|
p = FixedPatch2D.from_config(config(), (4, 4, 12, 12))
|
|
p.check_synchronization(0, 0, phase="postcollision_g")
|
|
p.check_synchronization(8, 16, phase="postcollision_g")
|
|
for c, f, phase in ((8, 16, "arrival_f"), (8, 15, "postcollision_g"),
|
|
(-1, -2, "postcollision_g"), (1., 2, "postcollision_g"),
|
|
(True, 2, "postcollision_g"), (1, np.bool_(True), "postcollision_g")):
|
|
with self.subTest(c=c, f=f, phase=phase), self.assertRaises(ValueError):
|
|
p.check_synchronization(c, f, phase=phase)
|
|
|
|
|
|
class FixedPatchSetTests(unittest.TestCase):
|
|
"""Several separated patches: level-wide re-addressing of FixedPatch2D tables."""
|
|
FACE = ((4, 4, 8, 12), (12, 4, 16, 12)) # 16x24 coarse, face gap 4
|
|
S1 = ((4, 4, 12, 12), (20, 4, 28, 12)) # 16x32 coarse, translation copy
|
|
|
|
def build(self, rects, ny=16, nx=24):
|
|
return refinement.FixedPatchSet2D.from_config(LBMConfig(nx=nx, ny=ny, viscosity=.1), rects)
|
|
|
|
def test_tables_are_shifted_single_patch_tables(self):
|
|
s = self.build(self.FACE)
|
|
fy, fx = s.fine_shape
|
|
self.assertEqual((fy, fx), (32, 48))
|
|
rows = cells = 0
|
|
for p in s.patches:
|
|
y0, x0 = p.fine_level_origin
|
|
ly, lx = p.fine_shape
|
|
np.testing.assert_array_equal(s.fine_kind[y0:y0 + ly, x0:x0 + lx], p.fine_kind)
|
|
gy, gx = np.divmod(p.ghost_fine.astype(np.int64), lx)
|
|
np.testing.assert_array_equal(s.ghost_fine[cells:cells + len(gy)],
|
|
(gy + y0) * fx + gx + x0)
|
|
np.testing.assert_array_equal(s.ghost_coarse[cells:cells + len(gy)], p.ghost_coarse)
|
|
q, c = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx)
|
|
cy, cx = np.divmod(c, lx)
|
|
np.testing.assert_array_equal(s.coalesce_fine[rows:rows + len(q)],
|
|
q * fy * fx + (cy + y0) * fx + cx + x0)
|
|
np.testing.assert_array_equal(s.coalesce_coarse[rows:rows + len(q)], p.coalesce_coarse)
|
|
rows, cells = rows + len(q), cells + len(gy)
|
|
self.assertEqual((rows, cells), (len(s.coalesce_coarse), len(s.ghost_fine)))
|
|
self.assertEqual(int((s.fine_kind == ACTIVE).sum()), 4 * int(s.covered.sum()))
|
|
|
|
def test_input_order_is_irrelevant(self):
|
|
a, b = self.build(self.FACE), self.build(self.FACE[::-1])
|
|
for name in ("fine_kind", "ghost_fine", "ghost_coarse", "coalesce_coarse", "coalesce_fine"):
|
|
np.testing.assert_array_equal(getattr(a, name), getattr(b, name))
|
|
self.assertEqual(a.bounds, b.bounds)
|
|
|
|
def test_one_patch_set_matches_fixed_patch(self):
|
|
p = FixedPatch2D.from_config(config(), (4, 4, 12, 12))
|
|
s = self.build(((4, 4, 12, 12),), nx=16)
|
|
np.testing.assert_array_equal(s.covered, p.covered)
|
|
self.assertEqual((s.tau_coarse, s.tau_fine), (p.tau_coarse, p.tau_fine))
|
|
np.testing.assert_array_equal(s.coalesce_coarse, p.coalesce_coarse)
|
|
np.testing.assert_array_equal(s.ghost_coarse, p.ghost_coarse)
|
|
|
|
def test_gap_and_overlap_rejected(self):
|
|
for gap in (0, 2, 3):
|
|
with self.subTest(gap=gap), self.assertRaisesRegex(ValueError, "closer"):
|
|
self.build(((4, 4, 8, 12), (8 + gap, 4, 12 + gap, 12)))
|
|
with self.assertRaises(ValueError):
|
|
self.build(())
|
|
with self.assertRaises(ValueError): # FixedPatch2D margin still applies
|
|
self.build(((1, 4, 8, 12), (12, 4, 16, 12)))
|
|
|
|
def test_sparse_fine_layout_stores_exactly_the_supports(self):
|
|
for rects, nx in ((self.FACE, 24), (self.S1, 32)):
|
|
s = self.build(rects, nx=nx)
|
|
layout = s.fine_layout((4, 4))
|
|
support = np.zeros(s.fine_shape, bool)
|
|
for p in s.patches:
|
|
y0, x0 = p.fine_level_origin
|
|
support[y0:y0 + p.fine_shape[0], x0:x0 + p.fine_shape[1]] = True
|
|
np.testing.assert_array_equal(layout.owned, support)
|
|
self.assertTrue(layout.owned.ravel()[s.ghost_fine.astype(np.intp)].all())
|
|
with self.assertRaisesRegex(ValueError, "aligned"):
|
|
self.build(self.FACE).fine_layout((5, 4))
|
|
for block in ((0, 4), (4, 0), (-4, 4), (4.0, 4), (True, 4), (4,), (4, 4, 4), (np.bool_(1), 4)):
|
|
with self.subTest(block=block), self.assertRaisesRegex(ValueError, "positive integers"):
|
|
self.build(self.FACE).fine_layout(block)
|
|
|
|
|
|
class PatchSetBoundaryTests(unittest.TestCase):
|
|
"""fine_objects/jets on a level-wide patch set (16x32, patches 4 cells apart)."""
|
|
PATCHES = ((2, 2, 14, 14), (18, 2, 30, 14))
|
|
|
|
@staticmethod
|
|
def body(x, y, obj_id=0):
|
|
from CelerisLab.body.geometry.circle import CircleGeometry
|
|
from CelerisLab.body.objects import SimObject
|
|
return SimObject(obj_id, geometry=CircleGeometry(x, y, 2.15), center=(x, y), radius=2.15)
|
|
|
|
def prepare(self, nx, rects, objects, jets=()):
|
|
cfg = LBMConfig(nx=nx, ny=16, viscosity=.1)
|
|
s = refinement.FixedPatchSet2D.from_config(cfg, rects)
|
|
return s, refinement.prepare_fine_boundary(cfg, s, objects, jets=jets)
|
|
|
|
def test_two_bodies_are_the_one_body_case_translated(self):
|
|
_, one = self.prepare(16, self.PATCHES[:1], [self.body(8., 8.)])
|
|
_, two = self.prepare(32, self.PATCHES, [self.body(8., 8., 7), self.body(24., 8., 7)])
|
|
f1, f2 = one.flags.reshape(32, 32), two.flags.reshape(32, 64)
|
|
np.testing.assert_array_equal(f2[:, :32], f1)
|
|
np.testing.assert_array_equal(f2[:, 32:], f1)
|
|
n = len(one.curved_host[0])
|
|
self.assertEqual(len(two.curved_host[0]), 2 * n)
|
|
body_id = two.curved_host[3]
|
|
self.assertEqual(sorted(set(body_id.tolist())), [0, 1]) # positions, not obj_id 7
|
|
# Translation must reproduce identical link geometry (prerequisite of the
|
|
# bitwise GPU tiling gate): dirs, q and fallback equal, receivers shifted.
|
|
fidx, dirs, q, fallback = (two.curved_host[i] for i in (0, 1, 2, 7))
|
|
for column in (dirs, q, fallback):
|
|
np.testing.assert_array_equal(column[:n], column[n:])
|
|
y, x = np.divmod(fidx.astype(np.int64), 64)
|
|
np.testing.assert_array_equal(y[:n], y[n:])
|
|
np.testing.assert_array_equal(x[:n] + 32, x[n:])
|
|
|
|
def test_body_outside_a_single_patch_rejected(self):
|
|
for x in (16., 13.): # in the gap / crossing the patch edge
|
|
with self.subTest(x=x), self.assertRaises(ValueError):
|
|
self.prepare(32, self.PATCHES, [self.body(x, 8.)])
|
|
|
|
def test_one_jet_over_two_bodies_rejected_even_at_zero_amplitude(self):
|
|
from CelerisLab.body.jets import Jet, Parallel, Tube
|
|
bodies = [self.body(8., 8.), self.body(24., 8.)]
|
|
span = Jet(Tube((5., 8.), 0., 22., .6), Parallel(0.), amplitude=0.)
|
|
with self.assertRaisesRegex(ValueError, "one jet"):
|
|
self.prepare(32, self.PATCHES, bodies, jets=[span])
|
|
own = [Jet(Tube((5., 8.), 0., 6., .6), Parallel(0.), amplitude=.002),
|
|
Jet(Tube((21., 8.), 0., 6., .6), Parallel(0.), amplitude=.002)]
|
|
self.prepare(32, self.PATCHES, bodies, jets=own) # one jet per body: accepted
|
|
|
|
|
|
class NestPatchTests(unittest.TestCase):
|
|
"""Three-level single chain: L2 bounds in L1 window cells (probe fixture)."""
|
|
def setUp(self):
|
|
self.cfg = config()
|
|
self.l1 = FixedPatch2D.from_config(self.cfg, (4, 4, 12, 12)) # 24x24 window
|
|
|
|
def test_probe_fixture_composes_kinds(self):
|
|
child, kind = refinement.nest_patch(self.cfg, self.l1, (8, 8, 16, 16))
|
|
self.assertEqual(child.fine_shape, (24, 24))
|
|
self.assertAlmostEqual(child.tau_coarse, self.l1.tau_fine)
|
|
base = self.l1.fine_kind
|
|
for role in (INNER_GHOST, OUTER_GHOST):
|
|
np.testing.assert_array_equal(kind[base == role], role)
|
|
self.assertTrue((kind[child.covered] == OUTER_GHOST).all())
|
|
self.assertEqual(int((kind == ACTIVE).sum()), 192)
|
|
self.assertEqual(int((kind == INNER_GHOST).sum()), 144)
|
|
self.assertEqual(int(child.covered.sum()), 64)
|
|
|
|
def test_invalid_nesting_rejected(self):
|
|
for bounds in ((2, 8, 16, 16), # covers inherited INNER
|
|
(4, 8, 16, 16), # explosion donors in INNER
|
|
(8, 8, 20, 16), # donors reach OUTER/INNER at the far side
|
|
(5, 8, 16, 16)): # one-cell outward: restriction footprint/donors
|
|
with self.subTest(bounds=bounds), self.assertRaises(ValueError):
|
|
refinement.nest_patch(self.cfg, self.l1, bounds)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
unittest.main(verbosity=2)
|