Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add host-only arbitrary-depth refinement hierarchy builder

Browse Source 
lbm/hierarchy.py loops the existing pair logic (FixedPatch2D/FixedPatchSet2D
for L0->L1, nest_patch/nest_patch_set for deeper pairs, _middle_config for
the per-level config chain) and returns validated per-level configs, pair
tables and composed kinds for any depth. Legacy scalar window forms are
accepted so existing public shapes are unchanged. Validates each config copy
because LBMConfig.omega is only computed by validate(). Not yet used by the
stepper.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
This commit is contained in:
Claude
2026-10-17 01:19:29 +00:00
co-authored by Claude Opus 5.5
parent 77c034a29a
commit eeb90ace9a
2 changed files with 191 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/CelerisLab/lbm/hierarchy.py
+97
View File
@@ -0,0 +1,97 @@
# CelerisLab/lbm/hierarchy.py
"""Host-only geometry of a fixed 2:1 refinement hierarchy of any depth (numpy, no CUDA).
Loops the existing adjacent-pair logic of refinement.py; it adds no transfer rule.
``refinements[l]`` describes level l+1 in level-l integer cells: a scalar
(x0, y0, x1, y1) is one FixedPatch2D window (its level keeps the window shape), a
sequence of rectangles is a FixedPatchSet2D on the level-wide grid (origin 0).
A scalar child is only accepted under a scalar parent. Every pair after the first
is built by nest_patch/nest_patch_set, whose _nested_kind validates coverage,
donors/receivers, surviving pulls and the parent's upward restriction footprint.
"""
from dataclasses import dataclass, replace
from numbers import Integral
import numpy as np
from .refinement import (ACTIVE, OUTER_GHOST, FixedPatch2D, FixedPatchSet2D,
_middle_config, nest_patch, nest_patch_set)
@dataclass(frozen=True, eq=False)
class Hierarchy:
"""configs[l]: validated LBMConfig of level l; pairs[l]: tables of (l, l+1);
kinds[l]: composed per-cell kind (child coverage OUTER on every non-finest level)."""
configs: tuple
pairs: tuple
kinds: tuple
@property
def depth(self):
return len(self.kinds)
def is_rectangle(value):
value = tuple(value)
return len(value) == 4 and all(isinstance(v, (Integral, np.integer))
and not isinstance(v, (bool, np.bool_)) for v in value)
def build_hierarchy(cfg, refinements, *, wall_aware=False):
"""Geometry for len(refinements) + 1 levels; wall_aware marks base rows 0 and -1 OUTER."""
configs, pairs, kinds = [_validated(cfg)], [], []
for level, entry in enumerate(refinements):
scalar = is_rectangle(entry)
if not scalar and not tuple(entry):
raise ValueError("every refinement level needs at least one rectangle")
if level == 0:
pair = (FixedPatch2D if scalar else FixedPatchSet2D).from_config(configs[0], entry)
kind = np.where(pair.covered, OUTER_GHOST, ACTIVE).astype(np.uint8)
if wall_aware:
kind[[0, -1]] = OUTER_GHOST
elif scalar:
if not isinstance(pairs[-1], FixedPatch2D):
raise ValueError("a scalar child requires a scalar parent; use a rectangle list")
pair, kind = nest_patch(configs[level - 1], pairs[-1], entry)
else:
pair, kind = nest_patch_set(configs[level - 1], pairs[-1], entry)
kinds.append(np.asarray(kind))
pairs.append(pair)
configs.append(_validated(_middle_config(configs[level], pair)))
_check_chain(configs[level], configs[level + 1], pair)
if pairs:
kinds.append(pairs[-1].fine_kind)
else:
kinds.append(np.full((cfg.ny, cfg.nx), ACTIVE, np.uint8))
return Hierarchy(tuple(configs), tuple(pairs), tuple(kinds))
def _validated(cfg):
"""LBMConfig.omega is computed only by validate(); replace() copies it stale."""
cfg = replace(cfg)
cfg.validate()
return cfg
def _check_chain(coarse, fine, pair):
if not (np.isclose(pair.tau_coarse, 1 / coarse.omega, rtol=1e-12, atol=0)
and np.isclose(pair.tau_fine, 1 / fine.omega, rtol=1e-12, atol=0)
and tuple(pair.fine_shape) == (fine.ny, fine.nx)):
raise RuntimeError("level config chain is inconsistent with the pair tables")
def leaf_partition(h):
"""Base-cell fraction published as ACTIVE over all levels (exact: powers of two).
Window (scalar) levels are offset by their fine_level_origin first.
"""
base = np.zeros(h.kinds[0].shape)
oy = ox = 0
for level, kind in enumerate(h.kinds):
if level:
pair = h.pairs[level - 1]
py, px = (pair.fine_level_origin if isinstance(pair, FixedPatch2D) else (0, 0))
oy, ox = 2 * oy + py, 2 * ox + px
y, x = np.nonzero(kind == ACTIVE)
np.add.at(base, ((y + oy) >> level, (x + ox) >> level), 0.25 ** level)
return base
tests/unit/test_hierarchy.py
+94
View File
@@ -0,0 +1,94 @@
"""CPU geometry tests for lbm.hierarchy (no CUDA). Geometry construction only, not LBM numerics."""
import unittest
import numpy as np
from CelerisLab.config import LBMConfig
from CelerisLab.lbm import refinement as r
from CelerisLab.lbm.hierarchy import build_hierarchy, leaf_partition
CFG = LBMConfig(nx=48, ny=48, viscosity=.02)
L1 = [(6, 6, 26, 26), (32, 30, 42, 42)] # L0 cells; L1 ACTIVE x/y 12..52 and x 64..84, y 60..84
L2 = [(20, 20, 44, 44), (72, 68, 76, 76)] # L1 cells; L2 ACTIVE 40..88 and x 144..152, y 136..152
L3 = [(56, 56, 72, 72)] # L2 cells
FOUR = [L1, L2, L3]
SMALL = LBMConfig(nx=16, ny=16, viscosity=.1)
class FourLevels(unittest.TestCase):
@classmethod
def setUpClass(cls):
cls.h = build_hierarchy(CFG, FOUR)
def test_config_chain(self):
self.assertEqual(self.h.depth, 4)
for level, c in enumerate(self.h.configs):
self.assertEqual((c.ny, c.nx), (48 << level, 48 << level))
self.assertAlmostEqual(c.viscosity, .02 * 2 ** level, places=12)
self.assertEqual(self.h.kinds[level].shape, (c.ny, c.nx))
for level, p in enumerate(self.h.pairs):
self.assertAlmostEqual(p.tau_coarse, 1 / self.h.configs[level].omega, places=12)
self.assertAlmostEqual(p.tau_fine, 1 / self.h.configs[level + 1].omega, places=12)
def test_partition_and_composition(self):
h = self.h
np.testing.assert_array_equal(leaf_partition(h), np.ones((48, 48)))
for level in (1, 2):
covered, raw = h.pairs[level].covered, h.pairs[level - 1].fine_kind
self.assertTrue((h.kinds[level][covered] == r.OUTER_GHOST).all())
np.testing.assert_array_equal(h.kinds[level][~covered], raw[~covered])
np.testing.assert_array_equal(h.kinds[3], h.pairs[2].fine_kind)
def test_depth_is_list_length(self):
for depth in range(1, 5):
h = build_hierarchy(CFG, FOUR[:depth - 1])
self.assertEqual(h.depth, depth)
np.testing.assert_array_equal(leaf_partition(h), np.ones((48, 48)))
for level in range(depth - 1):
np.testing.assert_array_equal(h.kinds[level], self.h.kinds[level])
def test_illegal_coverage_reports_first_guard(self):
cases = [([L1, [(10, 20, 30, 44)]], "inherited ACTIVE"),
([L1, L2[:1] + [(54, 54, 58, 58)]], "inherited ACTIVE"),
([L1, L2, [(40, 40, 60, 60)]], "donors/receivers"),
([L1, L2, [(41, 41, 60, 60)]], "donors/receivers"),
([L1, L2, [(56, 56, 64, 64), (66, 56, 72, 64)]], "closer than")]
for levels, reason in cases:
with self.subTest(levels=levels), self.assertRaisesRegex(ValueError, reason):
build_hierarchy(CFG, levels)
build_hierarchy(CFG, [L1, L2, [(42, 42, 60, 60)]]) # minimal legal margin
class LegacyForms(unittest.TestCase):
def test_scalar_parent_and_child_match_existing_calls(self):
local = r.FixedPatch2D.from_config(SMALL, (4, 4, 12, 12))
child, kind = r.nest_patch(SMALL, local, (8, 8, 16, 16))
h = build_hierarchy(SMALL, [(4, 4, 12, 12), (8, 8, 16, 16)])
np.testing.assert_array_equal(h.kinds[0], np.where(local.covered, r.OUTER_GHOST, r.ACTIVE))
np.testing.assert_array_equal(h.kinds[1], kind)
np.testing.assert_array_equal(h.kinds[2], child.fine_kind)
for name in ("ghost_fine", "ghost_coarse", "coalesce_coarse", "coalesce_fine"):
np.testing.assert_array_equal(getattr(h.pairs[1], name), getattr(child, name))
np.testing.assert_array_equal(leaf_partition(h), np.ones((16, 16)))
def test_set_child_under_scalar_parent(self):
local = r.FixedPatch2D.from_config(SMALL, (4, 4, 12, 12))
child, kind = r.nest_patch_set(SMALL, local, [(8, 8, 16, 16)])
h = build_hierarchy(SMALL, [(4, 4, 12, 12), [(8, 8, 16, 16)]])
np.testing.assert_array_equal(h.kinds[1], kind)
np.testing.assert_array_equal(h.pairs[1].coalesce_fine, child.coalesce_fine)
np.testing.assert_array_equal(leaf_partition(h), np.ones((16, 16)))
def test_scalar_child_under_set_parent_is_rejected(self):
with self.assertRaisesRegex(ValueError, "scalar parent"):
build_hierarchy(SMALL, [[(4, 4, 12, 12)], (12, 12, 20, 20)])
def test_wall_rows_and_single_level(self):
h = build_hierarchy(SMALL, [(4, 4, 12, 12)], wall_aware=True)
self.assertTrue((h.kinds[0][[0, -1]] == r.OUTER_GHOST).all())
one = build_hierarchy(SMALL, [])
self.assertEqual(one.depth, 1)
self.assertTrue((one.kinds[0] == r.ACTIVE).all())
if __name__ == "__main__":
unittest.main()