"""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"}, {"streaming": "esopull"}, {"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) ghosts = dict(zip(map(int, s.ghost_fine), map(int, s.ghost_coarse))) for f, c in zip((gy + y0) * fx + gx + x0, p.ghost_coarse): self.assertEqual(ghosts[int(f)], int(c)) q, c = np.divmod(p.coalesce_fine.astype(np.int64), ly * lx) cy, cx = np.divmod(c, lx) rows_by_target = dict(zip(map(int, s.coalesce_coarse), s.coalesce_fine)) for target, children in zip(p.coalesce_coarse, q * fy * fx + (cy + y0) * fx + cx + x0): np.testing.assert_array_equal(rows_by_target[int(target)], children) 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_use_dependency_checks(self): for gap in (-2, 0, 1, 2, 3): with self.subTest(gap=gap): pair = self.build(((4, 4, 8, 12), (8 + gap, 4, 12 + gap, 12))) refinement._transfer_footprint(pair) refinement._check_parent_pull_closure(pair) 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, "divide"): 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) class ChildBoundaryTests(unittest.TestCase): @staticmethod def body(center=(8., 8.), radius=.8, obj_id=99, spec=None): from CelerisLab.body.geometry.circle import CircleGeometry from CelerisLab.body.objects import SimObject return SimObject(obj_id, geometry=CircleGeometry(*center, radius), center=center, radius=radius, surface_actuation=spec) def prepare(self, objects, jets=(), large=False): from dataclasses import replace cfg = LBMConfig(nx=32 if large else 16, ny=32 if large else 16, viscosity=.1) parent = FixedPatch2D.from_config(cfg, (4, 4, 28, 28) if large else (4, 4, 12, 12)) child, _ = refinement.nest_patch(cfg, parent, (8, 8, 40, 40) if large else (8, 8, 16, 16)) middle = replace(cfg, nx=parent.fine_shape[1], ny=parent.fine_shape[0], viscosity=.2) return child, refinement.prepare_fine_boundary(middle, child, objects, jets=jets, parent_origin=parent.fine_origin, parent_spacing=.5) def test_root_frame_matches_independent_parent_frame_projection(self): from CelerisLab.body.jets import Jet, Parallel, Tube, jet_link_shapes jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002) original = self.body() child, b = self.prepare([original], [jet]) middle = LBMConfig(nx=24, ny=24, viscosity=.2) transformed = self.body((12., 12.), 1.6) local_jet = Jet(Tube((9., 12.), 0., 6., .6), Parallel(0.), amplitude=.002) other = refinement.prepare_fine_boundary(middle, child, [transformed], jets=[local_jet]) for a, c in zip(b.curved_host, other.curved_host): np.testing.assert_array_equal(a, c) np.testing.assert_array_equal(b.flags, other.flags) np.testing.assert_allclose(b.curved_hits, 2 + .5 * other.curved_hits, atol=0, rtol=0) for a, c in zip(jet_link_shapes([jet], b.curved_hits), jet_link_shapes([local_jet], other.curved_hits)): np.testing.assert_allclose(a, c, atol=1e-14, rtol=0) self.assertEqual(len(b.curved_host[0]), 60) self.assertAlmostEqual(child.tau_fine, .5 + 12 * .1) self.assertEqual(original.geometry.center, (8., 8.)) self.assertEqual(original.obj_id, 99) def test_invalid_body_and_jet_contracts(self): from CelerisLab.body.jets import Jet, Parallel, Tube from CelerisLab.body.objects import SurfaceActuation jet = Jet(Tube((6.5, 8.), 0., 3., .3), Parallel(0.), amplitude=.002) for body in (self.body((6., 8.)), self.body((6.9, 8.)), self.body(radius=2.1)): with self.assertRaises(ValueError): self.prepare([body]) moving = self.body() moving.state.vx = .01 with self.assertRaisesRegex(ValueError, "stationary"): self.prepare([moving]) with self.assertRaises(ValueError): self.prepare([self.body()], [jet, jet]) empty = Jet(Tube((0., 0.), 0., 1., .1), Parallel(0.), amplitude=0.) with self.assertRaises(ValueError): self.prepare([self.body()], [empty]) spec = SurfaceActuation(.7, .85, .002) with self.assertRaisesRegex(ValueError, "SurfaceActuation"): self.prepare([self.body(spec=spec)], [spec.as_jet((8., 8.), .8)]) def test_one_jet_two_real_bodies_same_caller_id_rejected(self): from CelerisLab.body.jets import Jet, Parallel, Tube bodies = [self.body((12., 14.), 1., 7), self.body((17., 14.), 1., 7)] _, boundary = self.prepare(bodies, large=True) # establish legal geometry first self.assertEqual(np.unique(boundary.curved_host[3]).tolist(), [0, 1]) for amplitude in (0., .002): jet = Jet(Tube((10., 14.), 0., 9., .4), Parallel(0.), amplitude=amplitude) with self.assertRaisesRegex(ValueError, "one jet may drive only one body"): self.prepare(bodies, [jet], large=True) if __name__ == "__main__": unittest.main(verbosity=2)