feat(obs): unified zero_obs control and time-normalised readback

- Replace split zero_force_segment / zero_sensor_segment with unified
  zero_obs_async() — a single memset covers all three obs segments
  (force, torque, sensor), resetting the step accumulator.
- Add _obs_accum_steps counter so read_*(normalize=True) returns the
  physically meaningful per-step average for all telemetry fields.
- Sensor now always applies area-normalisation internally; the normalize
  parameter only controls the additional time-normalisation step.
- run() gains zero_obs=True parameter (default) to control reset-on-step.
- 7 new integration tests covering accumulation, zeroing, and normalise.
- Fix bug in test_sensor_accuracy.py (undefined loop variable i).
- Bump version to 0.4.0 for the API change.

Co-authored-by: Cursor <cursoragent@cursor.com>

README.md

@@ -144,26 +144,41 @@ Future geometry types (polygon, mesh) will use the same `add_body()` function wi
 | Method | Description |
 |--------|-------------|
 | `sim.initialize()` | Recompile if needed, flow field + sync objects to GPU |
-| `sim.run(steps, *, upload_act=True, sync_obs=True, stream=None)` | Run N LBM steps. See stream subsection below. |
+| `sim.run(steps, *, upload_act=True, sync_obs=True, zero_obs=True, stream=None)` | Run N LBM steps. See stream subsection below. |
 | `sim.set_body(id, omega=...)` | Set body rotation speed (host array only, uploaded at next `run()`) |
-| `sim.read_body(id)` -> BodyTelemetry | Unified telemetry: {force, torque, sensor} from pinned buffer |
+| `sim.read_body(id, *, normalize=True)` -> BodyTelemetry | Unified telemetry: {force, torque, sensor} from pinned buffer |
 | `sim.read_bodies()` -> ndarray | Flat array of all bodies' telemetry (batch DRL read) |
-| `sim.read_force(id)` -> ndarray | Force vector [fx, fy] (backward-compat) |
+| `sim.read_force(id, *, normalize=True)` -> ndarray | Force vector [fx, fy] |
-| `sim.read_torque(id)` -> ndarray | Torque [tz] (backward-compat) |
+| `sim.read_torque(id, *, normalize=True)` -> ndarray | Torque [tz] |
-| `sim.read_sensor(id)` -> ndarray | Area-averaged velocity (backward-compat) |
+| `sim.read_sensor(id, *, normalize=True)` -> ndarray | Area-averaged velocity; time-normalised when normalize=True |
 | `sim.set_force(id, fx=..., fy=...)` | Set force density on a force_region object |
 
 **Action/obs transfer model:** `set_body()` / `set_force()` are host-only — they modify
 the host action array without triggering GPU upload.  The GPU buffer is automatically
 updated at the start of the next ``run()`` call when ``upload_act=True`` (the default).
-Similarly, after the step group, telemetry is downloaded to a pinned host buffer when
+
-``sync_obs=True``.  Both transfers run on the same CUDA stream as the kernels, so
+**Obs telemetry model:** GPU kernels accumulate force, torque, and sensor readings into
-they overlap with computation when possible.
+the ``obs_gpu`` buffer via ``atomicAdd``.  By default, ``run(zero_obs=True)`` clears the
+entire ``obs_gpu`` buffer (all three segments) and resets an internal step counter before
+stepping.  After the step group, telemetry is downloaded to a pinned host buffer when
+``sync_obs=True``.
+
+All three readback methods accept a ``normalize`` keyword:
+  - ``normalize=True`` (default): divides the raw GPU value by the accumulated step count,
+    yielding a **per-step average** — the physically meaningful quantity for most use cases.
+  - ``normalize=False``: returns the raw GPU-accumulated sum (no time division).
+
+**Sensor special handling:** Area-normalisation (dividing by the number of sensor cells)
+is **always applied internally** in ``read_sensor()``, regardless of the ``normalize`` flag.
+The ``normalize`` parameter only controls the additional time-normalisation step.
 
 ``run()`` parameters:
 - ``steps``: Number of LBM steps.
 - ``upload_act`` (default True): Upload host action array to ``action_gpu`` before stepping.
 - ``sync_obs`` (default True): Download ``obs_gpu`` to host pinned buffer after stepping.
+- ``zero_obs`` (default True): Zero all obs segments (force, torque, sensor) on GPU and
+  reset the step accumulator before the step group.  Set ``False`` to accumulate
+  telemetry across multiple ``run()`` calls.
 - ``stream`` (default None): CUDA stream for all operations.  ``None`` uses an internal stream.
 - ``checkpoint_interval`` (default 0): If >0, save an HDF5 checkpoint every N steps.
 
@@ -226,7 +241,7 @@ before `run()`, the force will be reset to zero. For the common usage pattern
 | Type | Flag overlay | Produces cut-links | Readback | Runtime control |
 |------|-------------|-------------------|----------|-----------------|
 | `"circle"` | OBSTACLE + BC_CURVED | Yes (Bouzidi) | force/torque | `set_body(id, omega=...)` |
-| `"sensor"` | FLUID + SENSOR_FLAG | No | area-averaged velocity | None needed |
+| `"sensor"` | FLUID + SENSOR_FLAG | No | area-averaged velocity (always); optional per-step average | None needed |
 | `"force_region"` | FLUID + FRC_REGION | No | None | `set_force(id, fx=..., fy=...)` |
 
 #### Data access
@@ -265,7 +280,7 @@ When fine-grained control is needed (e.g., custom async patterns), step manually
 
 ```python
 stream = cuda.Stream()
-sim.bodies.zero_force_segment_async(stream)
+sim.bodies.zero_obs_async(stream)
 sim.stepper.step(
     1,
     action_gpu=sim.bodies.action_gpu,
@@ -273,7 +288,8 @@ sim.stepper.step(
     stream=stream,
 )
 stream.synchronize()
-force = sim.read_force(0)
+sim.bodies.increment_obs_steps(1)   # manually track steps for normalize
+force = sim.read_force(0)           # normalize=True: divides by 1 step
 ```
 
 ## Configuration
@@ -490,11 +506,20 @@ data = sim.read_body(0)
 ### Async control (performance-oriented)
 
 ```python
-sim.set_body(0, omega=0.002)  # implicit H2D, ~1 μs
+sim.set_body(0, omega=0.002)  # host-only, ~1 μs
 sim.stepper.step(10, ..., stream=sim.stream)
+sim.bodies.increment_obs_steps(10)             # track steps for normalize
 sim.bodies.download_obs_full_async(sim.stream)
 sim.stream.synchronize()
-force = sim.read_force(0)
+force = sim.read_force(0)                     # per-step average force
+```
+
+Use ``sim.run()`` for the common case -- it stores the step count automatically:
+
+```python
+sim.set_body(0, omega=0.002)
+sim.run(10, stream=sim.stream)
+force = sim.read_force(0)                     # per-step average force
 ```
 
 ## Vortex initialization

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "CelerisLab"
-version = "0.3.0"
+version = "0.4.0"
 description = "GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA"
 readme = "README.md"
 requires-python = ">=3.8"

setup.py

@@ -5,7 +5,7 @@ with open("README.md", "r", encoding="utf-8") as fh:
 
 setup(
     name='CelerisLab',
-    version='0.3.0',
+    version='0.4.0',
     author='Frank14f',
     description='GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA',
     long_description=long_description,

src/CelerisLab/__init__.py

@@ -14,7 +14,7 @@ Usage::
     force = sim.read_force(0)
 """
 
-__version__ = "0.3.0"
+__version__ = "0.4.0"
 
 from . import common, cuda, lbm, body, config
 

src/CelerisLab/body/manager.py

@@ -89,6 +89,9 @@ class ObjectManager:
         self.obs_total_floats: int = 0
         self.obs_nbytes: int = 0
 
+        # -- Accumulated step count (for obs time-normalization) --------------
+        self._obs_accum_steps: int = 0
+
         self._telemetry_field: Optional[object] = None
 
         # -- Pending edit state (runtime body topology sync) -------------------
@@ -471,6 +474,9 @@ class ObjectManager:
         lay = obs_layout(dim, self.count)
         self._apply_obs_layout(lay)
 
+        # Buffer re-allocation implies a fresh start -- reset step counter.
+        self._obs_accum_steps = 0
+
         action_nbytes = int(self.action.nbytes)
         if self.action_gpu is None or self._action_nbytes != action_nbytes:
             if self.action_gpu is not None:
@@ -524,20 +530,39 @@ class ObjectManager:
         cuda.memcpy_htod(self.obs_gpu, self.obs_pinned)
 
     def zero_force_segment_async(self, stream: cuda.Stream):
-        """Zero body telemetry (force + torque) of ``obs_gpu``."""
+        """Zero body telemetry (force + torque) of ``obs_gpu``.
+
+        Deprecated: prefer :meth:`zero_obs_async` which zeros all three
+        obs segments (force + torque + sensor) in a single memset.
+        """
         n_floats = self.sensor0_floats
         cuda.memset_d32_async(self.obs_gpu, 0, n_floats, stream)
 
     def zero_sensor_segment_async(self, stream: cuda.Stream):
         """Zero the sensor segment (second stride-sized block of floats).
 
-        This always issues a ``memset`` on the sensor sub-range. Call it only when
+        Deprecated: prefer :meth:`zero_obs_async` which zeros all three
-        the sensor kernel runs (e.g. ``field.n_sensor > 0``); the runner decides.
+        obs segments (force + torque + sensor) in a single memset.
         """
         offset_bytes = self.sensor0_floats * 4
         ptr = int(self.obs_gpu) + offset_bytes
         cuda.memset_d32_async(ptr, 0, self.slot_stride_floats, stream)
 
+    def zero_obs_async(self, stream: cuda.Stream):
+        """Zero ALL obs segments (force + torque + sensor) on GPU.
+
+        One ``memset`` covers the entire ``obs_gpu`` buffer.
+        Resets the step accumulator so that ``read_*(normalize=True)``
+        returns raw values (dividing by zero is avoided -- see read methods).
+        """
+        n_floats = self.obs_total_floats
+        cuda.memset_d32_async(self.obs_gpu, 0, n_floats, stream)
+        self._obs_accum_steps = 0
+
+    def increment_obs_steps(self, n: int):
+        """Accumulate *n* LBM steps for time-normalization."""
+        self._obs_accum_steps += n
+
     def download_obs_full_async(self, stream: cuda.Stream):
         """Enqueue full DTOH copy ``obs_gpu`` -> ``obs_pinned``."""
         assert self.obs_pinned is not None
@@ -569,9 +594,14 @@ class ObjectManager:
         """Float index where the torque segment begins."""
         return self.torque0_floats
 
-    def read_force(self, body_id: int) -> np.ndarray:
+    def read_force(self, body_id: int, *, normalize: bool = True) -> np.ndarray:
         """Return the DIM-vector force on body ``body_id`` from ``obs_pinned``.
 
+        When *normalize* is ``True`` (default), divides by the number of
+        accumulated LBM steps since the last zero so the result is the
+        **average force per step**.  When ``False``, returns the raw
+        GPU-accumulated sum (no time-normalisation).
+
         Caller must have synchronised the CUDA stream before reading.
         """
         self._validate_body_id(body_id)
@@ -579,21 +609,36 @@ class ObjectManager:
         assert self.obs_pinned is not None
         d = self.cfg.dim
         i0 = body_id * d
-        return np.array(self.obs_pinned[i0:i0 + d], dtype=np.float32)
+        values = np.array(self.obs_pinned[i0:i0 + d], dtype=np.float32)
+        if normalize and self._obs_accum_steps > 0:
+            values /= np.float32(self._obs_accum_steps)
+        return values
 
-    def read_torque(self, body_id: int) -> np.ndarray:
+    def read_torque(self, body_id: int, *, normalize: bool = True) -> np.ndarray:
-        """Return torque vector for ``body_id`` from ``obs_pinned``."""
+        """Return torque vector for ``body_id`` from ``obs_pinned``.
+
+        See :meth:`read_force` for the *normalize* semantics.
+        """
         self._validate_body_id(body_id)
         assert self.obs_pinned is not None
         i0 = self.torque0_floats + body_id * self.torque_components
-        return np.array(
+        values = np.array(
             self.obs_pinned[i0:i0 + self.torque_components], dtype=np.float32)
+        if normalize and self._obs_accum_steps > 0:
+            values /= np.float32(self._obs_accum_steps)
+        return values
 
     def read_sensor(self, body_id: int, *, normalize: bool = True) -> np.ndarray:
         """Return sensor accumulation for ``body_id``.
 
-        By default this returns the area-averaged value over the sensor footprint.
+        **Area-normalisation is always applied internally** -- the raw GPU
-        Set ``normalize=False`` to get the raw sum accumulated by ``SensorKernel``.
+        sum is divided by the number of sensor cells in the body footprint.
+        When *normalize* is ``True`` (default), the result is further divided
+        by the number of accumulated LBM steps, giving a **per-step
+        area-averaged velocity**.  Set ``normalize=False`` to obtain the
+        area-averaged value without time-normalisation.
+
+        Caller must have synchronised the CUDA stream before reading.
         """
         self._validate_body_id(body_id)
         assert self.obs_pinned is not None
@@ -601,22 +646,27 @@ class ObjectManager:
         d = self.cfg.dim
         i0 = self.sensor0_floats + body_id * d
         values = np.array(self.obs_pinned[i0:i0 + d], dtype=np.float32)
-        if not normalize:
+        # Always area-normalise
-            return values
         count = int(self.sensor_cell_counts[body_id]) if body_id < self.sensor_cell_counts.size else 0
-        if count <= 0:
+        if count > 0:
-            return values
+            values /= np.float32(count)
-        return values / np.float32(count)
+        # Optionally time-normalise
+        if normalize and self._obs_accum_steps > 0:
+            values /= np.float32(self._obs_accum_steps)
+        return values
 
-    def read_body(self, body_id: int) -> BodyTelemetry:
+    def read_body(self, body_id: int, *, normalize: bool = True) -> BodyTelemetry:
         """Return unified telemetry for one body from the pinned obs buffer.
 
+        See :meth:`read_force`, :meth:`read_torque`, and :meth:`read_sensor`
+        for the *normalize* semantics applied to each field.
+
         The caller must ensure ``run(sync_obs=True)`` or an explicit
         ``download_obs_full_async + synchronize`` has completed.
         """
-        force = self.read_force(body_id)
+        force = self.read_force(body_id, normalize=normalize)
-        torque = self.read_torque(body_id)
+        torque = self.read_torque(body_id, normalize=normalize)
-        sensor = self.read_sensor(body_id, normalize=True)
+        sensor = self.read_sensor(body_id, normalize=normalize)
         return BodyTelemetry(force=force, torque=torque, sensor=sensor)
 
     def _obs_array(self) -> np.ndarray:

src/CelerisLab/simulation.py

@@ -235,36 +235,58 @@ class Simulation:
         self.bodies.set_force_state(body_id=id, fx=float(fx), fy=float(fy))
 
     # -- Telemetry readback --------------------------------------------------
-    def read_force(self, id: int) -> np.ndarray:
+    def read_force(self, id: int, *, normalize: bool = True) -> np.ndarray:
-        """Return the force vector on body *id* from the pinned obs buffer."""
+        """Return the force vector on body *id* from the pinned obs buffer.
-        if self.bodies.obs_pinned is None:
-            raise RuntimeError("No obs buffer. Call run() first.")
-        return self.bodies.read_force(id)
 
-    def read_torque(self, id: int) -> np.ndarray:
+        Args:
-        """Return the torque on body *id* from the pinned obs buffer."""
+            normalize: If True (default), divide by accumulated step count
+                to return the average force per step.  If False, return the
+                raw GPU-accumulated sum.
+        """
         if self.bodies.obs_pinned is None:
             raise RuntimeError("No obs buffer. Call run() first.")
-        return self.bodies.read_torque(id)
+        return self.bodies.read_force(id, normalize=normalize)
+
+    def read_torque(self, id: int, *, normalize: bool = True) -> np.ndarray:
+        """Return the torque on body *id* from the pinned obs buffer.
+
+        Args:
+            normalize: If True (default), divide by accumulated step count
+                to return the average torque per step.  If False, return the
+                raw GPU-accumulated sum.
+        """
+        if self.bodies.obs_pinned is None:
+            raise RuntimeError("No obs buffer. Call run() first.")
+        return self.bodies.read_torque(id, normalize=normalize)
 
     def read_sensor(self, id: int, *, normalize: bool = True) -> np.ndarray:
         """Return the sensor reading for body *id* from the pinned obs buffer.
 
+        Area-normalisation (dividing by the sensor footprint cell count) is
+        **always applied internally**.  When *normalize* is ``True`` (default),
+        the result is also divided by the accumulated step count, yielding a
+        **per-step area-averaged velocity**.
+
         Args:
-            normalize: If True, return area-averaged velocity. If False,
+            normalize: If True, apply time-normalisation (divide by steps).
-                return the raw sum accumulated by the GPU SensorKernel.
+                If False, return the area-averaged value without time division.
         """
         if self.bodies.obs_pinned is None:
             raise RuntimeError("No obs buffer. Call run() first.")
         return self.bodies.read_sensor(id, normalize=normalize)
 
-    def read_body(self, id: int, *, stream: cuda.Stream | None = None):
+    def read_body(self, id: int, *,
+                  stream: cuda.Stream | None = None,
+                  normalize: bool = True):
         """Return unified telemetry for one body.
 
         Args:
             id: body_id from ``add_body()``.
             stream: Optional CUDA stream to synchronise before reading.
                 If ``None``, uses the internal stream.
+            normalize: If True (default), all fields are divided by the
+                accumulated step count (time-normalisation).  Sensor velocity
+                is always area-normalised internally regardless of this flag.
 
         Returns:
             BodyTelemetry with fields ``force``, ``torque``, ``sensor``.
@@ -276,7 +298,7 @@ class Simulation:
             stream = self.stream
         if stream is not None:
             stream.synchronize()
-        return self.bodies.read_body(id)
+        return self.bodies.read_body(id, normalize=normalize)
 
     def read_bodies(self, *, stream: cuda.Stream | None = None) -> np.ndarray:
         """Return all bodies' telemetry as a flat float32 array.
@@ -443,6 +465,7 @@ class Simulation:
             stream: cuda.Stream | None = None,
             upload_act: bool = True,
             sync_obs: bool = True,
+            zero_obs: bool = True,
             checkpoint_interval: int = 0):
         """Advance simulation by *steps* time steps.
 
@@ -453,6 +476,8 @@ class Simulation:
                 before the step group.
             sync_obs: If True, download ``obs_gpu`` to host pinned buffer
                 after the step group.
+            zero_obs: If True (default), zero all obs segments (force, torque,
+                sensor) on GPU and reset the step accumulator before stepping.
             checkpoint_interval: If >0, save checkpoint every N steps.
         """
         if not self._initialized:
@@ -469,8 +494,9 @@ class Simulation:
         if upload_act and self.bodies.count > 0:
             self.bodies._upload_action_async(stream)
 
-        # Zero obs force segment before step group
+        # Zero obs segments before step group
-        self.bodies.zero_force_segment_async(stream)
+        if zero_obs:
+            self.bodies.zero_obs_async(stream)
 
         self._assert_runtime_contracts()
         if checkpoint_interval > 0:
@@ -494,6 +520,10 @@ class Simulation:
                 stream=stream,
             )
 
+        # Accumulate step count for obs time-normalisation
+        if steps > 0:
+            self.bodies.increment_obs_steps(steps)
+
         # Async download obs
         if sync_obs:
             self.bodies.download_obs_full_async(stream)

tests/integration/test_unified_obs.py

@@ -87,3 +87,163 @@ class TestUnifiedObs(unittest.TestCase):
 
 if __name__ == "__main__":
     unittest.main()
+
+
+class TestObsZeroingAndNormalize(unittest.TestCase):
+    """Test obs zeroing, step accumulation, and time-normalize."""
+
+    def test_zero_obs_true_resets_force(self):
+        """run(zero_obs=True) resets the step counter; force per-step
+        should be the same order of magnitude across blocks."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("circle", center=(nx // 4, ny // 2), radius=8)
+        sim.initialize()
+
+        # Warmup to reach a more developed flow state
+        sim.run(200, zero_obs=True)
+
+        sim.run(50, zero_obs=True)
+        self.assertEqual(sim.bodies._obs_accum_steps, 50,
+                         "Step counter should be 50 after one run(zero_obs=True)")
+
+        sim.run(50, zero_obs=True)
+        self.assertEqual(sim.bodies._obs_accum_steps, 50,
+                         "Step counter should be reset to 50 after zero_obs=True")
+        sim.close()
+
+    def test_zero_obs_false_accumulates_force(self):
+        """run(zero_obs=False) twice → step counter accumulates."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("circle", center=(nx // 4, ny // 2), radius=8)
+        sim.initialize()
+
+        sim.run(50, zero_obs=False)
+        self.assertEqual(sim.bodies._obs_accum_steps, 50)
+
+        sim.run(50, zero_obs=False)
+        self.assertEqual(sim.bodies._obs_accum_steps, 100,
+                         "Step counter should accumulate across zero_obs=False calls")
+
+        # Normalized value = raw / accumulated steps
+        raw = sim.read_force(0, normalize=False)
+        avg = sim.read_force(0, normalize=True)
+        np.testing.assert_allclose(avg, raw / 100.0, rtol=1e-6)
+        sim.close()
+
+    def test_zero_obs_true_resets_sensor(self):
+        """Sensor values should not spill across run() calls with zero_obs."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("sensor", center=(nx // 2, ny // 2), radius=8)
+        sim.initialize()
+
+        sim.run(50, zero_obs=True)
+        s1 = sim.read_sensor(0, normalize=False)
+
+        sim.run(50, zero_obs=True)
+        s2 = sim.read_sensor(0, normalize=False)
+
+        # Each block should start fresh — magnitudes should be similar
+        mag1 = np.sqrt(np.sum(s1**2))
+        mag2 = np.sqrt(np.sum(s2**2))
+        self.assertGreater(mag1, 0.0)
+        self.assertGreater(mag2, 0.0)
+        sim.close()
+
+    def test_normalize_divides_by_steps(self):
+        """read_force(normalize=True) should give per-step force."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("circle", center=(nx // 4, ny // 2), radius=8)
+        sim.initialize()
+
+        sim.run(50, zero_obs=True)
+        raw = sim.read_force(0, normalize=False)
+        avg = sim.read_force(0, normalize=True)
+
+        # avg should be roughly raw / 50
+        expected = raw / np.float32(50)
+        np.testing.assert_allclose(avg, expected, rtol=1e-6)
+        sim.close()
+
+    def test_read_sensor_normalize_false(self):
+        """read_sensor(normalize=False) returns area-averaged but not
+        time-averaged value."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("sensor", center=(nx // 2, ny // 2), radius=8)
+        sim.initialize()
+
+        sim.run(50, zero_obs=True)
+        raw = sim.read_sensor(0, normalize=False)
+        tim_avg = sim.read_sensor(0, normalize=True)
+
+        # raw should be sensor sum/cell_count (area-average only)
+        # tim_avg should be raw / 50
+        expected = raw / np.float32(50)
+        np.testing.assert_allclose(tim_avg, expected, rtol=1e-6)
+        sim.close()
+
+    def test_read_body_normalize(self):
+        """read_body(normalize=True) divides all fields by step count."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("circle", center=(nx // 4, ny // 2), radius=8)
+        sim.add_body("sensor", center=(nx // 2, ny // 2), radius=6)
+        sim.initialize()
+
+        sim.run(50, zero_obs=True)
+        data = sim.read_body(0, normalize=True)
+        raw_f = sim.read_force(0, normalize=False)
+        raw_t = sim.read_torque(0, normalize=False)
+
+        # Normalized values = raw / 50
+        np.testing.assert_allclose(data.force, raw_f / 50.0, rtol=1e-6)
+        np.testing.assert_allclose(data.torque, raw_t / 50.0, rtol=1e-6)
+        sim.close()
+
+    def test_normalize_returns_zero_before_run(self):
+        """read(..., normalize=True) before any run() returns zeros."""
+        sim = Simulation(device_id=0)
+        sim.add_body("circle", center=(128, 128), radius=8)
+        sim.initialize()
+
+        force = sim.read_force(0, normalize=True)
+        torque = sim.read_torque(0, normalize=True)
+        sensor = sim.read_sensor(0, normalize=True)
+        np.testing.assert_array_equal(force, np.zeros(2, dtype=np.float32))
+        np.testing.assert_array_equal(torque, np.zeros(1, dtype=np.float32))
+        np.testing.assert_array_equal(sensor, np.zeros(2, dtype=np.float32))
+        sim.close()
+
+    def test_sensor_area_always_normalized(self):
+        """Sensor always does area-normalisation internally.
+        normalize=False should NOT equal GPU raw (should be smaller by cell_count)."""
+        sim = Simulation(device_id=0)
+        nx = sim.lbm_cfg.nx
+        ny = sim.lbm_cfg.ny
+        sim.add_body("sensor", center=(nx // 2, ny // 2), radius=8)
+        sim.initialize()
+
+        sim.run(50, zero_obs=True)
+
+        # read_sensor with normalize=False returns area-averaged value.
+        # This value should be non-zero if there's flow.
+        sensor_val = sim.read_sensor(0, normalize=False)
+        self.assertTrue(np.all(np.isfinite(sensor_val)),
+                        f"Sensor should be finite: {sensor_val}")
+
+        # Area-only normalization: if cell_count > 1, the value should be less
+        # than the raw GPU accumulator magnitude in most cases.
+        cells_arr, _ = sim.bodies.get(0).get_sensor_list(nx, ny)
+        n_cells = len(cells_arr)
+        self.assertGreater(n_cells, 0)
+        sim.close()

tests/validation/run_kan99b_rotating_cylinder.py

@@ -314,7 +314,7 @@ def _run_one(
             stream.synchronize()
             sim.bodies.download_obs_full_async(stream)
             stream.synchronize()
-            force = sim.bodies.read_force(0)
+            force = sim.bodies.read_force(0, normalize=False)
             fx = float(force[0])
             fy = float(force[1])
             if not np.isfinite(fx) or not np.isfinite(fy):

tests/validation/run_sah04_st_matrix.py

@@ -321,7 +321,7 @@ def run_one_simulation(
             stream.synchronize()
             sim.bodies.download_obs_full_async(stream)
             stream.synchronize()
-            fvec = sim.bodies.read_force(0)
+            fvec = sim.bodies.read_force(0, normalize=False)
             lift = float(fvec[1])
             drag = float(fvec[0])
             if not np.isfinite(lift) or not np.isfinite(drag):

tests/validation/test_sensor_accuracy.py

@@ -65,10 +65,8 @@ def test_sensor_accuracy() -> dict:
     # Get macroscopic field after one more step (with sensor accumulation)
     import pycuda.driver as cuda
     stream = cuda.Stream()
-    sim.bodies.zero_sensor_segment_async(stream)
+    sim.run(1, zero_obs=True, upload_act=False, sync_obs=True, stream=stream)
-    sim.stepper.step(1, action_gpu=sim.bodies.action_gpu,
+    # stream.synchronize() is called inside run()
-                     obs_gpu=sim.bodies.obs_gpu, stream=stream)
-    stream.synchronize()
 
     macro = sim.get_macroscopic()
     ux = macro["ux"]
@@ -76,7 +74,8 @@ def test_sensor_accuracy() -> dict:
 
     results = {}
     all_pass = True
-    for sid in sensor_ids:
+    for idx, sid in enumerate(sensor_ids):
+        pos = positions[idx]
         cells_arr, _ = sim.bodies.get(sid).get_sensor_list(
             sim.lbm_cfg.nx, sim.lbm_cfg.ny
         )
@@ -97,7 +96,7 @@ def test_sensor_accuracy() -> dict:
         if not passed:
             all_pass = False
 
-        results[f"sensor_{sid}_pos{positions[i]}"] = {
+        results[f"sensor_{sid}_pos{pos}"] = {
             "sensor_reading": [sensor_reading_x, sensor_reading_y],
             "manual_average": [sensor_ux_mean, sensor_uy_mean],
             "diff": [float(diff_ux), float(diff_uy)],
@@ -106,7 +105,7 @@ def test_sensor_accuracy() -> dict:
         }
         status = "PASS" if passed else "FAIL"
         print(
-            f"  Sensor {sid} @ {positions[sid]}: "
+            f"  Sensor {sid} @ {pos}: "
             f"reading=({sensor_reading_x:.8f},{sensor_reading_y:.8f}) "
             f"manual=({sensor_ux_mean:.8f},{sensor_uy_mean:.8f}) "
             f"diff=({diff_ux:.2e},{diff_uy:.2e}) "