318 lines
15 KiB
Diff
318 lines
15 KiB
Diff
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)
|