diff --git a/src/CelerisLab/lbm/refinement.py b/src/CelerisLab/lbm/refinement.py new file mode 100644 index 0000000..0a6f563 --- /dev/null +++ b/src/CelerisLab/lbm/refinement.py @@ -0,0 +1,138 @@ +"""CPU preparation for one fixed D2Q9 2:1 patch; no coupled stepper yet. + +Cell indices use y*nx+x; population slots use q*(nx*ny)+cell, matching +core/layout.cuh. Bounds are half-open (x0, y0, x1, y1) in coarse cells. +The four fine halo cells implement the pinned CPUREF volumetric geometry. +Topology does not encode physical boundary conditions or obstacle safety. +""" +from dataclasses import dataclass, replace +from numbers import Integral + +import numpy as np + +from ..config import LBMConfig +from .descriptors import D2Q9_EX, D2Q9_EY + +ACTIVE, INNER_GHOST, OUTER_GHOST = 0, 1, 2 +HALO = 4 + + +def _readonly(values, dtype): + array = np.ascontiguousarray(values, dtype=dtype) + array.setflags(write=False) + return array + + +@dataclass(frozen=True, eq=False) +class FixedPatch2D: + """Layer metadata and compact exchange tables, constructed by from_config. + + Shapes are (ny, nx), including the fine halo. ``fine_origin`` is its lower + cell edge in coarse physical units, with coarse dx=dt=1 and fine dx=dt=.5. + ``ghost_fine`` and ``ghost_coarse`` pair fine cells with coarse donors. + ``coalesce_coarse`` contains incoming population slots; each corresponding + row of ``coalesce_fine`` contains four fine slots in (dy, dx) row order. + Weights are all 1/4. Closed corner links are included even when the coarse + donor is uncovered. These tables contain no DDF data and allocate no GPU. + + This is preparation, not permission to advance these levels with OneStep: + its channel boundaries and fused collision do not implement patch coupling. + """ + coarse_shape: tuple[int, int] + bounds: tuple[int, int, int, int] + tau_coarse: float + tau_fine: float + fine_kind: np.ndarray + ghost_fine: np.ndarray + ghost_coarse: np.ndarray + coalesce_coarse: np.ndarray + coalesce_fine: np.ndarray + + @classmethod + def from_config(cls, cfg: LBMConfig, bounds) -> "FixedPatch2D": + """Prepare maps without changing cfg or allocating coarse-sized tables. + + The initial backend contract is D2Q9/SRT/FP32/unshifted/double_buffer, + without LES. Two coarse cells of margin keep all halo donors in bounds; + this does NOT certify that those donors are fluid or periodic. + """ + checked = replace(cfg) + checked.validate() + if (not checked.is_d2q9 or checked.collision != "SRT" + or checked.streaming != "double_buffer" + or checked.store_precision != "FP32" + or checked.ddf_shifting or checked.les_enabled): + raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/unshifted/double_buffer without LES") + bounds = tuple(bounds) + if len(bounds) != 4 or any(isinstance(v, (bool, np.bool_)) + or not isinstance(v, Integral) for v in bounds): + raise ValueError("bounds must be four integers: (x0, y0, x1, y1)") + x0, y0, x1, y1 = bounds = tuple(map(int, bounds)) + nx, ny = checked.nx, checked.ny + if not (2 <= x0 < x1 <= nx - 2 and 2 <= y0 < y1 <= ny - 2): + raise ValueError("patch must be nonempty with a two-coarse-cell margin") + fx, fy = 2 * (x1 - x0) + 2 * HALO, 2 * (y1 - y0) + 2 * HALO + if max(nx * ny, fx * fy) > np.iinfo(np.uint32).max: + raise ValueError("layer cell counts exceed uint32 indexing") + tau_c = 3.0 * checked.viscosity + 0.5 + tau_f = 6.0 * checked.viscosity + 0.5 + if not checked.omega_min <= 1.0 / tau_f <= checked.omega_max: + raise ValueError("fine relaxation rate lies outside configured omega guards") + + kind = np.full((fy, fx), OUTER_GHOST, dtype=np.uint8) + kind[2:-2, 2:-2] = INNER_GHOST + kind[HALO:-HALO, HALO:-HALO] = ACTIVE + gy, gx = np.nonzero(kind != ACTIVE) + ghost_fine = gy * fx + gx + ghost_coarse = (y0 + (gy - HALO) // 2) * nx + x0 + (gx - HALO) // 2 + + targets, sources = [], [] + # Only the outside ring can receive fine-to-coarse populations. + # Integer doubled midpoints include all eight diagonal corner passes. + for y in range(y0 - 1, y1 + 1): + xs = range(x0 - 1, x1 + 1) if y in (y0 - 1, y1) else (x0 - 1, x1) + for x in xs: + for q in range(1, 9): + mx, my = 2 * x + 1 - D2Q9_EX[q], 2 * y + 1 - D2Q9_EY[q] + if 2 * x0 <= mx <= 2 * x1 and 2 * y0 <= my <= 2 * y1: + targets.append(q * nx * ny + y * nx + x) + a, b = HALO + 2 * (x - x0), HALO + 2 * (y - y0) + sources.append([q * fx * fy + (b + dy) * fx + a + dx + for dy in (0, 1) for dx in (0, 1)]) + return cls((ny, nx), bounds, tau_c, tau_f, + _readonly(kind, np.uint8), + _readonly(ghost_fine, np.uint32), + _readonly(ghost_coarse, np.uint32), + _readonly(targets, np.uint64), + _readonly(sources, np.uint64).reshape(-1, 4)) + + @property + def fine_shape(self) -> tuple[int, int]: + return self.fine_kind.shape + + @property + def fine_origin(self) -> tuple[float, float]: + return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2 + + @property + def table_bytes(self) -> int: + return sum(a.nbytes for a in (self.fine_kind, self.ghost_fine, + self.ghost_coarse, self.coalesce_coarse, + self.coalesce_fine)) + + def check_synchronization(self, coarse_steps: int, fine_steps: int, *, phase: str) -> None: + """Check a proposed backend handoff; does not convert or collide data. + + Production storage is g=C(f), while CPUREF returns arrival_f. A future + g-to-g coarse step must collide coarse cells once and fine cells twice; + fine halo ghosts never collide. Coalescence must precede coarse collision. + Converting reference f to g for comparison is an observation operation, + not an additional production collision. Counts here are completed steps. + """ + if phase != "postcollision_g": + raise ValueError("patch backend handoff requires postcollision_g") + for value in (coarse_steps, fine_steps): + if isinstance(value, (bool, np.bool_)) or not isinstance(value, Integral) or value < 0: + raise ValueError("step counts must be nonnegative integers") + if fine_steps != 2 * coarse_steps: + raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps") diff --git a/tests/unit/test_refinement.py b/tests/unit/test_refinement.py new file mode 100644 index 0000000..bfe8e74 --- /dev/null +++ b/tests/unit/test_refinement.py @@ -0,0 +1,167 @@ +"""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) + + +if __name__ == "__main__": + unittest.main(verbosity=2)