Files
CelerisLab/tests/unit/test_refinement.py

386 lines
20 KiB
Python

"""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_offset
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_offset
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)