From 999d3a0e9413f98de444c1c68c8d9f0937fdc80d Mon Sep 17 00:00:00 2001 From: TRON Date: Fri, 16 Oct 2026 10:07:34 -0700 Subject: [PATCH] Add synchronized CPU leaf population remap --- docs/synchronized_leaf_remap.md | 70 ++++++++ src/CelerisLab/lbm/synchronized_remap.py | 147 ++++++++++++++++ tests/unit/test_synchronized_remap.py | 205 +++++++++++++++++++++++ 3 files changed, 422 insertions(+) create mode 100644 docs/synchronized_leaf_remap.md create mode 100644 src/CelerisLab/lbm/synchronized_remap.py create mode 100644 tests/unit/test_synchronized_remap.py diff --git a/docs/synchronized_leaf_remap.md b/docs/synchronized_leaf_remap.md new file mode 100644 index 0000000..9d74d00 --- /dev/null +++ b/docs/synchronized_leaf_remap.md @@ -0,0 +1,70 @@ +# Synchronized leaf G remap (CPU only) + +This is a CPU state-transfer component, **not a runtime adaptive solver**. +It changes only a separated L2 rectangle set in a three-level, body-free, +periodic hierarchy. L0/L1 topology, relaxation parameters, direction order, +physical frame and completed time remain unchanged. Existing geometry checks +validate child support/nesting and separation. No CUDA, scheduler, allocator, +boundary handling, collision, ghost refill or runtime installation is added. + +```python +from CelerisLab.lbm.refinement import FixedPatchSet2D +from CelerisLab.lbm.synchronized_remap import RemapTopology, remap_synchronized + +# cfg: D2Q9/SRT/FP32/unshifted/double_buffer, nz=1, no LES +parent = FixedPatchSet2D.from_config(cfg, [(4, 4, 20, 28), (24, 4, 44, 28)]) +old = RemapTopology.from_config(cfg, parent, [(16, 16, 32, 40)]) +new = RemapTopology.from_config(cfg, parent, [(20, 16, 36, 40)]) +candidate = remap_synchronized( + old, new, current_g, base_count=7, level_steps=(7, 14, 28), + phase="postcollision_g", origin=(0., 0.), spacing=1., +) +``` + +The example has a 48×32 root; both child supports allow 4×4 blocks. +Child bounds use the parent's allocated L1 coordinates, not rectangle-local +indices. A scalar FixedPatch2D parent is also supported: its shifted L1 window +origin is retained. Physical lower edges are `origin + spacing*(x,y)` on L0, +`origin + spacing*(parent.fine_origin + (x,y)/2)` on L1 and +`origin + spacing*(parent.fine_origin + (x,y)/4)` on L2. + +Factory options `coarse_boundary` and `objects` must declare `"periodic"` and +an empty body collection (all levels). LBMConfig's channel settings do not +select this operator's boundary mode. Geometry metadata is factory-created; +manual mutation/forging of metadata is unsupported. Caller certifies the +actual state has no bodies/forcing, that all device work and old readers are +quiescent, and that inputs are the **current owned-leaf postcollision G**. +Count/phase stamps check consistency, not freshness or active-buffer identity. + +`current_g` is three NumPy Q-major arrays `(9,ny,nx)`. For blocked inputs, +pass `layouts=(base_layout, fine_layout, child_layout)`; each optional `None` +entry denotes dense input. Every leaf must have an owner. Only owner cores +are read, never redundant block halos. Mixed precision is rejected. FP64 is +accepted for CPU reference tests; this does not enable FP64 runtime storage. +Outputs `candidate.g` are separate dense arrays regardless of input layout. +Memory scales with public grid extents, not only sparse allocated blocks. + +Same-level leaf overlap is copied bitwise in physical coordinates. Newly +refined L2 copies all nine values from its valid old L1 parent. Newly exposed +L1 receives the complete four-child mean, including q0. Covered parents, +ghosts and holes never serve as donors; NaN/Inf outside leaves is ignored, +while nonfinite required leaves or nonpositive density reject the operation. +Restriction uses FP64 quarter-weight accumulation and rounds to the input +dtype. New nonleaves are NaN. Source arrays/metadata are not changed. + +`candidate` retains `base_count`, `level_steps`, `phase`, `origin`, `spacing` +and new `topology.leaf_masks`. Existing `FixedPatchStepper.initialize` resets +the counter: using it is **not** a continuation adapter. This component neither +publishes a candidate nor retires old state. Runtime table rebuild/halo refill +and correct first-step scheduling remain separate integration work. + +Tests accumulate leaf-area integrals in FP64 for each population, mass and +momentum. The declared absolute tolerance is `2e-13 * sum(abs(old population +integrals))` for FP64, `2e-7 * scale` for FP32, accounting for reduction and +restriction rounding (not a solver accuracy bound). Coarse→fine→coarse copy +round-trips preserve state; coarsening arbitrary fine structure loses it and +is not invertible. No kinetic accuracy or stability claim follows. + +```sh +PYTHONPATH=src python -m unittest discover -s tests/unit -p 'test_synchronized_remap.py' -v +``` diff --git a/src/CelerisLab/lbm/synchronized_remap.py b/src/CelerisLab/lbm/synchronized_remap.py new file mode 100644 index 0000000..91dbf4b --- /dev/null +++ b/src/CelerisLab/lbm/synchronized_remap.py @@ -0,0 +1,147 @@ +"""CPU-only whole-volume G transfer, not a runtime adaptive solver. + +Only body-free periodic three-level grids with a fixed L0/L1 topology are +supported. The caller certifies that current owned leaves are quiescent, +unshifted D2Q9 postcollision G at the supplied completed base interval. +No device buffers, boundary state, ghosts, collision or time advance occur. +""" +from dataclasses import dataclass +from numbers import Integral + +import numpy as np + +from .refinement import ACTIVE, FixedPatch2D, FixedPatchSet2D, nest_patch_set + + +@dataclass(frozen=True, eq=False) +class RemapTopology: + """Validated fixed parent and separated L2 set; use ``from_config``. + + Child rectangles are in the parent's allocated L1 coordinates, including + its window offset/halo when the parent is a scalar patch. Dense masks are + rebuilt from geometry, never inferred from finite values in input arrays. + """ + parent: object + child: FixedPatchSet2D + leaf_masks: tuple + + @classmethod + def from_config(cls, cfg, parent, child_rectangles, *, coarse_boundary="periodic", objects=()): + if coarse_boundary != "periodic" or tuple(objects) or cfg.nz != 1: + raise ValueError("remap supports only body-free periodic 2D state") + if type(parent) not in (FixedPatch2D, FixedPatchSet2D): + raise ValueError("parent must be a fixed patch or fixed patch set") + parent = type(parent).from_config(cfg, parent.bounds) + child, middle_kind = nest_patch_set(cfg, parent, child_rectangles) + masks = (~parent.covered, middle_kind == ACTIVE, child.fine_kind == ACTIVE) + for mask in masks: + mask.setflags(write=False) + return cls(parent, child, masks) + + +@dataclass(frozen=True, eq=False) +class RemapCandidate: + """Separate dense arrays plus explicit synchronization/frame metadata. + + Passing these arrays to the existing stepper.initialize resets its count; + this candidate is deliberately NOT an installed runtime continuation. + """ + topology: RemapTopology + g: tuple + base_count: int + level_steps: tuple + phase: str + origin: tuple + spacing: float + + +def _same_parent(a, b): + return (type(a) is type(b) and a.bounds == b.bounds + and a.coarse_shape == b.coarse_shape + and a.tau_coarse == b.tau_coarse and a.tau_fine == b.tau_fine) + + +def _read_leaves(topology, arrays, layouts): + """Unpack only owned cores; nonleaf values are not evidence of validity.""" + if len(arrays) != 3 or len(layouts) != 3: + raise ValueError("expected three G arrays and three layouts") + result = [] + for a, layout, leaf in zip(arrays, layouts, topology.leaf_masks): + if not isinstance(a, np.ndarray) or a.dtype not in (np.dtype('float32'), np.dtype('float64')): + raise ValueError("G arrays must be FP32 or FP64 NumPy arrays") + if layout is not None: + if layout.shape != leaf.shape or not layout.owned[leaf].all(): + raise ValueError("layout must own every required leaf donor") + a = layout.to_dense(a) + if a.shape != (9, *leaf.shape): + raise ValueError("G must have Q-major (9, ny, nx) shape") + values = a[:, leaf] + if not np.isfinite(values).all() or (values.sum(axis=0, dtype=np.float64) <= 0).any(): + raise ValueError("required leaf donors must be finite with positive density") + result.append(a) + if len({a.dtype for a in result}) != 1: + raise ValueError("all levels must use the same population dtype") + return result + + +def remap_synchronized(old, new, old_g, *, base_count, level_steps, phase, + origin=(0., 0.), spacing=1., layouts=None): + """Return a non-mutating candidate at the SAME physical time and frame. + + ``old``/``new`` are RemapTopology.from_config results. L0/L1 geometry and + relaxation parameters must match; only the L2 rectangle set may change. + ``old_g`` contains current owned-leaf G, dense or (with ``layouts``) blocked. + ``level_steps`` must equal (base_count, 2*base_count, 4*base_count), and + ``phase`` must be 'postcollision_g'. These stamps are caller assertions, + not proof of device quiescence or active-buffer freshness. + + Preserve same-level leaf overlap bitwise. Copy all nine parent G values + to newly refined children; average ALL four old children (including q0) + for newly exposed L1. Restriction accumulates in FP64, then rounds once to + the source dtype. No equilibrium reconstruction or nonequilibrium scaling. + All destination nonleaves are NaN. Inputs and their write flags are intact. + """ + if (not isinstance(old, RemapTopology) or not isinstance(new, RemapTopology) + or not _same_parent(old.parent, new.parent) + or old.child.tau_fine != new.child.tau_fine): + raise ValueError("remap requires unchanged L0/L1 geometry and relaxation") + if not isinstance(level_steps, (tuple, list)): + raise ValueError("level_steps must be three completed integer counts") + counts = (base_count, *level_steps) + if (len(counts) != 4 or any(isinstance(n, (bool, np.bool_)) + or not isinstance(n, Integral) or n < 0 for n in counts) + or tuple(level_steps) != (base_count, 2*base_count, 4*base_count) + or phase != 'postcollision_g'): + raise ValueError("requires synchronized completed counts and postcollision_g") + try: + frame = np.asarray(origin, dtype=np.float64) + spacing = float(spacing) + except (TypeError, ValueError): + raise ValueError("frame needs a finite (x,y) origin and positive spacing") from None + if (frame.shape != (2,) or not np.isfinite(frame).all() + or not np.isscalar(spacing) or not np.isfinite(spacing) or spacing <= 0): + raise ValueError("frame needs a finite (x,y) origin and positive spacing") + source = _read_leaves(old, old_g, (None,)*3 if layouts is None else layouts) + candidate = [np.full((9, *mask.shape), np.nan, dtype=source[0].dtype) + for mask in new.leaf_masks] + for level, (before, after) in enumerate(zip(old.leaf_masks, new.leaf_masks)): + overlap = before & after + candidate[level][:, overlap] = source[level][:, overlap] + + # A child set uses a level-wide grid in its parent's allocated coordinates. + # Thus child (y,x) belongs to parent (y//2,x//2), even for a shifted L1 window. + y, x = np.nonzero(new.leaf_masks[2] & ~old.leaf_masks[2]) + if not old.leaf_masks[1][y//2, x//2].all(): + raise ValueError("new L2 requires an old owned L1 leaf") + candidate[2][:, y, x] = source[1][:, y//2, x//2] + + y, x = np.nonzero(new.leaf_masks[1] & ~old.leaf_masks[1]) + children = [(2*y+dy, 2*x+dx) for dy in (0, 1) for dx in (0, 1)] + if any(not old.leaf_masks[2][cy, cx].all() for cy, cx in children): + raise ValueError("reactivated L1 requires all four old owned L2 leaves") + # Scale first to avoid an overflowing intermediate sum for finite FP64 G. + mean = sum(source[2][:, cy, cx].astype(np.float64) * .25 for cy, cx in children) + candidate[1][:, y, x] = mean + _read_leaves(new, candidate, (None,)*3) + return RemapCandidate(new, tuple(candidate), int(base_count), tuple(map(int, level_steps)), + phase, tuple(map(float, frame)), float(spacing)) diff --git a/tests/unit/test_synchronized_remap.py b/tests/unit/test_synchronized_remap.py new file mode 100644 index 0000000..0554a66 --- /dev/null +++ b/tests/unit/test_synchronized_remap.py @@ -0,0 +1,205 @@ +"""CPU remap tests: full-population conservation, never a runtime AMR claim.""" +from dataclasses import replace +import unittest + +import numpy as np + +from CelerisLab.config import LBMConfig +from CelerisLab.lbm.blocks import BlockLayout2D +from CelerisLab.lbm.descriptors import D2Q9_EX, D2Q9_EY +from CelerisLab.lbm.refinement import FixedPatch2D, FixedPatchSet2D +from CelerisLab.lbm.synchronized_remap import RemapTopology, remap_synchronized + + +CFG = LBMConfig(nx=48, ny=32, viscosity=.1) +PARENTS = [(4, 4, 20, 28), (24, 4, 44, 28)] +A, B = [(16, 16, 32, 40)], [(20, 16, 36, 40)] + + +def topology(rectangles, parent=None, cfg=CFG): + if parent is None: + parent = FixedPatchSet2D.from_config(cfg, PARENTS) + return RemapTopology.from_config(cfg, parent, rectangles) + + +def fields(mesh, dtype=np.float64): + rng = np.random.default_rng(8473) + result = [] + for leaf in mesh.leaf_masks: + # Direction-dependent and non-equilibrium, with independently varying + # children: q0 and momentum cannot be recovered from equilibrium rho/u. + g = rng.uniform(.03, .15, (9, *leaf.shape)).astype(dtype) + g += np.arange(9, dtype=dtype)[:, None, None] * .01 + g[:, ~leaf] = np.nan + result.append(g) + return result + + +def remap(old, new, g, **kwargs): + options = dict(base_count=7, level_steps=(7, 14, 28), phase='postcollision_g', + origin=(-3., 2.), spacing=1.25) + options.update(kwargs) + return remap_synchronized(old, new, g, **options) + + +def integrals(mesh, g, spacing=1.25): + populations = sum(a[:, mask].sum(axis=1, dtype=np.float64) * spacing**2 / 4**level + for level, (a, mask) in enumerate(zip(g, mesh.leaf_masks))) + return np.r_[populations, populations.sum(), populations @ D2Q9_EX, + populations @ D2Q9_EY] + + +class SynchronizedRemapTests(unittest.TestCase): + def setUp(self): + self.old, self.new = topology(A), topology(B) + self.g = fields(self.old) + + def test_identity_bitwise_and_no_mutation_or_count_advance(self): + copies = [a.copy() for a in self.g] + flags = [a.flags.writeable for a in self.g] + out = remap(self.old, self.old, self.g) + self.assertEqual((out.base_count, out.level_steps, out.phase, out.origin, out.spacing), + (7, (7, 14, 28), 'postcollision_g', (-3., 2.), 1.25)) + for a, b, saved, leaf, flag in zip(self.g, out.g, copies, self.old.leaf_masks, flags): + self.assertEqual(a[:, leaf].tobytes(), b[:, leaf].tobytes()) + np.testing.assert_array_equal(a, saved) + self.assertEqual(a.flags.writeable, flag) + self.assertFalse(np.shares_memory(a, b)) + self.assertTrue(np.isnan(b[:, ~leaf]).all()) + + def test_shift_preserves_physical_overlap_not_local_patch_index(self): + out = remap(self.old, self.new, self.g) + for level in range(3): + overlap = self.old.leaf_masks[level] & self.new.leaf_masks[level] + self.assertEqual(self.g[level][:, overlap].tobytes(), out.g[level][:, overlap].tobytes()) + y, x = np.nonzero(self.new.leaf_masks[2] & ~self.old.leaf_masks[2]) + np.testing.assert_array_equal(out.g[2][:, y, x], self.g[1][:, y//2, x//2]) + y, x = np.nonzero(self.new.leaf_masks[1] & ~self.old.leaf_masks[1]) + expected = sum(self.g[2][:, 2*y+dy, 2*x+dx] * .25 + for dy in (0, 1) for dx in (0, 1)) + np.testing.assert_array_equal(out.g[1][:, y, x], expected) + self.assertGreater(np.max(np.abs(expected[0] - self.g[2][0, 2*y, 2*x])), .001) + + def test_fp64_and_fp32_area_population_mass_momentum_conservation(self): + for dtype, tolerance in ((np.float64, 2e-13), (np.float32, 2e-7)): + g = fields(self.old, dtype) + out = remap(self.old, self.new, g) + before, after = integrals(self.old, g), integrals(self.new, out.g) + # FP64 accumulation; absolute scale handles near-cancelled momentum. + scale = np.sum(np.abs(before[:9])) + np.testing.assert_allclose(after, before, rtol=0, atol=tolerance * scale) + for array, leaf in zip(out.g, self.new.leaf_masks): + self.assertEqual(array.dtype, dtype) + self.assertTrue(np.isfinite(array[:, leaf]).all()) + self.assertTrue(np.isnan(array[:, ~leaf]).all()) + + def test_layout_and_rectangle_permutation_independence(self): + old = topology(A + [(56, 16, 80, 40)]) + new = topology(B + [(56, 16, 80, 40)]) + reordered = topology(list(reversed(B + [(56, 16, 80, 40)])), + FixedPatchSet2D.from_config(CFG, list(reversed(PARENTS)))) + g = fields(old) + expected = remap(old, new, g) + for block in ((4, 4), (2, 2)): + layouts = (BlockLayout2D.build(CFG.ny, CFG.nx, *block), + old.parent.fine_layout(block), old.child.fine_layout(block)) + blocked = [layout.to_blocked(a) for layout, a in zip(layouts, g)] + # Corrupt all block halos; only owner-core populations are admissible. + for layout, a in zip(layouts, blocked): + h = layout.halo + a[:, :, :h, :] = np.inf + a[:, :, -h:, :] = np.inf + a[:, :, :, :h] = np.inf + a[:, :, :, -h:] = np.inf + saved = [a.copy() for a in blocked] + actual = remap(old, reordered, blocked, layouts=layouts) + for a, b in zip(expected.g, actual.g): + np.testing.assert_array_equal(a, b) + for a, b in zip(saved, blocked): + np.testing.assert_array_equal(a, b) + + def test_shifted_parent_physical_frame(self): + # L1 window origin=(2,2) in root units; child coordinates remain parent-local. + cfg = LBMConfig(nx=32, ny=32, viscosity=.1) + parent = FixedPatch2D.from_config(cfg, (4, 4, 28, 28)) + old, new = topology(A, parent, cfg), topology(B, parent, cfg) + g = fields(old) + out = remap(old, new, g, origin=(10., -7.)) + self.assertEqual(out.topology.parent.fine_origin, (2., 2.)) + overlap = old.leaf_masks[2] & new.leaf_masks[2] + np.testing.assert_array_equal(out.g[2][:, overlap], g[2][:, overlap]) + np.testing.assert_allclose(integrals(old, g), integrals(new, out.g), rtol=2e-13, atol=2e-12) + + def test_outside_leaf_nan_inf_ignored_required_leaf_rejected(self): + for level in range(3): + dirty = [a.copy() for a in self.g] + dirty[level][:, ~self.old.leaf_masks[level]] = np.inf + actual = remap(self.old, self.new, dirty) + expected = remap(self.old, self.new, self.g) + for a, b in zip(expected.g, actual.g): + np.testing.assert_array_equal(a, b) + # Specifically test restriction children, prolongation donors, and overlap. + required = [(0, self.old.leaf_masks[0]), + (1, self.old.leaf_masks[1] & ~self.new.leaf_masks[1]), + (2, self.old.leaf_masks[2] & ~self.new.leaf_masks[2])] + for level, mask in required: + y, x = np.argwhere(mask)[0] + for bad in (np.nan, np.inf, -np.inf): + dirty = [a.copy() for a in self.g] + dirty[level][0, y, x] = bad + with self.assertRaisesRegex(ValueError, 'required leaf donors'): + remap(self.old, self.new, dirty) + dirty = [a.copy() for a in self.g] + dirty[0][:, 0, 0] = -1 + with self.assertRaisesRegex(ValueError, 'positive density'): + remap(self.old, self.new, dirty) + + def test_absent_required_block_rejected(self): + layouts = (None, None, BlockLayout2D.build(*self.old.leaf_masks[2].shape, 4, 4, + periodic=(False, False), origins=[(0, 0)])) + arrays = [*self.g[:2], layouts[2].to_blocked(self.g[2])] + with self.assertRaisesRegex(ValueError, 'own every required leaf'): + remap(self.old, self.new, arrays, layouts=layouts) + + def test_coarse_fine_coarse_identity_not_general_inverse(self): + # Grow then shrink: newly refined cells copy one parent value, so means + # recover it exactly. Shrink then grow loses arbitrary child variation. + wide = topology([(16, 16, 36, 40)]) + refined = remap(self.old, wide, self.g) + restored = remap(wide, self.old, refined.g) + for a, b in zip(self.g, restored.g): + np.testing.assert_array_equal(a, b) + shifted = remap(self.old, self.new, self.g) + returned = remap(self.new, self.old, shifted.g) + lost = self.old.leaf_masks[2] & ~self.new.leaf_masks[2] + self.assertFalse(np.array_equal(self.g[2][:, lost], returned.g[2][:, lost])) + + def test_invalid_geometry_counts_representation_frame(self): + for rectangles in ([], [(2, 2, 10, 10)], [(16, 16, 32, 40)] * 2, + [(16, 16, 33.5, 40)], [(16, 16, 32, 40), (33, 16, 36, 40)]): + with self.assertRaises(ValueError): + topology(rectangles) + parent = self.old.parent + for changes in (dict(lattice_model='D3Q19'), dict(ddf_shifting=True), + dict(collision='TRT'), dict(les_enabled=True), dict(nz=2)): + with self.assertRaises(ValueError): + RemapTopology.from_config(replace(CFG, **changes), parent, A) + for kwargs in (dict(coarse_boundary='slot_reservoir'), dict(objects=(object(),))): + with self.assertRaisesRegex(ValueError, 'body-free periodic'): + RemapTopology.from_config(CFG, parent, A, **kwargs) + different = topology(A, cfg=replace(CFG, viscosity=.12)) + with self.assertRaisesRegex(ValueError, 'unchanged'): + remap(self.old, different, self.g) + for kwargs in (dict(base_count=True), dict(base_count=-1), dict(level_steps=(7, 14, 27)), + dict(level_steps=(7, 14)), dict(level_steps=None), dict(phase='arrival_f'), + dict(origin=(np.nan, 0)), dict(spacing=0), dict(spacing='bad')): + with self.assertRaises(ValueError): + remap(self.old, self.new, self.g, **kwargs) + for g in (self.g[:2], [self.g[0].astype(np.float32), *self.g[1:]], + [self.g[0][0], *self.g[1:]]): + with self.assertRaises(ValueError): + remap(self.old, self.new, g) + + +if __name__ == '__main__': + unittest.main()