重构body api，性能分析，项目整理

src/CelerisLab/__init__.py

@@ -7,12 +7,11 @@ Usage::
     from CelerisLab import Simulation
 
     sim = Simulation("configs/config_lbm.json")
-    sim.add_cylinder((100, 50), radius=10)
+    sim.add_body("circle", center=(100, 50), radius=10)
     sim.initialize()
     sim.run(1000)
     data = sim.get_macroscopic()
-
-Legacy FlowField API is still importable from CelerisLab.lbm.driver.
+    force = sim.read_force(0)
 """
 
 __version__ = "0.3.0"

src/CelerisLab/body/__init__.py

@@ -2,7 +2,18 @@
 """
 Body / object management for immersed and rigid objects.
 """
-from .objects import SimObject, Cylinder, Sensor
+from .objects import SimObject, ObjectState, ObjectControl
 from .manager import ObjectManager
+from .registry import BodyRegistry
+from .geometry.base import CutLink, SensorCell, Geometry
+from .geometry.circle import CircleGeometry
+from .coupling.soa_packer import pack_cut_links_to_soa, pack_sensor_to_soa
 
-__all__ = ["SimObject", "Cylinder", "Sensor", "ObjectManager"]
+__all__ = [
+    "SimObject", "ObjectState", "ObjectControl",
+    "ObjectManager",
+    "BodyRegistry",
+    "CutLink", "SensorCell", "Geometry",
+    "CircleGeometry",
+    "pack_cut_links_to_soa", "pack_sensor_to_soa",
+]

src/CelerisLab/body/coupling/__init__.py

@@ -0,0 +1,4 @@
+# CelerisLab/body/coupling/__init__.py
+"""
+Body-fluid coupling: SoA packing, force/torque exchange.
+"""

src/CelerisLab/body/coupling/soa_packer.py

@@ -0,0 +1,132 @@
+# CelerisLab/body/coupling/soa_packer.py
+"""
+SoA packer -- transforms CutLink records into kernel-consumable SoA arrays.
+
+Responsibilities:
+    - Pack CutLink[] -> 8-column SoA for CurvedLinkSoA.assign_host()
+    - Pack SensorCell[] -> 2-column SoA for SensorSoA.assign_host()
+    - Document the host-kernel SoA contract in one concentrated location.
+
+Design:
+    This is the single host-side function that knows the kernel SoA layout.
+    Geometry builders (circle, polygon, mesh) produce CutLink[] and call this
+    function -- they do NOT know the SoA format.  If the kernel SoA layout
+    changes in the future, only this file needs updating.
+"""
+
+from __future__ import annotations
+
+from typing import List, Tuple
+
+import numpy as np
+
+from ..geometry.base import CutLink, SensorCell
+
+
+# Convenience aliases kept for backward compat during refactoring.
+# Geometry callers now produce CutLink[] which is packed here.
+# The constants remain for the fallback classification used during
+# cut-link enumeration in circle.py.
+FALLBACK_CLASS_BOUZIDI = np.uint8(0)
+FALLBACK_CLASS_HALFWAY_BOUNCEBACK = np.uint8(1)
+
+
+def pack_cut_links_to_soa(
+    links: List[CutLink],
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray,
+           np.ndarray, np.ndarray, np.ndarray,
+           np.ndarray, np.ndarray]:
+    """Pack CutLink list into 8-column SoA arrays for CurvedLinkSoA.
+
+    Output order matches CurvedLinkSoA.assign_host() signature::
+
+        fluid_idx (uint32), dir (uint8), q (float32),
+        body_id (int32), rx (float32), ry (float32),
+        rz (float32), fallback (uint8)
+
+    This is the single host-side SoA packing point.
+    All geometry builders (circle, polygon, mesh) call this same
+    function, isolating the kernel SoA layout from geometry logic.
+
+    Semantic contract for each column (reference for host and kernel readers):
+
+        fluid_idx
+            Linear index of the fluid cell whose lattice link crosses into solid.
+            This cell receives the Bouzidi-interpolated population on the next
+            pull step.
+        dir
+            Lattice direction from fluid_idx toward the solid neighbor (1..NQ-1).
+        q
+            Fractional distance fluid -> hit_point / fluid -> neighbor, in (0, 1].
+            Used by the kernel for Bouzidi interpolation weighting.
+        body_id
+            Object id used to look up wall velocity from the action buffer
+            (action[body_id * slot + slot - 1] = omega).
+        rx, ry, rz
+            Lever arm from body center to wall hit point (hit_point - center),
+            used by the kernel for torque accumulation:
+            torque_z = rx * fy - ry * fx (2D).
+        fallback
+            0 = Bouzidi linear interpolation is legal (donor cell is fluid and
+            in-domain).  1 = donor is invalid; use half-way bounce-back
+            fallback [Gin08b].
+
+    Args:
+        links: List of CutLink records from any geometry builder.
+
+    Returns:
+        8-element tuple of numpy arrays matching CurvedLinkSoA.assign_host().
+    """
+    n = len(links)
+    if n == 0:
+        return (
+            np.zeros(0, dtype=np.uint32),
+            np.zeros(0, dtype=np.uint8),
+            np.zeros(0, dtype=np.float32),
+            np.zeros(0, dtype=np.int32),
+            np.zeros(0, dtype=np.float32),
+            np.zeros(0, dtype=np.float32),
+            np.zeros(0, dtype=np.float32),
+            np.zeros(0, dtype=np.uint8),
+        )
+
+    fluid_idx = np.empty(n, dtype=np.uint32)
+    directions = np.empty(n, dtype=np.uint8)
+    q_vals = np.empty(n, dtype=np.float32)
+    body_ids = np.empty(n, dtype=np.int32)
+    rx_vals = np.empty(n, dtype=np.float32)
+    ry_vals = np.empty(n, dtype=np.float32)
+    rz_vals = np.empty(n, dtype=np.float32)
+    fallback_vals = np.empty(n, dtype=np.uint8)
+
+    for i, cl in enumerate(links):
+        fluid_idx[i] = np.uint32(cl.fluid_idx)
+        directions[i] = np.uint8(cl.dir)
+        q_vals[i] = np.float32(cl.q)
+        body_ids[i] = np.int32(cl.body_id)
+        rx_vals[i] = np.float32(cl.rx)
+        ry_vals[i] = np.float32(cl.ry)
+        rz_vals[i] = np.float32(cl.rz)
+        fallback_vals[i] = np.uint8(cl.fallback)
+
+    return (fluid_idx, directions, q_vals, body_ids,
+            rx_vals, ry_vals, rz_vals, fallback_vals)
+
+
+def pack_sensor_to_soa(
+    cells: List[SensorCell],
+) -> Tuple[np.ndarray, np.ndarray]:
+    """Pack SensorCell list into 2-column SoA arrays for SensorSoA.
+
+    Returns:
+        cells (uint32), obj_ids (int32) -- empty arrays if *cells* is empty.
+    """
+    n = len(cells)
+    if n == 0:
+        return np.zeros(0, dtype=np.uint32), np.zeros(0, dtype=np.int32)
+    idx_arr = np.empty(n, dtype=np.uint32)
+    id_arr = np.empty(n, dtype=np.int32)
+    for i, sc in enumerate(cells):
+        idx_arr[i] = np.uint32(sc.idx)
+        id_arr[i] = np.int32(sc.body_id)
+    return idx_arr, id_arr

src/CelerisLab/body/geometry/__init__.py

@@ -0,0 +1,4 @@
+# CelerisLab/body/geometry/__init__.py
+"""
+Geometry shape implementations for immersed / rigid bodies.
+"""

src/CelerisLab/body/geometry/base.py

@@ -0,0 +1,137 @@
+# CelerisLab/body/geometry/base.py
+"""
+Geometry shape interface: uniform intermediate representation for all shapes.
+
+Design:
+    Geometry is an abstract base that defines the interface for shape-specific
+    cut-link enumeration, sensor cell discovery, and flag mask generation.
+    The output is geometry-agnostic CutLink and SensorCell records -- no
+    boundary-method-specific fields leak into the geometry layer.
+"""
+
+from __future__ import annotations
+
+from abc import ABC, abstractmethod
+from dataclasses import dataclass
+from typing import List, Optional, Tuple
+
+import numpy as np
+
+
+@dataclass
+class CutLink:
+    """One cut link: fluid node -> wall -> solid neighbor.
+
+    Geometry-agnostic intermediate record.  Produced by circle/polygon/mesh
+    builders, consumed by the SoA packer.  No Bouzidi-specific fields.
+
+    Fields:
+        fluid_idx:  linear index of the fluid-side cell.
+        dir:        lattice direction towards wall (1..NQ-1).
+        q:          fractional distance along link, in (0, 1].
+        body_id:    owning body id.
+        hit_x:      wall hit-point x (lattice).
+        hit_y:      wall hit-point y.
+        hit_z:      wall hit-point z (0 for 2D).
+        rx:         hit_x - center_x  (torque lever arm).
+        ry:         hit_y - center_y.
+        rz:         hit_z - center_z (0 for 2D).
+        normal_x:   inward-facing surface normal x at hit point (0 for 2D stub).
+        normal_y:   inward-facing surface normal y at hit point (0 for 2D stub).
+        fallback:   0 = Bouzidi legal; 1 = half-way bounce-back required.
+        scheme_tag: reserved for future boundary-scheme selection (0 = Bouzidi).
+        motion_tag: reserved for future motion type (0 = stationary, 1 = rotating, 2 = translating).
+    """
+    fluid_idx: int
+    dir: int
+    q: float
+    body_id: int
+    hit_x: float = 0.0
+    hit_y: float = 0.0
+    hit_z: float = 0.0
+    rx: float = 0.0
+    ry: float = 0.0
+    rz: float = 0.0
+    normal_x: float = 0.0
+    normal_y: float = 0.0
+    fallback: int = 0
+    scheme_tag: int = 0
+    motion_tag: int = 0
+
+
+@dataclass
+class SensorCell:
+    """One sensor sampling cell."""
+    idx: int        # linear cell index
+    body_id: int    # owning body id
+
+
+class Geometry(ABC):
+    """Abstract geometry shape.
+
+    Subclasses implement shape-specific grid projection (cut links, sensor cells,
+    flag masks).  The output is a list of CutLink / SensorCell records --
+    geometry-agnostic, boundary-method-agnostic.
+    """
+
+    @property
+    @abstractmethod
+    def center(self) -> Tuple[float, ...]:
+        ...
+
+    @abstractmethod
+    def build_cut_links(
+        self, nx: int, ny: int,
+        domain_flags: Optional[np.ndarray] = None,
+    ) -> List[CutLink]:
+        """Enumerate cut links for this shape.
+
+        Args:
+            nx, ny: Grid dimensions.
+            domain_flags: Optional uint16 flags array for donor-cell
+                fluid/solid classification. None = all donors assumed fluid.
+
+        Returns:
+            List of CutLink records (empty if zero links found).
+        """
+        ...
+
+    @abstractmethod
+    def build_sensor_cells(self, nx: int, ny: int) -> List[SensorCell]:
+        """Enumerate sensor cells for this shape.
+
+        Args:
+            nx, ny: Grid dimensions.
+
+        Returns:
+            List of SensorCell records (empty if zero cells found).
+        """
+        ...
+
+    @abstractmethod
+    def build_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        """Return (nx*ny,) uint16 flag mask with bits set for this shape.
+
+        This variant marks obstacle cells (``SOLID|OBSTACLE|BC_CURVED``).
+
+        Args:
+            nx, ny: Grid dimensions.
+
+        Returns:
+            uint16 array of shape (nx*ny,) with flag bits set for obstacle cells.
+        """
+        ...
+
+    @abstractmethod
+    def build_sensor_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        """Return (nx*ny,) uint16 flag mask for sensor cells.
+
+        This variant marks sensor cells (``FLUID|SENSOR_FLAG``).
+
+        Args:
+            nx, ny: Grid dimensions.
+
+        Returns:
+            uint16 array of shape (nx*ny,) with flag bits set for sensor cells.
+        """
+        ...

src/CelerisLab/body/geometry/circle.py

@@ -0,0 +1,192 @@
+# CelerisLab/body/geometry/circle.py
+"""
+Circle geometry: cut-link enumeration, sensor cells, flag mask.
+
+2D only.  For 3D, a separate SphereGeometry would be added.
+"""
+
+from __future__ import annotations
+
+import math
+from typing import List, Optional, Tuple
+
+import numpy as np
+
+from .base import CutLink, Geometry, SensorCell
+from ...lbm.descriptors import (
+    FLUID, SOLID, OBSTACLE, BC_CURVED, SENSOR_FLAG,
+    D2Q9_EX, D2Q9_EY,
+)
+
+
+class CircleGeometry(Geometry):
+    """Circular geometry for 2D immersed bodies and sensors."""
+
+    def __init__(self, cx: float, cy: float, radius: float) -> None:
+        self._cx = cx
+        self._cy = cy
+        self._radius = radius
+
+    @property
+    def center(self) -> Tuple[float, float]:
+        return (self._cx, self._cy)
+
+    # ------------------------------------------------------------------
+    # Cut-link enumeration
+    # ------------------------------------------------------------------
+    def build_cut_links(
+        self, nx: int, ny: int,
+        domain_flags: Optional[np.ndarray] = None,
+    ) -> List[CutLink]:
+        """Enumerate cut links for this circle.
+
+        Iterates the AABB bounding box around the circle, walks D2Q9
+        directions from each fluid cell, records links that cross the
+        circle boundary into solid.
+
+        Args:
+            nx, ny: Grid dimensions.
+            domain_flags: Optional uint16 flags array for donor-cell
+                fluid/solid classification. None = all donors assumed fluid.
+
+        Returns:
+            List of CutLink records (empty if zero links found).
+        """
+        from ...common.preprocess import find_circle_intersection
+
+        cx, cy = self._cx, self._cy
+        r = self._radius
+        r_sq = r * r
+
+        # AABB clip
+        x0 = max(1, int(cx - r) - 2)
+        x1 = min(nx - 2, int(cx + r) + 2)
+        y0 = max(1, int(cy - r) - 2)
+        y1 = min(ny - 2, int(cy + r) + 2)
+
+        def _donor_is_fluid(xd: int, yd: int) -> bool:
+            if domain_flags is None:
+                return True
+            idx = xd + yd * nx
+            return (int(domain_flags[idx]) & FLUID) != 0
+
+        def _fallback_class(
+            q_val: float, donor_inside: bool,
+            donor_in_domain: bool, donor_is_fluid: bool,
+        ) -> int:
+            if not (0.0 < q_val <= 1.0):
+                raise ValueError(
+                    f"q must satisfy 0 < q <= 1, got {q_val}")
+            if q_val >= 0.5:
+                return 0
+            if (not donor_in_domain) or donor_inside or (not donor_is_fluid):
+                return 1
+            return 0
+
+        links: List[CutLink] = []
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - cx) ** 2 + (y - cy) ** 2 < r_sq:
+                    continue
+                for i in range(1, 9):  # D2Q9 non-zero directions
+                    xn, yn = x + D2Q9_EX[i], y + D2Q9_EY[i]
+                    if not (0 <= xn < nx and 0 <= yn < ny):
+                        continue
+                    if (xn - cx) ** 2 + (yn - cy) ** 2 >= r_sq:
+                        continue
+                    hit = find_circle_intersection(
+                        float(x), float(y),
+                        float(xn), float(yn),
+                        cx, cy, r,
+                    )
+                    if hit is None:
+                        continue
+                    dx = hit[0] - x
+                    dy = hit[1] - y
+                    norm = math.sqrt(
+                        float(D2Q9_EX[i]) ** 2 + float(D2Q9_EY[i]) ** 2)
+                    q = (math.sqrt(dx * dx + dy * dy) / norm
+                         if norm > 0 else 0.5)
+                    q = max(0.01, min(q, 1.0))
+                    xff = x - D2Q9_EX[i]
+                    yff = y - D2Q9_EY[i]
+                    donor_in_domain = (0 <= xff < nx and 0 <= yff < ny)
+                    donor_inside = False
+                    donor_is_fluid = False
+                    if donor_in_domain:
+                        donor_inside = (
+                            (xff - cx) ** 2 + (yff - cy) ** 2
+                        ) < r_sq
+                        donor_is_fluid = _donor_is_fluid(xff, yff)
+                    links.append(CutLink(
+                        fluid_idx=x + y * nx,
+                        dir=i,
+                        q=float(q),
+                        body_id=0,  # assigned during packing
+                        hit_x=float(hit[0]),
+                        hit_y=float(hit[1]),
+                        hit_z=0.0,
+                        rx=float(hit[0] - cx),
+                        ry=float(hit[1] - cy),
+                        rz=0.0,
+                        normal_x=0.0,
+                        normal_y=0.0,
+                        fallback=_fallback_class(
+                            q, donor_inside, donor_in_domain, donor_is_fluid),
+                        scheme_tag=0,
+                        motion_tag=1,
+                    ))
+
+        return links
+
+    # ------------------------------------------------------------------
+    # Sensor cells
+    # ------------------------------------------------------------------
+    def build_sensor_cells(self, nx: int, ny: int) -> List[SensorCell]:
+        """Enumerate sensor cells inside the circular footprint.
+
+        Returns:
+            List of SensorCell records.
+        """
+        cx, cy = self._cx, self._cy
+        r = self._radius
+        r_sq = r * r
+
+        x0 = max(1, int(cx - r) - 1)
+        x1 = min(nx - 2, int(cx + r) + 1)
+        y0 = max(1, int(cy - r) - 1)
+        y1 = min(ny - 2, int(cy + r) + 1)
+
+        cells: List[SensorCell] = []
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - cx) ** 2 + (y - cy) ** 2 < r_sq:
+                    cells.append(SensorCell(idx=x + y * nx, body_id=0))
+        return cells
+
+    # ------------------------------------------------------------------
+    # Flag mask
+    # ------------------------------------------------------------------
+    def build_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        """Return (nx*ny,) uint16 flag mask with OBSTACLE|SOLID|BC_CURVED."""
+        return self._build_flag_mask_inner(nx, ny, SOLID | OBSTACLE | BC_CURVED)
+
+    def build_sensor_flag_mask(self, nx: int, ny: int) -> np.ndarray:
+        """Return (nx*ny,) uint16 flag mask with FLUID|SENSOR_FLAG."""
+        return self._build_flag_mask_inner(nx, ny, FLUID | SENSOR_FLAG)
+
+    def _build_flag_mask_inner(self, nx: int, ny: int, flag_val: int) -> np.ndarray:
+        n = nx * ny
+        mask = np.zeros(n, dtype=np.uint16)
+        cx, cy = self._cx, self._cy
+        r = self._radius
+        r_sq = r * r
+        x0 = max(1, int(cx - r) - 1)
+        x1 = min(nx - 2, int(cx + r) + 1)
+        y0 = max(1, int(cy - r) - 1)
+        y1 = min(ny - 2, int(cy + r) + 1)
+        for x in range(x0, x1 + 1):
+            for y in range(y0, y1 + 1):
+                if (x - cx) ** 2 + (y - cy) ** 2 < r_sq:
+                    mask[x + y * nx] = np.uint16(flag_val)
+        return mask

src/CelerisLab/body/manager.py

@@ -1,9 +1,9 @@
 # CelerisLab/body/manager.py
 """
-ObjectManager — batch management of SimObjects.
+ObjectManager -- batch management of SimObjects.
 
 Responsibilities:
-  - Add / remove / query objects
+  - Add / remove / query objects (delegated to ``BodyRegistry``)
   - Build merged flag masks and compact cut-link / sensor lists
   - Allocate packed telemetry ``obs_gpu`` + pagelocked mirror ``obs_pinned``
   - Sync geometry to :class:`~CelerisLab.lbm.field.LBMField`
@@ -20,6 +20,8 @@ import pycuda.driver as cuda
 from typing import Dict, List, Optional
 
 from .objects import SimObject
+from .registry import BodyRegistry
+from .preprocess.flag_overlay import merge_flag_masks
 from ..cuda.compiler_v2 import obs_layout
 from ..cuda.obs_layout import ObsLayout
 from ..lbm.descriptors import (
@@ -46,8 +48,7 @@ class ObjectManager:
         self.nz = nz
         self.nq = nq
         self.cfg = cfg
-        self._objects: Dict[int, SimObject] = {}
-        self._next_id = 0
+        self._registry = BodyRegistry()
 
         self.action_gpu: Optional[cuda.DeviceAllocation] = None
         self.obs_gpu: Optional[cuda.DeviceAllocation] = None
@@ -73,29 +74,29 @@ class ObjectManager:
         # when bodies.count == 0 (matches _resize_buffers sizing for n=0).
         self._resize_buffers()
 
+    # -- Registry delegation -------------------------------------------------
     def add(self, obj: SimObject) -> int:
         """Register an object and return its id."""
-        obj.obj_id = self._next_id
-        self._objects[self._next_id] = obj
-        self._next_id += 1
+        body_id = self._registry.add(obj)
         self._resize_buffers()
-        return obj.obj_id
+        return body_id
 
     def remove(self, obj_id: int):
-        del self._objects[obj_id]
+        self._registry.remove(obj_id)
         self._resize_buffers()
 
     def get(self, obj_id: int) -> SimObject:
-        return self._objects[obj_id]
+        return self._registry.get(obj_id)
 
     @property
     def objects(self) -> List[SimObject]:
-        return list(self._objects.values())
+        return self._registry.objects
 
     @property
     def count(self) -> int:
-        return len(self._objects)
+        return self._registry.count
 
+    # -- Buffer sizing -------------------------------------------------------
     def _resize_buffers(self):
         n = self.count
         dim = self.cfg.dim if self.cfg else 2
@@ -112,6 +113,7 @@ class ObjectManager:
             new_sensor_counts[:copy_n] = self.sensor_cell_counts[:copy_n]
         self.sensor_cell_counts = new_sensor_counts
 
+    # -- Flag composition ----------------------------------------------------
     def build_flags(self, base_flags: np.ndarray) -> np.ndarray:
         """Merge object flag masks onto a clean domain base.
 
@@ -119,13 +121,13 @@ class ObjectManager:
         overlays stateless and prevents stale obstacle bits from surviving after
         geometry edits or future body motion.
         """
-        merged = np.array(base_flags, copy=True)
-        for obj in self._objects.values():
-            mask = obj.get_flag_mask(self.nx, self.ny, self.nz)
-            nonzero = mask != 0
-            merged[nonzero] = mask[nonzero]
-        return merged
+        masks = [
+            obj.get_flag_mask(self.nx, self.ny, self.nz)
+            for obj in self._registry.objects
+        ]
+        return merge_flag_masks(base_flags, masks)
 
+    # -- Compact list building -----------------------------------------------
     def build_compact_lists(self, domain_flags: np.ndarray | None = None):
         """Build cut-link SoA columns and sensor lists.
 
@@ -145,8 +147,8 @@ class ObjectManager:
         else:
             ex, ey = D2Q9_EX, D2Q9_EY
 
-        for obj in self._objects.values():
-            if hasattr(obj, 'get_curved_list'):
+        for obj in self._registry.objects:
+            if not obj.is_sensor and hasattr(obj, 'get_curved_list'):
                 (
                     fluid_idx, dirs, q_vals, body_ids,
                     rx_vals, ry_vals, rz_vals, fallback_vals
@@ -164,7 +166,7 @@ class ObjectManager:
                     cl_ry.append(ry_vals)
                     cl_rz.append(rz_vals)
                     cl_fallback.append(fallback_vals)
-            if hasattr(obj, 'get_sensor_list'):
+            if obj.is_sensor and hasattr(obj, 'get_sensor_list'):
                 cells, ids = obj.get_sensor_list(self.nx, self.ny, self.nz)
                 if len(cells) > 0:
                     s_cells.append(cells)
@@ -192,6 +194,7 @@ class ObjectManager:
             sensor_cells, sensor_obj_id,
         )
 
+    # -- GPU sync ------------------------------------------------------------
     def sync_to_gpu(self, field, *, rebuild_flags: bool = True):
         """Upload compact lists, action buffer, and packed ``obs`` GPU buffer.
 
@@ -301,7 +304,7 @@ class ObjectManager:
         cuda.memset_d32_async(ptr, 0, self.slot_stride_floats, stream)
 
     def download_obs_full_async(self, stream: cuda.Stream):
-        """Enqueue full DTOH copy ``obs_gpu`` → ``obs_pinned``."""
+        """Enqueue full DTOH copy ``obs_gpu`` -> ``obs_pinned``."""
         assert self.obs_pinned is not None
         cuda.memcpy_dtoh_async(self.obs_pinned, self.obs_gpu, stream)
 
@@ -376,7 +379,8 @@ class ObjectManager:
 
     def _validate_body_id(self, body_id: int) -> None:
         if body_id < 0 or body_id >= self.count:
-            raise IndexError(f"body_id {body_id} out of range for {self.count} objects")
+            raise IndexError(
+                f"body_id {body_id} out of range for {self.count} objects")
 
     def _refresh_action_from_objects(self) -> None:
         """Populate action slots from object state using the runtime contract."""
@@ -385,11 +389,12 @@ class ObjectManager:
         dim = self.cfg.dim if self.cfg else 2
         slot = 3 * dim
         self.action.fill(0)
-        for body_id, obj in self._objects.items():
-            base = body_id * slot
+        for obj in self._registry.objects:
+            base = obj.obj_id * slot
             if base + slot > self.action.size:
                 continue
-            self.action[base + slot - 1] = np.float32(getattr(obj.state, "omega", 0.0))
+            self.action[base + slot - 1] = np.float32(
+                getattr(obj.state, "omega", 0.0))
 
     def fallback_link_count(self) -> int:
         """Return number of curved links marked for fallback bounce-back."""
@@ -407,8 +412,10 @@ class ObjectManager:
         curved = self._telemetry_field.curved
         if curved.count == 0:
             return 0
-        return int(np.count_nonzero(curved.q[:curved.count] < np.float32(threshold)))
+        return int(np.count_nonzero(
+            curved.q[:curved.count] < np.float32(threshold)))
 
+    # -- Placeholders for future use -----------------------------------------
     def update_states(self, dt: float):
         """Integrate object motion (placeholder)."""
         pass

src/CelerisLab/body/objects.py

@@ -1,28 +1,30 @@
 # CelerisLab/body/objects.py
 """
-Lightweight flat object model for immersed / rigid bodies and sensors.
+SimObject -- container wrapping geometry, state, and control.
 
 Design:
-  - SimObject is a thin base with standard interface (state, control,
-    get_flag_mask, get_curved_list).
-  - Concrete types (Cylinder, Sensor, …) override geometry methods.
-  - No deep inheritance tree.  Users can subclass SimObject directly
-    for custom shapes — just implement get_flag_mask / get_curved_list.
+    SimObject is a thin container that holds a Geometry instance (shape),
+    an ObjectState (position/velocity), and an ObjectControl (control input).
+    Concrete geometry logic lives in ``body/geometry/*.py``.
+
+    The ``get_curved_list()`` / ``get_sensor_list()`` methods exist for
+    backward compatibility with ``ObjectManager.build_compact_lists()``.
+    They delegate to the Geometry's build methods and the SoA packer.
 """
 
-import math
-import numpy as np
+from __future__ import annotations
+
 from dataclasses import dataclass
-from typing import Tuple
+from typing import Optional, Tuple
 
-# Lattice constants — canonical Python source is lbm/descriptors.py
-from ..lbm.descriptors import (
-    FLUID, SOLID, GAS, INTERFACE,
-    BC_WALL, BC_INLET, BC_OUTLET, BC_CURVED, BC_PERIODIC, BC_MOVING,
-    SENSOR_FLAG, OBSTACLE,
-    D2Q9_EX, D2Q9_EY,
-)
+import numpy as np
 
+from .geometry.base import Geometry
+
+
+# ---------------------------------------------------------------------------
+# Constants
+# ---------------------------------------------------------------------------
 FALLBACK_CLASS_BOUZIDI = np.uint8(0)
 FALLBACK_CLASS_HALFWAY_BOUNCEBACK = np.uint8(1)
 
@@ -32,7 +34,7 @@ FALLBACK_CLASS_HALFWAY_BOUNCEBACK = np.uint8(1)
 # ---------------------------------------------------------------------------
 @dataclass
 class ObjectState:
-    """Position, velocity, orientation (z used for 3D layouts)."""
+    """Position, velocity, orientation (z reserved for 3D layouts)."""
     x: float = 0.0
     y: float = 0.0
     z: float = 0.0
@@ -52,305 +54,99 @@ class ObjectControl:
 
 
 # ---------------------------------------------------------------------------
-# Base class
+# SimObject container
 # ---------------------------------------------------------------------------
 class SimObject:
-    """Base for all simulation objects."""
+    """Container for a simulation object: geometry, state, control.
+
+    The Geometry instance handles shape-specific operations (cut-link
+    enumeration, flag mask generation).  SimObject wraps it with runtime
+    state and an id for bookkeeping.
+
+    The *is_sensor* flag selects the flag mask variant for sensors
+    (``FLUID|SENSOR_FLAG``) vs. obstacle bodies (``SOLID|OBSTACLE|BC_CURVED``).
+    """
 
     obj_type: str = "generic"
 
-    def __init__(self, obj_id: int, center: Tuple[float, ...],
-                 radius: float = 0.0):
+    def __init__(
+        self,
+        obj_id: int,
+        geometry: Optional[Geometry] = None,
+        center: Optional[Tuple[float, ...]] = None,
+        radius: float = 0.0,
+        is_sensor: bool = False,
+    ) -> None:
         self.obj_id = obj_id
-        self.center = center
+        self.center = center if center is not None else (0.0, 0.0)
         self.radius = radius
-        z0 = float(center[2]) if len(center) > 2 else 0.0
-        self.state = ObjectState(x=float(center[0]), y=float(center[1]), z=z0)
+        self._geometry = geometry
+        self._is_sensor = is_sensor
+
+        c = self.center
+        z0 = float(c[2]) if len(c) > 2 else 0.0
+        self.state = ObjectState(
+            x=float(c[0]), y=float(c[1]), z=z0,
+        )
         self.control = ObjectControl()
 
+    @property
+    def geometry(self) -> Geometry:
+        if self._geometry is None:
+            raise RuntimeError(
+                f"SimObject id={self.obj_id} has no geometry assigned.")
+        return self._geometry
+
+    @geometry.setter
+    def geometry(self, g: Geometry) -> None:
+        self._geometry = g
+
+    @property
+    def is_sensor(self) -> bool:
+        return self._is_sensor
+
+    # ------------------------------------------------------------------
+    # Methods consumed by ObjectManager.build_compact_lists() / build_flags()
+    # ------------------------------------------------------------------
     def get_flag_mask(self, nx: int, ny: int, nz: int = 1) -> np.ndarray:
-        """Return (n,) uint16 array with flag bits set for this object."""
-        raise NotImplementedError
+        """Return (nx*ny,) uint16 array with flag bits set for this object.
 
+        Sensors produce ``FLUID|SENSOR_FLAG``; bodies produce
+        ``SOLID|OBSTACLE|BC_CURVED``.
+        """
+        if self._is_sensor:
+            return self.geometry.build_sensor_flag_mask(nx, ny)
+        return self.geometry.build_flag_mask(nx, ny)
 
+    def get_curved_list(
+        self, nx: int, ny: int, nq: int,
+        nz: int = 1, ex=None, ey=None, ez=None,
+        domain_flags=None,
+    ):
+        """Return 8-column SoA arrays for curved-boundary kernel.
 
-
-# ---------------------------------------------------------------------------
-# Concrete types
-# ---------------------------------------------------------------------------
-class Cylinder(SimObject):
-    obj_type = "cylinder"
-
-    def get_flag_mask(self, nx: int, ny: int, nz: int = 1) -> np.ndarray:
-        n = nx * ny * nz
-        mask = np.zeros(n, dtype=np.uint16)
-        xc, yc, zc = self.state.x, self.state.y, self.state.z
-        r = self.radius
-        r_sq = r * r
-        x0 = max(1, int(xc - r) - 1)
-        x1 = min(nx - 2, int(xc + r) + 1)
-        y0 = max(1, int(yc - r) - 1)
-        y1 = min(ny - 2, int(yc + r) + 1)
-        z0 = max(1, int(zc - r) - 1) if nz > 1 else 0
-        z1 = min(nz - 2, int(zc + r) + 1) if nz > 1 else 0
-        if nz == 1:
-            for x in range(x0, x1 + 1):
-                for y in range(y0, y1 + 1):
-                    if (x - xc) ** 2 + (y - yc) ** 2 < r_sq:
-                        mask[x + y * nx] = SOLID | OBSTACLE | BC_CURVED
-        else:
-            for z in range(z0, z1 + 1):
-                for y in range(y0, y1 + 1):
-                    for x in range(x0, x1 + 1):
-                        if (x - xc) ** 2 + (y - yc) ** 2 + (z - zc) ** 2 < r_sq:
-                            k = x + y * nx + z * nx * ny
-                            mask[k] = SOLID | OBSTACLE | BC_CURVED
-        return mask
-
-    def get_curved_list(self, nx: int, ny: int, nq: int,
-                        nz: int = 1, ex=None, ey=None, ez=None,
-                        domain_flags=None):
-        """Return per-link curved boundary data (fluid-centric).
-
-        Each record represents one cut link: a fluid node with a lattice
-        direction that crosses the cylinder surface into solid.
-
-        Args:
-            nx, ny, nq: grid and lattice dimensions.
-            nz: z-depth; use 1 for 2D layouts (linear index x + y·nx).
-            ex, ey, ez: lattice vectors (length *nq*). *ez* required for a full
-                D3Q19 fluid→solid link walk when ``nz > 1``.
+        Delegates to the geometry object for cut-link enumeration and
+        to ``pack_cut_links_to_soa`` for SoA packing.
 
         Returns:
-            fluid_idx:  uint32  [n_links] fluid cell linear indices
-            directions: uint8   [n_links] direction from fluid toward wall
-            q_values:   float32 [n_links] fractional distance ∈ (0, 1]
-            body_ids:   int32   [n_links] object id for each link
-            rx_values:  float32 [n_links] x arm from center to wall hit point
-            ry_values:  float32 [n_links] y arm from center to wall hit point
-            rz_values:  float32 [n_links] z arm (zero for 2D)
-            fallback:   uint8   [n_links]
-                0=normal Bouzidi interpolation,
-                1=half-way moving bounce-back fallback
+            8-tuple matching CurvedLinkSoA.assign_host().
         """
-        from ..common.preprocess import (
-            find_circle_intersection,
-            find_sphere_ray_segment,
-        )
-        xc, yc, zc = self.state.x, self.state.y, self.state.z
-        r = self.radius
-        if ex is None or ey is None:
-            ex, ey = D2Q9_EX, D2Q9_EY
+        from .coupling.soa_packer import pack_cut_links_to_soa
 
-        idx_list = []
-        dir_list = []
-        q_list = []
-        rx_list = []
-        ry_list = []
-        rz_list = []
-        fallback_list = []
-
-        x0 = max(1, int(xc - r) - 2)
-        x1 = min(nx - 2, int(xc + r) + 2)
-        y0 = max(1, int(yc - r) - 2)
-        y1 = min(ny - 2, int(yc + r) + 2)
-        z0 = max(1, int(zc - r) - 2) if nz > 1 else 0
-        z1 = min(nz - 2, int(zc + r) + 2) if nz > 1 else 0
-
-        r_sq = r * r
-        use_3d = ez is not None
-
-        def _donor_is_fluid(xd: int, yd: int, zd: int = 0) -> bool:
-            if domain_flags is None:
-                return True
-            if nz > 1:
-                idx = xd + yd * nx + zd * nx * ny
-            else:
-                idx = xd + yd * nx
-            return (int(domain_flags[idx]) & FLUID) != 0
-
-        def _fallback_class(
-            q_val: float, donor_inside: bool, donor_in_domain: bool, donor_is_fluid: bool
-        ) -> np.uint8:
-            """Classify donor legality for Bouzidi interpolation fallback."""
-            if not (0.0 < q_val <= 1.0):
-                raise ValueError(f"q must satisfy 0 < q <= 1, got {q_val}")
-            if q_val >= 0.5:
-                return FALLBACK_CLASS_BOUZIDI
-            if (not donor_in_domain) or donor_inside or (not donor_is_fluid):
-                return FALLBACK_CLASS_HALFWAY_BOUNCEBACK
-            return FALLBACK_CLASS_BOUZIDI
-
-        if use_3d:
-            for z in range(z0, z1 + 1):
-                for y in range(y0, y1 + 1):
-                    for x in range(x0, x1 + 1):
-                        if (x - xc) ** 2 + (y - yc) ** 2 + (z - zc) ** 2 < r_sq:
-                            continue
-                        for i in range(1, nq):
-                            xn = x + ex[i]
-                            yn = y + ey[i]
-                            zn = z + ez[i]
-                            if not (0 <= xn < nx and 0 <= yn < ny
-                                    and 0 <= zn < nz):
-                                continue
-                            if (xn - xc) ** 2 + (yn - yc) ** 2 + (
-                                    zn - zc) ** 2 >= r_sq:
-                                continue
-                            hit = find_sphere_ray_segment(
-                                float(x), float(y), float(z),
-                                float(xn), float(yn), float(zn),
-                                xc, yc, zc, r,
-                            )
-                            if hit is None:
-                                continue
-                            dx = hit[0] - x
-                            dy = hit[1] - y
-                            dz = hit[2] - z
-                            norm = math.sqrt(
-                                ex[i] ** 2 + ey[i] ** 2 + ez[i] ** 2)
-                            q = (math.sqrt(dx * dx + dy * dy + dz * dz) / norm
-                                 if norm > 0 else 0.5)
-                            q = max(0.01, min(q, 1.0))
-                            xff = x - ex[i]
-                            yff = y - ey[i]
-                            zff = z - ez[i]
-                            donor_in_domain = (0 <= xff < nx and 0 <= yff < ny and 0 <= zff < nz)
-                            donor_inside = False
-                            donor_is_fluid = False
-                            if donor_in_domain:
-                                donor_inside = (
-                                    (xff - xc) ** 2 + (yff - yc) ** 2 + (zff - zc) ** 2
-                                ) < r_sq
-                                donor_is_fluid = _donor_is_fluid(xff, yff, zff)
-                            idx_list.append(
-                                np.uint32(x + y * nx + z * nx * ny))
-                            dir_list.append(np.uint8(i))
-                            q_list.append(np.float32(q))
-                            rx_list.append(np.float32(hit[0] - xc))
-                            ry_list.append(np.float32(hit[1] - yc))
-                            rz_list.append(np.float32(hit[2] - zc))
-                            fallback_list.append(
-                                _fallback_class(q, donor_inside, donor_in_domain, donor_is_fluid))
-        else:
-            for x in range(x0, x1 + 1):
-                for y in range(y0, y1 + 1):
-                    if (x - xc) ** 2 + (y - yc) ** 2 < r_sq:
-                        continue
-                    for i in range(1, nq):
-                        xn, yn = x + ex[i], y + ey[i]
-                        if not (0 <= xn < nx and 0 <= yn < ny):
-                            continue
-                        if (xn - xc) ** 2 + (yn - yc) ** 2 >= r_sq:
-                            continue
-                        hit = find_circle_intersection(
-                            float(x), float(y),
-                            float(xn), float(yn),
-                            xc, yc, r,
-                        )
-                        if hit is not None:
-                            dx = hit[0] - x
-                            dy = hit[1] - y
-                            norm = math.sqrt(ex[i] ** 2 + ey[i] ** 2)
-                            q = math.sqrt(dx * dx + dy * dy) / norm if norm > 0 else 0.5
-                            q = max(0.01, min(q, 1.0))
-                            xff = x - ex[i]
-                            yff = y - ey[i]
-                            donor_in_domain = (0 <= xff < nx and 0 <= yff < ny)
-                            donor_inside = False
-                            donor_is_fluid = False
-                            if donor_in_domain:
-                                donor_inside = ((xff - xc) ** 2 + (yff - yc) ** 2) < r_sq
-                                donor_is_fluid = _donor_is_fluid(xff, yff)
-                            idx_list.append(np.uint32(x + y * nx))
-                            dir_list.append(np.uint8(i))
-                            q_list.append(np.float32(q))
-                            rx_list.append(np.float32(hit[0] - xc))
-                            ry_list.append(np.float32(hit[1] - yc))
-                            rz_list.append(np.float32(0.0))
-                            fallback_list.append(
-                                _fallback_class(q, donor_inside, donor_in_domain, donor_is_fluid))
-
-        n = len(idx_list)
-        if n > 0:
-            return (np.array(idx_list, dtype=np.uint32),
-                    np.array(dir_list, dtype=np.uint8),
-                    np.array(q_list, dtype=np.float32),
-                    np.full(n, self.obj_id, dtype=np.int32),
-                    np.array(rx_list, dtype=np.float32),
-                    np.array(ry_list, dtype=np.float32),
-                    np.array(rz_list, dtype=np.float32),
-                    np.array(fallback_list, dtype=np.uint8))
-        return (np.zeros(0, dtype=np.uint32),
-                np.zeros(0, dtype=np.uint8),
-                np.zeros(0, dtype=np.float32),
-                np.zeros(0, dtype=np.int32),
-                np.zeros(0, dtype=np.float32),
-                np.zeros(0, dtype=np.float32),
-                np.zeros(0, dtype=np.float32),
-                np.zeros(0, dtype=np.uint8))
-
-
-class Sensor(SimObject):
-    obj_type = "sensor"
-
-    def get_flag_mask(self, nx: int, ny: int, nz: int = 1) -> np.ndarray:
-        n = nx * ny * nz
-        mask = np.zeros(n, dtype=np.uint16)
-        xc, yc, zc = self.state.x, self.state.y, self.state.z
-        r = self.radius
-        r_sq = r * r
-        x0 = max(1, int(xc - r) - 1)
-        x1 = min(nx - 2, int(xc + r) + 1)
-        y0 = max(1, int(yc - r) - 1)
-        y1 = min(ny - 2, int(yc + r) + 1)
-        z0 = max(1, int(zc - r) - 1) if nz > 1 else 0
-        z1 = min(nz - 2, int(zc + r) + 1) if nz > 1 else 0
-        if nz == 1:
-            for x in range(x0, x1 + 1):
-                for y in range(y0, y1 + 1):
-                    if (x - xc) ** 2 + (y - yc) ** 2 < r_sq:
-                        mask[x + y * nx] = FLUID | SENSOR_FLAG
-        else:
-            for z in range(z0, z1 + 1):
-                for y in range(y0, y1 + 1):
-                    for x in range(x0, x1 + 1):
-                        if (x - xc) ** 2 + (y - yc) ** 2 + (z - zc) ** 2 < r_sq:
-                            k = x + y * nx + z * nx * ny
-                            mask[k] = FLUID | SENSOR_FLAG
-        return mask
+        links = self.geometry.build_cut_links(nx, ny, domain_flags=domain_flags)
+        for cl in links:
+            cl.body_id = self.obj_id
+        return pack_cut_links_to_soa(links)
 
     def get_sensor_list(self, nx: int, ny: int, nz: int = 1):
-        """Return compact sensor list.
+        """Return compact sensor list (cells, obj_ids).
 
         Returns:
-            cells:   uint32 array of cell indices
-            obj_ids: int32 array of object ids
+            cells (uint32), obj_ids (int32).
         """
-        cell_list = []
-        xc, yc, zc = self.state.x, self.state.y, self.state.z
-        r = self.radius
-        r_sq = r * r
-        x0 = max(1, int(xc - r) - 1)
-        x1 = min(nx - 2, int(xc + r) + 1)
-        y0 = max(1, int(yc - r) - 1)
-        y1 = min(ny - 2, int(yc + r) + 1)
-        z0 = max(1, int(zc - r) - 1) if nz > 1 else 0
-        z1 = min(nz - 2, int(zc + r) + 1) if nz > 1 else 0
-        if nz == 1:
-            for x in range(x0, x1 + 1):
-                for y in range(y0, y1 + 1):
-                    if (x - xc) ** 2 + (y - yc) ** 2 < r_sq:
-                        cell_list.append(np.uint32(x + y * nx))
-        else:
-            for z in range(z0, z1 + 1):
-                for y in range(y0, y1 + 1):
-                    for x in range(x0, x1 + 1):
-                        if (x - xc) ** 2 + (y - yc) ** 2 + (z - zc) ** 2 < r_sq:
-                            cell_list.append(
-                                np.uint32(x + y * nx + z * nx * ny))
-        if cell_list:
-            cells = np.array(cell_list, dtype=np.uint32)
-            obj_ids = np.full(len(cell_list), self.obj_id, dtype=np.int32)
-            return cells, obj_ids
-        return np.zeros(0, dtype=np.uint32), np.zeros(0, dtype=np.int32)
+        from .coupling.soa_packer import pack_sensor_to_soa
+
+        cells = self.geometry.build_sensor_cells(nx, ny)
+        for sc in cells:
+            sc.body_id = self.obj_id
+        return pack_sensor_to_soa(cells)

src/CelerisLab/body/preprocess/__init__.py

@@ -0,0 +1,4 @@
+# CelerisLab/body/preprocess/__init__.py
+"""
+Grid preprocessing: flag overlay, cut-link packing, sensor mapping.
+"""

src/CelerisLab/body/preprocess/flag_overlay.py

@@ -0,0 +1,34 @@
+# CelerisLab/body/preprocess/flag_overlay.py
+"""
+Flag overlay: merge per-object flag masks onto a clean domain base.
+"""
+
+from __future__ import annotations
+
+from typing import List
+
+import numpy as np
+
+
+def merge_flag_masks(
+    base: np.ndarray,
+    masks: List[np.ndarray],
+) -> np.ndarray:
+    """Merge object flag masks onto a clean domain base.
+
+    For each mask in *masks*, non-zero entries overwrite the corresponding
+    positions in *base*.  Callers should pass a fresh channel-layout flag
+    array to prevent stale obstacle bits from surviving after geometry edits.
+
+    Args:
+        base: uint16 domain flag array (e.g. from build_channel_flags()).
+        masks: List of uint16 mask arrays from each object's get_flag_mask().
+
+    Returns:
+        New uint16 array with object overlays applied.
+    """
+    merged = np.array(base, copy=True)
+    for mask in masks:
+        nonzero = mask != 0
+        merged[nonzero] = mask[nonzero]
+    return merged

src/CelerisLab/body/registry.py

@@ -0,0 +1,49 @@
+# CelerisLab/body/registry.py
+"""
+BodyRegistry -- pure add/remove/query for simulation objects.
+No GPU or LBM knowledge.
+"""
+
+from __future__ import annotations
+
+from typing import Dict, List
+
+from .objects import SimObject
+
+
+class BodyRegistry:
+    """Registry for SimObject instances.  No GPU or LBM knowledge."""
+
+    def __init__(self) -> None:
+        self._objects: Dict[int, SimObject] = {}
+        self._next_id = 0
+
+    def add(self, obj: SimObject) -> int:
+        """Register a new object and return its assigned id."""
+        obj.obj_id = self._next_id
+        self._objects[self._next_id] = obj
+        self._next_id += 1
+        return obj.obj_id
+
+    def remove(self, obj_id: int) -> None:
+        """Unregister an object by id."""
+        del self._objects[obj_id]
+
+    def get(self, obj_id: int) -> SimObject:
+        """Return the registered object by id."""
+        return self._objects[obj_id]
+
+    def clear(self) -> None:
+        """Remove all objects and reset id counter."""
+        self._objects.clear()
+        self._next_id = 0
+
+    @property
+    def objects(self) -> List[SimObject]:
+        """List of all registered objects in insertion order."""
+        return list(self._objects.values())
+
+    @property
+    def count(self) -> int:
+        """Number of registered objects."""
+        return len(self._objects)

src/CelerisLab/common/__init__.py

@@ -3,7 +3,6 @@
 Common utilities and preprocessing functions.
 """
 
-from . import preprocess
-from . import streakline
+from . import preprocess, pathline, checkpoint, render, streakline
 
-__all__ = ["preprocess", "streakline"]
+__all__ = ["preprocess", "pathline", "checkpoint", "render", "streakline"]

src/CelerisLab/common/_types.py

@@ -0,0 +1,9 @@
+# CelerisLab/common/_types.py
+"""
+Shared type aliases for the ``common`` package.
+"""
+
+from typing import Tuple
+
+CylinderSpec = Tuple[Tuple[float, float], float]
+"""A cylinder defined by ``((cx, cy), radius)`` in lattice units."""

src/CelerisLab/common/pathline.py

@@ -0,0 +1,207 @@
+# CelerisLab/common/pathline.py
+"""
+Particle pathline history accumulation and rendering (debug tool).
+
+This is **not** a streakline — it accumulates each particle's full trajectory
+over time.  For streakline visualization use ``streakline.Streakline``.
+"""
+
+from __future__ import annotations
+
+from concurrent.futures import ThreadPoolExecutor
+from dataclasses import dataclass, field
+from typing import List, Optional, Sequence, Tuple
+
+import numpy as np
+
+from .streakline._render import gaussian_blur2d
+from ._types import CylinderSpec
+
+
+def configure_compute_threads(num_threads: int = 0) -> int:
+    """Set BLAS/OpenMP thread env vars for host-side numpy work."""
+    import os
+
+    if num_threads <= 0:
+        num_threads = min(32, os.cpu_count() or 4)
+    for var in (
+        "OMP_NUM_THREADS",
+        "OPENBLAS_NUM_THREADS",
+        "MKL_NUM_THREADS",
+        "NUMEXPR_NUM_THREADS",
+        "VECLIB_MAXIMUM_THREADS",
+    ):
+        os.environ[var] = str(int(num_threads))
+    return int(num_threads)
+
+
+@dataclass
+class ParticleTrailSet:
+    """Growing polylines per particle (pathline history, not streakline)."""
+
+    trails: List[np.ndarray] = field(default_factory=list)
+    active: List[bool] = field(default_factory=list)
+
+    def inject(self, points: np.ndarray) -> None:
+        pts = np.asarray(points, dtype=np.float64)
+        if pts.ndim != 2 or pts.shape[1] != 2:
+            raise ValueError("points must be [N,2].")
+        for p in pts:
+            self.trails.append(np.asarray([[p[0], p[1]]], dtype=np.float64))
+            self.active.append(True)
+
+    def keep_indices(self, survived_mask: np.ndarray) -> None:
+        mask = np.asarray(survived_mask, dtype=bool).ravel()
+        j = 0
+        for i, is_active in enumerate(self.active):
+            if not is_active:
+                continue
+            if j >= mask.size:
+                break
+            if not bool(mask[j]):
+                self.active[i] = False
+            j += 1
+
+    def append_positions(self, particles: np.ndarray) -> None:
+        pts = np.asarray(particles, dtype=np.float64)
+        j = 0
+        for i, is_active in enumerate(self.active):
+            if not is_active:
+                continue
+            if j >= pts.shape[0]:
+                break
+            self.trails[i] = np.vstack([self.trails[i], pts[j]])
+            j += 1
+
+
+def _rasterize_segment(
+    intensity: np.ndarray,
+    x0: float,
+    y0: float,
+    x1: float,
+    y1: float,
+    weight: float,
+) -> None:
+    """Stamp one segment into a scalar intensity map."""
+    ny, nx = intensity.shape
+    n = int(max(abs(x1 - x0), abs(y1 - y0), 1.0)) + 1
+    xs = np.linspace(x0, x1, n)
+    ys = np.linspace(y0, y1, n)
+    xi = np.clip(np.rint(xs).astype(np.int64), 0, nx - 1)
+    yi = np.clip(np.rint(ys).astype(np.int64), 0, ny - 1)
+    w = float(max(0.0, min(1.0, weight)))
+    np.add.at(intensity, (yi, xi), w)
+
+
+def _draw_trail(
+    intensity: np.ndarray,
+    trail: np.ndarray,
+    *,
+    fade_along_trail: bool,
+) -> None:
+    if trail.shape[0] < 2:
+        return
+    nseg = trail.shape[0] - 1
+    for i in range(nseg):
+        if fade_along_trail:
+            t = float(i + 1) / float(nseg)
+            weight = 0.15 + 0.85 * t
+        else:
+            weight = 1.0
+        _rasterize_segment(
+            intensity,
+            float(trail[i, 0]),
+            float(trail[i, 1]),
+            float(trail[i + 1, 0]),
+            float(trail[i + 1, 1]),
+            weight,
+        )
+
+
+def render_trails(
+    trail_set: ParticleTrailSet,
+    *,
+    nx: int,
+    ny: int,
+    out_path: str,
+    cylinders: Optional[Sequence[CylinderSpec]] = None,
+    blur_sigma: float = 0.0,
+    fade_along_trail: bool = True,
+    num_threads: int = 0,
+    background_color: Tuple[float, float, float] = (1.0, 1.0, 1.0),
+    streak_color: Tuple[float, float, float] = (1.0, 0.0, 0.0),
+) -> dict:
+    """Render accumulated path trails (pathline collage — debug only).
+
+    This is NOT a streakline renderer.  See ``streakline.Streakline``.
+    """
+    n_threads = configure_compute_threads(num_threads)
+    intensity = np.zeros((ny, nx), dtype=np.float64)
+
+    trails = trail_set.trails
+    if not trails:
+        density_max = 0.0
+    else:
+        chunk = max(1, (len(trails) + n_threads - 1) // n_threads)
+
+        def _worker(batch: List[np.ndarray]) -> np.ndarray:
+            local = np.zeros((ny, nx), dtype=np.float64)
+            for tr in batch:
+                _draw_trail(local, tr, fade_along_trail=fade_along_trail)
+            return local
+
+        batches = [trails[i : i + chunk] for i in range(0, len(trails), chunk)]
+        if len(batches) == 1:
+            intensity = _worker(batches[0])
+        else:
+            with ThreadPoolExecutor(max_workers=n_threads) as pool:
+                parts = list(pool.map(_worker, batches))
+            for part in parts:
+                intensity = np.maximum(intensity, part)
+
+        if np.any(intensity > 0):
+            vmax = float(np.percentile(intensity[intensity > 0], 99.0))
+            intensity = np.clip(intensity / max(vmax, 1e-12), 0.0, 1.0)
+        density_max = float(np.max(intensity))
+
+    rgb = np.ones((ny, nx, 3), dtype=np.float64)
+    rgb[:, :, 0] = background_color[0]
+    rgb[:, :, 1] = background_color[1] * (1.0 - intensity) + streak_color[1] * intensity
+    rgb[:, :, 2] = background_color[2] * (1.0 - intensity) + streak_color[2] * intensity
+    if streak_color[0] < 0.999:
+        rgb[:, :, 0] = background_color[0] * (1.0 - intensity) + streak_color[0] * intensity
+
+    if cylinders:
+        yy, xx = np.mgrid[0:ny, 0:nx]
+        for (cx, cy), radius in cylinders:
+            mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
+            rgb[mask] = 0.0
+
+    if float(blur_sigma) > 0.0:
+        for c in range(3):
+            rgb[:, :, c] = gaussian_blur2d(rgb[:, :, c], sigma=float(blur_sigma))
+
+    rgb = np.clip(rgb, 0.0, 1.0)
+    import os
+    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
+
+    try:
+        import matplotlib
+        matplotlib.use("Agg")
+        import matplotlib.pyplot as plt
+        fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
+        ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
+        ax.axis("off")
+        ax.imshow(rgb, origin="lower", interpolation="nearest", aspect="auto")
+        fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
+        plt.close(fig)
+    except ImportError:
+        from PIL import Image
+        Image.fromarray((np.clip(rgb, 0.0, 1.0) * 255.0).astype(np.uint8)).save(out_path)
+
+    return {
+        "image_path": out_path,
+        "trail_count": int(len(trails)),
+        "density_max": density_max,
+        "num_threads": int(n_threads),
+    }

src/CelerisLab/common/preprocess.py

@@ -2,7 +2,7 @@
 """Host-side geometric helpers for immersed-boundary preprocessing.
 
 Responsibilities:
-    Circle–ray intersection along fluid-to-neighbour links and simple
+    Circle-ray intersection along fluid-to-neighbour links and simple
     discrete sensor-area counting used when building IBM masks and probes.
 """
 
@@ -11,7 +11,7 @@ import numpy as np
 
 
 def find_circle_intersection(x, y, x_neb, y_neb, xc, yc, r0):
-    """Return (px, py) of the closest intersection along fluid→neighbour ray,
+    """Return (px, py) of the closest intersection along fluid->neighbour ray,
     or None if no valid root in [0, 1]."""
     dx, dy = x_neb - x, y_neb - y
     a = dx ** 2 + dy ** 2
@@ -32,30 +32,6 @@ def find_circle_intersection(x, y, x_neb, y_neb, xc, yc, r0):
     return None
 
 
-def find_sphere_ray_segment(px, py, pz, qx, qy, qz, cx, cy, cz, r0):
-    """Closest intersection of segment P→Q with sphere (center C, radius r0).
-
-    Returns:
-        (ix, iy, iz) on the segment with parameter in [0, 1], or None.
-    """
-    dx, dy, dz = qx - px, qy - py, qz - pz
-    a = dx * dx + dy * dy + dz * dz
-    if a == 0.0:
-        return None
-    b = 2.0 * (dx * (px - cx) + dy * (py - cy) + dz * (pz - cz))
-    c = (px - cx) ** 2 + (py - cy) ** 2 + (pz - cz) ** 2 - r0 ** 2
-    det = b * b - 4.0 * a * c
-    if det < 0.0:
-        return None
-    sqrt_det = math.sqrt(det)
-    t1 = (-b - sqrt_det) / (2.0 * a)
-    t2 = (-b + sqrt_det) / (2.0 * a)
-    for t in (t1, t2):
-        if 0.0 <= t <= 1.0:
-            return px + t * dx, py + t * dy, pz + t * dz
-    return None
-
-
 def find_sensor_area(radius):
     area = 0
     for i in range(int(np.floor(-radius)), int(np.ceil(radius))):
@@ -63,4 +39,13 @@ def find_sensor_area(radius):
             if i ** 2 + j ** 2 <= radius ** 2:
                 area += 1
     return area
-    
+
+
+def cylinders_from_triangle_layout(layout: dict) -> "list[tuple[tuple[float, float], float]]":
+    """Build cylinder list from triangle-layout dict used in experiment notebook."""
+    radius = float(layout["radius_lb"])
+    return [
+        ((float(layout["x_apex"]), float(layout["y_center"])), radius),
+        ((float(layout["x_rear"]), float(layout["y_lower"])), radius),
+        ((float(layout["x_rear"]), float(layout["y_upper"])), radius),
+    ]

src/CelerisLab/common/render.py

@@ -0,0 +1,92 @@
+# CelerisLab/common/render.py
+"""
+Field rendering utilities (vorticity computation and visualization).
+
+Separated from the streakline module: these are macroscopic field diagnostics,
+not particle-tracking operations.
+"""
+
+from __future__ import annotations
+
+import os
+from typing import Optional, Sequence
+
+import numpy as np
+
+from ._types import CylinderSpec
+
+
+def compute_vorticity(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
+    """Central-difference vorticity omega = d(uy)/dx - d(ux)/dy on interior."""
+    ux_arr = np.asarray(ux, dtype=np.float64)
+    uy_arr = np.asarray(uy, dtype=np.float64)
+    vort = np.zeros_like(ux_arr, dtype=np.float64)
+    dudy = 0.5 * (ux_arr[2:, 1:-1] - ux_arr[:-2, 1:-1])
+    dvdx = 0.5 * (uy_arr[1:-1, 2:] - uy_arr[1:-1, :-2])
+    vort[1:-1, 1:-1] = dvdx - dudy
+    return vort
+
+
+def render_vorticity_field(
+    vort: np.ndarray,
+    *,
+    nx: int,
+    ny: int,
+    out_path: str,
+    cylinders: Optional[Sequence[CylinderSpec]] = None,
+    vmin: Optional[float] = None,
+    vmax: Optional[float] = None,
+    percentile: float = 99.0,
+    minimal_axes: bool = True,
+) -> dict:
+    """Render vorticity with white at zero, green negative, purple positive."""
+    try:
+        import matplotlib
+
+        matplotlib.use("Agg")
+        import matplotlib.pyplot as plt
+        from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm
+    except ImportError as exc:
+        raise RuntimeError("render_vorticity_field requires matplotlib.") from exc
+
+    field = np.asarray(vort, dtype=np.float64)
+    vmax_field = float(np.nanmax(np.abs(field)))
+    abs_pct = float(np.nanpercentile(np.abs(field), percentile))
+    abs_max = max(abs_pct, 0.35 * vmax_field, 1e-6)
+    if vmin is None:
+        vmin = -abs_max
+    if vmax is None:
+        vmax = abs_max
+    if vmin >= 0.0 or vmax <= 0.0:
+        raise ValueError(f"vmin/vmax must bracket zero: vmin={vmin}, vmax={vmax}")
+
+    cmap = LinearSegmentedColormap.from_list(
+        "vort_green_white_purple",
+        ["#1a9850", "#f7f7f7", "#762a83"],
+        N=256,
+    )
+    norm = TwoSlopeNorm(vmin=float(vmin), vcenter=0.0, vmax=float(vmax))
+
+    rgba = cmap(norm(field))
+    if cylinders:
+        yy, xx = np.mgrid[0:ny, 0:nx]
+        for (cx, cy), radius in cylinders:
+            mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
+            rgba[mask] = (0.0, 0.0, 0.0, 1.0)
+
+    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
+    fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
+    ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
+    ax.imshow(rgba, origin="lower", interpolation="nearest", aspect="auto")
+    if minimal_axes:
+        ax.axis("off")
+    fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
+    plt.close(fig)
+
+    return {
+        "image_path": out_path,
+        "vmin": float(vmin),
+        "vmax": float(vmax),
+        "vort_min": float(np.nanmin(field)),
+        "vort_max": float(np.nanmax(field)),
+    }

src/CelerisLab/common/streakline.py

@@ -1,778 +0,0 @@
-# CelerisLab/common/streakline.py
-"""Streakline utilities for online/offline particle-based flow visualization.
-
-Responsibilities:
-    - Build dense release seeds (point/line/strip).
-    - Integrate particles with RK4 and space-time interpolated velocity fields.
-    - Support both online (in-memory frames from Simulation) and offline playback.
-    - Render density-like streakline images or accumulated path trails.
-"""
-
-from __future__ import annotations
-
-import math
-import os
-from concurrent.futures import ThreadPoolExecutor
-from dataclasses import dataclass, field
-from typing import Any, List, Optional, Sequence, Tuple
-
-import numpy as np
-
-CylinderSpec = Tuple[Tuple[float, float], float]
-
-
-@dataclass(frozen=True)
-class FlowFrame:
-    """One sampled velocity frame."""
-
-    step: int
-    ux: np.ndarray
-    uy: np.ndarray
-
-
-@dataclass(frozen=True)
-class ReleaseConfig:
-    """Release strategy configuration for dense streakline seeding."""
-
-    mode: str = "strip"  # point | line | strip
-    line_span: float = 12.0
-    line_count: int = 5
-    downstream_count: int = 4
-    downstream_spacing: float = 3.0
-    inject_per_seed: int = 2
-
-
-@dataclass(frozen=True)
-class IntegratorConfig:
-    """Particle integration and filtering settings."""
-
-    alpha_t: float = 0.2
-    alpha_x: float = 0.4
-    age_decay_steps: float = 20000.0
-    max_particle_age: Optional[float] = None
-    diffusion_coeff: float = 0.0
-    random_seed: int = 1234
-    num_threads: int = 0
-
-
-@dataclass
-class ParticleTrailSet:
-    """Growing polylines per particle (pathline history — not physical streakline)."""
-
-    trails: List[np.ndarray] = field(default_factory=list)
-    active: List[bool] = field(default_factory=list)
-
-    def inject(self, points: np.ndarray) -> None:
-        pts = np.asarray(points, dtype=np.float64)
-        if pts.ndim != 2 or pts.shape[1] != 2:
-            raise ValueError("points must be [N,2].")
-        for p in pts:
-            self.trails.append(np.asarray([[p[0], p[1]]], dtype=np.float64))
-            self.active.append(True)
-
-    def keep_indices(self, survived_mask: np.ndarray) -> None:
-        """Deactivate trails whose particle did not survive the last advance."""
-        mask = np.asarray(survived_mask, dtype=bool).ravel()
-        j = 0
-        for i, is_active in enumerate(self.active):
-            if not is_active:
-                continue
-            if j >= mask.size:
-                break
-            if not bool(mask[j]):
-                self.active[i] = False
-            j += 1
-
-    def keep_mask(self, survived_mask: np.ndarray) -> None:
-        """Alias for keep_indices (boolean mask over active particles)."""
-        self.keep_indices(survived_mask)
-
-    def append_positions(self, particles: np.ndarray) -> None:
-        pts = np.asarray(particles, dtype=np.float64)
-        j = 0
-        for i, is_active in enumerate(self.active):
-            if not is_active:
-                continue
-            if j >= pts.shape[0]:
-                break
-            self.trails[i] = np.vstack([self.trails[i], pts[j]])
-            j += 1
-
-
-def configure_compute_threads(num_threads: int = 0) -> int:
-    """Set BLAS/OpenMP thread env vars for host-side numpy work."""
-    if num_threads <= 0:
-        num_threads = min(32, os.cpu_count() or 4)
-    for var in (
-        "OMP_NUM_THREADS",
-        "OPENBLAS_NUM_THREADS",
-        "MKL_NUM_THREADS",
-        "NUMEXPR_NUM_THREADS",
-        "VECLIB_MAXIMUM_THREADS",
-    ):
-        os.environ[var] = str(int(num_threads))
-    return int(num_threads)
-
-
-def cylinders_from_triangle_layout(layout: dict) -> List[CylinderSpec]:
-    """Build cylinder list from triangle-layout dict used in experiment notebook."""
-    radius = float(layout["radius_lb"])
-    return [
-        ((float(layout["x_apex"]), float(layout["y_center"])), radius),
-        ((float(layout["x_rear"]), float(layout["y_lower"])), radius),
-        ((float(layout["x_rear"]), float(layout["y_upper"])), radius),
-    ]
-
-
-def estimate_sampling_plan(
-    *,
-    st_ref: float,
-    diameter: float,
-    u_ref: float,
-    snapshots_per_period: float = 24.0,
-    periods: int = 5,
-) -> dict:
-    """Estimate save interval and snapshot count from nominal shedding scale."""
-    period_steps = float(diameter) / (float(st_ref) * float(u_ref))
-    save_every = int(max(20, round(period_steps / snapshots_per_period / 10.0) * 10))
-    n_snapshots = int(max(20, round(float(periods) * float(snapshots_per_period))))
-    return {
-        "st_ref": float(st_ref),
-        "period_steps_est": float(period_steps),
-        "snapshots_per_period_target": float(snapshots_per_period),
-        "save_every_recommended": int(save_every),
-        "snapshot_count_recommended": int(n_snapshots),
-    }
-
-
-def build_release_points(
-    base_points: np.ndarray,
-    release_cfg: ReleaseConfig,
-) -> np.ndarray:
-    """Expand base release points to dense points by line/strip policy."""
-    pts = np.asarray(base_points, dtype=np.float64)
-    if pts.ndim != 2 or pts.shape[1] != 2:
-        raise ValueError("base_points must be [N,2].")
-    mode = str(release_cfg.mode).lower()
-    if mode not in ("point", "line", "strip"):
-        raise ValueError("release mode must be point|line|strip.")
-
-    line_count = max(1, int(release_cfg.line_count))
-    ds_count = max(1, int(release_cfg.downstream_count))
-    line_offsets = np.linspace(
-        -0.5 * float(release_cfg.line_span),
-        0.5 * float(release_cfg.line_span),
-        line_count,
-    )
-    ds_offsets = np.arange(ds_count, dtype=np.float64) * float(release_cfg.downstream_spacing)
-
-    out: list[list[float]] = []
-    for x0, y0 in pts:
-        line_iter = [0.0] if mode == "point" else line_offsets
-        ds_iter = [0.0] if mode in ("point", "line") else ds_offsets
-        for dy in line_iter:
-            for dx in ds_iter:
-                out.append([float(x0 + dx), float(y0 + dy)])
-    return np.asarray(out, dtype=np.float64)
-
-
-def _bilinear_sample(field: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
-    ny, nx = field.shape
-    x0 = np.floor(x).astype(np.int64)
-    y0 = np.floor(y).astype(np.int64)
-    x1 = x0 + 1
-    y1 = y0 + 1
-    inside = (x0 >= 0) & (x1 < nx) & (y0 >= 0) & (y1 < ny)
-    out = np.full_like(x, np.nan, dtype=np.float64)
-    if not np.any(inside):
-        return out
-
-    xi = x[inside] - x0[inside]
-    yi = y[inside] - y0[inside]
-    f00 = field[y0[inside], x0[inside]]
-    f10 = field[y0[inside], x1[inside]]
-    f01 = field[y1[inside], x0[inside]]
-    f11 = field[y1[inside], x1[inside]]
-    out[inside] = (
-        (1.0 - xi) * (1.0 - yi) * f00
-        + xi * (1.0 - yi) * f10
-        + (1.0 - xi) * yi * f01
-        + xi * yi * f11
-    )
-    return out
-
-
-def _velocity_at(pos: np.ndarray, theta: float, f0: FlowFrame, f1: FlowFrame) -> np.ndarray:
-    x = pos[:, 0]
-    y = pos[:, 1]
-    vx0 = _bilinear_sample(f0.ux, x, y)
-    vy0 = _bilinear_sample(f0.uy, x, y)
-    vx1 = _bilinear_sample(f1.ux, x, y)
-    vy1 = _bilinear_sample(f1.uy, x, y)
-    vx = (1.0 - theta) * vx0 + theta * vx1
-    vy = (1.0 - theta) * vy0 + theta * vy1
-    return np.stack([vx, vy], axis=1)
-
-
-def _particles_valid_mask(
-    p_new: np.ndarray,
-    *,
-    nx: int,
-    ny: int,
-    solid_center: Optional[Tuple[float, float]] = None,
-    solid_radius: Optional[float] = None,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-) -> np.ndarray:
-    valid = np.all(np.isfinite(p_new), axis=1)
-    valid &= (p_new[:, 0] >= 0.0) & (p_new[:, 0] <= nx - 1)
-    valid &= (p_new[:, 1] >= 0.0) & (p_new[:, 1] <= ny - 1)
-
-    specs: List[CylinderSpec] = []
-    if cylinders:
-        specs.extend(cylinders)
-    elif solid_center is not None and solid_radius is not None:
-        specs.append((solid_center, float(solid_radius)))
-
-    for (cx, cy), radius in specs:
-        dx = p_new[:, 0] - float(cx)
-        dy = p_new[:, 1] - float(cy)
-        valid &= (dx * dx + dy * dy) > (float(radius) * float(radius))
-    return valid
-
-
-def _inject_particles(
-    particles: np.ndarray,
-    ages: np.ndarray,
-    release_points: np.ndarray,
-    *,
-    inject_per_seed: int,
-) -> Tuple[np.ndarray, np.ndarray]:
-    new_pts = np.repeat(release_points, max(1, int(inject_per_seed)), axis=0)
-    if particles.size == 0:
-        return new_pts.copy(), np.zeros((new_pts.shape[0],), dtype=np.float64)
-    return (
-        np.vstack([particles, new_pts]),
-        np.concatenate([ages, np.zeros((new_pts.shape[0],), dtype=np.float64)]),
-    )
-
-
-def advance_particles_between_frames(
-    particles: np.ndarray,
-    ages: np.ndarray,
-    frame0: FlowFrame,
-    frame1: FlowFrame,
-    *,
-    nx: int,
-    ny: int,
-    cfg: IntegratorConfig,
-    solid_center: Optional[Tuple[float, float]] = None,
-    solid_radius: Optional[float] = None,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-    rng: Optional[np.random.Generator] = None,
-) -> Tuple[np.ndarray, np.ndarray, dict, np.ndarray]:
-    """Advance particles from frame0.step to frame1.step with RK4.
-
-    Returns:
-        particles, ages, diagnostics, survived_mask (bool, length = input count).
-    """
-    if particles.size == 0:
-        return particles, ages, {"n_substeps": 0, "max_speed": 0.0}, np.zeros((0,), dtype=bool)
-
-    dt_save = float(frame1.step - frame0.step)
-    max_speed = float(
-        max(
-            np.nanmax(np.sqrt(frame0.ux * frame0.ux + frame0.uy * frame0.uy)),
-            np.nanmax(np.sqrt(frame1.ux * frame1.ux + frame1.uy * frame1.uy)),
-        )
-    )
-    umax = max(max_speed, 1e-8)
-    dt_trace = min(float(cfg.alpha_t) * dt_save, float(cfg.alpha_x) / umax)
-    n_sub = max(1, int(math.ceil(dt_save / dt_trace)))
-    dt_sub = dt_save / float(n_sub)
-
-    p = np.asarray(particles, dtype=np.float64)
-    a = np.asarray(ages, dtype=np.float64)
-    track = np.ones((p.shape[0],), dtype=bool)
-    local_rng = rng or np.random.default_rng(int(cfg.random_seed))
-    diffusion = float(cfg.diffusion_coeff)
-
-    for sub_idx in range(n_sub):
-        if p.size == 0:
-            break
-        theta = float(sub_idx) / float(n_sub)
-        k1 = _velocity_at(p, theta, frame0, frame1)
-        k2 = _velocity_at(p + 0.5 * dt_sub * k1, min(1.0, theta + 0.5 / n_sub), frame0, frame1)
-        k3 = _velocity_at(p + 0.5 * dt_sub * k2, min(1.0, theta + 0.5 / n_sub), frame0, frame1)
-        k4 = _velocity_at(p + dt_sub * k3, min(1.0, theta + 1.0 / n_sub), frame0, frame1)
-        p_new = p + (k1 + 2.0 * k2 + 2.0 * k3 + k4) * (dt_sub / 6.0)
-
-        if diffusion > 0.0:
-            sigma = math.sqrt(max(0.0, 2.0 * diffusion * dt_sub))
-            p_new += local_rng.normal(0.0, sigma, size=p_new.shape)
-
-        valid = _particles_valid_mask(
-            p_new,
-            nx=nx,
-            ny=ny,
-            solid_center=solid_center,
-            solid_radius=solid_radius,
-            cylinders=cylinders,
-        )
-        track[track] = valid
-        p = p_new[valid]
-        a = a[valid] + dt_sub
-        if cfg.max_particle_age is not None:
-            age_keep = a <= float(cfg.max_particle_age)
-            track[track] = age_keep
-            p = p[age_keep]
-            a = a[age_keep]
-
-    survived_mask = track
-    return p, a, {"dt_trace": float(dt_sub), "n_substeps": int(n_sub), "max_speed": float(max_speed)}, survived_mask
-
-
-def run_streakline_offline(
-    frames: Sequence[FlowFrame],
-    *,
-    nx: int,
-    ny: int,
-    release_points: np.ndarray,
-    release_cfg: ReleaseConfig,
-    integrator_cfg: IntegratorConfig,
-    solid_center: Optional[Tuple[float, float]] = None,
-    solid_radius: Optional[float] = None,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-) -> Tuple[np.ndarray, np.ndarray, dict]:
-    """Reconstruct streakline from pre-collected velocity frames."""
-    if len(frames) < 2:
-        raise ValueError("At least two frames are required.")
-    pts = build_release_points(release_points, release_cfg)
-    particles = np.zeros((0, 2), dtype=np.float64)
-    ages = np.zeros((0,), dtype=np.float64)
-    rng = np.random.default_rng(int(integrator_cfg.random_seed))
-    max_speed = 0.0
-    last_step = int(frames[0].step)
-    total_substeps = 0
-
-    for i in range(len(frames) - 1):
-        particles, ages = _inject_particles(
-            particles,
-            ages,
-            pts,
-            inject_per_seed=release_cfg.inject_per_seed,
-        )
-        particles, ages, diag, _ = advance_particles_between_frames(
-            particles,
-            ages,
-            frames[i],
-            frames[i + 1],
-            nx=nx,
-            ny=ny,
-            cfg=integrator_cfg,
-            solid_center=solid_center,
-            solid_radius=solid_radius,
-            cylinders=cylinders,
-            rng=rng,
-        )
-        max_speed = max(max_speed, float(diag["max_speed"]))
-        total_substeps += int(diag["n_substeps"])
-        last_step = int(frames[i + 1].step)
-
-    return particles, ages, {
-        "frames_used": int(len(frames)),
-        "last_step": int(last_step),
-        "n_particles_alive": int(particles.shape[0]),
-        "max_speed_seen": float(max_speed),
-        "total_substeps": int(total_substeps),
-    }
-
-
-def run_streakline_online(
-    sim: Any,
-    *,
-    start_step: int,
-    sample_every: int,
-    n_samples: int,
-    release_points: np.ndarray,
-    release_cfg: ReleaseConfig,
-    integrator_cfg: IntegratorConfig,
-    solid_center: Optional[Tuple[float, float]] = None,
-    solid_radius: Optional[float] = None,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-) -> Tuple[np.ndarray, np.ndarray, dict]:
-    """Compute streakline online from a running Simulation without disk snapshots."""
-    if sample_every < 1 or n_samples < 2:
-        raise ValueError("sample_every must be >=1 and n_samples must be >=2.")
-    if not getattr(sim, "_initialized", False):
-        raise RuntimeError("Simulation must be initialized before online streakline run.")
-
-    nx = int(sim.lbm_cfg.nx)
-    ny = int(sim.lbm_cfg.ny)
-    frames: list[FlowFrame] = []
-    target_last = int(start_step + sample_every * (n_samples - 1))
-
-    while int(sim.stepper.step_count) < target_last:
-        sim.step(1)
-        step = int(sim.stepper.step_count)
-        if step < int(start_step):
-            continue
-        if (step - int(start_step)) % int(sample_every) != 0:
-            continue
-        macro = sim.get_macroscopic()
-        frames.append(
-            FlowFrame(
-                step=step,
-                ux=np.asarray(macro["ux"], dtype=np.float64),
-                uy=np.asarray(macro["uy"], dtype=np.float64),
-            )
-        )
-        if len(frames) >= int(n_samples):
-            break
-
-    if len(frames) < 2:
-        raise RuntimeError("Online run did not collect enough frames.")
-
-    particles, ages, diag = run_streakline_offline(
-        frames,
-        nx=nx,
-        ny=ny,
-        release_points=release_points,
-        release_cfg=release_cfg,
-        integrator_cfg=integrator_cfg,
-        solid_center=solid_center,
-        solid_radius=solid_radius,
-        cylinders=cylinders,
-    )
-    diag["mode"] = "online"
-    diag["sample_every"] = int(sample_every)
-    diag["start_step"] = int(start_step)
-    return particles, ages, diag
-
-
-def compute_vorticity(ux: np.ndarray, uy: np.ndarray) -> np.ndarray:
-    """Central-difference vorticity ω = ∂uy/∂x − ∂ux/∂y on interior cells."""
-    ux_arr = np.asarray(ux, dtype=np.float64)
-    uy_arr = np.asarray(uy, dtype=np.float64)
-    vort = np.zeros_like(ux_arr, dtype=np.float64)
-    dudy = 0.5 * (ux_arr[2:, 1:-1] - ux_arr[:-2, 1:-1])
-    dvdx = 0.5 * (uy_arr[1:-1, 2:] - uy_arr[1:-1, :-2])
-    vort[1:-1, 1:-1] = dvdx - dudy
-    return vort
-
-
-def render_vorticity_field(
-    vort: np.ndarray,
-    *,
-    nx: int,
-    ny: int,
-    out_path: str,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-    vmin: Optional[float] = None,
-    vmax: Optional[float] = None,
-    percentile: float = 99.0,
-    minimal_axes: bool = True,
-) -> dict:
-    """Render vorticity with white at zero, green negative, purple positive."""
-    try:
-        import matplotlib
-
-        matplotlib.use("Agg")
-        import matplotlib.pyplot as plt
-        from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm
-    except ImportError as exc:
-        raise RuntimeError("render_vorticity_field requires matplotlib.") from exc
-
-    field = np.asarray(vort, dtype=np.float64)
-    vmax_field = float(np.nanmax(np.abs(field)))
-    abs_pct = float(np.nanpercentile(np.abs(field), percentile))
-    # Freestream dominates |ω|; do not let a near-zero percentile wash out wake structure.
-    abs_max = max(abs_pct, 0.35 * vmax_field, 1e-6)
-    if vmin is None:
-        vmin = -abs_max
-    if vmax is None:
-        vmax = abs_max
-    if vmin >= 0.0 or vmax <= 0.0:
-        raise ValueError(f"vmin/vmax must bracket zero: vmin={vmin}, vmax={vmax}")
-
-    cmap = LinearSegmentedColormap.from_list(
-        "vort_green_white_purple",
-        ["#1a9850", "#f7f7f7", "#762a83"],
-        N=256,
-    )
-    norm = TwoSlopeNorm(vmin=float(vmin), vcenter=0.0, vmax=float(vmax))
-
-    rgba = cmap(norm(field))
-    if cylinders:
-        yy, xx = np.mgrid[0:ny, 0:nx]
-        for (cx, cy), radius in cylinders:
-            mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
-            rgba[mask] = (0.0, 0.0, 0.0, 1.0)
-
-    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
-    fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
-    ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
-    ax.imshow(rgba, origin="lower", interpolation="nearest", aspect="auto")
-    if minimal_axes:
-        ax.axis("off")
-    fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
-    plt.close(fig)
-
-    return {
-        "image_path": out_path,
-        "vmin": float(vmin),
-        "vmax": float(vmax),
-        "vort_min": float(np.nanmin(field)),
-        "vort_max": float(np.nanmax(field)),
-    }
-
-
-def gaussian_blur2d(img: np.ndarray, sigma: float = 1.0) -> np.ndarray:
-    """Separable Gaussian blur helper."""
-    if float(sigma) <= 0.0:
-        return np.asarray(img, dtype=np.float64)
-    radius = max(1, int(math.ceil(3.0 * float(sigma))))
-    x = np.arange(-radius, radius + 1, dtype=np.float64)
-    k = np.exp(-(x**2) / (2.0 * sigma * sigma))
-    k /= np.sum(k)
-    arr = np.asarray(img, dtype=np.float64)
-    arr = np.apply_along_axis(lambda m: np.convolve(m, k, mode="same"), 1, arr)
-    arr = np.apply_along_axis(lambda m: np.convolve(m, k, mode="same"), 0, arr)
-    return arr
-
-
-def _rasterize_segment(
-    intensity: np.ndarray,
-    x0: float,
-    y0: float,
-    x1: float,
-    y1: float,
-    weight: float,
-) -> None:
-    """Stamp one segment into a scalar intensity map (for white-bg red compose)."""
-    ny, nx = intensity.shape
-    n = int(max(abs(x1 - x0), abs(y1 - y0), 1.0)) + 1
-    xs = np.linspace(x0, x1, n)
-    ys = np.linspace(y0, y1, n)
-    xi = np.clip(np.rint(xs).astype(np.int64), 0, nx - 1)
-    yi = np.clip(np.rint(ys).astype(np.int64), 0, ny - 1)
-    w = float(max(0.0, min(1.0, weight)))
-    np.add.at(intensity, (yi, xi), w)
-
-
-def _draw_trail(
-    intensity: np.ndarray,
-    trail: np.ndarray,
-    *,
-    fade_along_trail: bool,
-) -> None:
-    if trail.shape[0] < 2:
-        return
-    nseg = trail.shape[0] - 1
-    for i in range(nseg):
-        if fade_along_trail:
-            t = float(i + 1) / float(nseg)
-            weight = 0.15 + 0.85 * t
-        else:
-            weight = 1.0
-        _rasterize_segment(
-            intensity,
-            float(trail[i, 0]),
-            float(trail[i, 1]),
-            float(trail[i + 1, 0]),
-            float(trail[i + 1, 1]),
-            weight,
-        )
-
-
-def _draw_cylinders(rgb: np.ndarray, cylinders: Sequence[CylinderSpec]) -> None:
-    ny, nx = rgb.shape[0], rgb.shape[1]
-    yy, xx = np.mgrid[0:ny, 0:nx]
-    for (cx, cy), radius in cylinders:
-        mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
-        rgb[mask] = 0.0
-
-
-def render_snapshot_trails(
-    trail_set: ParticleTrailSet,
-    *,
-    nx: int,
-    ny: int,
-    out_path: str,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-    blur_sigma: float = 0.0,
-    fade_along_trail: bool = True,
-    num_threads: int = 0,
-    background_color: Tuple[float, float, float] = (1.0, 1.0, 1.0),
-    streak_color: Tuple[float, float, float] = (1.0, 0.0, 0.0),
-) -> dict:
-    """Render accumulated path trails (pathline collage — debug only, not streakline)."""
-    n_threads = configure_compute_threads(num_threads)
-    intensity = np.zeros((ny, nx), dtype=np.float64)
-
-    trails = trail_set.trails
-    if not trails:
-        density_max = 0.0
-    else:
-        chunk = max(1, (len(trails) + n_threads - 1) // n_threads)
-
-        def _worker(batch: List[np.ndarray]) -> np.ndarray:
-            local = np.zeros((ny, nx), dtype=np.float64)
-            for tr in batch:
-                _draw_trail(local, tr, fade_along_trail=fade_along_trail)
-            return local
-
-        batches = [trails[i : i + chunk] for i in range(0, len(trails), chunk)]
-        if len(batches) == 1:
-            intensity = _worker(batches[0])
-        else:
-            with ThreadPoolExecutor(max_workers=n_threads) as pool:
-                parts = list(pool.map(_worker, batches))
-            for part in parts:
-                intensity = np.maximum(intensity, part)
-
-        if np.any(intensity > 0):
-            vmax = float(np.percentile(intensity[intensity > 0], 99.0))
-            intensity = np.clip(intensity / max(vmax, 1e-12), 0.0, 1.0)
-        density_max = float(np.max(intensity))
-
-    rgb = np.ones((ny, nx, 3), dtype=np.float64)
-    rgb[:, :, 0] = background_color[0]
-    rgb[:, :, 1] = background_color[1] * (1.0 - intensity) + streak_color[1] * intensity
-    rgb[:, :, 2] = background_color[2] * (1.0 - intensity) + streak_color[2] * intensity
-    if streak_color[0] < 0.999:
-        rgb[:, :, 0] = background_color[0] * (1.0 - intensity) + streak_color[0] * intensity
-
-    if cylinders:
-        _draw_cylinders(rgb, cylinders)
-
-    if float(blur_sigma) > 0.0:
-        for c in range(3):
-            rgb[:, :, c] = gaussian_blur2d(rgb[:, :, c], sigma=float(blur_sigma))
-
-    rgb = np.clip(rgb, 0.0, 1.0)
-    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
-    try:
-        import matplotlib
-
-        matplotlib.use("Agg")
-        import matplotlib.pyplot as plt
-
-        fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
-        ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
-        ax.axis("off")
-        ax.imshow(rgb, origin="lower", interpolation="nearest", aspect="auto")
-        fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
-        plt.close(fig)
-    except ImportError:
-        from PIL import Image
-
-        Image.fromarray((np.clip(rgb, 0.0, 1.0) * 255.0).astype(np.uint8)).save(out_path)
-
-    return {
-        "image_path": out_path,
-        "trail_count": int(len(trails)),
-        "density_max": density_max,
-        "num_threads": int(n_threads),
-    }
-
-
-def render_streakline_density(
-    positions: np.ndarray,
-    ages: np.ndarray,
-    *,
-    nx: int,
-    ny: int,
-    out_path: str,
-    release_points: Optional[np.ndarray] = None,
-    solid_center: Optional[Tuple[float, float]] = None,
-    solid_radius: Optional[float] = None,
-    cylinders: Optional[Sequence[CylinderSpec]] = None,
-    age_decay_steps: float = 20000.0,
-    blur_sigma: float = 1.2,
-    title: str = "Streakline",
-    minimal_axes: bool = False,
-    background_color: Tuple[float, float, float] = (0.0, 0.0, 0.0),
-    streak_color: Tuple[float, float, float] = (1.0, 0.5, 0.0),
-    show_release_points: bool = True,
-) -> dict:
-    """Render streakline as current particle positions (density cloud at observation time)."""
-    try:
-        import matplotlib
-
-        matplotlib.use("Agg")
-        import matplotlib.pyplot as plt
-        from matplotlib.colors import LinearSegmentedColormap
-    except ImportError as exc:
-        raise RuntimeError("render_streakline_density requires matplotlib.") from exc
-
-    if positions.size == 0:
-        density = np.zeros((ny, nx), dtype=np.float64)
-    else:
-        w = np.exp(-np.asarray(ages, dtype=np.float64) / max(1.0, float(age_decay_steps)))
-        y_idx = np.clip(np.floor(positions[:, 1]).astype(np.int64), 0, ny - 1)
-        x_idx = np.clip(np.floor(positions[:, 0]).astype(np.int64), 0, nx - 1)
-        density = np.zeros((ny, nx), dtype=np.float64)
-        np.add.at(density, (y_idx, x_idx), w)
-    density = gaussian_blur2d(density, sigma=float(blur_sigma))
-
-    if minimal_axes:
-        rgb = np.ones((ny, nx, 3), dtype=np.float64)
-        rgb[:, :, 0] = background_color[0]
-        rgb[:, :, 1] = background_color[1]
-        rgb[:, :, 2] = background_color[2]
-        if np.any(density > 0):
-            vmax = float(np.percentile(density[density > 0], 99.0))
-            alpha = np.clip(density / max(vmax, 1e-12), 0.0, 1.0)
-            for c, val in enumerate(streak_color):
-                rgb[:, :, c] = np.clip(rgb[:, :, c] * (1.0 - alpha) + val * alpha, 0.0, 1.0)
-        specs = list(cylinders or [])
-        if solid_center is not None and solid_radius is not None:
-            specs.append((solid_center, float(solid_radius)))
-        if specs:
-            _draw_cylinders(rgb, specs)
-        os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
-        fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
-        ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
-        ax.axis("off")
-        ax.imshow(rgb, origin="lower", interpolation="nearest", aspect="auto")
-        fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
-        plt.close(fig)
-        return {
-            "image_path": out_path,
-            "n_particles": int(positions.shape[0]) if positions.ndim == 2 else 0,
-            "density_max": float(np.max(density)) if density.size else 0.0,
-        }
-
-    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
-    fig, ax = plt.subplots(figsize=(12, 5))
-    vmax = np.percentile(density[density > 0], 99.0) if np.any(density > 0) else 1.0
-    cmap = LinearSegmentedColormap.from_list("streak", [(0.0, 0.0, 0.0), streak_color])
-    ax.imshow(
-        density,
-        origin="lower",
-        cmap=cmap,
-        extent=(0, nx - 1, 0, ny - 1),
-        vmin=0.0,
-        vmax=max(vmax, 1e-12),
-        aspect="equal",
-    )
-    if solid_center is not None and solid_radius is not None:
-        circ = plt.Circle(solid_center, radius=solid_radius, fill=False, color="cyan", linewidth=1.4)
-        ax.add_patch(circ)
-    if show_release_points and release_points is not None and np.asarray(release_points).size:
-        rp = np.asarray(release_points)
-        ax.scatter(rp[:, 0], rp[:, 1], c="white", s=10, marker="x")
-    ax.set_title(title)
-    ax.set_xlabel("x (lattice)")
-    ax.set_ylabel("y (lattice)")
-    ax.set_xlim(0, nx - 1)
-    ax.set_ylim(0, ny - 1)
-    fig.tight_layout()
-    fig.savefig(out_path, dpi=170, bbox_inches="tight")
-    plt.close(fig)
-
-    return {
-        "image_path": out_path,
-        "n_particles": int(positions.shape[0]) if positions.ndim == 2 else 0,
-        "density_max": float(np.max(density)) if density.size else 0.0,
-    }

src/CelerisLab/common/streakline/__init__.py

@@ -0,0 +1,49 @@
+# CelerisLab/common/streakline/__init__.py
+"""
+Streakline — continuous-release particle-based flow visualization.
+
+A post-processing sub-package under ``common/``.  The ``Streakline`` class is a
+passive consumer: users call ``observe()`` in their own simulation loop,
+controlling body actions, checkpointing, or other modules independently.
+
+Usage::
+
+    from CelerisLab.common.streakline import Streakline, ReleaseConfig, IntegratorConfig
+
+    streak = Streakline(release_points=..., nx=nx, ny=ny, ...)
+    streak.observe(ux=macro["ux"], uy=macro["uy"], step=step)
+    streak.render("output.png")
+"""
+
+from ._config import (
+    ReleaseConfig,
+    IntegratorConfig,
+    FlowFrame,
+    CylinderSpec,
+    build_release_points,
+)
+from ._streakline import Streakline
+from ._integrate import (
+    advance_between_frames as advance_particles_between_frames,
+    bilinear_sample,
+    inject_particles,
+    particles_valid_mask,
+)
+from ._render import render_density, gaussian_blur2d
+
+# configure_compute_threads is imported lazily (in _streakline.py) to avoid
+# circular dependency between common.pathline and common.streakline.
+
+# Backward-compat aliases
+render_streakline_density = render_density
+
+__all__ = [
+    "Streakline",
+    "ReleaseConfig",
+    "IntegratorConfig",
+    "FlowFrame",
+    "CylinderSpec",
+    "build_release_points",
+    "advance_particles_between_frames",
+    "render_streakline_density",
+]

src/CelerisLab/common/streakline/_config.py

@@ -0,0 +1,94 @@
+# CelerisLab/common/streakline/_config.py
+"""
+Configuration dataclasses for streakline particle integration.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Optional
+
+import numpy as np
+
+from .._types import CylinderSpec
+
+
+@dataclass(frozen=True)
+class FlowFrame:
+    """One sampled velocity frame (internal use only)."""
+    step: int
+    ux: np.ndarray
+    uy: np.ndarray
+
+
+@dataclass(frozen=True)
+class ReleaseConfig:
+    """Release strategy configuration for dense streakline seeding.
+
+    Args:
+        mode: ``"point"`` | ``"line"`` | ``"strip"``.
+        line_span: Total cross-stream span for ``line`` / ``strip`` modes.
+        line_count: Number of cross-stream release points.
+        downstream_count: Number of streamwise points (``strip`` mode only).
+        downstream_spacing: Spacing between downstream copies (lattice units).
+        inject_per_seed: Particles injected per release point each frame.
+    """
+    mode: str = "strip"
+    line_span: float = 12.0
+    line_count: int = 5
+    downstream_count: int = 4
+    downstream_spacing: float = 3.0
+    inject_per_seed: int = 2
+
+
+@dataclass(frozen=True)
+class IntegratorConfig:
+    """Particle integration and filtering settings.
+
+    Args:
+        alpha_t: Fraction of frame interval used as substep ceiling.
+        alpha_x: Fraction of lattice unit per substep (CFL-like).
+        age_decay_steps: Exponential decay scale for age-weight in rendering.
+        max_particle_age: Hard cutoff age (None = no cutoff).
+        diffusion_coeff: Random-walk diffusion for visual thickening (0 = off).
+        random_seed: RNG seed for diffusion.
+    """
+    alpha_t: float = 0.2
+    alpha_x: float = 0.4
+    age_decay_steps: float = 20000.0
+    max_particle_age: Optional[float] = None
+    diffusion_coeff: float = 0.0
+    random_seed: int = 1234
+
+
+def build_release_points(
+    base_points: np.ndarray,
+    release_cfg: ReleaseConfig,
+) -> np.ndarray:
+    """Expand *base_points* to a dense release grid per *release_cfg*."""
+    pts = np.asarray(base_points, dtype=np.float64)
+    if pts.ndim != 2 or pts.shape[1] != 2:
+        raise ValueError("base_points must be [N,2].")
+    mode = str(release_cfg.mode).lower()
+    if mode not in ("point", "line", "strip"):
+        raise ValueError("release mode must be point|line|strip.")
+
+    line_count = max(1, int(release_cfg.line_count))
+    ds_count = max(1, int(release_cfg.downstream_count))
+    line_offsets = np.linspace(
+        -0.5 * float(release_cfg.line_span),
+        0.5 * float(release_cfg.line_span),
+        line_count,
+    )
+    ds_offsets = np.arange(ds_count, dtype=np.float64) * float(
+        release_cfg.downstream_spacing
+    )
+
+    out: list[list[float]] = []
+    for x0, y0 in pts:
+        line_iter = [0.0] if mode == "point" else line_offsets
+        ds_iter = [0.0] if mode in ("point", "line") else ds_offsets
+        for dy in line_iter:
+            for dx in ds_iter:
+                out.append([float(x0 + dx), float(y0 + dy)])
+    return np.asarray(out, dtype=np.float64)

src/CelerisLab/common/streakline/_integrate.py

@@ -0,0 +1,180 @@
+# CelerisLab/common/streakline/_integrate.py
+"""
+Particle integration: bilinear interpolation, RK4 substep, validity filtering.
+"""
+
+from __future__ import annotations
+
+import math
+from typing import Optional, Sequence
+
+import numpy as np
+
+from ._config import CylinderSpec, FlowFrame, IntegratorConfig
+
+
+def bilinear_sample(field: np.ndarray, x: np.ndarray, y: np.ndarray) -> np.ndarray:
+    """Bilinear interpolation at floating-point positions ``(x, y)``.
+
+    Returns ``NaN`` for positions outside the domain.
+    """
+    ny, nx = field.shape
+    x0 = np.floor(x).astype(np.int64)
+    y0 = np.floor(y).astype(np.int64)
+    x1 = x0 + 1
+    y1 = y0 + 1
+    inside = (x0 >= 0) & (x1 < nx) & (y0 >= 0) & (y1 < ny)
+    out = np.full_like(x, np.nan, dtype=np.float64)
+    if not np.any(inside):
+        return out
+
+    xi = x[inside] - x0[inside]
+    yi = y[inside] - y0[inside]
+    f00 = field[y0[inside], x0[inside]]
+    f10 = field[y0[inside], x1[inside]]
+    f01 = field[y1[inside], x0[inside]]
+    f11 = field[y1[inside], x1[inside]]
+    out[inside] = (
+        (1.0 - xi) * (1.0 - yi) * f00
+        + xi * (1.0 - yi) * f10
+        + (1.0 - xi) * yi * f01
+        + xi * yi * f11
+    )
+    return out
+
+
+def velocity_at(
+    pos: np.ndarray, theta: float, f0: FlowFrame, f1: FlowFrame
+) -> np.ndarray:
+    """Space-time interpolated velocity at positions ``pos`` (N, 2)."""
+    x = pos[:, 0]
+    y = pos[:, 1]
+    vx0 = bilinear_sample(f0.ux, x, y)
+    vy0 = bilinear_sample(f0.uy, x, y)
+    vx1 = bilinear_sample(f1.ux, x, y)
+    vy1 = bilinear_sample(f1.uy, x, y)
+    vx = (1.0 - theta) * vx0 + theta * vx1
+    vy = (1.0 - theta) * vy0 + theta * vy1
+    return np.stack([vx, vy], axis=1)
+
+
+def particles_valid_mask(
+    p_new: np.ndarray,
+    *,
+    nx: int,
+    ny: int,
+    cylinders: Sequence[CylinderSpec],
+) -> np.ndarray:
+    """Boolean mask: particles within domain bounds and outside solids."""
+    valid = np.all(np.isfinite(p_new), axis=1)
+    valid &= (p_new[:, 0] >= 0.0) & (p_new[:, 0] <= nx - 1)
+    valid &= (p_new[:, 1] >= 0.0) & (p_new[:, 1] <= ny - 1)
+
+    for (cx, cy), radius in cylinders:
+        dx = p_new[:, 0] - float(cx)
+        dy = p_new[:, 1] - float(cy)
+        valid &= (dx * dx + dy * dy) > (float(radius) * float(radius))
+    return valid
+
+
+def inject_particles(
+    particles: np.ndarray,
+    ages: np.ndarray,
+    release_points: np.ndarray,
+    *,
+    inject_per_seed: int,
+) -> tuple[np.ndarray, np.ndarray]:
+    """Append new particles at *release_points* with age 0."""
+    new_pts = np.repeat(release_points, max(1, int(inject_per_seed)), axis=0)
+    if particles.size == 0:
+        return new_pts.copy(), np.zeros((new_pts.shape[0],), dtype=np.float64)
+    return (
+        np.vstack([particles, new_pts]),
+        np.concatenate([ages, np.zeros((new_pts.shape[0],), dtype=np.float64)]),
+    )
+
+
+def advance_between_frames(
+    particles: np.ndarray,
+    ages: np.ndarray,
+    frame0: FlowFrame,
+    frame1: FlowFrame,
+    *,
+    nx: int,
+    ny: int,
+    cfg: IntegratorConfig,
+    cylinders: Sequence[CylinderSpec],
+    rng: Optional[np.random.Generator] = None,
+) -> tuple[np.ndarray, np.ndarray, dict, np.ndarray]:
+    """Advance particles from *frame0* to *frame1* with RK4.
+
+    Returns:
+        particles, ages, diagnostics, survived_mask (input-length bool).
+    """
+    if particles.size == 0:
+        return (
+            particles,
+            ages,
+            {"n_substeps": 0, "max_speed": 0.0},
+            np.zeros((0,), dtype=bool),
+        )
+
+    dt_save = float(frame1.step - frame0.step)
+    max_speed = float(
+        max(
+            np.nanmax(np.sqrt(frame0.ux * frame0.ux + frame0.uy * frame0.uy)),
+            np.nanmax(np.sqrt(frame1.ux * frame1.ux + frame1.uy * frame1.uy)),
+        )
+    )
+    umax = max(max_speed, 1e-8)
+    dt_trace = min(float(cfg.alpha_t) * dt_save, float(cfg.alpha_x) / umax)
+    n_sub = max(1, int(math.ceil(dt_save / dt_trace)))
+    dt_sub = dt_save / float(n_sub)
+
+    p = np.asarray(particles, dtype=np.float64)
+    a = np.asarray(ages, dtype=np.float64)
+    track = np.ones((p.shape[0],), dtype=bool)
+    local_rng = rng or np.random.default_rng(int(cfg.random_seed))
+    diffusion = float(cfg.diffusion_coeff)
+
+    for sub_idx in range(n_sub):
+        if p.size == 0:
+            break
+        theta = float(sub_idx) / float(n_sub)
+        k1 = velocity_at(p, theta, frame0, frame1)
+        k2 = velocity_at(
+            p + 0.5 * dt_sub * k1, min(1.0, theta + 0.5 / n_sub), frame0, frame1
+        )
+        k3 = velocity_at(
+            p + 0.5 * dt_sub * k2, min(1.0, theta + 0.5 / n_sub), frame0, frame1
+        )
+        k4 = velocity_at(
+            p + dt_sub * k3, min(1.0, theta + 1.0 / n_sub), frame0, frame1
+        )
+        p_new = p + (k1 + 2.0 * k2 + 2.0 * k3 + k4) * (dt_sub / 6.0)
+
+        if diffusion > 0.0:
+            sigma = math.sqrt(max(0.0, 2.0 * diffusion * dt_sub))
+            p_new += local_rng.normal(0.0, sigma, size=p_new.shape)
+
+        valid = particles_valid_mask(p_new, nx=nx, ny=ny, cylinders=cylinders)
+        track[track] = valid
+        p = p_new[valid]
+        a = a[valid] + dt_sub
+        if cfg.max_particle_age is not None:
+            age_keep = a <= float(cfg.max_particle_age)
+            track[track] = age_keep
+            p = p[age_keep]
+            a = a[age_keep]
+
+    survived_mask = track
+    return (
+        p,
+        a,
+        {
+            "dt_trace": float(dt_sub),
+            "n_substeps": int(n_sub),
+            "max_speed": float(max_speed),
+        },
+        survived_mask,
+    )

src/CelerisLab/common/streakline/_render.py

@@ -0,0 +1,148 @@
+# CelerisLab/common/streakline/_render.py
+"""
+Render streakline particle clouds to density images.
+"""
+
+from __future__ import annotations
+
+import math
+import os
+from typing import Optional, Sequence, Tuple
+
+import numpy as np
+
+from ._config import CylinderSpec
+
+
+def gaussian_blur2d(img: np.ndarray, sigma: float = 1.0) -> np.ndarray:
+    """Separable Gaussian blur."""
+    if float(sigma) <= 0.0:
+        return np.asarray(img, dtype=np.float64)
+    radius = max(1, int(math.ceil(3.0 * float(sigma))))
+    x = np.arange(-radius, radius + 1, dtype=np.float64)
+    k = np.exp(-(x**2) / (2.0 * sigma * sigma))
+    k /= np.sum(k)
+    arr = np.asarray(img, dtype=np.float64)
+    for i in range(arr.shape[0]):
+        arr[i, :] = np.convolve(arr[i, :], k, mode="same")
+    for j in range(arr.shape[1]):
+        arr[:, j] = np.convolve(arr[:, j], k, mode="same")
+    return arr
+
+
+def _draw_cylinders(rgb: np.ndarray, cylinders: Sequence[CylinderSpec]) -> None:
+    ny, nx = rgb.shape[0], rgb.shape[1]
+    yy, xx = np.mgrid[0:ny, 0:nx]
+    for (cx, cy), radius in cylinders:
+        mask = (xx - float(cx)) ** 2 + (yy - float(cy)) ** 2 <= float(radius) ** 2
+        rgb[mask] = 0.0
+
+
+def render_density(
+    positions: np.ndarray,
+    ages: np.ndarray,
+    *,
+    nx: int,
+    ny: int,
+    out_path: str,
+    cylinders: Sequence[CylinderSpec] = (),
+    age_decay_steps: float = 20000.0,
+    blur_sigma: float = 1.2,
+    background_color: tuple[float, float, float] = (0.0, 0.0, 0.0),
+    streak_color: tuple[float, float, float] = (1.0, 0.5, 0.0),
+    show_release_points: Optional[np.ndarray] = None,
+) -> dict:
+    """Render streakline particle cloud as a density image.
+
+    Particles are binned onto the pixel grid with age-weighted exponential
+    decay, Gaussian-blurred, then composited onto *background_color* using
+    *streak_color*.
+    """
+    density = _build_density_map(
+        positions, ages, nx, ny, age_decay_steps=age_decay_steps
+    )
+    density = gaussian_blur2d(density, sigma=float(blur_sigma))
+    return _save_minimal_image(
+        density,
+        nx,
+        ny,
+        out_path,
+        cylinders,
+        background_color=background_color,
+        streak_color=streak_color,
+        release_points=show_release_points,
+    )
+
+
+def _build_density_map(
+    positions: np.ndarray,
+    ages: np.ndarray,
+    nx: int,
+    ny: int,
+    *,
+    age_decay_steps: float,
+) -> np.ndarray:
+    if positions.size == 0:
+        return np.zeros((ny, nx), dtype=np.float64)
+    w = np.exp(-np.asarray(ages, dtype=np.float64) / max(1.0, float(age_decay_steps)))
+    y_idx = np.clip(np.floor(positions[:, 1]).astype(np.int64), 0, ny - 1)
+    x_idx = np.clip(np.floor(positions[:, 0]).astype(np.int64), 0, nx - 1)
+    density = np.zeros((ny, nx), dtype=np.float64)
+    np.add.at(density, (y_idx, x_idx), w)
+    return density
+
+
+def _save_minimal_image(
+    density: np.ndarray,
+    nx: int,
+    ny: int,
+    out_path: str,
+    cylinders: Sequence[CylinderSpec],
+    *,
+    background_color: tuple[float, float, float],
+    streak_color: tuple[float, float, float],
+    release_points: Optional[np.ndarray] = None,
+) -> dict:
+    """Compose density into an RGB image and save (no matplotlib axes)."""
+    rgb = np.ones((ny, nx, 3), dtype=np.float64)
+    rgb[:, :, 0] = background_color[0]
+    rgb[:, :, 1] = background_color[1]
+    rgb[:, :, 2] = background_color[2]
+
+    if np.any(density > 0):
+        vmax = float(np.percentile(density[density > 0], 99.0))
+        alpha = np.clip(density / max(vmax, 1e-12), 0.0, 1.0)
+        for c, val in enumerate(streak_color):
+            rgb[:, :, c] = np.clip(
+                rgb[:, :, c] * (1.0 - alpha) + val * alpha, 0.0, 1.0
+            )
+
+    if cylinders:
+        _draw_cylinders(rgb, cylinders)
+
+    os.makedirs(os.path.dirname(os.path.abspath(out_path)), exist_ok=True)
+
+    try:
+        import matplotlib
+
+        matplotlib.use("Agg")
+        import matplotlib.pyplot as plt
+
+        fig = plt.figure(figsize=(nx / 100.0, ny / 100.0), dpi=100)
+        ax = fig.add_axes([0.0, 0.0, 1.0, 1.0])
+        ax.axis("off")
+        ax.imshow(rgb, origin="lower", interpolation="nearest", aspect="auto")
+        fig.savefig(out_path, dpi=100, bbox_inches="tight", pad_inches=0)
+        plt.close(fig)
+    except ImportError:
+        from PIL import Image
+
+        Image.fromarray((np.clip(rgb, 0.0, 1.0) * 255.0).astype(np.uint8)).save(
+            out_path
+        )
+
+    return {
+        "image_path": out_path,
+        "n_particles": int(density.sum()),
+        "density_max": float(np.max(density)) if density.size else 0.0,
+    }

src/CelerisLab/common/streakline/_streakline.py

@@ -0,0 +1,117 @@
+# CelerisLab/common/streakline/_streakline.py
+"""
+Streakline — passive-consumer particle-based flow visualization.
+"""
+
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+from typing import Optional, Sequence
+
+import numpy as np
+
+from ._config import (
+    CylinderSpec,
+    FlowFrame,
+    IntegratorConfig,
+    ReleaseConfig,
+    build_release_points,
+)
+from ._integrate import advance_between_frames, inject_particles
+from ._render import render_density
+
+
+@dataclass
+class Streakline:
+    """Passive streakline integrator.
+
+    Does **not** own a simulation or control the stepping loop.  Users call
+    ``observe()`` whenever they have a new velocity frame; body actions,
+    checkpointing, and other modules are managed independently.
+
+    Usage::
+
+        streak = Streakline(release_points=..., nx=nx, ny=ny, ...)
+        # In your own simulation loop:
+        macro = sim.get_macroscopic()
+        streak.observe(ux=macro["ux"], uy=macro["uy"], step=step)
+        streak.render("output.png")
+    """
+
+    release_points: np.ndarray
+    release_cfg: ReleaseConfig
+    integrator_cfg: IntegratorConfig
+    nx: int
+    ny: int
+    cylinders: Sequence[CylinderSpec] = field(default_factory=list)
+
+    _particles: np.ndarray = field(init=False, repr=False)
+    _ages: np.ndarray = field(init=False, repr=False)
+    _prev_frame: Optional[FlowFrame] = field(init=False, repr=False, default=None)
+
+    def __post_init__(self):
+        self._particles = np.zeros((0, 2), dtype=np.float64)
+        self._ages = np.zeros((0,), dtype=np.float64)
+        self._dense_pts = build_release_points(
+            np.asarray(self.release_points), self.release_cfg
+        )
+
+    def observe(self, ux: np.ndarray, uy: np.ndarray, step: int) -> dict:
+        """Feed one velocity snapshot.
+
+        First call caches the frame.  From the second call onward particles
+        are injected and advanced from the previous frame to this one.
+        """
+        frame = FlowFrame(
+            step=int(step),
+            ux=np.asarray(ux, dtype=np.float64),
+            uy=np.asarray(uy, dtype=np.float64),
+        )
+
+        if self._prev_frame is None:
+            self._prev_frame = frame
+            return {}
+
+        self._particles, self._ages = inject_particles(
+            self._particles,
+            self._ages,
+            self._dense_pts,
+            inject_per_seed=self.release_cfg.inject_per_seed,
+        )
+
+        self._particles, self._ages, diag, _ = advance_between_frames(
+            self._particles,
+            self._ages,
+            self._prev_frame,
+            frame,
+            nx=self.nx,
+            ny=self.ny,
+            cfg=self.integrator_cfg,
+            cylinders=self.cylinders,
+        )
+
+        self._prev_frame = frame
+        return diag
+
+    def render(self, out_path: str, **kwargs) -> dict:
+        """Render current particle cloud as a density image."""
+        # Default age_decay_steps from integrator config; **kwargs may override
+        kw = dict(kwargs, cylinders=self.cylinders)
+        kw.setdefault("age_decay_steps", self.integrator_cfg.age_decay_steps)
+        return render_density(
+            self._particles,
+            self._ages,
+            nx=self.nx,
+            ny=self.ny,
+            out_path=out_path,
+            **kw,
+        )
+
+    @property
+    def n_particles(self) -> int:
+        return self._particles.shape[0]
+
+    def reset(self):
+        self._particles = np.zeros((0, 2), dtype=np.float64)
+        self._ages = np.zeros((0,), dtype=np.float64)
+        self._prev_frame = None

src/CelerisLab/configs/CONFIG.md

@@ -45,7 +45,7 @@ Python `config.py` 只负责读取和校验，不是配置位置。
 | `outlet.srt_neq_damp` | float | 0.5 | [0, 1] | 仅当 `outlet.mode` 为 **`neq_extrap`** 且碰撞为 **SRT/TRT** 时：出口全分布 damped NEQ 的阻尼系数。`outlet.mode` 为 **`blended`** 时 SRT/TRT 走混合未知方向分支，不使用本键。MRT + `neq_extrap` 仍为少量方向的 NEQ 外推 |
 | `y_wall_bc` | string | `"bounce_back"` | `bounce_back`, `free_slip` | y 向上下壁面边界条件（仅壁面相邻流体行生效） |
 | `omega_guard.min` | float | 0.01 | | ω 下界（LES τ_eff 保护） |
-| `omega_guard.max` | float | 1.96 | | ω 上界，高 Re 建议 1.90-1.96 |
+| `omega_guard.max` | float | 1.99 | | ω 上界，高 Re 建议 1.90-1.99 |
 
 ## cuda — 编译
 
@@ -88,5 +88,5 @@ step/one_step_*.cu         → kernel 编排
 
 ### 稳定性默认窗口（Stability defaults）
 
-- 默认 `method.omega_guard = {min: 0.01, max: 1.96}`，用于约束碰撞频率与 LES 有效粘度路径。
-- 高 Re 调参建议将 `omega_guard.max` 保持在 `1.90-1.96` 区间；过高上界会增大接近 `omega -> 2` 的不稳定风险。
+- 默认 `method.omega_guard = {min: 0.01, max: 1.99}`，用于约束碰撞频率与 LES 有效粘度路径。
+- 高 Re 调参建议将 `omega_guard.max` 保持在 `1.90-1.99` 区间；过高上界会增大接近 `omega -> 2` 的不稳定风险。

src/CelerisLab/configs/config_lbm.json

@@ -53,8 +53,8 @@
         "_y_wall_bc": "bounce_back | free_slip",
         "omega_guard": {
             "min": 0.01,
-            "max": 1.96,
-            "_doc": "High-Re: keep max in 1.90-1.96"
+            "max": 1.99,
+            "_doc": "High-Re: keep max in 1.90-1.99"
         }
     },
     "cuda": {

src/CelerisLab/lbm/kernels/config/config_grid.h

@@ -6,8 +6,8 @@
 #define NT   256
 #define MULT_GPU 0
 
-#define NX  3000
-#define NY  300
+#define NX  384
+#define NY  192
 #define NZ  1
 
 // ---- Lattice model (single source of truth) ----

src/CelerisLab/lbm/kernels/config/config_method.h

@@ -20,7 +20,7 @@
 #define Y_WALL_BC            1
 
 #define OMEGA_COLLISION_MIN 0.01f
-#define OMEGA_COLLISION_MAX 1.960f
+#define OMEGA_COLLISION_MAX 1.990f
 #define TRT_MAGIC_PARAM     0.187500f
 
 // NEQ damping coefficients for inlet/outlet BC reconstruction.

src/CelerisLab/lbm/kernels/config/config_objects.h

@@ -3,6 +3,6 @@
 #ifndef CELERIS_CONFIG_OBJECTS_H
 #define CELERIS_CONFIG_OBJECTS_H
 
-#define N_OBJS 3
+#define N_OBJS 2
 
 #endif

src/CelerisLab/lbm/kernels/config/config_physics.h

@@ -4,9 +4,9 @@
 #define CELERIS_CONFIG_PHYSICS_H
 
 #define LBtype float
-#define VIS 0.0040000000
+#define VIS 0.0035000000
 #define RHO 1.0
-#define U0  0.04
+#define U0  0.03
 
 #define PI 3.141592653589793238
 

src/CelerisLab/lbm/kernels/operators/macro.cuh

@@ -155,7 +155,10 @@ __device__ __forceinline__ void compute_feq(float rho, float ux, float uy, float
 #endif  // LATTICE_MODEL (compute_feq)
 
 // ---------------------------------------------------------------------------
-// compute_pressure (diagnostic, from state equation p = cs² ρ)
+// --- Diagnostic only ---
+// ---------------------------------------------------------------------------
+// compute_pressure / compute_pressure_perturbation are diagnostic helpers
+// not consumed by any kernel on a hot path.
 // ---------------------------------------------------------------------------
 __device__ __forceinline__ float compute_pressure(float rho) {
     return CS2 * rho;

src/CelerisLab/lbm/kernels/step/one_step_double.cu

@@ -79,6 +79,8 @@ void OneStep(
     // has_bc(fl): BC subtype != NONE -> wall/inlet/outlet nodes only
     // Obstacle nodes (FLAG_OBSTACLE, BC subtype=NONE) are excluded here
     // and handled solely by the is_obstacle branch below.
+    // TODO: BC_MOVING and BC_PERIODIC are defined in flags.cuh but NOT
+    //       dispatched here — they silently fall through to bounce_back_swap.
     if (is_solid(fl) && has_bc(fl))
         apply_boundary_pull(f, fl, x, y, fi_in, k);
 

src/CelerisLab/lbm/kernels/step/one_step_esopull.cu

@@ -99,6 +99,22 @@ void EsoPullStep(
     if (is_obstacle(fl) && !is_curved(fl))
         bounce_back_swap(f);
 
+    // ---- y-wall adjacent row correction ----
+    // EsoPull's alternating read/write pattern means `fi` stores the previous
+    // step's post-collision value at this node's opposite directions — same
+    // semantic as `fi_in` in the double-buffer path.  The existing wall-BB
+    // helpers read from `fi` and work correctly.
+    //
+    // NOTE: This is a first-order correction.  A dedicated EsoPull bounce-back
+    //       scheme may be needed for higher accuracy (see FluidX3D).
+    if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) {
+#if Y_WALL_BC == 1
+        apply_wall_freeslip_d2q9_y_pull(y, f, fi, k);
+#else
+        apply_wall_bb_d2q9_y_pull(y, f, fi, k);
+#endif
+    }
+
     // Collision
     // Same donor/ghost warning as the double-buffer path: for the local Zou-He
     // inlet, the ghost-node state must be regularized before the next pull uses
@@ -138,6 +154,11 @@ void EsoPullStep(
     if (is_obstacle(fl) && !is_curved(fl))
         bounce_back_swap(f);
 
+    // ---- y-wall adjacent row correction (3D) ----
+    // Same logic as the D2Q9 path — see note above.
+    if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2)))
+        apply_wall_bb_d3q19_y_pull(y, f, fi, k);
+
     const bool collide_inlet_ghost = is_inlet(fl) && interior_y
                                    && inlet_scheme_uses_post_collision_ghost();
     if (is_fluid(fl) || collide_inlet_ghost) {

src/CelerisLab/simulation.py

@@ -1,17 +1,20 @@
 # CelerisLab/simulation.py
 """
-Top-level orchestrator — assembles LBM field, stepper, body manager,
+Top-level orchestrator -- assembles LBM field, stepper, body manager,
 and CUDA context into a single coherent simulation.
 
 Usage::
 
     sim = Simulation("configs/config_lbm.json")
-    sim.add_cylinder((100, 50), radius=10)
+    sim.add_body("circle", center=(50, 50), radius=10)
     sim.initialize()
     sim.run(1000)
     macro = sim.get_macroscopic()   # {"rho": ..., "ux": ..., "uy": ...}
+    force = sim.read_force(0)       # force on body 0
 """
 
+from __future__ import annotations
+
 from typing import Dict, Optional, Tuple, Any
 
 import os
@@ -24,7 +27,7 @@ from .cuda.context import CudaContext
 from .cuda import compiler_v2 as compiler
 from .lbm.field import LBMField
 from .lbm.stepper import LBMStepper
-from .body.objects import Cylinder, Sensor, SimObject
+from .body.objects import SimObject
 from .body.manager import ObjectManager
 
 
@@ -39,6 +42,9 @@ class Simulation:
                  lbm_config_path: Optional[str] = None,
                  body_config_path: Optional[str] = None,
                  device_id: int = 0):
+        # Internal stream for async DTOH coordination
+        self._stream: Optional[cuda.Stream] = None
+
         # Load configs
         self.lbm_cfg = load_lbm_config(lbm_config_path)
         self.body_cfg = load_body_config(body_config_path)
@@ -68,18 +74,64 @@ class Simulation:
             self.recompile()
         self._assert_object_count_contract(expected_count=self.bodies.count)
 
+    # -- Stream access -------------------------------------------------------
+    @property
+    def stream(self) -> cuda.Stream:
+        """Internal CUDA stream for async DTOH coordination."""
+        if self._stream is None:
+            self._stream = cuda.Stream()
+        return self._stream
+
     # -- Object management ---------------------------------------------------
+    def add_body(self, type: str = "circle", center=None,
+                 radius: float = 0.0, **kwargs) -> int:
+        """Add a simulation body. Returns body_id.
+
+        Args:
+            type: ``"circle"`` or ``"sensor"`` (future: ``"polygon"``, ``"mesh"``).
+            center: (x, y) center coordinates (lattice units).
+            radius: body radius (lattice units).
+            **kwargs: passed to the geometry constructor.
+
+        Returns:
+            int body_id.
+        """
+        type_lower = str(type).strip().lower()
+
+        if type_lower == "circle":
+            from .body.geometry.circle import CircleGeometry
+            cx, cy = float(center[0]), float(center[1])
+            geom = CircleGeometry(cx, cy, radius)
+            obj = SimObject(obj_id=-1, geometry=geom,
+                            center=center, radius=radius,
+                            is_sensor=False)
+        elif type_lower == "sensor":
+            from .body.geometry.circle import CircleGeometry
+            cx, cy = float(center[0]), float(center[1])
+            geom = CircleGeometry(cx, cy, radius)
+            obj = SimObject(obj_id=-1, geometry=geom,
+                            center=center, radius=radius,
+                            is_sensor=True)
+        else:
+            raise ValueError(
+                f"Unknown body type '{type}'. "
+                "Supported: 'circle', 'sensor'."
+            )
+        return self.bodies.add(obj)
+
+    # -- Legacy convenience (keep for backward compat) -----------------------
     def add_cylinder(self, center: Tuple[float, ...],
                      radius: float) -> int:
-        obj = Cylinder(obj_id=-1, center=center, radius=radius)
-        return self.bodies.add(obj)
+        """Add a cylinder body (convenience wrapper around add_body)."""
+        return self.add_body("circle", center=center, radius=radius)
 
     def add_sensor(self, center: Tuple[float, ...],
                    radius: float) -> int:
-        obj = Sensor(obj_id=-1, center=center, radius=radius)
-        return self.bodies.add(obj)
+        """Add a sensor body (convenience wrapper around add_body)."""
+        return self.add_body("sensor", center=center, radius=radius)
 
     def add_object(self, obj: SimObject) -> int:
+        """Add a pre-configured SimObject directly."""
         return self.bodies.add(obj)
 
     def _load_objects_from_body_config(self) -> None:
@@ -91,24 +143,75 @@ class Simulation:
             omega = spec.get("omega", None)
 
             if obj_type == "":
-                raise ValueError("Body config object is missing required key 'type'.")
+                raise ValueError(
+                    "Body config object is missing required key 'type'.")
             if len(center) not in (2, 3):
                 raise ValueError(
-                    f"Body config object '{obj_type}' requires center with 2 or 3 entries."
+                    f"Body config object '{obj_type}' requires center "
+                    "with 2 or 3 entries."
                 )
 
             if obj_type == "cylinder":
-                body_id = self.add_cylinder(center=center, radius=radius)
+                body_id = self.add_body("circle", center=center, radius=radius)
                 if omega is not None:
                     self.bodies.get(body_id).state.omega = float(omega)
             elif obj_type == "sensor":
-                self.add_sensor(center=center, radius=radius)
+                self.add_body("sensor", center=center, radius=radius)
             else:
                 raise ValueError(
                     f"Unsupported body config object type '{obj_type}'. "
                     "Supported types: cylinder, sensor."
                 )
 
+    # -- Runtime control -----------------------------------------------------
+    def set_body(self, id: int, omega: float = None,
+                 vx: float = None, vy: float = None) -> None:
+        """Set runtime body state. Updates are implicitly uploaded to GPU.
+
+        Args:
+            id: body_id from add_body().
+            omega: rotation speed (lattice units).  ``None`` = no change.
+            vx: translation velocity x (future).  ``None`` = no change.
+            vy: translation velocity y (future).  ``None`` = no change.
+        """
+        obj = self.bodies.get(id)
+        changed = False
+        if omega is not None:
+            obj.state.omega = float(omega)
+            changed = True
+        if vx is not None:
+            obj.state.vx = float(vx)
+            changed = True
+        if vy is not None:
+            obj.state.vy = float(vy)
+            changed = True
+        if changed:
+            self.bodies._refresh_action_from_objects()
+            if self.bodies.action_gpu is not None:
+                cuda.memcpy_htod(self.bodies.action_gpu, self.bodies.action)
+
+    # -- Telemetry readback --------------------------------------------------
+    def read_force(self, id: int) -> np.ndarray:
+        """Download and return the force vector on body *id* (async stream, synced)."""
+        self._stream_obs_download()
+        return self.bodies.read_force(id)
+
+    def read_torque(self, id: int) -> np.ndarray:
+        """Download and return the torque on body *id* (async stream, synced)."""
+        self._stream_obs_download()
+        return self.bodies.read_torque(id)
+
+    def read_sensor(self, id: int) -> np.ndarray:
+        """Download and return the sensor reading for body *id* (async stream, synced)."""
+        self._stream_obs_download()
+        return self.bodies.read_sensor(id)
+
+    def _stream_obs_download(self) -> None:
+        """Async obs download on internal stream, then synchronize."""
+        if self.bodies.obs_pinned is not None:
+            self.bodies.download_obs_full_async(self.stream)
+            self.stream.synchronize()
+
     # -- Compilation ---------------------------------------------------------
     def recompile(self):
         """Re-generate config headers and recompile kernel.
@@ -138,7 +241,8 @@ class Simulation:
     def initialize(self):
         """Initialize flow field and sync objects to GPU."""
         # Recompile if objects were added after construction.
-        if self.bodies.count > 0 and self._read_generated_n_objects() != self.bodies.count:
+        if self.bodies.count > 0 and \
+           self._read_generated_n_objects() != self.bodies.count:
             self.recompile()
 
         self.field.flag = self.field.build_channel_flags()
@@ -263,25 +367,33 @@ class Simulation:
                 line = line.strip()
                 if line.startswith("#define N_OBJS"):
                     return int(line.split()[-1])
-        raise RuntimeError(f"N_OBJS not found in generated header: {cfg_path}")
+        raise RuntimeError(
+            f"N_OBJS not found in generated header: {cfg_path}")
 
-    def _assert_object_count_contract(self, expected_count: int | None = None) -> None:
+    def _assert_object_count_contract(
+        self, expected_count: int | None = None
+    ) -> None:
         """Ensure generated compile-time N_OBJS matches runtime object count."""
-        runtime_count = self.bodies.count if expected_count is None else int(expected_count)
+        runtime_count = (
+            self.bodies.count if expected_count is None
+            else int(expected_count)
+        )
         compiled_count = self._read_generated_n_objects()
         if compiled_count != runtime_count:
             raise RuntimeError(
                 "Object count contract mismatch: generated N_OBJS="
-                f"{compiled_count}, runtime count={runtime_count}. Recompile required."
+                f"{compiled_count}, runtime count={runtime_count}. "
+                "Recompile required."
             )
 
     def _assert_obs_layout_contract(self) -> None:
-        """Ensure host obs layout matches runtime body count and generated contract."""
+        """Ensure host obs layout matches runtime body count and contract."""
         lay = compiler.obs_layout(self.lbm_cfg.dim, self.bodies.count)
         mgr = self.bodies
         if mgr.obs_n_slots != lay.n_slots:
             raise RuntimeError(
-                f"obs layout mismatch: n_slots={mgr.obs_n_slots}, expected={lay.n_slots}."
+                "obs layout mismatch: n_slots="
+                f"{mgr.obs_n_slots}, expected={lay.n_slots}."
             )
         if (
             mgr.slot_stride_floats != lay.slot_stride_floats
@@ -290,11 +402,12 @@ class Simulation:
             or mgr.obs_total_floats != lay.total_floats
         ):
             raise RuntimeError(
-                "obs layout mismatch between runtime buffers and expected layout "
-                f"for dim={self.lbm_cfg.dim}, n_objects={self.bodies.count}."
+                "obs layout mismatch between runtime buffers and expected "
+                f"layout for dim={self.lbm_cfg.dim}, "
+                f"n_objects={self.bodies.count}."
             )
 
     def _assert_runtime_contracts(self) -> None:
-        """Validate object-count and packed-observation contracts before stepping."""
+        """Validate object-count and obs layout contracts before stepping."""
         self._assert_object_count_contract()
-        self._assert_obs_layout_contract()
+        self._assert_obs_layout_contract()