diff --git a/src/CelerisLab/lbm/raw_exchange.py b/src/CelerisLab/lbm/raw_exchange.py index d19fbea..3c0e6e7 100644 --- a/src/CelerisLab/lbm/raw_exchange.py +++ b/src/CelerisLab/lbm/raw_exchange.py @@ -7,16 +7,26 @@ import numpy as np def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h, - time_end, window, reference_center): + time_end, window, reference_center, duration=None): """Reduce FP64 interval integrals or convert FP32 last-substep raw rows. dt_l=h_l under acoustic 2:1 refinement. Local momentum j -> base force h*j -> impulse h*dt*j; torque adds one h; mass is h**2*dm. + window="window" reduces an FP64 observation-window ledger spanning + ``duration`` (integer >= 1) base intervals; rates divide by that duration. """ - if window not in ("last_base_step", "last_substep"): + if window not in ("last_base_step", "last_substep", "window"): raise ValueError("unknown exchange window") + if window == "window": + if isinstance(duration, bool) or not isinstance(duration, (int, np.integer)) or duration < 1: + raise ValueError("window duration must be an integer number of base intervals >= 1") + elif duration is not None: + raise ValueError("duration applies only to window='window'") values = np.array(rows, dtype=np.float64, copy=True) - duration = 1.0 if window == "last_base_step" else h + if window == "window": + duration = float(duration) + else: + duration = 1.0 if window == "last_base_step" else h if window == "last_substep": values *= np.array([h*h, h*h, h*h*h, h*h]) selected = np.asarray(body_ids) == body_id @@ -31,6 +41,7 @@ def summarize_exchange(rows, body_ids, aperture, *, body_id, level, h, return dict(wall=reduce(selected & ~aperture), aperture=reduce(selected & aperture), total=reduce(selected), body_id=int(body_id), level=int(level), time_start=time_end-duration, time_end=time_end, dt=h, - sample_count=round(1/h) if window == "last_base_step" else 1, + sample_count=(round(duration/h) if window == "window" + else round(1/h) if window == "last_base_step" else 1), window=window, units="base_lattice_per_unit_depth", definition="raw_boundary_exchange", reference_center=tuple(reference_center), valid=True) diff --git a/src/CelerisLab/lbm/refinement_stepper.py b/src/CelerisLab/lbm/refinement_stepper.py index 143c805..883d59d 100644 --- a/src/CelerisLab/lbm/refinement_stepper.py +++ b/src/CelerisLab/lbm/refinement_stepper.py @@ -172,6 +172,9 @@ class FixedPatchStepper: base interval or last substep in base-lattice units per unit depth. This is boundary momentum reaction, not a validated nozzle/stress load. Dense body layouts are not supported by this new diagnostic in this candidate. + ``observe_window_exchange=True`` (requires observe_raw_exchange) adds a second + per-link FP64 ledger accumulated by the same kernel between + start_exchange_window() and read_exchange_window(); disabled by default. ``jets`` (body.jets.Jet) use the same edge frame as ``fine_objects`` and drive fine cut links selected at their true hits; ``set_jet_amplitude`` re-uploads the per-link velocity column only. @@ -216,7 +219,7 @@ class FixedPatchStepper: fine_block=None, base_block=None, jets=(), child=None, child_block=None, child_objects=(), child_jets=(), child_rectangles_in_parent=None, refinements=None, level_blocks=None, level_objects=None, level_jets=None, - observe_raw_exchange=False): + observe_raw_exchange=False, observe_window_exchange=False): if HALO != 4: raise RuntimeError("refinement.cu requires a four-cell fine halo") if coarse_boundary not in ("periodic", "slot_reservoir", "channel"): @@ -315,6 +318,13 @@ class FixedPatchStepper: for b, layout in zip(self._boundary_hosts, self._layouts)): raise ValueError("raw exchange currently requires blocked body levels") self._observe_raw_exchange = observe_raw_exchange + if not isinstance(observe_window_exchange, bool): + raise ValueError("observe_window_exchange must be bool") + if observe_window_exchange and not observe_raw_exchange: + raise ValueError("observe_window_exchange requires observe_raw_exchange") + self._observe_window = observe_window_exchange + self._window_valid = False + self._window_start = None self._exchange_valid = False self._exchange_spacing = tuple(spacings) self._exchange_centers = tuple(tuple(tuple(o.geometry.center) for o in objects) @@ -365,6 +375,8 @@ class FixedPatchStepper: count = len(b.tables["cell"]) b.last_raw = self._allocate(count * 4 * 4) b.interval_integral = self._allocate(count * 4 * 8) + if self._observe_window: + b.window_integral = self._allocate(count * 4 * 8) b.accumulate = level.module.get_function("BlockAccumulateExchange") b.kernel = level.module.get_function("BlockCurvedBoundary") if self._level_boundaries[1] is not None: @@ -559,6 +571,7 @@ class FixedPatchStepper: cuda.memset_d8(self._fine_curved.link_exchange_gpu, 0, self._fine_curved.count * 12) self._step_count = 0 self._exchange_valid = False + self._window_valid = False self._status = "ready" except BaseException: self._status = "failed" @@ -639,7 +652,8 @@ class FixedPatchStepper: fine_block=fine_block, child_block=child_block, child_rectangles_in_parent=new.child.bounds, child_objects=child_objects, child_jets=jets, - observe_raw_exchange=self._observe_raw_exchange) + observe_raw_exchange=self._observe_raw_exchange, + observe_window_exchange=self._observe_window) try: if self._ddf_shifting: candidate._initialize_raw(_StoredPopulations(transferred.g, transferred.shifted)) @@ -771,6 +785,10 @@ class FixedPatchStepper: b.accumulate(b.last_raw, b.interval_integral, np.uint32(count), np.float64(h), np.float64(h), grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) + if self._window_valid and b.window_integral is not None: + b.accumulate(b.last_raw, b.window_integral, np.uint32(count), + np.float64(h), np.float64(h), + grid=((count + self._block[0] - 1) // self._block[0], 1), **kw) flags = np.uintp(0) if level.boundary is None else level.boundary.flags level.kernel(src, dst, np.uint32(level.layout.halo), flags, level.kind, grid=level.grid, **kw) @@ -804,9 +822,57 @@ class FixedPatchStepper: self._exchange_valid = self._observe_raw_exchange self._status = "ready" except BaseException: + self._window_valid = False self._status = "failed" raise + def start_exchange_window(self): + """Zero the window ledgers and open a window at the current completed step count.""" + with self._lock, self.ctx.activate(): + self._require("ready") + if not self._observe_window: + raise RuntimeError("window exchange observation is disabled") + self._window_valid = False + try: + for b in self._level_boundaries: + if b is not None and b.window_integral is not None: + cuda.memset_d8(b.window_integral, 0, len(b.tables["cell"]) * 32) + except BaseException: + self._status = "failed" + raise + self._window_start = self._step_count + self._window_valid = True + + def read_exchange_window(self, body_id, *, level=1): + """Non-destructive read of [window start, current step count] in base-lattice units.""" + with self._lock, self.ctx.activate(): + self._require("ready") + if not self._observe_window: + raise RuntimeError("window exchange observation is disabled") + if not self._window_valid: + raise RuntimeError("no valid exchange window; call start_exchange_window()") + duration = self._step_count - self._window_start + if duration < 1: + raise RuntimeError("exchange window requires at least one completed base interval") + if (isinstance(level, bool) or not isinstance(level, Integral) + or not 0 <= level < len(self._levels)): + raise ValueError("level must name an existing level") + b = self._level_boundaries[level] + if (b is None or isinstance(body_id, bool) or not isinstance(body_id, Integral) + or not 0 <= body_id < len(b.host.objects)): + raise ValueError("body_id must name a body on this level") + rows = np.empty((len(b.tables["cell"]), 4), np.float64) + try: + cuda.memcpy_dtoh(rows, b.window_integral) + except BaseException: + self._status = "failed" + raise + return summarize_exchange(rows, b.tables["body_id"], b.tables["aperture"], + body_id=body_id, level=level, h=self._exchange_spacing[level], + time_end=self._step_count, window="window", + reference_center=self._exchange_centers[level][body_id], + duration=duration) + def read_raw_exchange(self, body_id, *, level=1, window="last_base_step"): """Read raw boundary reaction about its geometric center in the root edge frame. @@ -885,6 +951,7 @@ class FixedPatchStepper: if b is not None: b.flags = b.device = b.kernel = None b.last_raw = b.interval_integral = b.accumulate = None + b.window_integral = None self._levels.clear() self._pairs = [] if getattr(self, "_fine_curved", None) is not None: diff --git a/src/CelerisLab/lbm/refinement_tables.py b/src/CelerisLab/lbm/refinement_tables.py index b097f30..96423c4 100644 --- a/src/CelerisLab/lbm/refinement_tables.py +++ b/src/CelerisLab/lbm/refinement_tables.py @@ -85,6 +85,7 @@ class _BlockedBoundary: last_raw: object = None interval_integral: object = None accumulate: object = None + window_integral: object = None # optional observation-window FP64 ledger def _prepare_blocked_boundary(boundary, layout, cfg, jets, *, _shifted_blocked=False): diff --git a/tests/unit/test_raw_exchange_window.py b/tests/unit/test_raw_exchange_window.py new file mode 100644 index 0000000..d7f1a73 --- /dev/null +++ b/tests/unit/test_raw_exchange_window.py @@ -0,0 +1,45 @@ +"""Host reduction contract for the optional observation-window exchange ledger.""" +import unittest + +import numpy as np + +from CelerisLab.lbm.raw_exchange import summarize_exchange + + +def reduce(rows, window, duration=None, h=.5, time_end=10): + return summarize_exchange(rows, np.array([0, 0, 1]), np.array([False, True, False]), + body_id=0, level=1, h=h, time_end=time_end, window=window, + reference_center=(0., 0.), duration=duration) + + +class ExchangeWindowReduction(unittest.TestCase): + def test_window_duration_scaling(self): + rows = np.array([[1., 2., 3., 4.], [.5, -1., 0., 1.], [9., 9., 9., 9.]]) + out = reduce(rows, "window", duration=4) + np.testing.assert_array_equal(out["total"]["impulse"], [1.5, 1.]) + np.testing.assert_array_equal(out["total"]["force"], [1.5 / 4, 1. / 4]) + self.assertEqual(out["total"]["torque"], 3. / 4) + self.assertEqual(out["total"]["mass_rate"], 5. / 4) + np.testing.assert_array_equal(out["aperture"]["impulse"], [.5, -1.]) + self.assertEqual((out["time_start"], out["time_end"]), (6., 10)) + self.assertEqual(out["sample_count"], 8) # duration / h fine substeps + self.assertEqual(out["window"], "window") + + def test_existing_windows_unchanged(self): + rows = np.array([[1., 2., 3., 4.], [0., 0., 0., 0.], [0., 0., 0., 0.]]) + base = reduce(rows, "last_base_step") + self.assertEqual((base["time_start"], base["sample_count"]), (9., 2)) + np.testing.assert_array_equal(base["total"]["force"], [1., 2.]) + sub = reduce(rows, "last_substep") + self.assertEqual(sub["sample_count"], 1) + + def test_duration_validation(self): + rows = np.zeros((3, 4)) + for bad in (None, 0, -1, 1.5, True): + with self.subTest(duration=bad), self.assertRaises(ValueError): + reduce(rows, "window", duration=bad) + with self.assertRaises(ValueError): + reduce(rows, "last_base_step", duration=2) + + +if __name__ == "__main__": unittest.main() diff --git a/tests/validation/multigrid/test_exchange_window.py b/tests/validation/multigrid/test_exchange_window.py new file mode 100644 index 0000000..8776e01 --- /dev/null +++ b/tests/validation/multigrid/test_exchange_window.py @@ -0,0 +1,150 @@ +"""Optional observation-window exchange ledger: host contract (mocked driver) and GPU gate. + +Run CPU and GPU classes in separate processes; GPU callers select the device. +GPU tolerance: window rows are summed per link over time, the per-step reads over links +then over time. Reordering a sum of m FP64 terms changes it by at most +2*(m-1)*eps*sum|terms| (each order is within (m-1)*eps*sum|terms| of the exact sum), so the +bound is scaled by the absolute sum of the per-step row integrals and is robust to cancellation. +""" +from contextlib import ExitStack +import unittest +from unittest.mock import patch + +import numpy as np +import pycuda.driver as cuda + +from CelerisLab.lbm.refinement_stepper import FixedPatchStepper +from test_refinement_blocked_base import mocked_backend +from test_refinement_blocked_fine import BOUNDS, circle, config + +EPS = np.finfo(np.float64).eps + + +def build(m, window=True, raw=True): + return FixedPatchStepper(config("periodic"), BOUNDS, context=m["context"], base_block=(3, 4), + fine_block=(16, 20), fine_objects=[circle()], + observe_raw_exchange=raw, observe_window_exchange=window) + + +def accumulations(m): + return [c for c in m["calls"] if c[0] == "BlockAccumulateExchange"] + + +class ExchangeWindowCPU(unittest.TestCase): + def test_constructor_validation(self): + with mocked_backend() as m: + with self.assertRaisesRegex(ValueError, "requires observe_raw_exchange"): + build(m, window=True, raw=False) + with self.assertRaisesRegex(ValueError, "must be bool"): + build(m, window=1) + + def test_state_machine_and_launches(self): + with mocked_backend() as m, ExitStack() as stack: + stack.enter_context(patch.object(cuda, "memset_d8")) + stack.enter_context(patch.object(cuda, "memset_d8_async")) + stack.enter_context(patch.object(cuda, "memcpy_dtoh", + side_effect=lambda out, src: out.fill(0))) + backend = build(m) + b = backend._level_boundaries[1] + self.assertIsNotNone(b.window_integral) + with self.assertRaises(RuntimeError): + backend.start_exchange_window() # not initialized + arrays = [np.full(level.shape, 1 / 9, np.float32) for level in backend._levels] + backend.initialize(*arrays) + with self.assertRaisesRegex(RuntimeError, "no valid exchange window"): + backend.read_exchange_window(0) + backend.step(1) + self.assertEqual(len(accumulations(m)), 2) # closed window: interval ledger only + backend.start_exchange_window() + with self.assertRaisesRegex(RuntimeError, "at least one completed"): + backend.read_exchange_window(0) + m["calls"].clear() + backend.step(2) + calls = accumulations(m) + self.assertEqual(len(calls), 8) # two fine hops x two ledgers x two base steps + self.assertEqual({id(c[1][1]) for c in calls}, + {id(b.interval_integral), id(b.window_integral)}) + out = backend.read_exchange_window(0) + self.assertEqual((out["time_start"], out["time_end"], out["sample_count"]), (1., 3, 4)) + self.assertEqual(backend.read_exchange_window(0)["time_end"], 3) # non-destructive + backend.initialize(*arrays) # invalidates + with self.assertRaisesRegex(RuntimeError, "no valid exchange window"): + backend.read_exchange_window(0) + backend.start_exchange_window() + with patch.object(backend, "_hop", side_effect=RuntimeError("injected")): + with self.assertRaises(RuntimeError): + backend.step(1) + self.assertFalse(backend._window_valid) + with self.assertRaises(RuntimeError): + backend.read_exchange_window(0) # failed status + + def test_disabled_default_has_no_window_cost(self): + with mocked_backend() as m, ExitStack() as stack: + stack.enter_context(patch.object(cuda, "memset_d8_async")) + backend = build(m, window=False) + self.assertIsNone(backend._level_boundaries[1].window_integral) + backend.initialize(*[np.full(l.shape, 1 / 9, np.float32) for l in backend._levels]) + backend.step(1) + self.assertEqual(len(accumulations(m)), 2) + with self.assertRaisesRegex(RuntimeError, "disabled"): + backend.start_exchange_window() + + +class ExchangeWindowGPU(unittest.TestCase): + def test_window_matches_summed_base_reads_and_populations_unchanged(self): + from test_raw_exchange_gpu import BOUNDS as RB, CFG, circle as rc, jet, state + from CelerisLab.cuda.context import CudaContext + with CudaContext(0) as ctx: + kw = dict(context=ctx, fine_block=(16, 20), fine_objects=[rc()], jets=[jet()], + observe_raw_exchange=True) + with FixedPatchStepper(CFG, RB, observe_window_exchange=True, **kw) as windowed, \ + FixedPatchStepper(CFG, RB, **kw) as plain: + arrays = [state(level.shape) for level in windowed._levels] + windowed.initialize(*arrays) + plain.initialize(*arrays) + windowed.step(2) + plain.step(2) + windowed.start_exchange_window() + b = windowed._level_boundaries[1] + count = len(b.tables["cell"]) + keys = ("impulse", "angular_impulse", "mass_into_fluid") + summed = {part: [0., 0., 0., 0.] for part in ("total", "wall", "aperture")} + scale = np.zeros(4) + steps = 0 + for speed in (0., .002, .002, -.001, 0., 0.): # amplitude updates stay in-window + windowed.set_jet_amplitude(0, speed) + plain.set_jet_amplitude(0, speed) + windowed.step(1) + plain.step(1) + steps += 1 + obs = windowed.read_raw_exchange(0) + for part in summed: + o = obs[part] + summed[part] = [summed[part][0] + o["impulse"][0], summed[part][1] + o["impulse"][1], + summed[part][2] + o["angular_impulse"], + summed[part][3] + o["mass_into_fluid"]] + rows = np.empty((count, 4), np.float64) + cuda.memcpy_dtoh(rows, b.interval_integral) + scale += np.abs(rows).sum(axis=0) + for a, e in zip(windowed._download_raw().values, plain._download_raw().values): + np.testing.assert_array_equal(np.isnan(a), np.isnan(e)) + np.testing.assert_array_equal(a[np.isfinite(a)], e[np.isfinite(e)]) + out = windowed.read_exchange_window(0) + self.assertEqual((out["time_end"] - out["time_start"], out["sample_count"]), + (steps, 2 * steps)) + bound = 2 * (count * steps) * EPS * scale + np.finfo(np.float64).tiny + for part in summed: + o = out[part] + got = np.array([o["impulse"][0], o["impulse"][1], o["angular_impulse"], + o["mass_into_fluid"]]) + diff = np.abs(got - np.array(summed[part])) + self.assertTrue((diff <= bound).all(), f"{part}: {diff} > {bound}") + np.testing.assert_array_equal(o["force"], o["impulse"] / steps) + again = windowed.read_exchange_window(0) + np.testing.assert_array_equal(again["total"]["impulse"], out["total"]["impulse"]) + windowed.initialize(*arrays) + with self.assertRaises(RuntimeError): + windowed.read_exchange_window(0) + + +if __name__ == "__main__": unittest.main()