From 1e8e2c9a69ad991a2ac3983fe69a0094a03b04a1 Mon Sep 17 00:00:00 2001 From: TRON Date: Fri, 16 Oct 2026 18:58:30 -0700 Subject: [PATCH] Clarify raw exchange metadata and streamline host tests --- docs/development/raw_exchange.md | 8 + src/CelerisLab/lbm/refinement_stepper.py | 7 +- src/CelerisLab/lbm/refinement_tables.py | 11 +- tests/__init__.py | 1 + tests/unit/test_raw_exchange.py | 245 +++++++++++------------ 5 files changed, 142 insertions(+), 130 deletions(-) create mode 100644 tests/__init__.py diff --git a/docs/development/raw_exchange.md b/docs/development/raw_exchange.md index 169016b..8c7a9a4 100644 --- a/docs/development/raw_exchange.md +++ b/docs/development/raw_exchange.md @@ -20,6 +20,14 @@ their geometric classification when their prescribed velocity is zero. All results are base-lattice units per unit out-of-plane depth. For local spacing and timestep h=dt, raw `(jx,jy,tz,dm)` converts to impulse `(h*dt*jx,h*dt*jy,h*h*dt*tz)` and mass `h*h*dm`. +The accumulator intentionally accepts h and dt separately to make the spatial +and temporal conversion explicit; the current acoustic hierarchy supplies dt=h. +The stored rx/ry arms run from geometry.center to the actual wall hit, not to +the fluid node. For this raw exchange j=c_i(a+b), moving the arm from the fluid +node to the hit adds q*c_i, whose cross product with j is zero. Thus either arm +gives the same raw torque (up to floating-point rounding); this equivalence +must not be assumed for a future nonparallel stress/advective correction. + Force/torque are impulse/angular impulse divided by the reported window length. `dm=b-a` is positive into the exterior fluid, while j is the opposite of the fluid's momentum change across that cut link. Interpolated BFL can introduce diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 1c8155a..44e95a1 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -21,7 +21,8 @@ from .block_stepper import MAX_GRID_DIM from .hierarchy import build_hierarchy, is_rectangle from .refinement import ACTIVE, HALO, FixedPatch2D, FixedPatchSet2D, prepare_fine_boundary from .refinement_tables import (blocked_base_tables, blocked_fine_tables, # noqa: F401 (re-export) - blocked_kind, _prepare_blocked_boundary) + blocked_kind, _prepare_blocked_boundary, + BLOCKED_BOUNDARY_DEVICE_COLUMNS) from .descriptors import SOLID from .raw_exchange import summarize_exchange @@ -335,8 +336,8 @@ class FixedPatchStepper: b = level.boundary if b is not None: b.flags = self._upload(b.tables["flags"]) - columns = ("block", "cell", "dir", "q", "fallback", "surface_velocity") - b.device = {k: self._upload(v) for k, v in b.tables.items() if k in columns} + b.device = {k: self._upload(v) for k, v in b.tables.items() + if k in BLOCKED_BOUNDARY_DEVICE_COLUMNS} if self._observe_raw_exchange: b.device.update(rx=self._upload(b.tables["rx"]), ry=self._upload(b.tables["ry"])) diff --git a/src/CelerisLab/lbm/refinement_tables.py b/src/CelerisLab/lbm/refinement_tables.py index 0f8d5f3..65acd3a 100644 --- a/src/CelerisLab/lbm/refinement_tables.py +++ b/src/CelerisLab/lbm/refinement_tables.py @@ -66,6 +66,11 @@ def blocked_base_tables(patch, layout, order, wall_aware=False, *, composed_kind blocked_kind(layout, dense_kind)) +# Required device columns plus optional actuation; observation arms are uploaded +# separately only when enabled. body_id/aperture remain host reduction metadata. +BLOCKED_BOUNDARY_DEVICE_COLUMNS = ("block", "cell", "dir", "q", "fallback", "surface_velocity") + + @dataclass class _BlockedBoundary: """One level's stationary BFL data and mutable Jet velocity column.""" @@ -87,10 +92,8 @@ def _prepare_blocked_boundary(boundary, layout, cfg, jets): return None fluid_idx, dirs, q, body_id, rx, ry, _, fallback = boundary.curved_host tables, _ = blocked_link_tables(layout, boundary.flags, fluid_idx, dirs, q, fallback) - # Physical row identity excludes body_id: duplicate owners are always wrong. - keys = np.column_stack((fluid_idx, dirs)) - if len(np.unique(keys, axis=0)) != len(keys): - raise ValueError("duplicate physical boundary cut-link") + # blocked_link_tables/map_cut_links already rejects repeated (fluid_idx, dir), + # independently of body_id, and repeated local writers. tables.update(body_id=np.asarray(body_id, np.int32), rx=np.asarray(rx, np.float32), ry=np.asarray(ry, np.float32), aperture=np.zeros(len(dirs), bool)) actuated = bool(jets) or any(o.surface_actuation is not None for o in boundary.objects) diff --git a/tests/__init__.py b/tests/__init__.py new file mode 100644 index 0000000..ff4f8f9 --- /dev/null +++ b/tests/__init__.py @@ -0,0 +1 @@ +"""CelerisLab test package; keeps shared test helpers independent of discovery paths.""" diff --git a/tests/unit/test_raw_exchange.py b/tests/unit/test_raw_exchange.py index 28f8e70..63b3fdc 100644 --- a/tests/unit/test_raw_exchange.py +++ b/tests/unit/test_raw_exchange.py @@ -1,11 +1,9 @@ """CPU reduction and real host scheduling with CUDA doubles, not GPU physics.""" from contextlib import nullcontext from dataclasses import replace -from unittest.mock import Mock, patch - +from unittest.mock import Mock import numpy as np import unittest - from CelerisLab.body.geometry.circle import CircleGeometry from CelerisLab.body.objects import SimObject, SurfaceActuation from CelerisLab.config import LBMConfig @@ -13,7 +11,7 @@ from CelerisLab.lbm.blocks import BlockLayout2D from CelerisLab.lbm.refinement import FixedPatch2D, prepare_fine_boundary from CelerisLab.lbm.refinement_tables import _prepare_blocked_boundary from CelerisLab.lbm.raw_exchange import summarize_exchange -from test_refinement_three_level_multipatch import host_module +from tests.unit.test_refinement_three_level_multipatch import host_module def body(center=(16., 16.), radius=.5, speed=0.): @@ -21,49 +19,6 @@ def body(center=(16., 16.), radius=.5, speed=0.): surface_actuation=SurfaceActuation(.7, .85, speed)) -def check_weights_split_and_window(level): - h = 2. ** -level - rows = np.array([[2, -3, 5, 7], [-1, 4, -2, -6], [99, 99, 99, 99]], np.float32) - kw = dict(body_id=0, level=level, h=h, time_end=3, reference_center=(7, 8)) - last = summarize_exchange(rows, [0, 0, 1], [False, True, False], window='last_substep', **kw) - integral = rows.astype(float)*[h*h, h*h, h**3, h*h]/h - full = summarize_exchange(integral, [0, 0, 1], [False, True, False], window='last_base_step', **kw) - np.testing.assert_array_equal(last['total']['force'], np.array([1, 1])*h) - assert last['total']['torque'] == 3*h*h - assert last['total']['mass_into_fluid'] == h*h - for name in ('force', 'torque', 'mass_rate'): - np.testing.assert_array_equal(last['total'][name], full['total'][name]) - for name in ('force', 'impulse', 'torque', 'angular_impulse', 'mass_into_fluid'): - np.testing.assert_array_equal(last['total'][name], - np.asarray(last['wall'][name])+last['aperture'][name]) - assert last['time_start'] == 3-h and full['time_start'] == 2 - assert last['sample_count'] == 1 and full['sample_count'] == 2**level - assert full['reference_center'] == (7, 8) - - -def test_final_substep_is_not_interval_mean(): - kw = dict(body_id=0, level=1, h=.5, time_end=1, reference_center=(0, 0)) - last = summarize_exchange([[4, 0, 0, 0]], [0], [False], window='last_substep', **kw) - full = summarize_exchange([[1.5, 0, 0, 0]], [0], [False], window='last_base_step', **kw) - assert last['total']['force'][0] == 2 and full['total']['force'][0] == 1.5 - - -def test_host_metadata_zero_jet_and_rows(): - cfg = LBMConfig(nx=32, ny=32, viscosity=.1) - pair = FixedPatch2D.from_config(cfg, (4, 4, 28, 28)) - boundary = prepare_fine_boundary(cfg, pair, [body(radius=2.15)]) - local = replace(cfg, nx=pair.fine_shape[1], ny=pair.fine_shape[0], viscosity=.2) - local.validate() - layout = BlockLayout2D.build(*pair.fine_shape, 8, 8, periodic=(False, False)) - b = _prepare_blocked_boundary(boundary, layout, local, ()) - for key, idx in [('body_id', 3), ('rx', 4), ('ry', 5)]: - np.testing.assert_array_equal(b.tables[key], boundary.curved_host[idx]) - assert b.tables['aperture'].any() and (~b.tables['aperture']).any() - assert not b.tables['surface_velocity'].any() - # Classification does not disappear when prescribed velocity is zero. - assert len(b.tables['cell']) == len(boundary.curved_host[0]) - - def backend(observe=True, obj=None): host = host_module() calls = [] @@ -83,89 +38,133 @@ def backend(observe=True, obj=None): return run, host, calls, arrays -def test_four_level_schedule_publication_and_fixed_storage(): - run, host, calls, arrays = backend() - b = run._level_boundaries[3] - rows = len(b.tables['cell']) - assert b.last_raw.size == rows*16 and b.interval_integral.size == rows*32 - with unittest.TestCase().assertRaisesRegex(RuntimeError, 'completed'): - run.read_raw_exchange(0, level=3) - run.step() - names = [name for name, _, _ in calls] - assert names.count('BlockAccumulateExchange') == 8 - for i, (name, args, kw) in enumerate(calls): - if name == 'BlockAccumulateExchange': - assert names[i-1] == 'BlockCurvedBoundary' and names[i+1] == 'BlockStep' - assert args[0] is b.last_raw and args[1] is b.interval_integral - assert args[3:] == (.125, .125) - assert kw['stream'] is run._stream - assert host.cuda.memset_d8_async.call_count == 1 - value = run.read_raw_exchange(0, level=3) - assert value['sample_count'] == 8 and value['time_start'] == 0 - run.step(0) - assert run.read_raw_exchange(0, level=3)['time_end'] == 1 - assert run.read_raw_exchange(0, level=3, window='last_substep')['time_start'] == .875 - run.step(2) - assert run.read_raw_exchange(0, level=3)['time_start'] == 2 - run.initialize(*arrays) - with unittest.TestCase().assertRaisesRegex(RuntimeError, 'completed'): - run.read_raw_exchange(0, level=3) - run.close() +class RawExchange(unittest.TestCase): + def _check_weights_split_and_window(self, level): + h = 2. ** -level + rows = np.array([[2, -3, 5, 7], [-1, 4, -2, -6], [99, 99, 99, 99]], np.float32) + kw = dict(body_id=0, level=level, h=h, time_end=3, reference_center=(7, 8)) + last = summarize_exchange(rows, [0, 0, 1], [False, True, False], window='last_substep', **kw) + integral = rows.astype(float)*[h*h, h*h, h**3, h*h]/h + full = summarize_exchange(integral, [0, 0, 1], [False, True, False], window='last_base_step', **kw) + np.testing.assert_array_equal(last['total']['force'], np.array([1, 1])*h) + self.assertEqual(last['total']['torque'], 3 * h * h) + self.assertEqual(last['total']['mass_into_fluid'], h * h) + for name in ('force', 'torque', 'mass_rate'): + np.testing.assert_array_equal(last['total'][name], full['total'][name]) + for name in ('force', 'impulse', 'torque', 'angular_impulse', 'mass_into_fluid'): + np.testing.assert_array_equal(last['total'][name], + np.asarray(last['wall'][name])+last['aperture'][name]) + self.assertEqual(last['time_start'], 3 - h) + self.assertEqual(full['time_start'], 2) + self.assertEqual(last['sample_count'], 1) + self.assertEqual(full['sample_count'], 2 ** level) + self.assertEqual(full['reference_center'], (7, 8)) + def test_final_substep_is_not_interval_mean(self): + kw = dict(body_id=0, level=1, h=.5, time_end=1, reference_center=(0, 0)) + last = summarize_exchange([[4, 0, 0, 0]], [0], [False], window='last_substep', **kw) + full = summarize_exchange([[1.5, 0, 0, 0]], [0], [False], window='last_base_step', **kw) + self.assertEqual(last['total']['force'][0], 2) + self.assertEqual(full['total']['force'][0], 1.5) + def test_host_metadata_zero_jet_and_rows(self): + cfg = LBMConfig(nx=32, ny=32, viscosity=.1) + pair = FixedPatch2D.from_config(cfg, (4, 4, 28, 28)) + boundary = prepare_fine_boundary(cfg, pair, [body(radius=2.15)]) + local = replace(cfg, nx=pair.fine_shape[1], ny=pair.fine_shape[0], viscosity=.2) + local.validate() + layout = BlockLayout2D.build(*pair.fine_shape, 8, 8, periodic=(False, False)) + b = _prepare_blocked_boundary(boundary, layout, local, ()) + for key, idx in [('body_id', 3), ('rx', 4), ('ry', 5)]: + np.testing.assert_array_equal(b.tables[key], boundary.curved_host[idx]) + self.assertTrue(b.tables['aperture'].any()) + self.assertTrue((~b.tables['aperture']).any()) + self.assertFalse(b.tables['surface_velocity'].any()) + # Classification does not disappear when prescribed velocity is zero. + self.assertEqual(len(b.tables['cell']), len(boundary.curved_host[0])) -def test_torque_reference_uses_geometry_center_not_optional_object_center(): - for center in (None, (3., 4.)): - obj = SimObject(91, geometry=CircleGeometry(16., 16., .5), center=center, radius=.5) - run, _, _, _ = backend(obj=obj) + def test_four_level_schedule_publication_and_fixed_storage(self): + run, host, calls, arrays = backend() + b = run._level_boundaries[3] + rows = len(b.tables['cell']) + self.assertEqual(b.last_raw.size, rows * 16) + self.assertEqual(b.interval_integral.size, rows * 32) + with self.assertRaisesRegex(RuntimeError, 'completed'): + run.read_raw_exchange(0, level=3) run.step() - result = run.read_raw_exchange(0, level=3) - assert result['reference_center'] == (16., 16.) - assert result['reference_center'] != tuple(obj.center) + names = [name for name, _, _ in calls] + self.assertEqual(names.count('BlockAccumulateExchange'), 8) + for i, (name, args, kw) in enumerate(calls): + if name == 'BlockAccumulateExchange': + self.assertEqual(names[i - 1], 'BlockCurvedBoundary') + self.assertEqual(names[i + 1], 'BlockStep') + self.assertIs(args[0], b.last_raw) + self.assertIs(args[1], b.interval_integral) + self.assertEqual(args[3:], (0.125, 0.125)) + self.assertIs(kw['stream'], run._stream) + self.assertEqual(host.cuda.memset_d8_async.call_count, 1) + value = run.read_raw_exchange(0, level=3) + self.assertEqual(value['sample_count'], 8) + self.assertEqual(value['time_start'], 0) + run.step(0) + self.assertEqual(run.read_raw_exchange(0, level=3)['time_end'], 1) + self.assertEqual(run.read_raw_exchange(0, level=3, window='last_substep')['time_start'], 0.875) + run.step(2) + self.assertEqual(run.read_raw_exchange(0, level=3)['time_start'], 2) + run.initialize(*arrays) + with self.assertRaisesRegex(RuntimeError, 'completed'): + run.read_raw_exchange(0, level=3) run.close() -def test_disabled_has_nulls_and_no_observation_allocation_or_launch(): - run, host, calls, _ = backend(False) - b = run._level_boundaries[3] - assert b.last_raw is b.interval_integral is None - run.step() - for name, args, _ in calls: - if name == 'BlockCurvedBoundary': - assert tuple(args[8:]) == (0, 0, 0) - assert all(n != 'BlockAccumulateExchange' for n, _, _ in calls) - host.cuda.memset_d8_async.assert_not_called() - with unittest.TestCase().assertRaisesRegex(RuntimeError, 'disabled'): - run.read_raw_exchange(0, level=3) - run.close() + def test_torque_reference_uses_geometry_center_not_optional_object_center(self): + for center in (None, (3., 4.)): + obj = SimObject(91, geometry=CircleGeometry(16., 16., .5), center=center, radius=.5) + run, _, _, _ = backend(obj=obj) + run.step() + result = run.read_raw_exchange(0, level=3) + self.assertEqual(result['reference_center'], (16.0, 16.0)) + self.assertNotEqual(result['reference_center'], tuple(obj.center)) + run.close() - -def test_failed_sync_never_publishes_partial_interval(): - run, host, calls, _ = backend() - run.step() - run._stream.synchronize.side_effect = RuntimeError('device failure') - with unittest.TestCase().assertRaisesRegex(RuntimeError, 'device failure'): + def test_disabled_has_nulls_and_no_observation_allocation_or_launch(self): + run, host, calls, _ = backend(False) + b = run._level_boundaries[3] + self.assertIs(b.last_raw, b.interval_integral) + self.assertIs(b.interval_integral, None) run.step() - assert not run._exchange_valid - with unittest.TestCase().assertRaisesRegex(RuntimeError, 'failed'): - run.read_raw_exchange(0, level=3) - run._stream.synchronize.side_effect = None - run.close() + for name, args, _ in calls: + if name == 'BlockCurvedBoundary': + self.assertEqual(tuple(args[8:]), (0, 0, 0)) + self.assertTrue(all((n != 'BlockAccumulateExchange' for n, _, _ in calls))) + host.cuda.memset_d8_async.assert_not_called() + with self.assertRaisesRegex(RuntimeError, 'disabled'): + run.read_raw_exchange(0, level=3) + run.close() + + def test_failed_sync_never_publishes_partial_interval(self): + run, host, calls, _ = backend() + run.step() + run._stream.synchronize.side_effect = RuntimeError('device failure') + with self.assertRaisesRegex(RuntimeError, 'device failure'): + run.step() + self.assertFalse(run._exchange_valid) + with self.assertRaisesRegex(RuntimeError, 'failed'): + run.read_raw_exchange(0, level=3) + run._stream.synchronize.side_effect = None + run.close() + + def test_dense_observation_rejected_before_device_access(self): + host = host_module() + ctx = Mock() + with self.assertRaisesRegex(ValueError, 'blocked body'): + host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28), + context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True) + ctx.activate.assert_not_called() + + def test_all_depth_weights(self): + for level in (1, 2, 3, 7): + self._check_weights_split_and_window(level) -def test_dense_observation_rejected_before_device_access(): - host = host_module() - ctx = Mock() - with unittest.TestCase().assertRaisesRegex(ValueError, 'blocked body'): - host.FixedPatchStepper(LBMConfig(nx=32, ny=32, viscosity=.1), (4, 4, 28, 28), - context=ctx, fine_objects=[body(radius=2.15)], observe_raw_exchange=True) - ctx.activate.assert_not_called() - - -def test_all_depth_weights(): - for level in (1, 2, 3, 7): - check_weights_split_and_window(level) - - -def load_tests(loader, tests, pattern): - return unittest.TestSuite(unittest.FunctionTestCase(value) for name, value in globals().items() - if name.startswith("test_") and callable(value)) +if __name__ == "__main__": + unittest.main()