feat: optional observation-window raw exchange ledger

FixedPatchStepper(observe_window_exchange=True), which requires
observe_raw_exchange, allocates a second per-link FP64 ledger per body
level. While a window is open, the unchanged BlockAccumulateExchange
kernel also accumulates into it after each fine hop. Small API:
start_exchange_window() zeroes and opens at the current completed step;
read_exchange_window(body_id, level=) is non-destructive and reports
impulse with rates over the integer window duration. initialize/
_initialize_raw and step failure invalidate; remap builds a new instance
with no open window; Jet amplitude updates may stay in-window. Default
disabled: no allocation and no extra launch; last_substep/last_base_step
unchanged. summarize_exchange gains window='window' with duration.
Host contract tests (mocked driver) and a GPU gate against summed
per-base reads with a reordering-bound FP64 tolerance.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/CelerisLab/lbm/raw_exchange.py

@@ -7,15 +7,25 @@ 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)
+    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])
@@ -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)

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:

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):

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()

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()