Preserve D2 repairs and add bounded jet and fixed-patch kernels

pyproject.toml

@@ -7,7 +7,7 @@ name = "CelerisLab"
 version = "0.5.1"
 description = "GPU-accelerated Lattice Boltzmann Method (LBM) CFD solver using CUDA"
 readme = "README.md"
-requires-python = ">=3.8"
+requires-python = ">=3.10"
 license = {text = "MIT"}
 authors = [
     {name = "Frank14f"}
@@ -19,8 +19,6 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Physics",
     "License :: OSI Approved :: MIT License",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
-    "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
 ]

setup.py

@@ -26,15 +26,13 @@ setup(
         'numpy>=1.19.0',
         'scipy>=1.5.0',
     ],
-    python_requires='>=3.8',
+    python_requires='>=3.10',
     classifiers=[
         'Development Status :: 3 - Alpha',
         'Intended Audience :: Science/Research',
         'Topic :: Scientific/Engineering :: Physics',
         'License :: OSI Approved :: MIT License',
         'Programming Language :: Python :: 3',
-        'Programming Language :: Python :: 3.8',
-        'Programming Language :: Python :: 3.9',
         'Programming Language :: Python :: 3.10',
         'Programming Language :: Python :: 3.11',
     ],

src/CelerisLab/__init__.py

@@ -1,28 +1,16 @@
-# CelerisLab/__init__.py
+"""GPU LBM library. CPU utilities do not initialize or import CUDA eagerly."""
-"""
+from importlib import import_module
-CelerisLab: GPU-Accelerated Computational Physics Simulation Library
 
-Usage::
+__version__ = "0.5.1"
+__all__ = ["Simulation", "lbm", "body", "common", "cuda", "config", "__version__"]
 
-    from CelerisLab import Simulation
 
-    sim = Simulation("configs/config_lbm.json")
+def __getattr__(name):
-    sim.add_body("circle", center=(100, 50), radius=10)
+    if name == "Simulation":
-    sim.initialize()
+        value = import_module(".simulation", __name__).Simulation
-    sim.run(1000)
+    elif name in ("lbm", "body", "common", "cuda", "config"):
-    data = sim.get_macroscopic()
+        value = import_module("." + name, __name__)
-    force = sim.read_force(0)
+    else:
-"""
+        raise AttributeError(name)
-
+    globals()[name] = value
-__version__ = "0.5.0"
+    return value
-
-from . import common, cuda, lbm, body, config
-
-try:
-    from .simulation import Simulation
-    __all__ = ["Simulation", "lbm", "body", "common", "cuda", "config",
-               "__version__"]
-except ImportError as e:
-    import warnings
-    warnings.warn(f"Simulation not available: {e}", ImportWarning)
-    __all__ = ["lbm", "body", "common", "cuda", "config", "__version__"]

src/CelerisLab/body/__init__.py

@@ -2,8 +2,8 @@
 """
 Body / object management for immersed and rigid objects.
 """
-from .objects import SimObject, ObjectState, ObjectControl
+from .objects import SimObject, ObjectState, ObjectControl, SurfaceActuation
-from .manager import ObjectManager
+from importlib import import_module
 from .registry import BodyRegistry
 from .sync_plan import BodySyncPlan
 from .geometry.base import CutLink, SensorCell, Geometry
@@ -11,7 +11,7 @@ from .geometry.circle import CircleGeometry
 from .coupling.soa_packer import pack_cut_links_to_soa, pack_sensor_to_soa
 
 __all__ = [
-    "SimObject", "ObjectState", "ObjectControl",
+    "SimObject", "ObjectState", "ObjectControl", "SurfaceActuation",
     "ObjectManager",
     "BodyRegistry",
     "BodySyncPlan",
@@ -19,3 +19,11 @@ __all__ = [
     "CircleGeometry",
     "pack_cut_links_to_soa", "pack_sensor_to_soa",
 ]
+
+
+def __getattr__(name):
+    if name != "ObjectManager":
+        raise AttributeError(name)
+    value = import_module(".manager", __name__).ObjectManager
+    globals()[name] = value
+    return value

src/CelerisLab/body/ddf_patch.py

@@ -56,7 +56,7 @@ def patch_ddf_for_body_sync(
 
     Args:
         field: LBMField instance (``field.ddf`` is modified in-place).
-        old_ddf: Physical DDF snapshot before topology change (n*NQ, float32).
+        old_ddf: Q-major physical DDF snapshot (NQ*n, float32).
         old_flags: Flag array before topology change.
         new_flags: Flag array after topology change.
         added_solid_mask: Bool array -- True where fluid became solid.
@@ -78,31 +78,34 @@ def upload_patched_ddf(
 ) -> None:
     """Upload only the DDF cells that changed during patching.
 
-    Uses per-cell ``cuda.memcpy_htod`` for both ``ddf_gpu`` and
+    Copies Q-major contiguous cell runs to both ``ddf_gpu`` and
-    ``temp_gpu``.  Leaves all other cells (especially boundary cells
+    ``temp_gpu`` using the field storage encoding.  Leaves all other cells (especially boundary cells
     with NaN) untouched on the GPU, preserving their streaming-
     compatible values from before the sync.
     """
     changed_indices: set[int] = set()
     changed_indices.update(np.where(added_solid_mask)[0].tolist())
     changed_indices.update(np.where(released_fluid_mask)[0].tolist())
-    if curved_fluid_indices is not None:
+    if curved_fluid_indices is not None and np.any(added_solid_mask):
         changed_indices.update(curved_fluid_indices.tolist())
     if not changed_indices:
         field.invalidate_host_ddf_cache()
         field._host_ddf_step = None
         return
 
-    nq = field.nq
+    # Bulk fields are Q-major; the compact EsoPull scratch is cell-major.
-    cell_bytes = nq * 4
+    # Encode only selected values with the same shift/precision path as a
-    ddf_ptr = int(field.ddf_gpu)
+    # whole-field upload, then copy contiguous cell runs in each direction.
-    temp_ptr = int(field.temp_gpu)
+    indices = np.array(sorted(changed_indices), dtype=np.intp)
-
+    encoded = field._encode_ddf(field.ddf.reshape(field.nq, field.n)[:, indices])
-    for idx in sorted(changed_indices):
+    breaks = np.flatnonzero(np.diff(indices) != 1) + 1
-        offset = idx * cell_bytes
+    for positions in np.split(np.arange(indices.size), breaks):
-        cell_data = field.ddf[idx * nq:(idx + 1) * nq].tobytes()
+        start = int(indices[positions[0]])
-        cuda.memcpy_htod(ddf_ptr + offset, cell_data)
+        for direction in range(field.nq):
-        cuda.memcpy_htod(temp_ptr + offset, cell_data)
+            offset = (direction * field.n + start) * field.store_bytes
+            values = np.ascontiguousarray(encoded[direction, positions])
+            cuda.memcpy_htod(int(field.ddf_gpu) + offset, values)
+            cuda.memcpy_htod(int(field.temp_gpu) + offset, values)
 
     field.invalidate_host_ddf_cache()
     field._host_ddf_step = None
@@ -119,7 +122,7 @@ def _patch_added_solid(field, added_solid_mask: np.ndarray) -> None:
     indices = np.where(added_solid_mask)[0]
     nq = field.nq
     for idx in indices:
-        field.ddf[idx * nq:(idx + 1) * nq] = feq
+        field.ddf.reshape(nq, field.n)[:, idx] = feq
 
 
 def _patch_boundary_fluid(
@@ -135,13 +138,13 @@ def _patch_boundary_fluid(
     nq = field.nq
     feq_fallback = compute_feq_d2q9(field.cfg.rho, 0.0, 0.0)
     for idx in curved_fluid_indices:
-        f = field.ddf[idx * nq:(idx + 1) * nq]
+        f = field.ddf.reshape(nq, field.n)[:, idx]
         if np.any(np.isnan(f)):
-            field.ddf[idx * nq:(idx + 1) * nq] = feq_fallback
+            field.ddf.reshape(nq, field.n)[:, idx] = feq_fallback
             continue
         rho, ux, uy = compute_macro_from_ddf(f)
         feq = compute_feq_d2q9(float(rho), float(ux), float(uy))
-        field.ddf[idx * nq:(idx + 1) * nq] = feq
+        field.ddf.reshape(nq, field.n)[:, idx] = feq
 
 
 def _patch_released_fluid(
@@ -173,7 +176,7 @@ def _patch_released_fluid(
     old_macro: dict[int, tuple[float, float, float]] = {}
     fluid_indices = np.where(old_is_fluid)[0]
     for idx in fluid_indices:
-        f = old_ddf[idx * nq:(idx + 1) * nq]
+        f = old_ddf.reshape(nq, field.n)[:, idx]
         if np.any(np.isnan(f)):
             continue
         rho, ux, uy = compute_macro_from_ddf(f)
@@ -199,7 +202,7 @@ def _patch_released_fluid(
         ux_avg = np.mean([old_macro[n][1] for n in fluid_nbrs])
         uy_avg = np.mean([old_macro[n][2] for n in fluid_nbrs])
         feq = compute_feq_d2q9(float(rho_avg), float(ux_avg), float(uy_avg))
-        field.ddf[idx * nq:(idx + 1) * nq] = feq
+        field.ddf.reshape(nq, field.n)[:, idx] = feq
         filled[idx] = True
         prev_layer_filled[idx] = (
             float(rho_avg), float(ux_avg), float(uy_avg))
@@ -221,7 +224,7 @@ def _patch_released_fluid(
             uy_avg = np.mean([prev_layer_filled[n][2] for n in filled_nbrs])
             feq = compute_feq_d2q9(
                 float(rho_avg), float(ux_avg), float(uy_avg))
-            field.ddf[idx * nq:(idx + 1) * nq] = feq
+            field.ddf.reshape(nq, field.n)[:, idx] = feq
             filled[idx] = True
             current_layer[idx] = (
                 float(rho_avg), float(ux_avg), float(uy_avg))
@@ -236,4 +239,4 @@ def _patch_released_fluid(
             u_inlet * 4.0 * y * (ny - 1 - y) / ((ny - 1) ** 2)
             if ny > 2 else u_inlet)
         feq = compute_feq_d2q9(field.cfg.rho, float(ux_fallback), 0.0)
-        field.ddf[idx * nq:(idx + 1) * nq] = feq
+        field.ddf.reshape(nq, field.n)[:, idx] = feq

src/CelerisLab/body/manager.py

@@ -31,6 +31,7 @@ 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 (
+    FLUID, SOLID,
     D2Q9_EX,
     D2Q9_EY,
     D3Q19_EX,
@@ -76,6 +77,7 @@ class ObjectManager:
         self._obs_alloc_nbytes: int = 0
 
         self.action = np.zeros(0, dtype=np.float32)
+        self._action_pinned: Optional[np.ndarray] = None
         self.obs_pinned: Optional[np.ndarray] = None
         self.obs_force_view: Optional[np.ndarray] = None
         self.obs_torque_view: Optional[np.ndarray] = None
@@ -106,13 +108,23 @@ class ObjectManager:
     # -- Registry delegation -------------------------------------------------
     def add(self, obj: SimObject) -> int:
         """Register an object and return its id."""
+        obj.validate_surface_actuation(self.cfg)
         body_id = self._registry.add(obj)
         self._resize_buffers()
         return body_id
 
     def remove(self, obj_id: int):
+        """Remove immediately and renumber remaining IDs, as runtime sync does.
+
+        Previously returned numeric IDs are invalidated by this operation.
+        """
         self._registry.remove(obj_id)
+        remaining = self.objects
+        self._registry.clear()
+        for obj in remaining:
+            self._registry.add(obj)
         self._resize_buffers()
+        self._refresh_action_from_objects()
 
     def get(self, obj_id: int) -> SimObject:
         return self._registry.get(obj_id)
@@ -149,6 +161,7 @@ class ObjectManager:
         The object is not added to the formal registry until ``commit_pending()``
         is called.  This implicitly activates the edit window.
         """
+        obj.validate_surface_actuation(self.cfg)
         self._edit_active = True
         self._pending_add.append(obj)
 
@@ -190,7 +203,8 @@ class ObjectManager:
             new = SimObject(obj_id=oid, geometry=obj.geometry,
                             center=obj.center, radius=obj.radius,
                             is_sensor=obj.is_sensor,
-                            is_force_region=obj.is_force_region)
+                            is_force_region=obj.is_force_region,
+                            surface_actuation=obj.surface_actuation)
             new.state = obj.state
             new.control = obj.control
             return new
@@ -295,6 +309,7 @@ class ObjectManager:
 
         cl_host = plan.curved_host
         field.curved.assign_host(*cl_host)
+        field.curved.assign_surface(plan.next_objects, self.cfg, plan.next_flags)
         s_cells, s_ids = plan.sensor_host
         field.sensors.assign_host(s_cells, s_ids)
         fr_cells, fr_ids = plan.force_region_host
@@ -363,6 +378,13 @@ class ObjectManager:
         cl_rx, cl_ry, cl_rz, cl_fallback = [], [], [], []
         s_cells, s_ids = [], []
         fr_cells, fr_ids = [], []
+        if domain_flags is None:
+            domain_flags = self.build_flags_for(
+                objects, np.full(self.nx * self.ny * self.nz, FLUID, dtype=np.uint16),
+                nx=self.nx, ny=self.ny, nz=self.nz,
+            )
+        active_fluid = ((domain_flags & FLUID) != 0) & ((domain_flags & SOLID) == 0)
+        force_owners: set[int] = set()
 
         if sensor_cell_counts_out is not None:
             sensor_cell_counts_out.fill(0)
@@ -374,6 +396,7 @@ class ObjectManager:
             ex, ey = D2Q9_EX, D2Q9_EY
 
         for obj in objects:
+            obj.validate_surface_actuation(self.cfg)
             if (not obj.is_sensor) and (not obj.is_force_region) and hasattr(obj, 'get_curved_list'):
                 (
                     fluid_idx, dirs, q_vals, body_ids,
@@ -394,6 +417,8 @@ class ObjectManager:
                     cl_fallback.append(fallback_vals)
             if obj.is_sensor and hasattr(obj, 'get_sensor_list'):
                 cells, ids = obj.get_sensor_list(self.nx, self.ny, self.nz)
+                keep = active_fluid[cells]
+                cells, ids = cells[keep], ids[keep]
                 if len(cells) > 0:
                     s_cells.append(cells)
                     s_ids.append(ids)
@@ -402,7 +427,12 @@ class ObjectManager:
                         sensor_cell_counts_out[obj.obj_id] = int(len(cells))
             if obj.is_force_region and hasattr(obj, 'get_sensor_list'):
                 cells, ids = obj.get_sensor_list(self.nx, self.ny, self.nz)
+                keep = active_fluid[cells]
+                cells, ids = cells[keep], ids[keep]
                 if len(cells) > 0:
+                    if force_owners.intersection(map(int, cells)):
+                        raise ValueError("Overlapping force_region fluid cells are unsupported")
+                    force_owners.update(map(int, cells))
                     fr_cells.append(cells)
                     fr_ids.append(ids)
 
@@ -422,6 +452,13 @@ class ObjectManager:
         r_fr_cells = np.concatenate(fr_cells) if fr_cells else np.zeros(0, dtype=np.uint32)
         r_fr_ids = np.concatenate(fr_ids) if fr_ids else np.zeros(0, dtype=np.int32)
 
+        if any(obj.surface_actuation is not None for obj in objects):
+            # Reject unresolved apertures/ownership before sync or restart mutates state.
+            from ..lbm.curved_links import CurvedLinkSoA
+            links = CurvedLinkSoA()
+            links.assign_host(r_idx, r_dir, r_q, r_bid, r_rx, r_ry, r_rz, r_fallback)
+            links.assign_surface(objects, self.cfg, domain_flags)
+
         return (
             r_idx, r_dir, r_q, r_bid,
             r_rx, r_ry, r_rz, r_fallback,
@@ -441,21 +478,24 @@ class ObjectManager:
         """
         self._telemetry_field = field
 
-        if rebuild_flags:
+        next_flags = self.build_flags(field.build_channel_flags()) if rebuild_flags else field.flag
-            field.flag = self.build_flags(field.build_channel_flags())
-            field.upload_flags()
 
         (
             cl_fluid_idx, cl_dir, cl_q, cl_body_id,
             cl_rx, cl_ry, cl_rz, cl_fallback_class,
             sensor_cells, sensor_obj_id,
             fr_cells, fr_obj_id,
-        ) = self.build_compact_lists(domain_flags=field.flag)
+        ) = self.build_compact_lists(domain_flags=next_flags)
+
+        if rebuild_flags:
+            field.flag = next_flags
+            field.upload_flags()
 
         field.curved.assign_host(
             cl_fluid_idx, cl_dir, cl_q, cl_body_id,
             cl_rx, cl_ry, cl_rz, cl_fallback_class,
         )
+        field.curved.assign_surface(self.objects, self.cfg, field.flag)
         field.sensors.assign_host(sensor_cells, sensor_obj_id)
         field.force_regions.assign_host(fr_cells, fr_obj_id)
 
@@ -570,9 +610,16 @@ class ObjectManager:
 
     def _upload_action_async(self, stream: cuda.Stream) -> None:
         """Async H2D of host action array to action_gpu."""
+        for obj in self.objects:
+            obj.validate_surface_actuation(self.cfg)
+            if obj.surface_actuation is not None and self.action[(obj.obj_id+1)*3*self.cfg.dim-1] != 0:
+                raise ValueError("Surface actuation cannot share the rotation action")
         if self.action_gpu is None or self.action.size == 0:
             return
-        cuda.memcpy_htod_async(self.action_gpu, self.action, stream)
+        if self._action_pinned is None or self._action_pinned.size != self.action.size:
+            self._action_pinned = cuda.pagelocked_empty(self.action.size, np.float32)
+        self._action_pinned[:] = self.action
+        cuda.memcpy_htod_async(self.action_gpu, self._action_pinned, stream)
 
     @property
     def obs_n_slots(self) -> int:
@@ -684,9 +731,12 @@ class ObjectManager:
         The GPU action buffer is updated at the next ``run()`` call.
         """
         self._validate_body_id(body_id)
+        obj = self.get(body_id)
+        obj.validate_surface_actuation(self.cfg, dict(vars(obj.state), omega=omega))
         dim = self.cfg.dim if self.cfg else 2
         slot = 3 * dim
         base = body_id * slot
+        self.get(body_id).state.omega = float(omega)
         self.action[base + slot - 1] = np.float32(omega)
 
     def set_force_state(self, body_id: int, fx: float = 0.0, fy: float = 0.0) -> None:
@@ -698,6 +748,14 @@ class ObjectManager:
         dim = self.cfg.dim if self.cfg else 2
         slot = 3 * dim
         base = body_id * slot
+        obj = self.get(body_id)
+        if not obj.is_force_region:
+            raise ValueError("set_force_state requires a force_region object")
+        # vx/vy are the existing force-region action aliases. Keep their
+        # authoritative host state in sync so edits to other bodies cannot
+        # discard this control during a full action refresh or topology sync.
+        obj.state.vx = float(fx)
+        obj.state.vy = float(fy)
         self.action[base + 0] = np.float32(fx)
         self.action[base + 1] = np.float32(fy)
 

src/CelerisLab/body/objects.py

@@ -15,6 +15,7 @@ Design:
 from __future__ import annotations
 
 from dataclasses import dataclass
+import math
 from typing import Optional, Tuple
 
 import numpy as np
@@ -53,6 +54,37 @@ class ObjectControl:
     mode: int = 0  # 0=force, 1=velocity, 2=displacement
 
 
+@dataclass(frozen=True)
+class SurfaceActuation:
+    """Fixed circular aperture: angles in radians, velocity in lattice units.
+
+    ``angle`` is measured counterclockwise from +x; ``half_width`` is the arc
+    half-angle. Positive normal_velocity blows from solid into fluid. The C1
+    profile (1-s*s)^2 uses the true hit angle, not a fluid-node aperture mask.
+    This immutable specification is set before initialization; no moving nozzle.
+    """
+    angle: float
+    half_width: float
+    normal_velocity: float
+
+    def __post_init__(self):
+        if not all(math.isfinite(v) for v in (self.angle, self.half_width, self.normal_velocity)):
+            raise ValueError("Surface actuation values must be finite")
+        if not 0 < self.half_width <= math.pi or abs(self.normal_velocity) >= 1:
+            raise ValueError("Require 0 < half_width <= pi and abs(normal_velocity) < 1")
+
+    def project(self, rx, ry):
+        """Project stored true hit-center vectors; returns velocity and aperture mask."""
+        radius = np.hypot(rx, ry)
+        if np.any(radius <= 0):
+            raise ValueError("Surface normal requires nonzero hit-center distance")
+        delta = np.arctan2(np.sin(np.arctan2(ry, rx)-self.angle),
+                           np.cos(np.arctan2(ry, rx)-self.angle))
+        aperture = np.abs(delta) < self.half_width
+        speed = self.normal_velocity*np.maximum(0., 1-(delta/self.half_width)**2)**2
+        return np.stack((speed*rx/radius, speed*ry/radius), axis=-1).astype(np.float32), aperture
+
+
 # ---------------------------------------------------------------------------
 # SimObject container
 # ---------------------------------------------------------------------------
@@ -77,6 +109,7 @@ class SimObject:
         radius: float = 0.0,
         is_sensor: bool = False,
         is_force_region: bool = False,
+        surface_actuation: Optional[SurfaceActuation] = None,
     ) -> None:
         self.obj_id = obj_id
         self.center = center if center is not None else (0.0, 0.0)
@@ -84,6 +117,10 @@ class SimObject:
         self._geometry = geometry
         self._is_sensor = is_sensor
         self._is_force_region = is_force_region
+        self._surface_actuation = (SurfaceActuation(**surface_actuation)
+                                  if isinstance(surface_actuation, dict) else surface_actuation)
+        if self._surface_actuation is not None and not isinstance(self._surface_actuation, SurfaceActuation):
+            raise TypeError("surface_actuation must be a SurfaceActuation or dict")
 
         c = self.center
         z0 = float(c[2]) if len(c) > 2 else 0.0
@@ -92,6 +129,22 @@ class SimObject:
         )
         self.control = ObjectControl()
 
+    @property
+    def surface_actuation(self):
+        return self._surface_actuation
+
+    def validate_surface_actuation(self, cfg, state=None):
+        if self.surface_actuation is None:
+            return
+        from .geometry.circle import CircleGeometry
+        if (type(self.geometry) is not CircleGeometry or self.is_sensor or self.is_force_region
+                or cfg is None or not cfg.is_d2q9 or cfg.store_precision != "FP32"
+                or cfg.streaming != "double_buffer" or cfg.ddf_shifting):
+            raise ValueError("Surface actuation requires a circle and D2Q9/FP32/unshifted/double_buffer")
+        values = vars(self.state) if state is None else state
+        if any(values.get(k, 0.) != 0 for k in ("vx", "vy", "omega")):
+            raise ValueError("Surface actuation requires a stationary circle")
+
     @property
     def geometry(self) -> Geometry:
         if self._geometry is None:

src/CelerisLab/body/preprocess/flag_overlay.py

@@ -9,6 +9,8 @@ from typing import List
 
 import numpy as np
 
+from ...lbm.descriptors import FLUID, SOLID, SENSOR_FLAG, FRC_REGION
+
 
 def merge_flag_masks(
     base: np.ndarray,
@@ -16,9 +18,9 @@ def merge_flag_masks(
 ) -> np.ndarray:
     """Merge object flag masks onto a clean domain base.
 
-    For each mask in *masks*, non-zero entries overwrite the corresponding
+    Apply solid occupancy first, then sensor/force-region annotations only
-    positions in *base*.  Callers should pass a fresh channel-layout flag
+    to fluid cells. Object order cannot turn an obstacle or channel wall into
-    array to prevent stale obstacle bits from surviving after geometry edits.
+    fluid. Callers supply a fresh channel-layout base after geometry edits.
 
     Args:
         base: uint16 domain flag array (e.g. from build_channel_flags()).
@@ -29,6 +31,10 @@ def merge_flag_masks(
     """
     merged = np.array(base, copy=True)
     for mask in masks:
-        nonzero = mask != 0
+        solid = (mask & SOLID) != 0
-        merged[nonzero] = mask[nonzero]
+        merged[solid] = mask[solid]
+    fluid = ((merged & FLUID) != 0) & ((merged & SOLID) == 0)
+    for mask in masks:
+        annotations = mask & (SENSOR_FLAG | FRC_REGION)
+        merged[fluid] |= annotations[fluid]
     return merged

src/CelerisLab/common/__init__.py

@@ -1,8 +1,12 @@
-# CelerisLab/common/__init__.py
+"""Shared host utilities; optional GPU/plotting dependencies load on demand."""
-"""
+from importlib import import_module
-Common utilities and preprocessing functions.
-"""
-
-from . import preprocess, pathline, checkpoint, render, streakline
 
 __all__ = ["preprocess", "pathline", "checkpoint", "render", "streakline"]
+
+
+def __getattr__(name):
+    if name not in __all__:
+        raise AttributeError(name)
+    value = import_module("." + name, __name__)
+    globals()[name] = value
+    return value

src/CelerisLab/common/checkpoint.py

@@ -12,15 +12,19 @@ Usage::
     ``lattice_model``, ``nx``, ``ny``, ``nz``, and ``store_precision``.
 """
 
+import hashlib
 import json
+from dataclasses import asdict
 import os
 import time
 
 import numpy as np
 import pycuda.driver as cuda
 
+from ..lbm.descriptors import FLUID, SOLID, OBSTACLE, SENSOR_FLAG, FRC_REGION
+
 # File format version — bump on incompatible layout changes.
-CKPT_VERSION = 1
+CKPT_VERSION = 4
 _SUPPORTED_STORE_PRECISIONS = ("FP32", "FP16S")
 
 
@@ -33,6 +37,66 @@ def _require_supported_precision(store_precision: str, source: str) -> None:
         )
 
 
+def _domain_flags(field, bodies, flags):
+    """Keep actual domain boundaries, removing only registered object overlays.
+
+    A solid object's hidden underlay is unknowable from merged flags: use FLUID
+    there canonically. This base is for fixed-geometry restart, not body removal.
+    """
+    base = flags & np.uint16(0xffff ^ (SENSOR_FLAG | FRC_REGION))
+    for obj in bodies.objects:
+        mask = obj.get_flag_mask(field.nx, field.ny, field.nz)
+        base[(mask & SOLID) != 0] = FLUID
+    if np.any(base & OBSTACLE):
+        raise ValueError("Checkpoint flags contain an unregistered obstacle")
+    return base
+
+
+def _geometry_records(field, bodies, base_flags=None):
+    """Describe actual geometry/projection, not just the number of objects."""
+    records = []
+    flags = bodies.build_flags(field.build_channel_flags() if base_flags is None else base_flags)
+    for obj in bodies.objects:
+        digest = hashlib.sha256()
+        digest.update(obj.get_flag_mask(field.nx, field.ny, field.nz).tobytes())
+        if obj.is_sensor or obj.is_force_region:
+            arrays = obj.get_sensor_list(field.nx, field.ny, field.nz)
+        else:
+            arrays = obj.get_curved_list(field.nx, field.ny, field.nq,
+                                         nz=field.nz, domain_flags=flags)
+        for array in arrays:
+            digest.update(np.ascontiguousarray(array).tobytes())
+        geometry = obj.geometry
+        records.append({
+            "obj_id": obj.obj_id,
+            "geometry": type(geometry).__module__ + "." + type(geometry).__qualname__,
+            "center": list(obj.center), "radius": float(obj.radius),
+            "is_sensor": obj.is_sensor, "is_force_region": obj.is_force_region,
+            "projection_sha256": digest.hexdigest(),
+        })
+    return records, flags
+
+
+def _validate_legacy_objects(obj_data, bodies):
+    """v1 can safely restart built-in circles/sensors with unchanged geometry.
+
+    It did not save authoritative force-region controls. Refuse those rather
+    than silently restoring zero force. Unknown/custom geometry is not inferred.
+    """
+    from ..body.geometry.circle import CircleGeometry
+
+    for od, obj in zip(obj_data, bodies.objects):
+        if obj.is_force_region:
+            raise ValueError("Legacy v1 checkpoint cannot restore force_region controls")
+        if type(obj.geometry) is not CircleGeometry:
+            raise ValueError("Legacy v1 checkpoint only supports built-in circle/sensor geometry")
+        expected = (float(od["x"]), float(od["y"]), float(od.get("z", 0.0)))
+        actual = (*obj.geometry.center, float(obj.state.z))
+        if (expected != actual or float(od["radius"]) != obj.geometry._radius
+                or od.get("type") != obj.obj_type):
+            raise ValueError("Checkpoint object geometry/type mismatch")
+
+
 def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None):
     """Save full simulation state to an HDF5 file.
 
@@ -61,10 +125,25 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None):
     cuda.memcpy_dtoh(tmp_buf, field.temp_gpu)
     flag_buf = np.empty(field.n, dtype=np.uint16)
     cuda.memcpy_dtoh(flag_buf, field.flag_gpu)
+    domain_flags = _domain_flags(field, bodies, flag_buf)
+    geometry, expected_flags = _geometry_records(field, bodies, domain_flags)
+    if not np.array_equal(flag_buf, expected_flags):
+        raise ValueError("Checkpoint flags do not match current object/domain geometry")
+    bodies.build_compact_lists_for(bodies.objects, domain_flags=expected_flags)
+    applied_action = None
+    if bodies.action_gpu is not None:
+        applied_action = np.empty_like(bodies.action)
+        cuda.memcpy_dtoh(applied_action, bodies.action_gpu)
+        for obj in bodies.objects:
+            if (obj.surface_actuation is not None
+                    and applied_action[(obj.obj_id + 1) * 3 * lbm_cfg.dim - 1] != 0):
+                raise ValueError("Surface actuation requires zero applied omega")
 
     with h5py.File(path, "w") as hf:
         # Version and metadata
-        hf.attrs["ckpt_version"] = CKPT_VERSION
+        # Older readers must not silently restart a driven surface as a wall.
+        hf.attrs["ckpt_version"] = (CKPT_VERSION if any(
+            obj.surface_actuation is not None for obj in bodies.objects) else 3)
         hf.attrs["celerislab_version"] = __import__(
             "CelerisLab").__version__
         hf.attrs["timestamp"] = time.strftime("%Y-%m-%dT%H:%M:%S")
@@ -74,25 +153,7 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None):
         hf.attrs["n_objects"] = bodies.count
 
         # Config as JSON string (human-readable in HDF5 viewers)
-        cfg = lbm_cfg
+        hf.attrs["config_json"] = json.dumps(asdict(lbm_cfg))
-        config_dict = {
-            "lattice_model": cfg.lattice_model,
-            "nx": cfg.nx, "ny": cfg.ny, "nz": cfg.nz,
-            "data_type": cfg.data_type,
-            "viscosity": cfg.viscosity, "velocity": cfg.velocity,
-            "rho": cfg.rho,
-            "collision": cfg.collision, "streaming": cfg.streaming,
-            "store_precision": cfg.store_precision,
-            "ddf_shifting": cfg.ddf_shifting,
-            "les_enabled": cfg.les_enabled, "les_cs": cfg.les_cs,
-            "les_closed_form": cfg.les_closed_form,
-            "trt_magic_param": cfg.trt_magic_param,
-            "inlet_profile": cfg.inlet_profile,
-            "inlet_scheme": cfg.inlet_scheme,
-            "outlet_mode": cfg.outlet_mode,
-            "omega_min": cfg.omega_min, "omega_max": cfg.omega_max,
-        }
-        hf.attrs["config_json"] = json.dumps(config_dict)
 
         # GPU arrays — gzip compression for portable, compact files
         hf.create_dataset("ddf", data=ddf_buf, compression="gzip",
@@ -101,8 +162,24 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None):
                           compression_opts=4)
         hf.create_dataset("flag", data=flag_buf, compression="gzip",
                           compression_opts=4)
+        hf.create_dataset("domain_flag", data=domain_flags, compression="gzip",
+                          compression_opts=4)
 
-        # Object states
+        hf.attrs["geometry_json"] = json.dumps(geometry)
+        hf.attrs["surface_actuation_json"] = json.dumps([
+            asdict(obj.surface_actuation) if obj.surface_actuation is not None else None
+            for obj in bodies.objects
+        ])
+        hf.attrs["runtime_params_json"] = json.dumps(field._runtime_params)
+        hf.attrs["obs_accum_steps"] = bodies._obs_accum_steps
+        if applied_action is not None:
+            hf.create_dataset("applied_action", data=applied_action)
+        if bodies.obs_gpu is not None:
+            obs = np.empty(bodies.obs_total_floats, dtype=np.float32)
+            cuda.memcpy_dtoh(obs, bodies.obs_gpu)
+            hf.create_dataset("obs", data=obs)
+
+        # Object states and host controls (next run uploads the restored action).
         obj_states = []
         for obj in bodies.objects:
             s = obj.state
@@ -112,6 +189,7 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None):
                 "x": s.x, "y": s.y, "z": s.z, "theta": s.theta,
                 "vx": s.vx, "vy": s.vy, "omega": s.omega,
                 "radius": obj.radius,
+                "control": asdict(obj.control),
             })
         hf.attrs["objects_json"] = json.dumps(obj_states)
 
@@ -124,6 +202,9 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies):
     The current config's lattice_model, nx, ny, nz, and
     store_precision must match the checkpoint.  Other parameters
     (collision, LES, …) may differ for parameter-study restarts.
+    v3 restores the saved domain boundaries; v1/v2 retain their default-channel
+    validation. v4 also requires matching static surface actuation. Object
+    geometry and roles must match in every version.
 
     Args:
         path: HDF5 checkpoint file path.
@@ -138,7 +219,7 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies):
     with h5py.File(path, "r") as hf:
         # Version check
         ver = hf.attrs.get("ckpt_version", 0)
-        if ver != CKPT_VERSION:
+        if ver not in (1, 2, 3, CKPT_VERSION):
             raise ValueError(
                 f"Checkpoint version mismatch: file={ver}, "
                 f"expected={CKPT_VERSION}. "
@@ -146,6 +227,16 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies):
 
         # Config compatibility check
         saved_cfg = json.loads(hf.attrs["config_json"])
+        if "on_site_velocity" in (saved_cfg.get("y_wall_bc"), lbm_cfg.y_wall_bc):
+            for key in ("y_wall_bc", "velocity", "wall_normal_velocity"):
+                if saved_cfg.get(key) != getattr(lbm_cfg, key):
+                    raise ValueError(f"Checkpoint on-site boundary mismatch on '{key}'")
+            if saved_cfg.get("wall_slot_width", 0) != lbm_cfg.wall_slot_width:
+                raise ValueError("Checkpoint on-site boundary mismatch on 'wall_slot_width'")
+            if saved_cfg.get("top_wall_bc", "velocity") != lbm_cfg.top_wall_bc:
+                raise ValueError("Checkpoint on-site boundary mismatch on 'top_wall_bc'")
+            if lbm_cfg.top_wall_bc == "reservoir" and saved_cfg.get("rho") != lbm_cfg.rho:
+                raise ValueError("Checkpoint reservoir boundary mismatch on 'rho'")
         for key in ("lattice_model", "nx", "ny", "nz", "store_precision",
                      "streaming", "ddf_shifting"):
             saved_val = saved_cfg.get(key)
@@ -170,31 +261,112 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies):
                 "the object count before loading, or rebuild the "
                 "simulation with the correct number of objects.")
 
-        # Restore GPU arrays
+        # Read and validate everything before modifying host/device state.
+        obj_data = json.loads(hf.attrs.get("objects_json", "[]"))
+        if (len(obj_data) != bodies.count
+                or [od.get("obj_id") for od in obj_data] != list(range(bodies.count))):
+            raise ValueError("Checkpoint object IDs must match the dense current registry")
+        surface_specs = [asdict(obj.surface_actuation) if obj.surface_actuation is not None else None
+                         for obj in bodies.objects]
+        if ver >= 4 and "surface_actuation_json" not in hf.attrs:
+            raise ValueError("Checkpoint is missing surface actuation specification")
+        saved_surface_specs = json.loads(hf.attrs.get(
+            "surface_actuation_json", json.dumps([None] * bodies.count)))
+        if saved_surface_specs != surface_specs:
+            raise ValueError("Checkpoint surface actuation mismatch")
+        domain_flags = None
+        if ver >= 3:
+            if "domain_flag" not in hf:
+                raise ValueError("Checkpoint is missing domain_flag")
+            domain_flags = hf["domain_flag"][:]
+            if domain_flags.shape != (field.n,) or domain_flags.dtype != np.uint16:
+                raise ValueError("Invalid checkpoint domain_flag shape/dtype")
+        geometry, expected_flags = _geometry_records(field, bodies, domain_flags)
+        if ver == 1:
+            _validate_legacy_objects(obj_data, bodies)
+        elif json.loads(hf.attrs["geometry_json"]) != geometry:
+            raise ValueError("Checkpoint object geometry/role mismatch; recreate the matching case")
+
+        dtype = np.dtype(np.float16 if field.store_bytes == 2 else np.float32)
         ddf_buf = hf["ddf"][:]
         tmp_buf = hf["temp"][:]
         flag_buf = hf["flag"][:]
+        for name, array, shape, expected_dtype in (
+            ("ddf", ddf_buf, (field.n * field.nq,), dtype),
+            ("temp", tmp_buf, (field.n * field.nq,), dtype),
+            ("flag", flag_buf, (field.n,), np.dtype(np.uint16)),
+        ):
+            if array.shape != shape or array.dtype != expected_dtype:
+                raise ValueError(f"Invalid checkpoint {name} shape/dtype")
+        if not np.array_equal(flag_buf, expected_flags):
+            raise ValueError("Checkpoint flags do not match current object/domain geometry")
+        if ver >= 3 and not np.array_equal(domain_flags, _domain_flags(field, bodies, flag_buf)):
+            raise ValueError("Checkpoint domain_flag is not a canonical object-free domain")
+        restored_step = int(hf.attrs["step_count"])
+        if restored_step < 0:
+            raise ValueError("Checkpoint step_count must be nonnegative")
+        raw_params = json.loads(hf.attrs.get("runtime_params_json", "{}"))
+        params = {key: float(raw_params[key]) for key in ("omega", "fx", "fy", "fz")
+                  if key in raw_params}
+        if not all(np.isfinite(value) for value in params.values()):
+            raise ValueError("Checkpoint runtime parameters must be finite")
+        if lbm_cfg.y_wall_bc == "on_site_velocity" and (
+                bodies.count or any(params.get(k, 0.0) != 0 for k in ("fx", "fy", "fz"))):
+            raise ValueError("on_site_velocity does not support body forces or objects")
+        obs_accum_steps = int(hf.attrs.get("obs_accum_steps", 0))
+        if obs_accum_steps < 0:
+            raise ValueError("Checkpoint obs_accum_steps must be nonnegative")
+        # Validate compact ownership before restoring any state or GPU arrays.
+        bodies.build_compact_lists_for(bodies.objects, domain_flags=expected_flags)
+        # A changed viscosity is an intentional parameter-study restart.
+        if saved_cfg.get("viscosity") != lbm_cfg.viscosity:
+            params["omega"] = float(lbm_cfg.omega)
+        params["n_objects"] = bodies.count
+        applied_action = hf["applied_action"][:] if "applied_action" in hf and ver >= 2 else None
+        if applied_action is not None and (
+                applied_action.shape != bodies.action.shape or applied_action.dtype != np.float32):
+            raise ValueError("Invalid checkpoint applied_action shape/dtype")
+        obs = hf["obs"][:] if "obs" in hf else None
+        if obs is not None:
+            from ..cuda.obs_layout import obs_layout
+            size = obs_layout(lbm_cfg.dim, bodies.count).total_floats
+            if obs.shape != (size,) or obs.dtype != np.float32:
+                raise ValueError("Invalid checkpoint obs shape/dtype")
+        state_fields = ("x", "y", "z", "theta", "vx", "vy", "omega")
+        states = [{key: float(od.get(key, 0.0)) for key in state_fields}
+                  for od in obj_data]
+        controls = [{key: float(od.get("control", {}).get(key, 0.0))
+                     for key in ("ax", "ay", "alpha")} for od in obj_data]
+        modes = [int(od.get("control", {}).get("mode", 0)) for od in obj_data]
+        if not all(np.isfinite(value) for record in states + controls for value in record.values()):
+            raise ValueError("Checkpoint object state/control must be finite")
+        for obj, state in zip(bodies.objects, states):
+            obj.validate_surface_actuation(lbm_cfg, state)
+            if (obj.surface_actuation is not None and applied_action is not None
+                    and applied_action[(obj.obj_id + 1) * 3 * lbm_cfg.dim - 1] != 0):
+                raise ValueError("Checkpoint surface actuation requires zero applied omega")
 
+        for obj, state, control, mode in zip(bodies.objects, states, controls, modes):
+            for key, value in state.items():
+                setattr(obj.state, key, value)
+            for key, value in control.items():
+                setattr(obj.control, key, value)
+            obj.control.mode = mode
+
+        field.flag = flag_buf.copy()
+        # Rebuild compact lists and action from the validated matching geometry.
+        # This also supports loading into a constructed but uninitialized case.
+        bodies.sync_to_gpu(field, rebuild_flags=False)
+        cuda.memcpy_htod(bodies.action_gpu,
+                         bodies.action if applied_action is None else applied_action)
         cuda.memcpy_htod(field.ddf_gpu, ddf_buf)
         cuda.memcpy_htod(field.temp_gpu, tmp_buf)
         cuda.memcpy_htod(field.flag_gpu, flag_buf)
-        field.flag = flag_buf.copy()
+        field.update_params(**params)
-
+        stepper._step_count = restored_step
-        # Restore step count
+        field.invalidate_host_ddf_cache()
-        stepper._step_count = int(hf.attrs["step_count"])
+        if obs is not None:
-
+            bodies.obs_pinned[:] = obs
-        # Restore object states
+            cuda.memcpy_htod(bodies.obs_gpu, obs)
-        obj_data = json.loads(hf.attrs.get("objects_json", "[]"))
+            bodies._obs_accum_steps = obs_accum_steps
-        for od in obj_data:
+        # v1 lacks observations; sync_to_gpu has explicitly reset the interval.
-            oid = od["obj_id"]
-            try:
-                obj = bodies.get(oid)
-                obj.state.x = od["x"]
-                obj.state.y = od["y"]
-                obj.state.z = float(od.get("z", 0.0))
-                obj.state.theta = od["theta"]
-                obj.state.vx = od["vx"]
-                obj.state.vy = od["vy"]
-                obj.state.omega = od["omega"]
-            except KeyError:
-                pass  # object not in current sim — skip

src/CelerisLab/config.py

@@ -12,6 +12,7 @@ macro dicts consumed by `cuda/compiler_v2.py` for header generation.
 
 import json
 import os
+import math
 from dataclasses import dataclass, field
 from typing import Any, Dict, List, Optional
 
@@ -33,10 +34,16 @@ INLET_SCHEME_MAP = {
     "regularized": 3,
 }
 OUTLET_MAP = {"neq_extrap": 0, "zero_gradient": 1, "blended": 2}
-Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1}
+Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1, "on_site_velocity": 2}
+TOP_WALL_BC_MAP = {"velocity": 0, "reservoir": 1}
 DTYPE_MAP = {"FP32": "float", "FP64": "double"}
 
 
+def _float_literal(value: float) -> str:
+    """Emit a round-trip numeric literal; display precision must not alter physics."""
+    return repr(float(value)) + "f"
+
+
 # ---------------------------------------------------------------------------
 # LBM config
 # ---------------------------------------------------------------------------
@@ -80,6 +87,9 @@ class LBMConfig:
     inlet_regularized_neq_damp: float = 0.5
     outlet_srt_neq_damp: float = 0.5
     y_wall_bc: str = "bounce_back"
+    wall_normal_velocity: float = 0.0
+    wall_slot_width: int = 0
+    top_wall_bc: str = "velocity"
     omega_min: float = 0.01
     omega_max: float = 1.99
 
@@ -101,21 +111,38 @@ class LBMConfig:
         return self.lattice_model == "D3Q19"
 
     def validate(self):
-        assert self.lattice_model in LATTICE_MODEL_MAP, (
+        for name, allowed in (
-            f"Unknown lattice_model: {self.lattice_model}. "
+            ("lattice_model", LATTICE_MODEL_MAP), ("collision", COLLISION_MAP),
-            f"Supported: {list(LATTICE_MODEL_MAP.keys())}")
+            ("streaming", STREAMING_MAP), ("y_wall_bc", Y_WALL_BC_MAP),
+            ("top_wall_bc", TOP_WALL_BC_MAP),
+            ("data_type", DTYPE_MAP), ("store_precision", PRECISION_MAP),
+            ("inlet_profile", INLET_MAP), ("inlet_scheme", INLET_SCHEME_MAP),
+            ("outlet_mode", OUTLET_MAP),
+        ):
+            if getattr(self, name) not in allowed:
+                raise ValueError(f"Unknown {name}: {getattr(self, name)!r}")
         lm = LATTICE_MODEL_MAP[self.lattice_model]
-        self.dim = lm["dim"]
+        self.dim, self.nq, self.lattice_model_id = lm["dim"], lm["nq"], lm["id"]
-        self.nq = lm["nq"]
+        for name in ("nx", "ny", "nz", "threads_per_block"):
-        self.lattice_model_id = lm["id"]
+            value = getattr(self, name)
-        assert self.collision in COLLISION_MAP, f"Unknown collision: {self.collision}"
+            if isinstance(value, bool) or not isinstance(value, int) or value <= 0:
-        assert self.streaming in STREAMING_MAP, f"Unknown streaming: {self.streaming}"
+                raise ValueError(f"{name} must be a positive integer.")
-        assert self.y_wall_bc in Y_WALL_BC_MAP, f"Unknown y_wall_bc: {self.y_wall_bc}"
+        if self.threads_per_block > 1024:
-        assert self.data_type in DTYPE_MAP, f"Unknown data_type: {self.data_type}"
+            raise ValueError("threads_per_block must not exceed 1024.")
-        assert self.store_precision in PRECISION_MAP, (
+        for name in ("viscosity", "velocity", "rho", "les_cs", "trt_magic_param",
-            f"Unknown store_precision: {self.store_precision}")
+                     "outlet_blend_alpha", "inlet_trt_neq_damp",
-        assert self.inlet_scheme in INLET_SCHEME_MAP, (
+                     "inlet_regularized_neq_damp", "outlet_srt_neq_damp",
-            f"Unknown inlet_scheme: {self.inlet_scheme}")
+                     "omega_min", "omega_max", "wall_normal_velocity"):
+            if not math.isfinite(getattr(self, name)):
+                raise ValueError(f"{name} must be finite.")
+        if self.viscosity <= 0 or self.rho <= 0:
+            raise ValueError("viscosity and rho must be positive.")
+        if self.les_cs < 0 or self.trt_magic_param <= 0:
+            raise ValueError("les_cs must be nonnegative and trt_magic_param positive.")
+        if not 0 <= self.outlet_blend_alpha <= 1:
+            raise ValueError("outlet_blend_alpha must lie in [0, 1].")
+        if self.omega_min <= 0:
+            raise ValueError("omega_min must be positive.")
         if self.store_precision == "FP16C":
             raise ValueError(
                 "store_precision='FP16C' is not supported in the current runtime path. "
@@ -126,9 +153,22 @@ class LBMConfig:
                 "Only data_type='FP32' is supported: CUDA kernels and LBMField "
                 "paths use float32; LBtype is not wired through device code yet."
             )
-        assert self.nx > 0 and self.ny > 0 and self.nz > 0
+        if self.dim == 2 and self.nz != 1:
-        if self.dim == 2:
+            raise ValueError("nz must be 1 for 2D")
-            assert self.nz == 1, "nz must be 1 for 2D"
+        if (isinstance(self.wall_slot_width, bool) or not isinstance(self.wall_slot_width, int)
+                or self.wall_slot_width != 0 and not 2 <= self.wall_slot_width <= self.nx):
+            raise ValueError("wall_slot_width must be 0 (disabled) or an integer in [2, nx]")
+        if self.wall_slot_width and (self.y_wall_bc != "on_site_velocity" or self.velocity != 0):
+            raise ValueError("wall_slot_width requires on_site_velocity and velocity=0")
+        if self.top_wall_bc == "reservoir" and (self.y_wall_bc != "on_site_velocity" or self.velocity != 0):
+            raise ValueError("reservoir requires on_site_velocity and velocity=0")
+        if self.y_wall_bc == "on_site_velocity":
+            if (not self.is_d2q9 or self.collision != "SRT"
+                    or self.streaming != "double_buffer" or self.store_precision != "FP32"
+                    or self.ddf_shifting or self.les_enabled):
+                raise ValueError("on_site_velocity requires D2Q9/SRT/FP32/unshifted/double_buffer without LES")
+            if abs(self.wall_normal_velocity) >= 1:
+                raise ValueError("on_site_velocity requires abs(wall_normal_velocity) < 1")
         if not (0.0 <= self.inlet_trt_neq_damp <= 1.0):
             raise ValueError("method.inlet.trt_neq_damp must lie in [0, 1].")
         if not (0.0 <= self.inlet_regularized_neq_damp <= 1.0):
@@ -150,6 +190,36 @@ class LBMConfig:
                 f"[{self.omega_min:.6f}, {self.omega_max:.6f}]."
             )
 
+    def validate_runtime(self):
+        """Validate the currently complete Simulation path, not local operators."""
+        self.validate()
+        if self.is_d3q19:
+            raise NotImplementedError(
+                "D3Q19 Simulation is not supported yet: boundary/auxiliary kernels "
+                "and object projection are not dimension-complete. D3Q19 local "
+                "descriptors and collision operators remain available for research."
+            )
+        if self.nx < 3 or self.ny < 4:
+            raise ValueError("Simulation requires nx >= 3 and ny >= 4: x ghost layers and distinct wall-adjacent rows.")
+
+    def allocation_signature(self) -> tuple:
+        """Parameters that cannot change while field allocations/state are retained."""
+        self.validate()
+        return (self.lattice_model, self.nx, self.ny, self.nz, self.data_type,
+                self.store_precision, bool(self.ddf_shifting), self.streaming,
+                (self.velocity, self.wall_normal_velocity, self.wall_slot_width,
+                 self.top_wall_bc, self.rho if self.top_wall_bc == "reservoir" else None)
+                if self.y_wall_bc == "on_site_velocity" else None)
+
+    @property
+    def top_wall_normal_velocity(self) -> float:
+        """Volume-balanced top velocity target, or reservoir initialization guess."""
+        if self.wall_slot_width in (0, self.nx):
+            return self.wall_normal_velocity
+        weights = (max(0.0, 1.0-((2*x-(self.nx-1))/self.wall_slot_width)**2)**2
+                   for x in range(self.nx))
+        return self.wall_normal_velocity * math.fsum(weights) / self.nx
+
     def to_macros(self) -> Dict[str, Any]:
         """Return flat dict of macro_name → value for header generation."""
         self.validate()
@@ -161,28 +231,32 @@ class LBMConfig:
             "NZ": self.nz,
             "LATTICE_MODEL": self.lattice_model_id,
             "LBtype": DTYPE_MAP[self.data_type],
-            "VIS": f"{self.viscosity:.10f}",
+            "VIS": repr(float(self.viscosity)),
-            "RHO": f"{self.rho:.1f}",
+            "RHO": repr(float(self.rho)),
-            "U0": f"{self.velocity}",
+            "U0": repr(float(self.velocity)),
             "COLLISION_MODEL": COLLISION_MAP[self.collision],
             "STREAMING_MODEL": STREAMING_MAP[self.streaming],
             "STORE_PRECISION": PRECISION_MAP[self.store_precision],
             "USE_DDF_SHIFTING": int(self.ddf_shifting),
             "USE_LES": int(self.les_enabled),
-            "LES_CS": f"{self.les_cs:.6f}f",
+            "LES_CS": _float_literal(self.les_cs),
             "LES_CLOSED_FORM": int(self.les_closed_form),
             "INLET_PROFILE": INLET_MAP[self.inlet_profile],
             "INLET_SCHEME": INLET_SCHEME_MAP[self.inlet_scheme],
             "OUTLET_MODE": OUTLET_MAP[self.outlet_mode],
-            "OUTLET_BLEND_ALPHA": f"{self.outlet_blend_alpha:.3f}f",
+            "OUTLET_BLEND_ALPHA": _float_literal(self.outlet_blend_alpha),
             "OUTLET_BACKFLOW_CLAMP": int(self.outlet_backflow_clamp),
             "Y_WALL_BC": Y_WALL_BC_MAP[self.y_wall_bc],
-            "OMEGA_COLLISION_MIN": f"{self.omega_min:.2f}f",
+            "WALL_NORMAL_VELOCITY": _float_literal(self.wall_normal_velocity),
-            "OMEGA_COLLISION_MAX": f"{self.omega_max:.3f}f",
+            "WALL_SLOT_WIDTH": self.wall_slot_width,
-            "TRT_MAGIC_PARAM": f"{self.trt_magic_param:.6f}f",
+            "TOP_WALL_NORMAL_VELOCITY": _float_literal(self.top_wall_normal_velocity),
-            "INLET_TRT_NEQ_DAMP": f"{self.inlet_trt_neq_damp:.4f}f",
+            "TOP_WALL_BC": TOP_WALL_BC_MAP[self.top_wall_bc],
-            "INLET_REGULARIZED_NEQ_DAMP": f"{self.inlet_regularized_neq_damp:.4f}f",
+            "OMEGA_COLLISION_MIN": _float_literal(self.omega_min),
-            "OUTLET_SRT_NEQ_DAMP": f"{self.outlet_srt_neq_damp:.4f}f",
+            "OMEGA_COLLISION_MAX": _float_literal(self.omega_max),
+            "TRT_MAGIC_PARAM": _float_literal(self.trt_magic_param),
+            "INLET_TRT_NEQ_DAMP": _float_literal(self.inlet_trt_neq_damp),
+            "INLET_REGULARIZED_NEQ_DAMP": _float_literal(self.inlet_regularized_neq_damp),
+            "OUTLET_SRT_NEQ_DAMP": _float_literal(self.outlet_srt_neq_damp),
         }
 
 
@@ -200,8 +274,11 @@ class BodyConfig:
 def _find_config(filename: str, explicit_path: Optional[str] = None) -> str:
     """Search for config file in standard locations."""
     candidates = []
-    if explicit_path:
+    if explicit_path is not None:
-        candidates.append(explicit_path)
+        explicit_path = os.fspath(explicit_path)
+        if not os.path.isfile(explicit_path):
+            raise FileNotFoundError(f"Explicit config path does not exist: {explicit_path}")
+        return os.path.abspath(explicit_path)
     env_dir = os.environ.get("CELERISLAB_CONFIG_DIR")
     if env_dir:
         candidates.append(os.path.join(env_dir, filename))
@@ -242,6 +319,9 @@ def load_lbm_config(path: Optional[str] = None) -> LBMConfig:
         inlet_regularized_neq_damp=m["inlet"].get("regularized_neq_damp", 0.5),
         outlet_srt_neq_damp=m["outlet"].get("srt_neq_damp", 0.5),
         y_wall_bc=m.get("y_wall_bc", "bounce_back"),
+        wall_normal_velocity=m.get("wall_normal_velocity", 0.0),
+        wall_slot_width=m.get("wall_slot_width", 0),
+        top_wall_bc=m.get("top_wall_bc", "velocity"),
         omega_min=m["omega_guard"]["min"],
         omega_max=m["omega_guard"]["max"],
         threads_per_block=c["threads_per_block"],
@@ -255,4 +335,4 @@ def load_body_config(path: Optional[str] = None) -> BodyConfig:
     fpath = _find_config("config_body.json", path)
     with open(fpath) as f:
         d = json.load(f)
-    return BodyConfig(objects=d.get("objects", []))
+    return BodyConfig(objects=d.get("objects", []))

src/CelerisLab/configs/CONFIG.md

@@ -1,7 +1,7 @@
 # CelerisLab 配置说明
 
 **唯一配置入口：`config_lbm.json`**  
-Python `config.py` 只负责读取和校验，不是配置位置。
+Python `config.py` 只负责读取和校验，不是配置位置。显式配置路径不存在时会报错，不会退回默认配置；仅未指定路径时按环境目录、工作目录、包内配置顺序查找。
 
 ---
 
@@ -9,7 +9,7 @@ Python `config.py` 只负责读取和校验，不是配置位置。
 
 | 字段 | 类型 | 默认 | 允许值 | 说明 |
 |------|------|------|--------|------|
-| `lattice_model` | string | `"D2Q9"` | `D2Q9`, `D3Q19` | 格子模型，决定维度和速度数量 |
+| `lattice_model` | string | `"D2Q9"` | `D2Q9`, `D3Q19` | 格子模型，决定维度和速度数量；`Simulation` 当前仅支持 D2Q9。D3Q19 保留供 descriptor / 局部算子研究，完整运行时在创建 GPU 上下文前明确拒绝 |
 | `nx` | int | 512 | >0 | x 方向格点数（流向），建议整除 `threads_per_block` 以减少尾块浪费 |
 | `ny` | int | 256 | >0 | y 方向格点数 |
 | `nz` | int | 1 | >0 | z 方向格点数，D2Q9 时必须为 1 |
@@ -58,16 +58,16 @@ Python `config.py` 只负责读取和校验，不是配置位置。
 
 ## 配置 → 代码映射
 
-```
 **`lbm/kernels/config/*.h` 由编译器根据 `LBMConfig` 自动生成，请勿手改。**
 
-仓库中可能检入过历史版本的 `config/*.h`，与当前默认 JSON **不一定一致**；若以 Python 入口运行仿真，`Simulation` / `compiler_v2.generate_config` 会在编译前重写这些头文件。若只用 `nvcc` 单独编译 `kernel_v2.cu` 且跳过 `generate_config`，则可能读到陈旧宏，属不受支持用法。
+仓库中检入的 `config/*.h` 可能是历史版本。`compiler_v2.generate_config` 会复制 kernel 源码到独立临时构建目录，并在副本中生成配置；不会改写包内头文件或 PTX。返回目录须显式传给 `compile_kernel(build_dir=...)`，编译结果与配置摘要均保存在该目录。`Simulation` 负责目录清理；直接使用编译器的调用者须调用 `cleanup_build_dir()`。跳过生成配置直接编译仓库中的 `kernel_v2.cu` 属不受支持用法。
 
+```
 config_lbm.json
   ↓  load_lbm_config()
 config.py:LBMConfig       (读取 + 校验)
   ↓  to_macros()
-cuda/compiler_v2.py       (生成 config/*.h)
+cuda/compiler_v2.py       (独立构建目录中生成 config/*.h)
   ↓
 config_grid.h             → NX, NY, NZ, NT, LATTICE_MODEL；
                             DIM / NQ 由 LATTICE_MODEL 推导写入（与 grid.lattice_model 一致，非 JSON 顶层字段）
@@ -100,4 +100,66 @@ step/one_step_*.cu         → kernel 编排
 
 D2Q9/D3Q19 的 flag 扩展位 `0x0800` 定义为 `FRC_REGION`，用于标记 force_region 类型物体的格子。
 - ForceRegionKernel 不依赖 flag 位执行，flag 只用于后处理可视化。
-- `merge_flag_masks` 会正确合并 force_region 的 `FLUID | FRC_REGION` 标记。
+- `merge_flag_masks` 会正确合并 force_region 的 `FLUID | FRC_REGION` 标记。
+# J1 on-site velocity planes (limited reference case)
+
+Set `method.y_wall_bc` to `"on_site_velocity"` and
+`method.wall_normal_velocity` to V (default 0). This selects an x-periodic
+domain with physical planes through fluid rows y=1 and y=NY−2, H=NY−3.
+The bottom velocity is (0,V); the top is (`physics.velocity`,V).
+These are constant **Cartesian** velocities: positive V injects through the
+bottom and leaves through the top. The selected rows do not also receive
+halfway reflection. Initialization uses a linear tangential profile and V.
+
+Supported only for D2Q9/SRT/FP32/unshifted/double_buffer without LES, objects
+or body forces. Boundary targets and the mode are fixed at construction;
+checkpoint restart requires matching targets. With the slot disabled this is an
+entire-plane Couette-throughflow reference. The default
+`bounce_back` mode and existing configurations retain their previous path.
+
+### Static centered bottom slot
+
+Optional `method.wall_slot_width` is an integer W: 0 (default) keeps the above
+uniform planes; 2≤W≤NX enables a centered bottom slot and requires
+`physics.velocity=0`. The same support restrictions apply. Width and amplitude
+are fixed at construction and must match on checkpoint restart; older v3
+checkpoints without this field mean W=0.
+
+For W<NX, let c=(NX−1)/2, s=2(x−c)/W at integer x=0..NX−1. The bottom target is
+`v=A*max(0,1-s*s)^2`, where A=`wall_normal_velocity` is the continuous peak,
+not a renormalized discrete peak. This profile and its first derivative vanish
+at the edges. W=NX explicitly means uniform v=A, not this shaped profile.
+With equal unit area per node, requested Q is the sum of these bottom targets;
+the top target is uniform Q/NX (`cfg.top_wall_normal_velocity`). The sum is
+computed in host double precision; device targets are FP32. This balances
+requested volume, not actual density-weighted mass exactly.
+
+Positive A injects from the bottom and exits through the top; negative A reverses
+both roles. Both are Cartesian y velocities; inward-normal velocity has the
+opposite sign at the top. The slot and its zero-velocity flanks share the same
+on-site y=1 plane. Initialization sets rho0, u=0 and linearly interpolates v
+between the two planes; it is not an assumed steady solution. This mode is
+static, single-level and flat; it does not support moving/curved slots.
+The finite-slot mode is experimental: the bounded W=16, A=±0.002 checks in
+`research/j2-slot-2026-10-12/REVIEW.md` exceeded the fixed late-time mass-drift
+tolerance. It is not yet validated for steady mass balance or long runs.
+
+### Opt-in top reservoir
+
+Set `method.top_wall_bc="reservoir"` with `on_site_velocity` and
+`physics.velocity=0` to fix the top density to `physics.rho` and tangential
+velocity to zero. The top Cartesian normal velocity is solved locally from
+known poststream populations; positive means outflow, negative means inflow.
+It may vary along x and is **not** constrained to the bottom requested Q/NX.
+The bottom slot/flanks and the supported lattice/collision/storage scope stay
+the same. This changes the physical top boundary condition.
+
+The initializer still uses uniform Q/NX at the top as an initial guess, then
+the density reconstruction applies from the first step. Mode, reservoir density
+and tangential velocity are fixed at construction and checked on restart.
+Missing `top_wall_bc` in older configs/checkpoints means `"velocity"`, which
+remains the default and preserves the previous two-velocity-boundary path.
+This reservoir mode is experimental. The bounded positive-slot test in
+`research/j3-reservoir-2026-10-12/REVIEW.md` passed its mass thresholds but failed
+the frozen full-field velocity convergence criterion; negative/zero cases and
+GPU restart validation were therefore not run. Long-run stability is unverified.

src/CelerisLab/cuda/__init__.py

@@ -1,8 +1,4 @@
-# CelerisLab/cuda/__init__.py
+"""CUDA build helpers (compilation itself does not need a live context)."""
-"""
-CUDA compiler and kernel management utilities.
-"""
-
 from . import compiler_v2 as compiler
 
-__all__ = ['compiler']
+__all__ = ["compiler"]

src/CelerisLab/cuda/compiler_v2.py

@@ -3,17 +3,24 @@
 Kernel configuration, compilation, and module loading.
 
 Workflow:
-  1. generate_config(lbm_cfg, n_objects) — write config/*.h from LBMConfig
+  1. build_dir = generate_config(lbm_cfg, n_objects) — isolated source snapshot
-  2. compile_kernel(arch)               — nvcc -ptx kernel_v2.cu
+  2. compile_kernel(arch, build_dir=build_dir)       — nvcc -ptx kernel_v2.cu
   3. load_module(ptx_path)              — PyCUDA module from PTX
 """
 
+from __future__ import annotations
+
 import os
 import subprocess
+import shutil
+import tempfile
+import hashlib
+import json
+from contextvars import ContextVar
+from dataclasses import asdict
+from pathlib import Path
 from typing import Optional
 
-import pycuda.driver as cuda
-
 from ..config import LBMConfig
 from .obs_layout import obs_layout  # re-exported for body.manager import path
 
@@ -26,6 +33,10 @@ _KERNEL_DIR = os.path.join(
 )
 
 
+_last_build_dir = ContextVar("celeris_last_build_dir", default=None)
+_last_ptx_path = ContextVar("celeris_last_ptx_path", default=None)
+
+
 def kernel_path(filename: str) -> str:
     return os.path.join(_KERNEL_DIR, filename)
 
@@ -37,18 +48,60 @@ _HEADER = "// AUTO-GENERATED by CelerisLab compiler – DO NOT EDIT MANUALLY\n"
 
 
 def _write(path: str, content: str):
-    os.makedirs(os.path.dirname(path), exist_ok=True)
+    # Headers and metadata belong to this caller's private build directory.
-    tmp_path = f"{path}.tmp"
+    Path(path).parent.mkdir(parents=True, exist_ok=True)
-    with open(tmp_path, "w") as f:
+    Path(path).write_text(content, encoding="utf-8")
-        f.write(content)
-    os.replace(tmp_path, path)
 
 
-def generate_config(cfg: LBMConfig, n_objects: int = 0):
+def generate_config(cfg: LBMConfig, n_objects: int = 0,
-    """Render all config/*.h files from *cfg*."""
+                    *, build_root: Optional[str] = None) -> str:
-    m = cfg.to_macros()
+    """Snapshot kernels and render headers in a fresh private directory.
-    cdir = os.path.join(_KERNEL_DIR, "config")
 
+    Returns its path. Pass it explicitly to compile_kernel for interleaved
+    builds. Legacy generate_config(); compile_kernel() remains context-local.
+    The caller owns cleanup_build_dir() once the loaded module is no longer
+    needed for recompilation/debugging. No package files are changed.
+    """
+    macros = cfg.to_macros()
+    if isinstance(n_objects, bool) or not isinstance(n_objects, int) or n_objects < 0:
+        raise ValueError("n_objects must be a nonnegative integer")
+    if cfg.y_wall_bc == "on_site_velocity" and n_objects:
+        raise ValueError("on_site_velocity currently requires an object-free domain")
+    build_dir = tempfile.mkdtemp(prefix="celeris-build-", dir=build_root)
+    try:
+        shutil.copytree(_KERNEL_DIR, build_dir, dirs_exist_ok=True,
+                        ignore=shutil.ignore_patterns("*.ptx", "*.cubin", "*.tmp"))
+        _render_config(macros, n_objects, os.path.join(build_dir, "config"))
+        digest = hashlib.sha256()
+        for path in sorted(Path(build_dir).rglob("*")):
+            if path.is_file():
+                digest.update(str(path.relative_to(build_dir)).encode())
+                digest.update(path.read_bytes())
+        metadata = {"source_sha256": digest.hexdigest(), "config": asdict(cfg),
+                    "n_objects": n_objects}
+        _write(os.path.join(build_dir, ".celeris-build"), json.dumps(metadata, sort_keys=True))
+    except BaseException:
+        shutil.rmtree(build_dir)
+        raise
+    _last_build_dir.set(build_dir)
+    _last_ptx_path.set(None)
+    return build_dir
+
+
+def cleanup_build_dir(build_dir: Optional[str]) -> None:
+    """Remove only a directory created by generate_config; safe to repeat."""
+    if build_dir is None or not os.path.exists(build_dir):
+        return
+    path = Path(build_dir).resolve()
+    if path == Path(_KERNEL_DIR).resolve() or not (path / ".celeris-build").is_file():
+        raise ValueError("Refusing to remove a directory not owned by the compiler")
+    shutil.rmtree(path)
+    if _last_build_dir.get() == str(build_dir):
+        _last_build_dir.set(None)
+        _last_ptx_path.set(None)
+
+
+def _render_config(m: dict, n_objects: int, cdir: str):
     _write(os.path.join(cdir, "config_grid.h"), f"""{_HEADER}\
 // Layer: Global / Grid
 #ifndef CELERIS_CONFIG_GRID_H
@@ -122,6 +175,10 @@ def generate_config(cfg: LBMConfig, n_objects: int = 0):
 #define OUTLET_BLEND_ALPHA   {m['OUTLET_BLEND_ALPHA']}
 #define OUTLET_BACKFLOW_CLAMP {m['OUTLET_BACKFLOW_CLAMP']}
 #define Y_WALL_BC            {m['Y_WALL_BC']}
+#define WALL_NORMAL_VELOCITY {m['WALL_NORMAL_VELOCITY']}
+#define WALL_SLOT_WIDTH      {m['WALL_SLOT_WIDTH']}
+#define TOP_WALL_NORMAL_VELOCITY {m['TOP_WALL_NORMAL_VELOCITY']}
+#define TOP_WALL_BC          {m['TOP_WALL_BC']}
 
 #define OMEGA_COLLISION_MIN {m['OMEGA_COLLISION_MIN']}
 #define OMEGA_COLLISION_MAX {m['OMEGA_COLLISION_MAX']}
@@ -184,36 +241,47 @@ def generate_config(cfg: LBMConfig, n_objects: int = 0):
 # ---------------------------------------------------------------------------
 # Compilation
 # ---------------------------------------------------------------------------
-def compile_kernel(arch: str = "sm_70",
+def compile_kernel(arch: str = "sm_70", source: str = "kernel_v2.cu",
-                   source: str = "kernel_v2.cu",
+                   output: Optional[str] = None, *,
-                   output: Optional[str] = None) -> str:
+                   build_dir: Optional[str] = None) -> str:
-    """Compile *source* to PTX.  Returns path to the .ptx file.
+    """Compile a generated private source snapshot; never use package headers."""
-
+    # A failed attempt must not make implicit load_module() reuse older code.
-    The translation unit includes ``config/*.h``. Those headers must match the
+    _last_ptx_path.set(None)
-    intended ``LBMConfig`` — call :func:`generate_config` before compiling.
+    build_dir = build_dir or _last_build_dir.get()
-    Invoking ``nvcc`` on ``kernel_v2.cu`` alone can pick up stale checked-in
+    if not build_dir or not (Path(build_dir) / ".celeris-build").is_file():
-    headers from ``lbm/kernels/config/``.
+        raise RuntimeError("Call generate_config() and pass its build_dir before compiling")
-    """
+    root = Path(build_dir).resolve()
-    src = kernel_path(source)
+    src = (root / source).resolve()
-    if output is None:
+    dst = (root / (output or source.replace(".cu", ".ptx"))).resolve()
-        output = source.replace(".cu", ".ptx")
+    if root not in src.parents or root not in dst.parents:
-    dst = kernel_path(output)
+        raise ValueError("Compiler source and output must stay inside the build directory")
-    import shutil
+    if dst == src or dst.suffix != ".ptx":
+        raise ValueError("Compiler output must be a separate .ptx file")
+    dst.parent.mkdir(parents=True, exist_ok=True)
     nvcc = shutil.which("nvcc") or "/usr/local/cuda/bin/nvcc"
-    subprocess.run(
+    fd, tmp = tempfile.mkstemp(prefix=dst.stem + "-", suffix=".ptx", dir=dst.parent)
-        [nvcc, "-ptx", f"-arch={arch}", "--use_fast_math", src, "-o", dst],
+    os.close(fd)
-        check=True,
+    command = [nvcc, "-ptx", f"-arch={arch}", "--use_fast_math", str(src), "-o", tmp]
-    )
+    try:
-    return dst
+        subprocess.run(command, check=True)
+        os.replace(tmp, dst)
+        metadata = json.loads((root / ".celeris-build").read_text())
+        metadata.update(arch=arch, flags=["--use_fast_math"],
+                        ptx_sha256=hashlib.sha256(dst.read_bytes()).hexdigest())
+        _write(str(dst) + ".json", json.dumps(metadata, sort_keys=True))
+    finally:
+        if os.path.exists(tmp):
+            os.unlink(tmp)
+    _last_ptx_path.set(str(dst))
+    return str(dst)
 
 
-# ---------------------------------------------------------------------------
+def load_module(ptx_path: Optional[str] = None):
-# Module loading
+    """Load explicit PTX or the last PTX compiled in this Python context."""
-# ---------------------------------------------------------------------------
+    ptx_path = ptx_path or _last_ptx_path.get()
-def load_module(ptx_path: Optional[str] = None) -> cuda.Module:
+    if not ptx_path:
-    """Load a compiled PTX as a PyCUDA module."""
+        raise RuntimeError("No compiled PTX; call compile_kernel() first")
-    if ptx_path is None:
+    import pycuda.driver as cuda
-        ptx_path = kernel_path("kernel_v2.ptx")
     return cuda.module_from_file(ptx_path)
 
 

src/CelerisLab/cuda/context.py

@@ -1,41 +1,77 @@
-# CelerisLab/cuda/context.py
+"""CUDA context ownership with balanced, operation-scoped activation."""
-"""
+from __future__ import annotations
-CUDA device selection, context lifecycle, and device info queries.
-"""
 
-import subprocess
+from contextlib import contextmanager
-from typing import Optional
+from functools import wraps
+from threading import RLock
 
 import pycuda.driver as cuda
 
 
-class CudaContext:
+def context_bound(function):
-    """Manages a single CUDA device context.
+    """Activate an object's ctx for one GPU operation (nested calls are safe).
 
-    Usage::
+    Low-level fields constructed under an explicitly current context may omit
-
+    ctx. Simulation attaches its owner before exposing the field to callers.
-        ctx = CudaContext(device_id=0)
-        # ... use PyCUDA ...
-        ctx.close()
-
-    Or as a context-manager::
-
-        with CudaContext(0) as ctx:
-            ...
     """
+    @wraps(function)
+    def wrapped(self, *args, **kwargs):
+        ctx = getattr(self, "ctx", None)
+        if ctx is None:
+            return function(self, *args, **kwargs)
+        with ctx.activate():
+            return function(self, *args, **kwargs)
+    return wrapped
 
+
+class CudaContext:
+    """Own a context without leaving it on the caller's context stack.
+
+    Use ``with ctx.activate():`` for GPU work. Independent instances may be
+    closed in any order, provided no activation of the closing owner is active.
+    ``with CudaContext(device_id) as ctx:`` activates and closes one owner.
+    """
     def __init__(self, device_id: int = 0):
+        self._ctx = None
+        self._lock = RLock()
+        self._active = 0
+        self._scopes = []
         cuda.init()
         if device_id < 0 or device_id >= cuda.Device.count():
-            raise ValueError(
+            raise ValueError(f"device_id {device_id} invalid; {cuda.Device.count()} device(s) available.")
-                f"device_id {device_id} invalid; "
-                f"{cuda.Device.count()} device(s) available."
-            )
         self.device_id = device_id
         self.device = cuda.Device(device_id)
         self._ctx = self.device.make_context()
+        # make_context pushes immediately. Only balanced activate scopes may
+        # leave this owner current; close must never pop another instance.
+        self._ctx.pop()
+
+    @contextmanager
+    def activate(self):
+        with self._lock:
+            if self._ctx is None:
+                raise RuntimeError("CUDA context is closed")
+            context = self._ctx
+            context.push()
+            self._active += 1
+        try:
+            yield self
+        finally:
+            try:
+                context.pop()
+            finally:
+                with self._lock:
+                    self._active -= 1
+
+    @property
+    def is_active(self) -> bool:
+        with self._lock:
+            return self._active > 0
+
+    @property
+    def is_closed(self) -> bool:
+        return self._ctx is None
 
-    # -- properties ----------------------------------------------------------
     @property
     def name(self) -> str:
         return self.device.name()
@@ -47,7 +83,6 @@ class CudaContext:
 
     @property
     def sm_arch(self) -> str:
-        """Return 'sm_XX' string for nvcc -arch."""
         cc = self.device.compute_capability()
         return f"sm_{cc[0]}{cc[1]}"
 
@@ -55,25 +90,36 @@ class CudaContext:
     def total_memory(self) -> int:
         return self.device.total_memory()
 
-    # -- lifecycle -----------------------------------------------------------
     def close(self):
-        if self._ctx is not None:
+        with self._lock:
-            self._ctx.pop()
+            if self._ctx is None:
+                return
+            if self._active:
+                raise RuntimeError("Cannot close a CUDA context during an active operation")
+            self._ctx.detach()
             self._ctx = None
 
     def __enter__(self):
+        scope = self.activate()
+        scope.__enter__()
+        self._scopes.append(scope)
         return self
 
     def __exit__(self, *exc):
-        self.close()
+        self._scopes.pop().__exit__(*exc)
+        if not self._scopes:
+            self.close()
 
     def __del__(self):
-        self.close()
+        try:
+            self.close()
+        except Exception:
+            # Explicit close/context-manager exit is the supported lifecycle.
+            pass
 
 
 def detect_compute_capability(device_id: int = 0) -> str:
-    """Return 'sm_XX' for the given device without creating a persistent context."""
+    """Query architecture without creating a persistent current context."""
     cuda.init()
-    dev = cuda.Device(device_id)
+    cc = cuda.Device(device_id).compute_capability()
-    cc = dev.compute_capability()
     return f"sm_{cc[0]}{cc[1]}"

src/CelerisLab/lbm/__init__.py

@@ -1,16 +1,17 @@
-# CelerisLab/lbm/__init__.py
+"""Lattice Boltzmann operators and lazily loaded GPU resources."""
-"""
+from importlib import import_module
-Lattice Boltzmann Method module for fluid simulation.
+from .descriptors import *  # preserve existing flag/lattice exports
-"""
 
-try:
+__all__ = ["LBMField", "LBMStepper", "add_vortex"]
-    from .descriptors import *  # noqa: flags, lattice vectors
+_LAZY = {"LBMField": ("field", "LBMField"),
-    from .field import LBMField
+         "LBMStepper": ("stepper", "LBMStepper"),
-    from .stepper import LBMStepper
+         "add_vortex": ("initializers", "add_vortex")}
-    from .initializers import add_vortex
 
-    __all__ = ["LBMField", "LBMStepper", "add_vortex"]
+
-except ImportError as e:
+def __getattr__(name):
-    import warnings
+    if name not in _LAZY:
-    warnings.warn(f"LBM module not fully available: {e}", ImportWarning)
+        raise AttributeError(name)
-    __all__ = []
+    module, symbol = _LAZY[name]
+    value = getattr(import_module("." + module, __name__), symbol)
+    globals()[name] = value
+    return value

src/CelerisLab/lbm/curved_links.py

@@ -23,6 +23,25 @@ def _free_alloc(ptr: Optional[cuda.DeviceAllocation]) -> None:
         ptr.free()
 
 
+def _upload_columns(columns, stream=None):
+    """Copy low-frequency topology data; keep pinned staging alive to completion."""
+    if stream is None:
+        for device, host in columns:
+            cuda.memcpy_htod(device, host)
+        return
+    staging = []
+    try:
+        for device, host in columns:
+            pinned = cuda.pagelocked_empty(host.shape, host.dtype)
+            pinned[:] = host
+            staging.append(pinned)
+            cuda.memcpy_htod_async(device, pinned, stream)
+    finally:
+        # Returning pageable/temporary storage while DMA is active is invalid.
+        # Topology uploads are infrequent and intentionally complete on return.
+        stream.synchronize()
+
+
 @dataclass
 class CurvedLinkSoA:
     """Cut-link columns on host with matching ``*_gpu`` device buffers."""
@@ -56,6 +75,48 @@ class CurvedLinkSoA:
     count: int = 0
     capacity: int = 0
 
+    # Optional fixed surface velocity derived from true hit-center columns.
+    surface_velocity: Optional[np.ndarray] = None
+    aperture: Optional[np.ndarray] = None
+    surface_velocity_gpu: Optional[cuda.DeviceAllocation] = None
+    link_exchange_gpu: Optional[cuda.DeviceAllocation] = None
+
+    def assign_surface(self, objects, cfg, flags):
+        """Keep the original eight-column geometry/checkpoint contract unchanged."""
+        enabled = [obj for obj in objects if obj.surface_actuation is not None]
+        self.surface_velocity = self.aperture = None
+        if not enabled:
+            return
+        from .descriptors import FLUID, SOLID
+        if np.any((flags[self.fluid_idx] & FLUID) == 0) or np.any(flags[self.fluid_idx] & SOLID):
+            raise ValueError("Surface cut links must have active fluid receivers")
+        if len(np.unique(self.fluid_idx.astype(np.uint64)*cfg.nq+self.dir)) != len(self.dir):
+            raise ValueError("Surface boundary requires unique cut-link writers")
+        self.surface_velocity = np.zeros((len(self.dir), 2), np.float32)
+        self.aperture = np.zeros(len(self.dir), bool)
+        for obj in enabled:
+            obj.validate_surface_actuation(cfg)
+            selected = self.body_id == obj.obj_id
+            velocity, aperture = obj.surface_actuation.project(self.rx[selected], self.ry[selected])
+            if not aperture.any():
+                raise ValueError("Surface aperture hits no resolved cut links")
+            self.surface_velocity[selected] = velocity
+            self.aperture[selected] = aperture
+
+    def surface_exchange(self, body_id):
+        """Latest boundary step's (Fx,Fy,Tz), split by aperture; not nozzle thrust.
+
+        Caller synchronizes the CUDA stream/context. Rebuild/restart clears this
+        diagnostic; it is not the time-accumulated rigid-body observation.
+        """
+        if self.link_exchange_gpu is None:
+            raise RuntimeError("No surface actuation diagnostic is allocated")
+        values = np.empty((self.count, 3), np.float32)
+        cuda.memcpy_dtoh(values, self.link_exchange_gpu)
+        body = self.body_id == body_id
+        return {name: values[body & mask].sum(axis=0, dtype=np.float64)
+                for name, mask in (("wall", ~self.aperture), ("aperture", self.aperture))}
+
     def assign_host(
         self,
         fluid_idx: np.ndarray,
@@ -125,16 +186,21 @@ class CurvedLinkSoA:
                 raise ValueError(
                     "CurvedLinkSoA column 'fallback_class' must use classes {0, 1}"
                 )
+        self.surface_velocity = self.aperture = None
 
     def upload(self, stream: Optional[cuda.Stream] = None) -> None:
         """Upload host columns to GPU with capacity-based allocation reuse.
 
         Args:
-            stream: If set, use async host-to-device copies on this stream.
+            stream: If set, use pinned copies ordered on this stream and wait
+                for completion before releasing staging memory.
         """
         n = int(len(self.fluid_idx))
         self.count = n
         if n == 0:
+            _free_alloc(self.surface_velocity_gpu)
+            _free_alloc(self.link_exchange_gpu)
+            self.surface_velocity_gpu = self.link_exchange_gpu = None
             return
 
         if n > self.capacity or self.fluid_idx_gpu is None:
@@ -160,25 +226,24 @@ class CurvedLinkSoA:
         assert self.rz_gpu is not None
         assert self.fallback_class_gpu is not None
 
-        if stream is not None:
+        _upload_columns([
-            cuda.memcpy_htod_async(self.fluid_idx_gpu, self.fluid_idx, stream)
+            (self.fluid_idx_gpu, self.fluid_idx),
-            cuda.memcpy_htod_async(self.dir_gpu, self.dir, stream)
+            (self.dir_gpu, self.dir),
-            cuda.memcpy_htod_async(self.q_gpu, self.q, stream)
+            (self.q_gpu, self.q),
-            cuda.memcpy_htod_async(self.body_id_gpu, self.body_id, stream)
+            (self.body_id_gpu, self.body_id),
-            cuda.memcpy_htod_async(self.rx_gpu, self.rx, stream)
+            (self.rx_gpu, self.rx),
-            cuda.memcpy_htod_async(self.ry_gpu, self.ry, stream)
+            (self.ry_gpu, self.ry),
-            cuda.memcpy_htod_async(self.rz_gpu, self.rz, stream)
+            (self.rz_gpu, self.rz),
-            cuda.memcpy_htod_async(
+            (self.fallback_class_gpu, self.fallback_class),
-                self.fallback_class_gpu, self.fallback_class, stream)
+        ], stream)
-        else:
+        _free_alloc(self.surface_velocity_gpu)
-            cuda.memcpy_htod(self.fluid_idx_gpu, self.fluid_idx)
+        _free_alloc(self.link_exchange_gpu)
-            cuda.memcpy_htod(self.dir_gpu, self.dir)
+        self.surface_velocity_gpu = self.link_exchange_gpu = None
-            cuda.memcpy_htod(self.q_gpu, self.q)
+        if self.surface_velocity is not None:
-            cuda.memcpy_htod(self.body_id_gpu, self.body_id)
+            self.surface_velocity_gpu = cuda.mem_alloc(self.surface_velocity.nbytes)
-            cuda.memcpy_htod(self.rx_gpu, self.rx)
+            self.link_exchange_gpu = cuda.mem_alloc(n*3*4)
-            cuda.memcpy_htod(self.ry_gpu, self.ry)
+            _upload_columns([(self.surface_velocity_gpu, self.surface_velocity),
-            cuda.memcpy_htod(self.rz_gpu, self.rz)
+                             (self.link_exchange_gpu, np.zeros((n, 3), np.float32))], stream)
-            cuda.memcpy_htod(self.fallback_class_gpu, self.fallback_class)
 
     def free_gpu_columns(self) -> None:
         """Release device columns only (host arrays unchanged)."""
@@ -190,6 +255,9 @@ class CurvedLinkSoA:
         _free_alloc(self.ry_gpu)
         _free_alloc(self.rz_gpu)
         _free_alloc(self.fallback_class_gpu)
+        _free_alloc(self.surface_velocity_gpu)
+        _free_alloc(self.link_exchange_gpu)
+        self.surface_velocity_gpu = self.link_exchange_gpu = None
         self.fluid_idx_gpu = None
         self.dir_gpu = None
         self.q_gpu = None
@@ -229,7 +297,8 @@ class SensorSoA:
         """Reallocate GPU columns and upload.
 
         Args:
-            stream: If set, use async host-to-device copies on this stream.
+            stream: If set, use pinned copies ordered on this stream and wait
+                for completion before releasing staging memory.
         """
         self.free_gpu_columns()
         n = int(len(self.cells))
@@ -239,12 +308,10 @@ class SensorSoA:
 
         self.cells_gpu = cuda.mem_alloc(int(self.cells.nbytes))
         self.obj_id_gpu = cuda.mem_alloc(int(self.obj_id.nbytes))
-        if stream is not None:
+        _upload_columns([
-            cuda.memcpy_htod_async(self.cells_gpu, self.cells, stream)
+            (self.cells_gpu, self.cells),
-            cuda.memcpy_htod_async(self.obj_id_gpu, self.obj_id, stream)
+            (self.obj_id_gpu, self.obj_id),
-        else:
+        ], stream)
-            cuda.memcpy_htod(self.cells_gpu, self.cells)
-            cuda.memcpy_htod(self.obj_id_gpu, self.obj_id)
 
     def free_gpu_columns(self) -> None:
         _free_alloc(self.cells_gpu)
@@ -281,7 +348,8 @@ class ForceRegionSoA:
         """Reallocate GPU columns and upload.
 
         Args:
-            stream: If set, use async host-to-device copies on this stream.
+            stream: If set, use pinned copies ordered on this stream and wait
+                for completion before releasing staging memory.
         """
         self.free_gpu_columns()
         n = int(len(self.cells))
@@ -291,12 +359,10 @@ class ForceRegionSoA:
 
         self.cells_gpu = cuda.mem_alloc(int(self.cells.nbytes))
         self.obj_id_gpu = cuda.mem_alloc(int(self.obj_id.nbytes))
-        if stream is not None:
+        _upload_columns([
-            cuda.memcpy_htod_async(self.cells_gpu, self.cells, stream)
+            (self.cells_gpu, self.cells),
-            cuda.memcpy_htod_async(self.obj_id_gpu, self.obj_id, stream)
+            (self.obj_id_gpu, self.obj_id),
-        else:
+        ], stream)
-            cuda.memcpy_htod(self.cells_gpu, self.cells)
-            cuda.memcpy_htod(self.obj_id_gpu, self.obj_id)
 
     def free_gpu_columns(self) -> None:
         _free_alloc(self.cells_gpu)
@@ -341,7 +407,8 @@ class EsoPullCurvedBuffer:
         """Reallocate GPU columns and upload.
 
         Args:
-            stream: If set, use async host-to-device copies on this stream.
+            stream: If set, use pinned copies ordered on this stream and wait
+                for completion before releasing staging memory.
         """
         self.free_gpu_columns()
         n = int(len(self.corr_value))
@@ -352,14 +419,11 @@ class EsoPullCurvedBuffer:
         self.corr_value_gpu = cuda.mem_alloc(int(self.corr_value.nbytes))
         self.target_cell_gpu = cuda.mem_alloc(int(self.target_cell.nbytes))
         self.target_slot_gpu = cuda.mem_alloc(int(self.target_slot.nbytes))
-        if stream is not None:
+        _upload_columns([
-            cuda.memcpy_htod_async(self.corr_value_gpu, self.corr_value, stream)
+            (self.corr_value_gpu, self.corr_value),
-            cuda.memcpy_htod_async(self.target_cell_gpu, self.target_cell, stream)
+            (self.target_cell_gpu, self.target_cell),
-            cuda.memcpy_htod_async(self.target_slot_gpu, self.target_slot, stream)
+            (self.target_slot_gpu, self.target_slot),
-        else:
+        ], stream)
-            cuda.memcpy_htod(self.corr_value_gpu, self.corr_value)
-            cuda.memcpy_htod(self.target_cell_gpu, self.target_cell)
-            cuda.memcpy_htod(self.target_slot_gpu, self.target_slot)
 
     def free_gpu_columns(self) -> None:
         _free_alloc(self.corr_value_gpu)

src/CelerisLab/lbm/field.py

@@ -21,6 +21,7 @@ import numpy as np
 import pycuda.driver as cuda
 
 from ..config import LBMConfig
+from ..cuda.context import context_bound
 from .curved_links import CurvedLinkSoA, SensorSoA, ForceRegionSoA, EsoPullCurvedBuffer
 from .descriptors import FLUID, SOLID, BC_WALL, BC_INLET, BC_OUTLET
 
@@ -122,6 +123,8 @@ class LBMField:
         if self.cfg.is_d2q9:
             flag[0:self.nx] = SOLID | BC_WALL
             flag[(self.ny - 1) * self.nx:self.ny * self.nx] = SOLID | BC_WALL
+            if self.cfg.y_wall_bc == "on_site_velocity":
+                return flag  # x periodic; physical y planes are the adjacent fluid rows.
             flag[::self.nx] = SOLID | BC_INLET
             flag[self.nx - 1::self.nx] = SOLID | BC_OUTLET
             flag[0] = SOLID | BC_WALL
@@ -154,6 +157,7 @@ class LBMField:
         return self.force_regions.count
 
     # -- DDF download with precision decode ----------------------------------
+    @context_bound
     def download_ddf(self, *, step_id: int | None = None, force: bool = False):
         """Copy DDF from GPU to host self.ddf (float32), decoding FP16S if needed.
         If DDF shifting is enabled, un-shifts by adding back lattice weights."""
@@ -190,6 +194,7 @@ class LBMField:
         # See load_physical_node_esopull / store_physical_node_esopull.
         self._host_ddf_step = int(step_id) if step_id is not None else None
 
+    @context_bound
     def _read_lattice_weights(self):
         """Read d_w[] from CUDA __constant__ memory (single source of truth in descriptors.cuh)."""
         if self._cached_w is not None:
@@ -201,6 +206,7 @@ class LBMField:
         return self._cached_w
 
     # -- d_params upload -----------------------------------------------------
+    @context_bound
     def _upload_params(self):
         """Pack LBMParams struct and upload to __constant__ d_params.
 
@@ -221,11 +227,16 @@ class LBMField:
         ptr, size = self.module.get_global("d_params")
         cuda.memcpy_htod(ptr, data)
 
+    @context_bound
     def update_params(self, **kwargs):
         """Re-upload d_params after changing runtime-adjustable values.
 
         Accepted keys: omega, fx, fy, fz, n_objects.
         """
+        if self.cfg.y_wall_bc == "on_site_velocity" and (
+                any(float(kwargs.get(k, 0.0)) != 0 for k in ("fx", "fy", "fz"))
+                or int(kwargs.get("n_objects", 0)) != 0):
+            raise ValueError("on_site_velocity does not support body forces or objects")
         self._runtime_params.update({
             "omega": float(kwargs.get("omega", self._runtime_params["omega"])),
             "fx": float(kwargs.get("fx", self._runtime_params["fx"])),
@@ -234,46 +245,34 @@ class LBMField:
             "n_objects": int(kwargs.get("n_objects", self._runtime_params["n_objects"])),
         })
 
-        p = self._runtime_params
+        self._upload_params()
-        fmt = "=f fff I"
-        data = struct.pack(
-            fmt,
-            p["omega"],
-            p["fx"], p["fy"], p["fz"],
-            p["n_objects"],
-        )
-        ptr, _ = self.module.get_global("d_params")
-        cuda.memcpy_htod(ptr, data)
 
     # -- Host ↔ GPU transfers ------------------------------------------------
-    def upload_ddf(self):
+    def _encode_ddf(self, physical: np.ndarray) -> np.ndarray:
-        upload_data = self.ddf
+        """Encode a Q-major physical array (full field or selected cells)."""
+        encoded = np.asarray(physical, dtype=np.float32).copy()
         if getattr(self.cfg, 'ddf_shifting', False):
-            # Host stores physical DDF; GPU stores shifted DDF (f - w).
+            encoded.reshape(self.nq, -1)[:] -= self._read_lattice_weights()[:, None]
-            w = self._read_lattice_weights()
-            upload_data = upload_data.copy()
-            ddf_2d = upload_data.reshape(self.nq, -1)
-            for i in range(self.nq):
-                ddf_2d[i] -= w[i]
-
         if self.store_bytes == 2:
             if self.store_dtype != np.float16:
-                raise RuntimeError(
+                raise RuntimeError("Invalid 2-byte DDF storage dtype; expected float16 for FP16S.")
-                    "Invalid 2-byte DDF storage dtype; expected float16 for FP16S."
+            encoded *= np.float32(32768.0)
-                )
+        return np.ascontiguousarray(encoded, dtype=self.store_dtype)
-            buf = (upload_data * 32768.0).astype(self.store_dtype)
+
-            cuda.memcpy_htod(self.ddf_gpu, buf)
+    @context_bound
-            cuda.memcpy_htod(self.temp_gpu, buf)
+    def upload_ddf(self):
-        else:
+        upload_data = self._encode_ddf(self.ddf)
-            cuda.memcpy_htod(self.ddf_gpu, upload_data)
+        cuda.memcpy_htod(self.ddf_gpu, upload_data)
-            cuda.memcpy_htod(self.temp_gpu, upload_data)
+        cuda.memcpy_htod(self.temp_gpu, upload_data)
         self.invalidate_host_ddf_cache()
 
     # The canonical download_ddf is defined above in the class body.
 
+    @context_bound
     def upload_flags(self):
         cuda.memcpy_htod(self.flag_gpu, self.flag)
 
+    @context_bound
     def upload_compact_lists(self, stream: Optional[cuda.Stream] = None) -> None:
         """Upload compact lists (curved, sensor, force_region) to GPU.
 
@@ -294,6 +293,7 @@ class LBMField:
         self.esopull_curved.upload(stream=stream)
 
     # -- Read lattice descriptors from CUDA module ---------------------------
+    @context_bound
     def _read_lattice_vectors(self):
         """Read d_cx, d_cy [, d_cz] from CUDA __constant__ memory."""
         if self.cfg.is_d2q9 and self._cached_vectors_2d is not None:
@@ -318,17 +318,20 @@ class LBMField:
         return self._cached_vectors_2d
 
     # -- Macroscopic field extraction ----------------------------------------
+    @context_bound
     def get_macroscopic(self, *, step_id: int | None = None):
         """Download DDF and compute rho, ux, uy [, uz].
 
-        For double_buffer streaming: host-side decode from raw DDF.
+        Double-buffer returns moments of the current spatial postcollision
-        For esopull streaming: GPU-side MacroscopicEsoPullKernel.
+        field. EsoPull returns moments of the completed-count next-pull
+        (precollision) state via MacroscopicEsoPullKernel. These are different
+        numerical phases, so pointwise equality across streaming modes is not
+        a general contract.
 
         .. note::
-            In EsoPull mode, ``get_ddf()`` returns backing-layout raw storage,
+            In EsoPull mode, ``get_ddf()`` returns backing-layout values after
-            NOT physical DDF.  Only ``get_macroscopic()`` is semantically correct.
+            precision and shifting decode, not spatial physical DDF. Use this
-            Verified: Kan99b K2 (200k steps, MRT) -- bit-identical to
+            macroscopic path for the documented next-pull interpretation.
-            double-buffer for all metrics (St, mean_CL, mean_CD, amp_CL, amp_CD).
         """
         if getattr(self.cfg, 'streaming', 'double_buffer') == 'esopull':
             return self.get_macroscopic_esopull(step_id=step_id)
@@ -358,6 +361,7 @@ class LBMField:
 
         raise ValueError(f"Unsupported lattice_model: {self.cfg.lattice_model}")
 
+    @context_bound
     def get_macroscopic_esopull(self, *, step_id: int | None = None) -> dict:
         """Compute macroscopic field from EsoPull backing layout on GPU.
 
@@ -366,11 +370,10 @@ class LBMField:
 
         .. note::
             Only supports D2Q9.  D3Q19 raises NotImplementedError.
-            This is the ONLY correct path for macroscopic diagnostics in EsoPull
+            At the completed step count this reads the next-pull/precollision
-            mode.  Host-side raw DDF download (get_ddf) is NOT physical DDF
+            state. Host get_ddf decodes precision and shifting but retains the
-            when streaming is esopull -- the backing layout scatters paired
+            scattered backing layout; it is not spatial physical DDF. This
-            directions across neighbor cells.
+            phase differs from the double-buffer postcollision diagnostic.
-            Verified: Kan99b K2 -- bit-identical to double-buffer.
 
         Args:
             step_id: Current step count used for EsoPull parity decoding.
@@ -408,6 +411,7 @@ class LBMField:
     # LBMField.
 
     # -- EsoPull sync decode/encode (compact list) --------------------------
+    @context_bound
     def esopull_sync_decode(self, patch_indices, step_count, stream=None):
         """Decode a subset of cells from EsoPull backing to physical DDF.
 
@@ -420,9 +424,11 @@ class LBMField:
             (scratch_gpu, scratch_host, idx_gpu): scratch buffer handles.
         """
         n = len(patch_indices)
+        if n == 0:
+            return None, np.empty(0, dtype=np.float32), None
         nq = self.nq
         scratch_gpu = cuda.mem_alloc(n * nq * 4)
-        scratch_host = np.empty(n * nq, dtype=np.float32)
+        scratch_host = cuda.pagelocked_empty(n * nq, np.float32)
         idx_gpu = cuda.mem_alloc(n * 4)
         cuda.memcpy_htod(idx_gpu, np.asarray(patch_indices, dtype=np.uint32))
 
@@ -433,20 +439,27 @@ class LBMField:
            block=(tpb, 1, 1), grid=((n + tpb - 1) // tpb, 1, 1),
            stream=stream)
 
-        cuda.memcpy_dtoh_async(scratch_host, scratch_gpu, stream)
+        if stream is None:
-        stream.synchronize()
+            cuda.memcpy_dtoh(scratch_host, scratch_gpu)
+        else:
+            cuda.memcpy_dtoh_async(scratch_host, scratch_gpu, stream)
+            stream.synchronize()
         return scratch_gpu, scratch_host, idx_gpu
 
+    @context_bound
     def esopull_sync_encode(self, scratch_gpu, scratch_host, idx_gpu,
                             n, step_count, stream=None):
-        """Upload patched physical DDF, collide + encode to EsoPull layout.
+        """Restore the same-count next-pull physical state, without collision.
 
-        The backing layout stores post-collision values.  This method
+        Scratch is cell-major physical FP32. The kernel applies shifting and
-        applies collision to the patched physical DDF before encoding,
+        storage encoding and is the exact inverse of same-count decode.
-        matching the store_f_esopull convention at parity *step_count*.
-        Frees scratch_gpu and idx_gpu after encoding.
         """
-        cuda.memcpy_htod_async(scratch_gpu, scratch_host, stream)
+        if n == 0:
+            return
+        if stream is None:
+            cuda.memcpy_htod(scratch_gpu, scratch_host)
+        else:
+            cuda.memcpy_htod_async(scratch_gpu, scratch_host, stream)
 
         fn = self.module.get_function("EncodePhysicalToCells")
         tpb = self.cfg.threads_per_block
@@ -455,15 +468,45 @@ class LBMField:
            block=(tpb, 1, 1), grid=((n + tpb - 1) // tpb, 1, 1),
            stream=stream)
 
+        # Keep host/device scratch alive until the async upload and kernel end.
+        if stream is None:
+            cuda.Context.synchronize()
+        else:
+            stream.synchronize()
         scratch_gpu.free()
         idx_gpu.free()
+        self.invalidate_host_ddf_cache()
 
-    def snapshot(self):
+    @context_bound
-        self.download_ddf(force=True)
+    def snapshot(self, *, step_id: int | None = None):
-        self._ddf_snap = self.ddf.copy()
+        """Save a field-only snapshot; simulation time and controls are unchanged."""
+        if self.cfg.streaming == "esopull":
+            if step_id is None:
+                raise ValueError("EsoPull field snapshots require the current step_id")
+            indices = np.arange(self.n, dtype=np.uint32)
+            scratch, physical, idx_gpu = self.esopull_sync_decode(indices, step_id)
+            self._ddf_snap = physical.reshape(self.n, self.nq).T.copy().reshape(-1)
+            scratch.free()
+            idx_gpu.free()
+        else:
+            self.download_ddf(force=True)
+            self._ddf_snap = self.ddf.copy()
 
-    def restore(self):
+    @context_bound
+    def restore(self, *, step_id: int | None = None):
+        """Restore the saved field at the current clock, including EsoPull parity."""
         if self._ddf_snap is None:
             raise RuntimeError("No snapshot to restore")
+        if self.cfg.streaming == "esopull" and step_id is None:
+            raise ValueError("EsoPull field snapshots require the current step_id")
         self.ddf = self._ddf_snap.copy()
-        self.upload_ddf()
+        if self.cfg.streaming == "esopull":
+            indices = np.arange(self.n, dtype=np.uint32)
+            scratch_gpu = cuda.mem_alloc(self.n * self.nq * 4)
+            physical = cuda.pagelocked_empty(self.n * self.nq, np.float32)
+            physical[:] = self.ddf.reshape(self.nq, self.n).T.reshape(-1)
+            idx_gpu = cuda.mem_alloc(indices.nbytes)
+            cuda.memcpy_htod(idx_gpu, indices)
+            self.esopull_sync_encode(scratch_gpu, physical, idx_gpu, self.n, step_id)
+        else:
+            self.upload_ddf()

src/CelerisLab/lbm/initializers.py

@@ -10,7 +10,7 @@ import numpy as np
 from scipy.special import jv, expi
 from typing import Tuple
 
-import pycuda.driver as cuda
+from .equilibrium import compute_feq_d2q9
 
 
 def add_vortex(field, center: Tuple[float, float],
@@ -22,6 +22,12 @@ def add_vortex(field, center: Tuple[float, float],
     """
     if not field.cfg.is_d2q9:
         raise NotImplementedError("Vortex init only for D2Q9")
+    if field.cfg.streaming != "double_buffer":
+        raise NotImplementedError(
+            "Vortex initialization requires double_buffer physical DDF; "
+            "raw EsoPull storage cannot be edited without its time/parity.")
+    if not np.isfinite(radius) or radius <= 0:
+        raise ValueError("Vortex radius must be finite and positive")
     if vortex_type not in ("lamb", "oseen", "taylor"):
         raise ValueError(f"Unknown vortex type: {vortex_type}")
 
@@ -49,7 +55,6 @@ def add_vortex(field, center: Tuple[float, float],
 
     # Read lattice descriptors from CUDA module (single source of truth)
     cx, cy = field._read_lattice_vectors()
-    w = field._read_lattice_weights()
 
     # Compute current macroscopic from DDF
     rho_old = np.sum(f, axis=0)
@@ -68,10 +73,7 @@ def add_vortex(field, center: Tuple[float, float],
     for j in range(1, ny - 1):
         for i in range(1, nx - 1):
             u, v, p = float(u_new[j, i]), float(v_new[j, i]), float(p_new[j, i])
-            for e in range(nq):
+            f[:, j, i] = compute_feq_d2q9(3*p, u, v)
-                eu = cx[e] * u + cy[e] * v
-                u2 = u * u + v * v
-                f[e, j, i] = w[e] * (3*p + 3*eu + 4.5*eu*eu - 1.5*u2)
 
     field.ddf = f.reshape(-1)
     field.upload_ddf()
@@ -103,16 +105,23 @@ def _lamb_vortex(r, theta, radius, strength):
 
 def _oseen_vortex(r, theta, radius, strength):
     kappa = 2 * np.pi * radius**2 * strength
-    r_safe = np.where(r > 1e-12, r, 1e-12)
-    exp_term = 1 - np.exp(-4 * r**2 / radius**2)
-    u_vor = -kappa / (2*np.pi*r_safe) * exp_term * np.sin(theta)
-    v_vor = kappa / (2*np.pi*r_safe) * exp_term * np.cos(theta)
     zeta = 4 * r**2 / radius**2
-    p_vor = (
+    exp_term = -np.expm1(-zeta)
-        -kappa**2 / (8*np.pi**2 * r_safe**2)
+    speed = kappa / (2*np.pi) * np.divide(
-        * (-2*zeta*(expi(-zeta) - expi(-2*zeta))
+        exp_term, r, out=np.zeros_like(r, dtype=float), where=r != 0)
-           + (1 - np.exp(-zeta))**2)
+    u_vor = -speed * np.sin(theta)
-    )
+    v_vor = speed * np.cos(theta)
+    # Factor out zeta/r**2 before evaluating the pressure. At the center,
+    # Ei(-zeta)-Ei(-2*zeta) tends to -log(2), not NaN. The expansion also
+    # avoids cancellation for near-center samples and preserves the same model.
+    pressure_shape = np.empty_like(zeta, dtype=float)
+    small = zeta < 1e-5
+    z = zeta[small]
+    pressure_shape[small] = 2*np.log(2.0) - z + z*z/2 - 7*z*z*z/36
+    z = zeta[~small]
+    pressure_shape[~small] = (
+        -2*(expi(-z) - expi(-2*z)) + np.expm1(-z)**2 / z)
+    p_vor = -kappa**2 / (2*np.pi**2*radius**2) * pressure_shape
     return u_vor, v_vor, p_vor
 
 
@@ -127,4 +136,4 @@ def _taylor_vortex(r, theta, radius, strength, nu):
         * np.exp(-2 * r**2 / radius**2)
         / (np.pi**2 * radius**6)
     )
-    return u_vor, v_vor, p_vor
+    return u_vor, v_vor, p_vor

src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh

@@ -5,10 +5,7 @@
 // We replace those directions with the opposite-direction population at the
 // SAME node from the INPUT buffer (half-way BB).
 //
-// In EsoPull mode (t > 0), reads via load_physical_dir_esopull to resolve
+// EsoPull stages its pre-streaming donors in PrepareEsoPullBoundaries.
-// the opposite direction's value from the EsoPull single-buffer layout.
-// For Kan99b (cylinder far from walls), this is correct in practice.
-// A more rigorous pre_streaming fix is pending further investigation.
 // ============================================================================
 
 #ifndef CELERIS_BOUNDARY_BOUNCE_BACK_CUH
@@ -21,33 +18,19 @@
 __device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y,
                                                   float* __restrict__ f,
                                                   const fpxx* __restrict__ fi_in,
-                                                  unsigned long k,
+                                                  unsigned long k)
-                                                  const unsigned long* j = nullptr,
-                                                  unsigned long t = 0ul)
 {
     if (y == 1) {
         // Directions sourced from y=0 wall in pull step: +y, +x+y, -x+y.
-        if (t == 0ul) {
+        f[3] = load_ddf(fi_in, index_f(k, 4u));
-            f[3] = load_ddf(fi_in, index_f(k, 4u));
+        f[5] = load_ddf(fi_in, index_f(k, 6u));
-            f[5] = load_ddf(fi_in, index_f(k, 6u));
+        f[8] = load_ddf(fi_in, index_f(k, 7u));
-            f[8] = load_ddf(fi_in, index_f(k, 7u));
-        } else {
-            f[3] = load_physical_dir_esopull(k, 4u, fi_in, j, t);  // N ← S
-            f[5] = load_physical_dir_esopull(k, 6u, fi_in, j, t);  // NE ← SW
-            f[8] = load_physical_dir_esopull(k, 7u, fi_in, j, t);  // NW ← SE
-        }
     }
     else if (y == (unsigned int)(NY - 2)) {
         // Directions sourced from y=NY-1 wall in pull step: -y, -x-y, +x-y.
-        if (t == 0ul) {
+        f[4] = load_ddf(fi_in, index_f(k, 3u));
-            f[4] = load_ddf(fi_in, index_f(k, 3u));
+        f[6] = load_ddf(fi_in, index_f(k, 5u));
-            f[6] = load_ddf(fi_in, index_f(k, 5u));
+        f[7] = load_ddf(fi_in, index_f(k, 8u));
-            f[7] = load_ddf(fi_in, index_f(k, 8u));
-        } else {
-            f[4] = load_physical_dir_esopull(k, 3u, fi_in, j, t);  // S ← N
-            f[6] = load_physical_dir_esopull(k, 5u, fi_in, j, t);  // SW ← NE
-            f[7] = load_physical_dir_esopull(k, 8u, fi_in, j, t);  // SE ← NW
-        }
     }
 }
 
@@ -62,31 +45,17 @@ __device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y,
 __device__ inline void apply_wall_freeslip_d2q9_y_pull(unsigned int y,
                                                         float* __restrict__ f,
                                                         const fpxx* __restrict__ fi_in,
-                                                        unsigned long k,
+                                                        unsigned long k)
-                                                        const unsigned long* j = nullptr,
-                                                        unsigned long t = 0ul)
 {
     if (y == 1) {
-        if (t == 0ul) {
+        f[3] = load_ddf(fi_in, index_f(k, 4u));
-            f[3] = load_ddf(fi_in, index_f(k, 4u));
+        f[5] = load_ddf(fi_in, index_f(k, 7u));
-            f[5] = load_ddf(fi_in, index_f(k, 7u));
+        f[8] = load_ddf(fi_in, index_f(k, 6u));
-            f[8] = load_ddf(fi_in, index_f(k, 6u));
-        } else {
-            f[3] = load_physical_dir_esopull(k, 4u, fi_in, j, t);  // N ← S  (normal)
-            f[5] = load_physical_dir_esopull(k, 7u, fi_in, j, t);  // NE ← SE (tangential mirror)
-            f[8] = load_physical_dir_esopull(k, 6u, fi_in, j, t);  // NW ← SW (tangential mirror)
-        }
     }
     else if (y == (unsigned int)(NY - 2)) {
-        if (t == 0ul) {
+        f[4] = load_ddf(fi_in, index_f(k, 3u));
-            f[4] = load_ddf(fi_in, index_f(k, 3u));
+        f[6] = load_ddf(fi_in, index_f(k, 8u));
-            f[6] = load_ddf(fi_in, index_f(k, 8u));
+        f[7] = load_ddf(fi_in, index_f(k, 5u));
-            f[7] = load_ddf(fi_in, index_f(k, 5u));
-        } else {
-            f[4] = load_physical_dir_esopull(k, 3u, fi_in, j, t);  // S ← N  (normal)
-            f[6] = load_physical_dir_esopull(k, 8u, fi_in, j, t);  // SW ← NW (tangential mirror)
-            f[7] = load_physical_dir_esopull(k, 5u, fi_in, j, t);  // SE ← NE (tangential mirror)
-        }
     }
 }
 

src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh

@@ -48,26 +48,27 @@ __device__ __forceinline__ float bouzidi_linear_moving_correction(
 {
     if (fallback_class != CURVED_FALLBACK_BOUZIDI) {
         // Fallback replaces the q < 1/2 branch when its donor is illegal.
-        // For D2Q9, alpha_i = 3 w_i, so the moving-wall addend is
+        // Here dir points from fluid to solid and alpha_i = 3 w_i rho.
-        // +2 * alpha_i * (c_i . u_w) in the q < 1/2 branch.
+        // The opposite-direction reflected population receives
-        return 2.0f * alpha_ci_dot_uw;
+        // -2 * alpha_i * (c_i . u_w) in the q < 1/2 branch.
+        return -2.0f * alpha_ci_dot_uw;
     }
     if (q < 0.5f) {
-        return 2.0f * alpha_ci_dot_uw;
+        return -2.0f * alpha_ci_dot_uw;
     }
-    return (alpha_ci_dot_uw / q);
+    return -(alpha_ci_dot_uw / q);
 }
 
 __device__ __forceinline__ float compute_moving_wall_addend(
     unsigned int dir, float q, unsigned char fallback_class,
-    float Uw, float Vw)
+    float Uw, float Vw, float rho = 1.0f)
 {
 #if DIM == 2
     const float ci_dot_uw = (float)d_cx[dir] * Uw + (float)d_cy[dir] * Vw;
 #elif DIM == 3
     const float ci_dot_uw = (float)d_cx[dir] * Uw + (float)d_cy[dir] * Vw + (float)d_cz[dir] * 0.0f;
 #endif
-    const float alpha_ci_dot_uw = 3.0f * d_w[dir] * ci_dot_uw;
+    const float alpha_ci_dot_uw = 3.0f * d_w[dir] * rho * ci_dot_uw;
     return bouzidi_linear_moving_correction(q, fallback_class, alpha_ci_dot_uw);
 }
 
@@ -81,7 +82,7 @@ __device__ inline void apply_bouzidi_link(
     float Uw, float Vw,
     float rx, float ry,
     unsigned char fallback_class,
-    float* __restrict__ obs, int body_id)
+    float* __restrict__ obs, int body_id, float* link_exchange = nullptr)
 {
     const unsigned int dir_opp = (unsigned int)opp_dir((int)dir);
 
@@ -108,13 +109,18 @@ __device__ inline void apply_bouzidi_link(
     } else {
         f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback_class);
     }
-    f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw);
+    float rho_wall = USE_DDF_SHIFTING ? 1.0f : 0.0f;
+    if (Uw != 0.0f || Vw != 0.0f) {
+        for (unsigned int i = 0; i < NQ; ++i)
+            rho_wall += load_ddf(fi, index_f(k_f, i));
+    }
+    f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw, rho_wall);
 
     // Write to the solid source node so the subsequent pull step loads the
     // corrected incoming population into the adjacent fluid node.
     store_ddf(fi, index_f(k_s, dir_opp), f_reflected);
 
-    if (obs != nullptr) {
+    if (obs != nullptr || link_exchange != nullptr) {
 #if USE_DDF_SHIFTING
         const float w_i = (float)d_w[dir];
         const float fx = (float)d_cx[dir] * ((f_toward + w_i) + (f_reflected + w_i));
@@ -123,10 +129,15 @@ __device__ inline void apply_bouzidi_link(
         const float fx = (float)d_cx[dir] * (f_toward + f_reflected);
         const float fy = (float)d_cy[dir] * (f_toward + f_reflected);
 #endif
-        atomicAdd(&obs[obs_force_index(body_id, 0)], fx);
-        atomicAdd(&obs[obs_force_index(body_id, 1)], fy);
         const float tz = rx * fy - ry * fx;
-        atomicAdd(&obs[obs_torque_index(body_id, 0)], tz);
+        if (link_exchange != nullptr) {
+            link_exchange[0] = fx; link_exchange[1] = fy; link_exchange[2] = tz;
+        }
+        if (obs != nullptr) {
+            atomicAdd(&obs[obs_force_index(body_id, 0)], fx);
+            atomicAdd(&obs[obs_force_index(body_id, 1)], fy);
+            atomicAdd(&obs[obs_torque_index(body_id, 0)], tz);
+        }
     }
 }
 #endif // DIM == 2
@@ -170,7 +181,12 @@ __device__ inline void apply_bouzidi_link(
     } else {
         f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback_class);
     }
-    f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw);
+    float rho_wall = USE_DDF_SHIFTING ? 1.0f : 0.0f;
+    if (Uw != 0.0f || Vw != 0.0f) {
+        for (unsigned int i = 0; i < NQ; ++i)
+            rho_wall += load_ddf(fi, index_f(k_f, i));
+    }
+    f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw, rho_wall);
 
     store_ddf(fi, index_f(k_s, dir_opp), f_reflected);
 

src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh

@@ -46,11 +46,8 @@ __device__ __forceinline__ float west_velocity_rho_closure_d2q9(
 {
     float rho_num = f[0] + f[3] + f[4] + 2.0f * (f[2] + f[6] + f[8]);
 #if USE_DDF_SHIFTING
-    // f[i] = physical_f[i] - w[i]; restore the sum of weights for
+    // Sum of the descriptor weights in this closure (including doubles).
-    // the directions appearing in the numerator (w[0]=4/9, w[3]=w[4]=1/9,
+    rho_num += 1.0f;
-    // w[2]=w[6]=w[8]=1/36, with f[2],f[6],f[8] double-counted):
-    // sum(w) = 4/9 + 1/9 + 1/9 + 2*(1/36+1/36+1/36) = 6/9 + 6/36 = 5/6
-    rho_num += 5.0f / 6.0f;
 #endif
     return rho_num / (1.0f - ux_target);
 }
@@ -64,11 +61,8 @@ __device__ __forceinline__ float west_velocity_rho_closure_d3q19(
     float rho_num = f[0] + f[3] + f[4] + f[5] + f[6] + f[11] + f[12] + f[17] + f[18]
                   + 2.0f * (f[2] + f[8] + f[10] + f[14] + f[16]);
 #if USE_DDF_SHIFTING
-    // f[i] = physical_f[i] - w[i]; restore the sum of weights for
+    // Sum of the descriptor weights in this closure (including doubles).
-    // the directions appearing in the numerator (w[0]=1/3, w[3..6,11,12,17,18]=1/18,
+    rho_num += 1.0f;
-    // w[2,8,10,14,16]=1/36 double-counted):
-    // sum(w) = 1/3 + 8*(1/18) + 5*(2/36) = 12/36 + 16/36 + 10/36 = 38/36 = 19/18
-    rho_num += 19.0f / 18.0f;
 #endif
     return rho_num / (1.0f - ux_target);
 }

src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh

@@ -21,34 +21,32 @@
 // west_velocity_rho_closure_d2q9() need clean known-direction values.  The
 // free-slip wall interferes with these at the top/bottom inlet corners.
 // The legacy name is kept for now (not renamed yet).
-__device__ inline void repair_zou_he_west_knowns_d2q9(
+__device__ inline void repair_west_knowns_from_donor_d2q9(
-    float* __restrict__ f,
+    float* __restrict__ f, const float* __restrict__ donor,
-    const fpxx* __restrict__ fi_in,
+    unsigned int x, unsigned int y)
-    unsigned int x,
-    unsigned int y,
-    const unsigned long* j_neb = nullptr,
-    unsigned long t = 0ul)
 {
     if (x != 0u) return;
-
-    const unsigned long k_int = linear_index(x + 1u, y);
-
     if (y == 1u) {
-        if (t == 0ul || j_neb == nullptr) {
+        f[3] = donor[3]; f[8] = donor[8];
-            f[3] = load_ddf(fi_in, index_f(k_int, 3u));
-            f[8] = load_ddf(fi_in, index_f(k_int, 8u));
-        } else {
-            f[3] = load_physical_dir_pre_streaming(k_int, 3u, fi_in, j_neb, t);
-            f[8] = load_physical_dir_pre_streaming(k_int, 8u, fi_in, j_neb, t);
-        }
     } else if (y == (unsigned int)(NY - 2)) {
-        if (t == 0ul || j_neb == nullptr) {
+        f[4] = donor[4]; f[6] = donor[6];
-            f[4] = load_ddf(fi_in, index_f(k_int, 4u));
+    }
-            f[6] = load_ddf(fi_in, index_f(k_int, 6u));
+}
-        } else {
+
-            f[4] = load_physical_dir_pre_streaming(k_int, 4u, fi_in, j_neb, t);
+// Double-buffer donors are immutable during the pull kernel. EsoPull calls
-            f[6] = load_physical_dir_pre_streaming(k_int, 6u, fi_in, j_neb, t);
+// the helper above with the same spatial state staged before in-place writes.
-        }
+__device__ inline void repair_zou_he_west_knowns_d2q9(
+    float* __restrict__ f, const fpxx* __restrict__ fi_in,
+    unsigned int x, unsigned int y)
+{
+    if (x != 0u || (y != 1u && y != (unsigned int)(NY - 2))) return;
+    const unsigned long k_int = linear_index(x + 1u, y);
+    if (y == 1u) {
+        f[3] = load_ddf(fi_in, index_f(k_int, 3u));
+        f[8] = load_ddf(fi_in, index_f(k_int, 8u));
+    } else {
+        f[4] = load_ddf(fi_in, index_f(k_int, 4u));
+        f[6] = load_ddf(fi_in, index_f(k_int, 6u));
     }
 }
 

src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh

@@ -131,29 +131,19 @@ __device__ __forceinline__ void apply_inlet_esopull_d2q9(
     float* __restrict__ f,
     unsigned int x,
     unsigned int y,
-    const fpxx* __restrict__ fi,
+    const float* __restrict__ f_neb)
-    unsigned long t)
 {
     const float u_target = inlet_target_u((float)y);
     const float v_target = 0.0f;
 
     // Free-slip y-walls: repair west inlet ghost-node known directions.
 #if Y_WALL_BC == 1
-    {
+    repair_west_knowns_from_donor_d2q9(f, f_neb, x, y);
-        unsigned long j_neb[NQ];
-        compute_neighbors(x + 1u, y, j_neb);
-        repair_zou_he_west_knowns_d2q9(f, fi, x, y, j_neb, t);
-    }
 #endif
 
 #if INLET_SCHEME == 0
     apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target);
 #elif INLET_SCHEME == 1 || INLET_SCHEME == 3
-    const unsigned long k_neb = linear_index(x + 1u, y);
-    unsigned long j_neb[NQ];
-    compute_neighbors(x + 1u, y, j_neb);
-    float f_neb[NQ];
-    load_f_esopull(k_neb, f_neb, fi, j_neb, t);
   #if INLET_SCHEME == 1
     apply_channel_stabilized_inlet_d2q9(f, f_neb, (float)y);
   #else
@@ -167,17 +157,8 @@ __device__ __forceinline__ void apply_inlet_esopull_d2q9(
 }
 
 __device__ __forceinline__ void apply_outlet_esopull_d2q9(
-    float* __restrict__ f,
+    float* __restrict__ f, const float* __restrict__ f_neb)
-    unsigned int x,
-    unsigned int y,
-    const fpxx* __restrict__ fi,
-    unsigned long t)
 {
-    const unsigned long k_neb = linear_index(x - 1u, y);
-    unsigned long j_neb[NQ];
-    compute_neighbors(x - 1u, y, j_neb);
-    float f_neb[NQ];
-    load_f_esopull(k_neb, f_neb, fi, j_neb, t);
     apply_pressure_outlet_d2q9(f, f_neb);
 }
 #endif

src/CelerisLab/lbm/kernels/boundary/planar_velocity.cuh

@@ -0,0 +1,44 @@
+// On-site D2Q9 Zou-He velocity closure, paired lattice ordering.
+// Input is poststream f at the arrival time. Only unknown populations change;
+// the caller owns collision/store. No force, shifting or off-lattice walls.
+#ifndef CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH
+#define CELERIS_BOUNDARY_PLANAR_VELOCITY_CUH
+
+__device__ __forceinline__ float bottom_wall_normal_velocity(unsigned int x)
+{
+#if WALL_SLOT_WIDTH > 0 && WALL_SLOT_WIDTH < NX
+    // C1 compact profile centered at (NX-1)/2; amplitude is its continuous peak.
+    const float s = (2.0f*(float)x-(float)(NX-1))/(float)WALL_SLOT_WIDTH;
+    const float w = fmaxf(0.0f, 1.0f-s*s);
+    return WALL_NORMAL_VELOCITY*w*w;
+#else
+    return WALL_NORMAL_VELOCITY; // Disabled or explicitly uniform full width.
+#endif
+}
+
+__device__ __forceinline__ void reconstruct_bottom_velocity(float* f, float u, float v)
+{
+    const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[4] + f[6] + f[7])) / (1.0f-v);
+    f[3] = f[4] + (2.0f/3.0f)*rho*v;
+    f[5] = f[6] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u + (1.0f/6.0f)*rho*v;
+    f[8] = f[7] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u + (1.0f/6.0f)*rho*v;
+}
+
+__device__ __forceinline__ void reconstruct_top_velocity(float* f, float u, float v)
+{
+    const float rho = (f[0] + f[1] + f[2] + 2.0f*(f[3] + f[5] + f[8])) / (1.0f+v);
+    f[4] = f[3] - (2.0f/3.0f)*rho*v;
+    f[6] = f[5] + 0.5f*(f[1]-f[2]) - 0.5f*rho*u - (1.0f/6.0f)*rho*v;
+    f[7] = f[8] + 0.5f*(f[2]-f[1]) + 0.5f*rho*u - (1.0f/6.0f)*rho*v;
+}
+
+// Fixed top density and tangential u=0; positive Cartesian v exits upward.
+// J_y=rho*v follows from known poststream populations (ordinary D2Q9).
+__device__ __forceinline__ void reconstruct_top_reservoir(float* f, float rho)
+{
+    const float jy = f[0]+f[1]+f[2]+2.0f*(f[3]+f[5]+f[8])-rho;
+    f[4] = f[3] - (2.0f/3.0f)*jy;
+    f[6] = f[5] + 0.5f*(f[1]-f[2]) - (1.0f/6.0f)*jy;
+    f[7] = f[8] + 0.5f*(f[2]-f[1]) - (1.0f/6.0f)*jy;
+}
+#endif

src/CelerisLab/lbm/kernels/kernel_v2.cu

@@ -56,6 +56,9 @@
 #include "boundary/bounce_back.cuh"
 #include "boundary/curved_boundary.cuh"
 #include "boundary/inlet_outlet.cuh"
+#if Y_WALL_BC == 2
+#include "boundary/planar_velocity.cuh"
+#endif
 #include "boundary/ibm_kernels.cuh"
 
 // ============================================================================

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

@@ -16,21 +16,13 @@
 //   m[7] = pxx   (stress,     s₇ = s_nu = ω = 1/(3ν + 0.5))
 //   m[8] = pxy   (stress,     s₈ = s_nu)
 //
-// NOTE on ddf_shifting (USE_DDF_SHIFTING):
+// D2Q9 transforms the stored populations directly. For shifted storage,
-//   MRT is NOT shift-invariant in moment space because the moment transform
+// M*w = (1,-2,1,0,0,0,0,0,0) is removed from the equilibrium moments.
-//   of the weight vector w produces non-zero entries in m[1] (energy) and
+// This avoids rounding small populations through an O(1) weight background.
-//   m[2] (energy^2).  Unlike SRT/TRT where (f-feq) cancels weights pairwise,
+// D3Q19 projects g-feq; both arrays use the same storage representation.
-//   the MRT moment construction and inverse transform do not preserve the
-//   shifted-space semantics.
 //
-//   Therefore when USE_DDF_SHIFTING=1, this function temporarily unshifts
+// Fin is the raw Guo source. Applying half before and half after linear
-//   the distributions to physical space before collision, then shifts them
+// relaxation gives the model-specific M^-1 (I-S/2) M source correction.
-//   back after.  This ensures the moment construction, relaxation, and
-//   inverse transform all operate on physically meaningful distributions.
-//
-// Forcing contract:
-//   Fin[i] is the raw Guo source term without the (1 - ω/2) prefactor.
-//   This operator applies the prefactor internally, matching collide_srt().
 // ============================================================================
 
 #ifndef CELERIS_OPERATORS_COLLISION_MRT_CUH
@@ -43,14 +35,6 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
                                             const float* __restrict__ Fin,
                                             float omega)
 {
-    // MRT is not shift-invariant in moment space.  When USE_DDF_SHIFTING=1,
-    // temporarily unshift to physical space before collision.
-#if USE_DDF_SHIFTING
-    g[0] += (float)d_w[0]; g[1] += (float)d_w[1]; g[2] += (float)d_w[2];
-    g[3] += (float)d_w[3]; g[4] += (float)d_w[4]; g[5] += (float)d_w[5];
-    g[6] += (float)d_w[6]; g[7] += (float)d_w[7]; g[8] += (float)d_w[8];
-#endif
-
     const float s_rho = 0.0f;
     const float s_e   = 1.2f;
     const float s_eps = 1.2f;
@@ -58,9 +42,10 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
     const float s_q   = 1.2f;
     const float s_jy  = 0.0f;
     const float s_nu  = omega;
-    const float c_tau = 1.0f - 0.5f * omega;
 
-    const float p = rho * (1.0f / 3.0f);
+    // Trapezoidal source: collide (g + Fin/2), then add Fin/2.
+    // Equivalent to M^-1 (I-S/2) M Fin for this relaxation matrix.
+    for (int i = 0; i < 9; ++i) g[i] += 0.5f * Fin[i];
 
     float m[9];
     m[0] =    g[0] +   g[1] +   g[2] +   g[3] +   g[4] +   g[5] +   g[6] +   g[7] +   g[8];
@@ -75,9 +60,17 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
 
     const float u2 = ux * ux + uy * uy;
     float meq[9];
+#if USE_DDF_SHIFTING
+    const float drho = rho - 1.0f;
+    meq[0] = drho;
+    meq[1] = -2.0f * drho + 3.0f * rho * u2;
+    meq[2] = drho - 3.0f * rho * u2;
+#else
+    const float p = rho * (1.0f / 3.0f);
     meq[0] = 3.0f * p;
     meq[1] = -6.0f * p + 3.0f * rho * u2;
     meq[2] =  3.0f * p - 3.0f * rho * u2;
+#endif
     meq[3] = rho * ux;
     meq[4] = -rho * ux;
     meq[5] = rho * uy;
@@ -108,14 +101,9 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
 
     #pragma unroll
     for (int i = 0; i < 9; i++) {
-        g[i] += c_tau * Fin[i];
+        g[i] += 0.5f * Fin[i];
     }
 
-#if USE_DDF_SHIFTING
-    g[0] -= (float)d_w[0]; g[1] -= (float)d_w[1]; g[2] -= (float)d_w[2];
-    g[3] -= (float)d_w[3]; g[4] -= (float)d_w[4]; g[5] -= (float)d_w[5];
-    g[6] -= (float)d_w[6]; g[7] -= (float)d_w[7]; g[8] -= (float)d_w[8];
-#endif
 }
 
 __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
@@ -133,11 +121,7 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
                                             const float* __restrict__ Fin,
                                             float omega)
 {
-    // MRT is not shift-invariant.  Temporarily unshift to physical space.
+    // g and compute_feq use the same shifted/unshifted representation.
-#if USE_DDF_SHIFTING
-    #pragma unroll
-    for (int i = 0; i < 19; i++) g[i] += (float)d_w[i];
-#endif
 
     float feq[19];
     compute_feq(rho, ux, uy, uz, feq);
@@ -145,7 +129,7 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
     const float s_nu   = omega;
     const float s_bulk = 1.2f;
     const float s_high = 1.4f;
-    const float c_tau  = 1.0f - 0.5f * omega;
+
 
     const float one_minus_s_nu   = 1.0f - s_nu;
     const float one_minus_s_bulk = 1.0f - s_bulk;
@@ -154,7 +138,8 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
     float neq[19];
     #pragma unroll
     for (int i = 0; i < 19; i++) {
-        neq[i] = g[i] - feq[i];
+        // Half-source before relaxation; the other half is added below.
+        neq[i] = g[i] - feq[i] + 0.5f * Fin[i];
     }
 
     float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
@@ -214,13 +199,9 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
              + one_minus_s_nu   * neq_dev
              + one_minus_s_bulk * neq_iso
              + one_minus_s_high * neq_h
-             + c_tau * Fin[i];
+             + 0.5f * Fin[i];
     }
 
-#if USE_DDF_SHIFTING
-    #pragma unroll
-    for (int i = 0; i < 19; i++) g[i] -= (float)d_w[i];
-#endif
 }
 
 #endif  // NQ

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

@@ -11,7 +11,7 @@
 //
 // Forcing contract:
 //   Fin[i] is the raw Guo source term without the (1 - ω/2) prefactor.
-//   This operator applies the prefactor internally, matching collide_srt().
+//   Even/odd source components use their own (1 - omega_plus/minus / 2).
 // ============================================================================
 
 #ifndef CELERIS_OPERATORS_COLLISION_TRT_CUH
@@ -33,6 +33,7 @@ __device__ __forceinline__ void collide_trt(float* __restrict__ f,
     const float wp = omega;
     const float wm = compute_omega_minus(wp);
     const float c_tau = 1.0f - 0.5f * wp;
+    const float c_odd = 1.0f - 0.5f * wm;
 
     f[0] = f[0] - wp * (f[0] - feq[0]) + c_tau * Fin[0];
 
@@ -49,8 +50,10 @@ __device__ __forceinline__ void collide_trt(float* __restrict__ f,
         const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb));
         const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb));
 
-        f[i]  = fi  - wp * delta_s - wm * delta_a + c_tau * Fin[i];
+        const float force_s = c_tau * 0.5f * (Fin[i] + Fin[ib]);
-        f[ib] = fb  - wp * delta_s + wm * delta_a + c_tau * Fin[ib];
+        const float force_a = c_odd * 0.5f * (Fin[i] - Fin[ib]);
+        f[i]  = fi  - wp * delta_s - wm * delta_a + force_s + force_a;
+        f[ib] = fb  - wp * delta_s + wm * delta_a + force_s - force_a;
     }
 }
 

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

@@ -3,7 +3,8 @@
 //
 // Velocity correction:   u* = u + F·dt / (2ρ)
 // Forcing term:
-//   Fin[i] = (1 - ω/2) · w_i · [ (c_i - u*)/cs² + (c_i · u*)/cs⁴ · c_i ] · F
+//   Fin[i] = w_i · [ (c_i - u*)/cs² + (c_i · u*)/cs⁴ · c_i ] · F
+//   The collision operator applies its relaxation-dependent prefactor.
 //   simplified: Fin[i] = 9·w_i · [ (c_i·F)(c_i·u* + 1/3) - u*·F/3 ]
 // ============================================================================
 

src/CelerisLab/lbm/kernels/refinement.cu

@@ -0,0 +1,120 @@
+// Fixed, unforced internal 2:1 patch. Separate compilation entry;
+// default OneStep/channel kernels are deliberately not included.
+#include <stdint.h>
+#include <math.h>
+#include "config.h"
+#include "core/precision.cuh"
+#include "core/descriptors.cuh"
+#include "core/layout.cuh"
+#include "core/params.cuh"
+#include "core/flags.cuh"
+#include "core/obs.cuh"
+#include "operators/macro.cuh"
+#include "operators/collision_srt.cuh"
+#include "streaming/pull_double_buffer.cuh"
+#include "boundary/planar_velocity.cuh"
+#include "boundary/curved_boundary.cuh"
+
+extern "C" {
+#include "step/curved_double.cu"
+}
+
+#if DIM != 2 || NQ != 9 || COLLISION_MODEL != 0 || STORE_PRECISION != 0 || USE_DDF_SHIFTING || USE_LES || STREAMING_MODEL != 0
+#error "Refinement requires unshifted FP32 D2Q9 SRT double buffers without LES"
+#endif
+
+// Enabled only in the validation build; no diagnostic arrays in normal builds.
+#ifdef CELERIS_REFINEMENT_AUDIT
+extern "C" {
+__device__ unsigned int collision_counts[TOTAL_CELLS];
+__device__ unsigned int explosion_counts[TOTAL_CELLS];
+__device__ unsigned int incoming_counts[NQ * TOTAL_CELLS];
+}
+#endif
+
+__device__ __forceinline__ void patch_collide(float* f, unsigned long k) {
+    float rho, ux, uy, feq[NQ];
+    compute_rho_u(f, rho, ux, uy);
+    compute_feq(rho, ux, uy, feq);
+    collide_srt_no_force(f, feq, d_params.omega);
+#ifdef CELERIS_REFINEMENT_AUDIT
+    ++collision_counts[k];
+#endif
+}
+
+extern "C" __global__ void PatchExplode(
+    const float* coarse, float* fine, const uint32_t* targets,
+    const uint32_t* donors, unsigned int count, unsigned long coarse_cells)
+{
+    unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
+    if (row >= count) return;
+    for (unsigned int q = 0; q < NQ; ++q)
+        fine[index_f(targets[row], q)] = coarse[q * coarse_cells + donors[row]];
+#ifdef CELERIS_REFINEMENT_AUDIT
+    ++explosion_counts[targets[row]];
+#endif
+}
+
+extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint16_t* flags) {
+    unsigned int x, y; unsigned long k;
+    index_from_thread(x, y, k);
+    if (x >= NX || y >= NY) return;
+    float f[NQ];
+    // Outer halo is invalid after EACH hop, including when reusing the buffer
+    // that held explosion data. Inner halo only streams; it never collides.
+    if (x < 2 || y < 2 || x >= NX - 2 || y >= NY - 2
+            || (flags != nullptr && is_solid(flags[k]))) {
+        for (int q = 0; q < NQ; ++q) f[q] = NAN;
+    } else {
+        unsigned long j[NQ];
+        compute_neighbors(x, y, j);  // These interior targets never wrap.
+        stream_pull_load(k, f, in, j);
+        if (x >= 4 && y >= 4 && x < NX - 4 && y < NY - 4)
+            patch_collide(f, k);
+    }
+    stream_pull_store(k, f, out);
+}
+
+extern "C" __global__ void PatchCoarseStep(
+    const float* in, float* out, const float* fine,
+    const uint64_t* targets, const uint64_t* children, unsigned int count,
+    unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1,
+    unsigned int wall_aware)
+{
+    unsigned int x, y; unsigned long k;
+    index_from_thread(x, y, k);
+    if (x >= NX || y >= NY) return;
+    // Both buffers were poisoned at initialization; covered cells stay invalid.
+    if (x >= x0 && x < x1 && y >= y0 && y < y1) return;
+    if (wall_aware && (y == 0 || y == NY - 1)) return;
+    unsigned long j[NQ];
+    compute_neighbors(x, y, j);
+    float f[NQ];
+    stream_pull_load(k, f, in, j);
+    if (x >= x0 - 1 && x <= x1 && y >= y0 - 1 && y <= y1) {
+        for (unsigned int q = 1; q < NQ; ++q) {
+            uint64_t slot = index_f(k, q);
+            unsigned int lo = 0, hi = count;
+            while (lo < hi) {
+                unsigned int mid = lo + (hi - lo) / 2;
+                if (targets[mid] < slot) lo = mid + 1;
+                else hi = mid;
+            }
+            if (lo < count && targets[lo] == slot) {
+                const uint64_t* s = children + 4ul * lo;
+                f[q] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]);
+#ifdef CELERIS_REFINEMENT_AUDIT
+                ++incoming_counts[slot];
+#endif
+            }
+        }
+    }
+    // Includes the eight grazing-corner links with UNCOVERED coarse donors.
+    // All replacement arrivals must be in place before computing rho/u.
+    if (wall_aware) {
+        if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x));
+        else if (y == NY - 2) reconstruct_top_reservoir(f, RHO);
+    }
+    patch_collide(f, k);
+    stream_pull_store(k, f, out);
+}

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

@@ -15,15 +15,7 @@
 //     (pre-streaming semantics = matching double-buffer fi_in[k, dir]).
 //   SensorKernel / ForceRegionKernel run after EsoPullStep, using
 //     t = step_count+1 for post-streaming physical-value semantics.
-//   Wall boundary functions (bounce_back.cuh) use load_physical_dir_esopull
+//   Y-wall donors are staged by PrepareEsoPullBoundaries before any backing writes.
-//     inside EsoPullStep (post-streaming, matching load_f_esopull output).
-//
-// NOTE: The wall boundary helpers use post-streaming reads (unlike Bouzidi)
-//       because they operate inside EsoPullStep on f[] already reconstructed
-//       by load_f_esopull.  Replacing post-streaming f[i] with post-streaming
-//       f[opp(i)] is the correct half-way bounce-back for EsoPull.
-//
-// Verified: Kan99b K2 (60k steps, MRT) — bit-identical to double-buffer.
 // ============================================================================
 
 #ifndef CELERIS_STEP_AUX_KERNELS_CU
@@ -32,54 +24,7 @@
 #include "../streaming/esopull_semantic_helpers.cuh"
 #include "../boundary/esopull_curved_helpers.cuh"
 
-__device__ __forceinline__ float action_omega(const float* action, int body_id)
+#include "curved_double.cu"
-{
-    if (action == nullptr || body_id < 0) return 0.0f;
-    unsigned int base = (unsigned int)body_id * (unsigned int)OBS_BODY_SLOT_FLOATS;
-    return action[base + (unsigned int)(OBS_BODY_SLOT_FLOATS - 1)];
-}
-
-__global__ void CurvedBoundaryKernel(
-    fpxx* fi,
-    const unsigned int* cl_fluid_idx,
-    const unsigned char* cl_dir,
-    const float* cl_q,
-    const float* cl_rx,
-    const float* cl_ry,
-    const float* cl_rz,
-    const unsigned char* cl_fallback_class,
-    const int* cl_body_id,
-    const float* action,
-    float* obs,
-    unsigned int n_curved)
-{
-    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
-    if (tid >= n_curved) return;
-
-    const unsigned long k_f = (unsigned long)cl_fluid_idx[tid];
-    const unsigned int dir = (unsigned int)cl_dir[tid];
-    const float q = cl_q[tid];
-    const int bid = cl_body_id[tid];
-    const float rx = cl_rx[tid];
-    const float ry = cl_ry[tid];
-    const unsigned char fallback_class = cl_fallback_class[tid];
-    const float omega = action_omega(action, bid);
-
-#if DIM == 2
-    const float Uw = -omega * ry;
-    const float Vw =  omega * rx;
-    (void)cl_rz;
-    apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid);
-#elif DIM == 3
-    const float rz = cl_rz[tid];
-    const float Uw = -omega * ry;
-    const float Vw =  omega * rx;
-    const float Ww = 0.0f;
-    apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, Ww,
-                       rx, ry, rz, fallback_class, obs, bid);
-#endif
-}
-// End of CurvedBoundaryKernel
 
 // ---------------------------------------------------------------------------
 // PrepareEsoPullCurvedKernel -- read-only curved prep for EsoPull
@@ -146,7 +91,12 @@ __global__ void PrepareEsoPullCurvedKernel(
     } else {
         f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback);
     }
-    f_reflected += compute_moving_wall_addend(dir, q, fallback, Uw, Vw);
+    float rho_wall = USE_DDF_SHIFTING ? 1.0f : 0.0f;
+    if (Uw != 0.0f || Vw != 0.0f) {
+        for (unsigned int i = 0; i < NQ; ++i)
+            rho_wall += load_physical_dir_pre_streaming(k_f, i, fi, j_f, t);
+    }
+    f_reflected += compute_moving_wall_addend(dir, q, fallback, Uw, Vw, rho_wall);
 
     unsigned long dst_cell;
     unsigned int dst_slot;
@@ -249,14 +199,13 @@ __global__ void SensorKernel(
 }
 
 // ---------------------------------------------------------------------------
-// ForceRegionKernel — operator-split local Guo forcing correction
+// ForceRegionKernel — operator-split exact-difference force-density impulse
 //
-// Runs after the main collide+stream step.  Applies a Guo forcing term to
+// Runs after the main step. The raw source evaluated at u + F/(2 rho)
-// cells within a force_region footprint.  This is an operator-split correction,
+// equals feq(rho, u + F/rho) - feq(rho, u) for the quadratic equilibrium.
-// not inline in OneStep's collide_dispatch.
+// No collision prefactor: collision has already happened. This preserves
-//
+// density and adds exactly F momentum, independently of viscosity. It is
-// The Guo formula and c_tau = 1 - omega/2 conventions match collide_srt,
+// a split impulse, not the collision-integrated global Guo discretization.
-// collide_mrt, and collide_trt exactly.
 // ---------------------------------------------------------------------------
 __device__ __forceinline__ float action_force_x(const float* action, int body_id)
 {
@@ -319,11 +268,9 @@ __global__ void ForceRegionKernel(
     float Fin[NQ];
     compute_guo_forcing(ux, uy, fx, fy, Fin);
 
-    float omega = d_params.omega;
-    float c_tau = 1.0f - 0.5f * omega;
     #pragma unroll
     for (int i = 0; i < NQ; i++)
-        f[i] += c_tau * Fin[i];
+        f[i] += Fin[i];
 #elif DIM == 3
     float rho_n, ux, uy, uz;
     compute_rho_u(f, rho_n, ux, uy, uz);
@@ -334,11 +281,9 @@ __global__ void ForceRegionKernel(
     float Fin[NQ];
     compute_guo_forcing(ux, uy, uz, fx, fy, 0.0f, Fin);
 
-    float omega = d_params.omega;
-    float c_tau = 1.0f - 0.5f * omega;
     #pragma unroll
     for (int i = 0; i < NQ; i++)
-        f[i] += c_tau * Fin[i];
+        f[i] += Fin[i];
 #endif
 
     // Write back with correct semantics

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

@@ -0,0 +1,61 @@
+// Shared compact double-buffer curved boundary; caller supplies core/BFL headers.
+#ifndef CELERIS_CURVED_DOUBLE_CU
+#define CELERIS_CURVED_DOUBLE_CU
+
+__device__ __forceinline__ float action_omega(const float* action, int body_id)
+{
+    if (action == nullptr || body_id < 0) return 0.0f;
+    unsigned int base = (unsigned int)body_id * (unsigned int)OBS_BODY_SLOT_FLOATS;
+    return action[base + (unsigned int)(OBS_BODY_SLOT_FLOATS - 1)];
+}
+
+__global__ void CurvedBoundaryKernel(
+    fpxx* fi,
+    const unsigned int* cl_fluid_idx,
+    const unsigned char* cl_dir,
+    const float* cl_q,
+    const float* cl_rx,
+    const float* cl_ry,
+    const float* cl_rz,
+    const unsigned char* cl_fallback_class,
+    const int* cl_body_id,
+    const float* action,
+    float* obs,
+    unsigned int n_curved,
+    const float* surface_velocity = nullptr, float* link_exchange = nullptr)
+{
+    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_curved) return;
+
+    const unsigned long k_f = (unsigned long)cl_fluid_idx[tid];
+    const unsigned int dir = (unsigned int)cl_dir[tid];
+    const float q = cl_q[tid];
+    const int bid = cl_body_id[tid];
+    const float rx = cl_rx[tid];
+    const float ry = cl_ry[tid];
+    const unsigned char fallback_class = cl_fallback_class[tid];
+    const float omega = action_omega(action, bid);
+
+#if DIM == 2
+    float Uw = -omega * ry;
+    float Vw =  omega * rx;
+    if (surface_velocity != nullptr) {
+        // Independent prescribed fluid velocity on fixed geometry.
+        const float u = surface_velocity[2*tid], v = surface_velocity[2*tid+1];
+        if (u != 0.0f || v != 0.0f) { Uw += u; Vw += v; }
+    }
+    (void)cl_rz;
+    apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, rx, ry, fallback_class, obs, bid,
+                       link_exchange == nullptr ? nullptr : link_exchange + 3*tid);
+#elif DIM == 3
+    const float rz = cl_rz[tid];
+    const float Uw = -omega * ry;
+    const float Vw =  omega * rx;
+    const float Ww = 0.0f;
+    apply_bouzidi_link(dir, q, k_f, fi, Uw, Vw, Ww,
+                       rx, ry, rz, fallback_class, obs, bid);
+#endif
+}
+// End of CurvedBoundaryKernel
+
+#endif

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

@@ -47,10 +47,8 @@ __global__ void DecodeCellsToPhysical(
 // ---------------------------------------------------------------------------
 // EncodePhysicalToCells:  physical DDF -> backing layout (compact list)
 //
-// Reads physical (pre-collision) DDF from *fi_phys*, applies collision
+// Inverse of DecodeCellsToPhysical at the SAME completed step count.
-// to produce post-collision DDF, then encodes into the EsoPull backing
+// Conversion only: no collision or time advancement, including unchanged halo.
-// layout on *fi*.  The backing layout stores post-collision values,
-// matching what store_f_esopull would produce at step *t*.
 // ---------------------------------------------------------------------------
 __global__ void EncodePhysicalToCells(
     fpxx* fi,
@@ -70,10 +68,9 @@ __global__ void EncodePhysicalToCells(
     for (int i = 0; i < NQ; i++)
         f[i] = fi_phys[tid * NQ + i];
 
-    // Collide: physical (pre-collision) -> post-collision
+#if USE_DDF_SHIFTING
-    float rho_n, ux, uy;
+    for (int i = 0; i < NQ; ++i) f[i] -= (float)d_w[i];
-    compute_rho_u(f, rho_n, ux, uy);
+#endif
-    collide_dispatch(f, rho_n, ux, uy);
 
     store_physical_node_esopull(k, f, fi, j, t);
 }

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

@@ -6,98 +6,41 @@
 #ifndef CELERIS_STEP_INIT_FLOW_CU
 #define CELERIS_STEP_INIT_FLOW_CU
 
-// ---- Helper: write equilibrium for a given velocity ----
+// Reuse the canonical equilibrium, including its cancellation-free shifted
+// expression. Recomputing physical feq then subtracting w erases tiny speeds.
 __device__ __forceinline__ void write_equilibrium(
     fpxx* fi, unsigned long k, float u_init)
 {
-    for (int i = 0; i < NQ; i++) {
+    float f[NQ];
-        float cu = (float)d_cx[i] * u_init;
+#if DIM == 2
-        float val = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu
+    compute_feq(RHO, u_init, 0.0f, f);
-                                      - 1.5f*u_init*u_init);
+#else
-#if USE_DDF_SHIFTING
+    compute_feq(RHO, u_init, 0.0f, 0.0f, f);
-        val -= d_w[i];
 #endif
-        store_ddf(fi, index_f(k, (unsigned int)i), val);
+    for (int i = 0; i < NQ; ++i)
-    }
+        store_ddf(fi, index_f(k, (unsigned int)i), f[i]);
 }
 
-__device__ __forceinline__ void write_rest_equilibrium(
+__device__ __forceinline__ void write_rest_equilibrium(fpxx* fi, unsigned long k)
-    fpxx* fi, unsigned long k)
 {
-    for (int i = 0; i < NQ; i++) {
+    write_equilibrium(fi, k, 0.0f);
-        float val = d_w[i] * RHO;
-#if USE_DDF_SHIFTING
-        val -= d_w[i];
-#endif
-        store_ddf(fi, index_f(k, (unsigned int)i), val);
-    }
 }
 
-// ---- EsoPull layout helpers: equilibrium in EsoPull backing layout at t=0 ----
+// Init is load-compatible at t=0. Same-time encoding is the inverse of that
-// In EsoPull t=0 (even step), store_f_esopull writes:
+// read, unlike store_f_esopull(0), which advances the streaming interpretation.
-//   f[i_odd]  → fi[ j[i_odd],  i_odd      ]   (neighbor, forward slot)
-//   f[i_even] → fi[ k,         i_even      ]   (local, reverse slot)
-// These helpers write equilibrium values directly in the correct EsoPull
-// layout so that load_f_esopull(t=0) reconstructs the correct f_eq[].
 #if DIM == 2
 __device__ __forceinline__ void write_equilibrium_esopull(
     fpxx* fi, unsigned long k, unsigned int x, unsigned int y, float u_init)
 {
-    unsigned long j[NQ];
+    unsigned long j[NQ]; compute_neighbors(x, y, j);
-    compute_neighbors(x, y, j);
+    float f[NQ]; compute_feq(RHO, u_init, 0.0f, f);
-
+    store_physical_node_esopull(k, f, fi, j, 0ul);
-    // --- rest particle: always local, slot 0 ---
-    float val0 = d_w[0] * RHO * (1.0f - 1.5f * u_init * u_init);
-#if USE_DDF_SHIFTING
-    val0 -= d_w[0];
-#endif
-    store_ddf(fi, index_f(k, 0u), val0);
-
-    // --- paired directions (i=1,3,5,7): {E,W}, {N,S}, {NE,SW}, {SE,NW} ---
-    for (int i = 1; i < NQ; i += 2) {
-        float cu_odd  = (float)d_cx[i] * u_init;
-        float val_odd = d_w[i] * RHO
-            * (1.0f + 3.0f * cu_odd + 4.5f * cu_odd * cu_odd
-               - 1.5f * u_init * u_init);
-        float cu_even = (float)d_cx[i + 1] * u_init;
-        float val_even = d_w[i + 1] * RHO
-            * (1.0f + 3.0f * cu_even + 4.5f * cu_even * cu_even
-               - 1.5f * u_init * u_init);
-#if USE_DDF_SHIFTING
-        val_odd  -= d_w[i];
-        val_even -= d_w[i + 1];
-#endif
-        // t=0 even step (load_f_esopull): f[i] reads from local slot i+1,
-        // f[i+1] reads from neighbor slot i.
-        // Therefore: fi[n, i+1] = f_eq[i]  (odd goes to local reverse),
-        //            fi[j[i], i] = f_eq[i+1] (even goes to neighbor forward)
-        store_ddf(fi, index_f(k,    (unsigned int)(i + 1)), val_odd);
-        store_ddf(fi, index_f(j[i], (unsigned int)i),          val_even);
-    }
 }
 
 __device__ __forceinline__ void write_rest_equilibrium_esopull(
     fpxx* fi, unsigned long k, unsigned int x, unsigned int y)
 {
-    unsigned long j[NQ];
+    write_equilibrium_esopull(fi, k, x, y, 0.0f);
-    compute_neighbors(x, y, j);
-
-    float val0 = d_w[0] * RHO;
-#if USE_DDF_SHIFTING
-    val0 -= d_w[0];
-#endif
-    store_ddf(fi, index_f(k, 0u), val0);
-
-    for (int i = 1; i < NQ; i += 2) {
-        float val_odd  = d_w[i] * RHO;
-        float val_even = d_w[i + 1] * RHO;
-#if USE_DDF_SHIFTING
-        val_odd  -= d_w[i];
-        val_even -= d_w[i + 1];
-#endif
-        store_ddf(fi, index_f(k,    (unsigned int)(i + 1)), val_odd);
-        store_ddf(fi, index_f(j[i], (unsigned int)i),          val_even);
-    }
 }
 #endif // DIM == 2
 
@@ -106,10 +49,12 @@ __device__ __forceinline__ uint16_t channel_flag_from_coords(
 {
     if (y == 0u || y == (unsigned int)(NY - 1))
         return (uint16_t)(FLAG_SOLID | FLAG_BC_WALL);
+#if Y_WALL_BC != 2
     if (x == 0u)
         return (uint16_t)(FLAG_SOLID | FLAG_BC_INLET);
     if (x == (unsigned int)(NX - 1))
         return (uint16_t)(FLAG_SOLID | FLAG_BC_OUTLET);
+#endif
     return FLAG_FLUID;
 }
 
@@ -142,7 +87,21 @@ __global__ void InitTubeFlow_v2(uint16_t* flag, fpxx* fi)
     if (is_solid(fl)) {
         write_rest_equilibrium(fi, k);
     } else {
+#if Y_WALL_BC == 2
+        // Initial g(0)=feq: both boundary velocities are already satisfied.
+        // Start from linear Couette, not the sought exponential steady profile.
+        float f[NQ];
+#if WALL_SLOT_WIDTH > 0 && WALL_SLOT_WIDTH < NX
+        const float t = (float)(y-1u)/(float)(NY-3);
+        const float v = (1.0f-t)*bottom_wall_normal_velocity(x)+t*TOP_WALL_NORMAL_VELOCITY;
+        compute_feq(RHO, 0.0f, v, f);
+#else
+        compute_feq(RHO, U0*(float)(y-1u)/(float)(NY-3), WALL_NORMAL_VELOCITY, f);
+#endif
+        stream_pull_store(k, f, fi);
+#else
         write_equilibrium(fi, k, inlet_target_u((float)y));
+#endif
     }
 #elif DIM == 3
     unsigned int x, y, z; unsigned long k;
@@ -195,4 +154,4 @@ __global__ void InitEsoPull(uint16_t* flag, fpxx* fi)
 #endif
 }
 
-#endif // CELERIS_STEP_INIT_FLOW_CU
+#endif // CELERIS_STEP_INIT_FLOW_CU

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

@@ -88,7 +88,16 @@ void OneStep(
         bounce_back_swap(f);
 
     if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) {
-#if Y_WALL_BC == 1
+#if Y_WALL_BC == 2
+        // Reconstruct f(t+dt) on the physical plane, then collide below.
+        // This replaces, rather than compounds, the halfway wall correction.
+        if (y == 1u) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x));
+#if TOP_WALL_BC == 1
+        else reconstruct_top_reservoir(f, RHO);
+#else
+        else reconstruct_top_velocity(f, U0, TOP_WALL_NORMAL_VELOCITY);
+#endif
+#elif Y_WALL_BC == 1
         apply_wall_freeslip_d2q9_y_pull(y, f, fi_in, k);
 #else
         apply_wall_bb_d2q9_y_pull(y, f, fi_in, k);

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

@@ -13,7 +13,7 @@
 #if DIM == 2
 __device__ __forceinline__ void apply_boundary_esopull(
     float* f, uint16_t fl, unsigned int x, unsigned int y,
-    const fpxx* fi, unsigned long t)
+    const float* boundary_values)
 {
     // Curved-BC nodes have been pre-processed by PrepareEsoPullCurvedKernel
     // + ApplyEsoPullCurvedKernel, so is_curved() cells here have valid backed
@@ -24,10 +24,10 @@ __device__ __forceinline__ void apply_boundary_esopull(
     bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
 
     if (is_inlet(fl) && interior_y) {
-        apply_inlet_esopull_d2q9(f, x, y, fi, t);
+        apply_inlet_esopull_d2q9(f, x, y, boundary_values + 6u * NX + y * NQ);
     }
     else if (is_outlet(fl) && interior_y) {
-        apply_outlet_esopull_d2q9(f, x, y, fi, t);
+        apply_outlet_esopull_d2q9(f, boundary_values + 6u * NX + (NY + y) * NQ);
     }
     else {
         bounce_back_swap(f);
@@ -56,6 +56,43 @@ __device__ __forceinline__ void apply_boundary_esopull_3d(
 }
 #endif
 
+// One read-only preparation for all D2 planar boundary donors, before curved
+// or in-place writes. [6*NX wall values][2*NY*NQ x-boundary donor values].
+// t=0 uses the initialized physical state; later t uses spatial postcollision
+// values, matching the double-buffer donor reads. Slots retain DDF shifting.
+__global__ void PrepareEsoPullBoundaries(const fpxx* fi, float* boundary_values,
+                                        unsigned long t)
+{
+#if DIM == 2
+    const unsigned int n = threadIdx.x + blockIdx.x * blockDim.x;
+    const unsigned int side = blockIdx.y;
+    if (side >= 2u) return;
+    if (n < NX) {
+        const unsigned int y = side == 0u ? 1u : NY - 2u;
+        const unsigned long k = linear_index(n, y);
+        unsigned long j[NQ]; compute_neighbors(n, y, j);
+        const unsigned int lower[3] = {4u, 6u, 7u};
+        const unsigned int upper[3] = {3u, 5u, 8u};
+        for (unsigned int a = 0; a < 3u; ++a) {
+            unsigned int donor = side == 0u ? lower[a] : upper[a];
+#if Y_WALL_BC == 1
+            if (a > 0u) donor = side == 0u ? lower[3u-a] : upper[3u-a];
+#endif
+            boundary_values[(side * NX + n) * 3u + a] =
+                load_physical_dir_pre_streaming(k, donor, fi, j, t);
+        }
+    }
+    if (n < NY) {
+        const unsigned int x = side == 0u ? 1u : NX - 2u;
+        const unsigned long k = linear_index(x, n);
+        unsigned long j[NQ]; compute_neighbors(x, n, j);
+        float* donor = boundary_values + 6u * NX + (side * NY + n) * NQ;
+        for (unsigned int i = 0; i < NQ; ++i)
+            donor[i] = load_physical_dir_pre_streaming(k, i, fi, j, t);
+    }
+#endif
+}
+
 // ============================================================================
 // EsoPullStep — Esoteric-Pull single-buffer step kernel
 // ============================================================================
@@ -63,7 +100,7 @@ __global__ __launch_bounds__(NT, 2)
 void EsoPullStep(
     fpxx* fi, uint16_t* flag,
     float* action, float* obs,
-    unsigned long t)
+    const float* boundary_values, unsigned long t)
 {
 #if DIM == 2
     unsigned int x, y; unsigned long k;
@@ -87,18 +124,17 @@ void EsoPullStep(
 
     // Boundary nodes
     if (is_solid(fl))
-        apply_boundary_esopull(f, fl, x, y, fi, t);
+        apply_boundary_esopull(f, fl, x, y, boundary_values);
 
     if (is_obstacle(fl) && !is_curved(fl))
         bounce_back_swap(f);
 
     // ---- y-wall adjacent row correction ----
     if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) {
-#if Y_WALL_BC == 1
+        const unsigned int side = y == 1u ? 0u : 1u;
-        apply_wall_freeslip_d2q9_y_pull(y, f, fi, k, j, t);
+        const float* staged = boundary_values + (side * NX + x) * 3u;
-#else
+        if (side == 0u) { f[3] = staged[0]; f[5] = staged[1]; f[8] = staged[2]; }
-        apply_wall_bb_d2q9_y_pull(y, f, fi, k, j, t);
+        else { f[4] = staged[0]; f[6] = staged[1]; f[7] = staged[2]; }
-#endif
     }
 
     // Collision

src/CelerisLab/lbm/kernels/streaming/esopull_semantic_helpers.cuh

@@ -14,14 +14,14 @@
 //   ForceRegionKernel(t):              reads/writes fi with post-streaming semantics
 //   SensorKernel(t):                   reads fi with post-streaming semantics
 //   MacroscopicEsoPullKernel(t):       reads fi with post-streaming semantics
-//   Wall boundary functions(t):        reads fi with pre-streaming semantics
+//   PrepareEsoPullBoundaries(t):             stages pre-streaming donors before writes
 //
 // The "t" parameter always represents the completed step count corresponding
 // to the current physical interpretation of the backing layout.
 // Double-buffer mode: t = 0, raw fi[k,i] = physical direction i.
 //
-// Verified: Kan99b K2 (200k steps, MRT) — bit-identical to double-buffer for
+// The same-time read/write mapping is checked by actual-header regression.
-// Strouhal, CD, CL, C'L, C'D metrics (FP32, D=20/30).
+// A physical node here is the next pull state, not the spatial postcollision node.
 // ============================================================================
 
 #ifndef CELERIS_STREAMING_ESOPULL_SEMANTIC_HELPERS_CUH
@@ -130,14 +130,8 @@ __device__ __forceinline__ void store_physical_node_esopull(
     fpxx* fi, const unsigned long* j, unsigned long t)
 {
     store_ddf(fi, index_f(k, 0u), f[0]);
-    for (unsigned int i = 1; i < NQ; i += 2) {
+    for (unsigned int i = 1; i < NQ; ++i) {
-        if (t & 1ul) {
+        store_ddf(fi, index_f(esopull_src_cell(k, i, j), esopull_slot(i, t)), f[i]);
-            store_ddf(fi, index_f(j[i], i + 1), f[i]);
-            store_ddf(fi, index_f(k,    i),     f[i + 1]);
-        } else {
-            store_ddf(fi, index_f(j[i], i),     f[i]);
-            store_ddf(fi, index_f(k,    i + 1), f[i + 1]);
-        }
     }
 }
 

src/CelerisLab/lbm/refinement.py

@@ -0,0 +1,292 @@
+"""CPU preparation for one fixed D2Q9 2:1 patch and optional interior circles.
+
+Cell indices use y*nx+x; population slots use q*(nx*ny)+cell, matching
+core/layout.cuh. Bounds are half-open (x0, y0, x1, y1) in coarse cells.
+The four fine halo cells implement the pinned CPUREF volumetric geometry.
+Base topology is independent of boundaries; prepare_fine_boundary checks them.
+"""
+from dataclasses import dataclass, replace
+from numbers import Integral
+
+import numpy as np
+
+from ..config import LBMConfig
+from .descriptors import D2Q9_EX, D2Q9_EY, FLUID, SOLID
+
+ACTIVE, INNER_GHOST, OUTER_GHOST = 0, 1, 2
+HALO = 4
+
+
+def _readonly(values, dtype):
+    array = np.ascontiguousarray(values, dtype=dtype)
+    array.setflags(write=False)
+    return array
+
+
+@dataclass(frozen=True, eq=False)
+class FixedPatch2D:
+    """Layer metadata and compact exchange tables, constructed by from_config.
+
+    Shapes are (ny, nx), including the fine halo. ``fine_origin`` is its lower
+    cell edge in coarse physical units, with coarse dx=dt=1 and fine dx=dt=.5.
+    ``ghost_fine`` and ``ghost_coarse`` pair fine cells with coarse donors.
+    ``coalesce_coarse`` contains incoming population slots; each corresponding
+    row of ``coalesce_fine`` contains four fine slots in (dy, dx) row order.
+    Weights are all 1/4. Closed corner links are included even when the coarse
+    donor is uncovered. These tables contain no DDF data and allocate no GPU.
+
+    This is preparation, not permission to advance these levels with OneStep:
+    its channel boundaries and fused collision do not implement patch coupling.
+    """
+    coarse_shape: tuple[int, int]
+    bounds: tuple[int, int, int, int]
+    tau_coarse: float
+    tau_fine: float
+    fine_kind: np.ndarray
+    ghost_fine: np.ndarray
+    ghost_coarse: np.ndarray
+    coalesce_coarse: np.ndarray
+    coalesce_fine: np.ndarray
+
+    @classmethod
+    def from_config(cls, cfg: LBMConfig, bounds) -> "FixedPatch2D":
+        """Prepare maps without changing cfg or allocating coarse-sized tables.
+
+        The initial backend contract is D2Q9/SRT/FP32/unshifted/double_buffer,
+        without LES. Two coarse cells of margin keep all halo donors in bounds;
+        this does NOT certify that those donors are fluid or periodic.
+        """
+        checked = replace(cfg)
+        checked.validate()
+        if (not checked.is_d2q9 or checked.collision != "SRT"
+                or checked.streaming != "double_buffer"
+                or checked.store_precision != "FP32"
+                or checked.ddf_shifting or checked.les_enabled):
+            raise ValueError("FixedPatch2D requires D2Q9/SRT/FP32/unshifted/double_buffer without LES")
+        bounds = tuple(bounds)
+        if len(bounds) != 4 or any(isinstance(v, (bool, np.bool_))
+                                   or not isinstance(v, Integral) for v in bounds):
+            raise ValueError("bounds must be four integers: (x0, y0, x1, y1)")
+        x0, y0, x1, y1 = bounds = tuple(map(int, bounds))
+        nx, ny = checked.nx, checked.ny
+        if not (2 <= x0 < x1 <= nx - 2 and 2 <= y0 < y1 <= ny - 2):
+            raise ValueError("patch must be nonempty with a two-coarse-cell margin")
+        fx, fy = 2 * (x1 - x0) + 2 * HALO, 2 * (y1 - y0) + 2 * HALO
+        if max(nx * ny, fx * fy) > np.iinfo(np.uint32).max:
+            raise ValueError("layer cell counts exceed uint32 indexing")
+        tau_c = 3.0 * checked.viscosity + 0.5
+        tau_f = 6.0 * checked.viscosity + 0.5
+        if not checked.omega_min <= 1.0 / tau_f <= checked.omega_max:
+            raise ValueError("fine relaxation rate lies outside configured omega guards")
+
+        kind = np.full((fy, fx), OUTER_GHOST, dtype=np.uint8)
+        kind[2:-2, 2:-2] = INNER_GHOST
+        kind[HALO:-HALO, HALO:-HALO] = ACTIVE
+        gy, gx = np.nonzero(kind != ACTIVE)
+        ghost_fine = gy * fx + gx
+        ghost_coarse = (y0 + (gy - HALO) // 2) * nx + x0 + (gx - HALO) // 2
+
+        targets, sources = [], []
+        # Only the outside ring can receive fine-to-coarse populations.
+        # Integer doubled midpoints include all eight diagonal corner passes.
+        for y in range(y0 - 1, y1 + 1):
+            xs = range(x0 - 1, x1 + 1) if y in (y0 - 1, y1) else (x0 - 1, x1)
+            for x in xs:
+                for q in range(1, 9):
+                    mx, my = 2 * x + 1 - D2Q9_EX[q], 2 * y + 1 - D2Q9_EY[q]
+                    if 2 * x0 <= mx <= 2 * x1 and 2 * y0 <= my <= 2 * y1:
+                        targets.append(q * nx * ny + y * nx + x)
+                        a, b = HALO + 2 * (x - x0), HALO + 2 * (y - y0)
+                        sources.append([q * fx * fy + (b + dy) * fx + a + dx
+                                        for dy in (0, 1) for dx in (0, 1)])
+        return cls((ny, nx), bounds, tau_c, tau_f,
+                   _readonly(kind, np.uint8),
+                   _readonly(ghost_fine, np.uint32),
+                   _readonly(ghost_coarse, np.uint32),
+                   _readonly(targets, np.uint64),
+                   _readonly(sources, np.uint64).reshape(-1, 4))
+
+    @property
+    def fine_shape(self) -> tuple[int, int]:
+        return self.fine_kind.shape
+
+    @property
+    def fine_origin(self) -> tuple[float, float]:
+        return self.bounds[0] - HALO / 2, self.bounds[1] - HALO / 2
+
+    @property
+    def table_bytes(self) -> int:
+        return sum(a.nbytes for a in (self.fine_kind, self.ghost_fine,
+                                     self.ghost_coarse, self.coalesce_coarse,
+                                     self.coalesce_fine))
+
+    def check_synchronization(self, coarse_steps: int, fine_steps: int, *, phase: str) -> None:
+        """Check a proposed backend handoff; does not convert or collide data.
+
+        Production storage is g=C(f), while CPUREF returns arrival_f. A future
+        g-to-g coarse step must collide coarse cells once and fine cells twice;
+        fine halo ghosts never collide. Coalescence must precede coarse collision.
+        Converting reference f to g for comparison is an observation operation,
+        not an additional production collision. Counts here are completed steps.
+        """
+        if phase != "postcollision_g":
+            raise ValueError("patch backend handoff requires postcollision_g")
+        for value in (coarse_steps, fine_steps):
+            if isinstance(value, (bool, np.bool_)) or not isinstance(value, Integral) or value < 0:
+                raise ValueError("step counts must be nonnegative integers")
+        if fine_steps != 2 * coarse_steps:
+            raise ValueError("levels are not synchronized: fine_steps must equal 2*coarse_steps")
+
+
+@dataclass(frozen=True, eq=False)
+class FineBoundary2D:
+    """Fixed local projection plus auditable, finite transfer dependency sets.
+
+    transfer_reads rows are (hop, direction, cell), traced from coalescence at
+    hop 2 back to exploded G at hop 0. No previous coarse interval is traversed.
+    Cell footprints also include the whole explosion write/donor tables.
+    """
+    objects: tuple
+    flags: np.ndarray
+    curved_host: tuple
+    boundary_cells: np.ndarray
+    transfer_fine_cells: np.ndarray
+    transfer_coarse_cells: np.ndarray
+    transfer_reads: np.ndarray
+
+
+def _transfer_footprint(patch):
+    ny, nx = patch.fine_shape
+    cells = nx * ny
+    pending = [(2, *divmod(int(slot), cells)) for slot in patch.coalesce_fine.ravel()]
+    reads = set()
+    while pending:
+        hop, q, k = pending.pop()
+        if (hop, q, k) in reads:
+            continue
+        reads.add((hop, q, k))
+        if hop == 0:
+            continue
+        y, x = divmod(k, nx)
+        kind = patch.fine_kind[y, x]
+        if kind == OUTER_GHOST:
+            raise ValueError("Transfer reads an invalid post-hop outer ghost")
+        # A collided output depends on all nine arrivals; ghosts only stream.
+        for direction in (range(9) if kind == ACTIVE else (q,)):
+            sx, sy = x - D2Q9_EX[direction], y - D2Q9_EY[direction]
+            if not (0 <= sx < nx and 0 <= sy < ny):
+                raise ValueError("Transfer dependency leaves the fine allocation")
+            pending.append((hop - 1, direction, sy * nx + sx))
+    fine = set(map(int, patch.ghost_fine)) | {k for _, _, k in reads}
+    coarse = set(map(int, patch.ghost_coarse)) | set(map(
+        int, patch.coalesce_coarse % np.prod(patch.coarse_shape)))
+    return fine, coarse, _readonly(sorted(reads), np.int64)
+
+
+def prepare_fine_boundary(cfg, patch, objects, *, wall_aware=False):
+    """Project fixed circles and reject boundary/transfer overlap before CUDA.
+
+    Input centers/radii are in coarse physical coordinates (dx=dt=1).
+    SurfaceActuation retains its lattice-velocity meaning: input coarse units
+    are converted into independent fine-local objects with dx=dt=.5.
+    Original objects are not retained or mutated. No runtime geometry updates.
+    """
+    if not objects:
+        return None
+    checked = replace(cfg)
+    checked.validate()
+    if (patch.coarse_shape != (cfg.ny, cfg.nx) or not checked.is_d2q9
+            or checked.collision != "SRT" or checked.streaming != "double_buffer"
+            or checked.store_precision != "FP32" or checked.ddf_shifting or checked.les_enabled
+            or patch.tau_coarse != 3*checked.viscosity+.5):
+        raise ValueError("Fine boundary requires the matching fixed-patch configuration")
+    from ..body.geometry.circle import CircleGeometry
+    from ..body.objects import SimObject
+    from ..body.coupling.soa_packer import pack_cut_links_to_soa
+
+    ny, nx = patch.fine_shape
+    fine_cfg = replace(cfg, nx=nx, ny=ny, viscosity=2 * cfg.viscosity)
+    flags = np.full(nx * ny, FLUID, np.uint16)
+    active = patch.fine_kind.ravel() == ACTIVE
+    local, links = [], []
+    x0, y0, x1, y1 = patch.bounds
+    dx = dt = .5
+    for bid, obj in enumerate(objects):
+        if (type(obj.geometry) is not CircleGeometry or obj.is_sensor or obj.is_force_region
+                or any(getattr(obj.state, name) != 0 for name in ("vx", "vy", "omega"))
+                or any(vars(obj.control).values())):
+            raise ValueError("Fine boundaries require stationary circles without other roles")
+        cx, cy = obj.geometry.center
+        radius = obj.geometry._radius
+        if (not np.isfinite([cx, cy, radius]).all() or radius <= 0
+                or cx-radius < x0 or cx+radius > x1 or cy-radius < y0 or cy+radius > y1):
+            raise ValueError("Entire fine circle must lie inside the physical patch")
+        center = tuple((v-o)/dx-.5 for v, o in zip((cx, cy), patch.fine_origin))
+        spec = obj.surface_actuation
+        if spec is not None:
+            spec = replace(spec, normal_velocity=spec.normal_velocity * dt / dx)
+        item = SimObject(bid, geometry=CircleGeometry(*center, radius/dx),
+                         center=center, radius=radius/dx, surface_actuation=spec)
+        item.validate_surface_actuation(fine_cfg)
+        mask = item.get_flag_mask(nx, ny)
+        solid = (mask & SOLID) != 0
+        if not solid.any() or np.any(solid & ~active) or np.any(flags[solid] & SOLID):
+            raise ValueError("Fine solids must be resolved, disjoint and entirely active")
+        flags[solid] = mask[solid]
+        local.append(item)
+    for item in local:
+        cut = item.geometry.build_cut_links(nx, ny, domain_flags=flags)
+        if not cut:
+            raise ValueError("Fine circle has no resolved cut links")
+        if item.surface_actuation is not None:
+            _, aperture = item.surface_actuation.project(
+                np.array([link.rx for link in cut], np.float32),
+                np.array([link.ry for link in cut], np.float32))
+            if not aperture.any():
+                raise ValueError("Surface aperture hits no resolved cut links")
+        for link in cut:
+            link.body_id = item.obj_id
+        links.extend(cut)
+
+    boundary, writers, receivers = set(), set(), set()
+    for link in links:
+        r, q = link.fluid_idx, link.dir
+        y, x = divmod(r, nx)
+        sx, sy = x + D2Q9_EX[q], y + D2Q9_EY[q]
+        s = sy * nx + sx
+        if not (active[r] and flags[r] & FLUID and active[s] and flags[s] & SOLID):
+            raise ValueError("Cut-link receiver/source must be active fluid/solid")
+        writer = (s, q+1 if q % 2 else q-1)
+        if writer in writers:
+            raise ValueError("Duplicate fine boundary writer")
+        writers.add(writer)
+        receivers.add((r, q))
+        boundary.update((r, s))
+        if link.q < .5:
+            dx_, dy_ = x - D2Q9_EX[q], y - D2Q9_EY[q]
+            if not (0 <= dx_ < nx and 0 <= dy_ < ny and active[dy_*nx+dx_]):
+                raise ValueError("q<.5 donor crosses the fine active interface; fallback cannot hide it")
+            boundary.add(dy_*nx+dx_)
+    # Reverse coverage: every fluid pull across a solid neighbor needs a writer.
+    required = set()
+    for r in np.flatnonzero(active & ((flags & FLUID) != 0)):
+        for q in range(1, 9):
+            if flags[r + D2Q9_EX[q] + nx*D2Q9_EY[q]] & SOLID:
+                required.add((int(r), q))
+    if required != receivers:
+        raise ValueError("Fine solid pull coverage does not match the cut links")
+    fine, coarse, reads = _transfer_footprint(patch)
+    if fine & (boundary | set(map(int, np.flatnonzero(flags & SOLID)))):
+        raise ValueError("Boundary stencil overlaps the two-hop transfer footprint")
+    for k in coarse:
+        y, x = divmod(k, cfg.nx)
+        if x0 <= x < x1 and y0 <= y < y1:
+            raise ValueError("Transfer coarse donor/target is not an uncovered bulk cell")
+        # Coarse reconstruction targets y=1,ny-2 pull from y=0..2,ny-3..ny-1.
+        if wall_aware and (y <= 2 or y >= cfg.ny-3):
+            raise ValueError("Transfer footprint touches the coarse wall reconstruction stencil")
+    return FineBoundary2D(tuple(local), _readonly(flags, np.uint16),
+                          tuple(_readonly(a, a.dtype) for a in pack_cut_links_to_soa(links)),
+                          _readonly(sorted(boundary), np.uint32),
+                          _readonly(sorted(fine), np.uint32),
+                          _readonly(sorted(coarse), np.uint32), reads)

src/CelerisLab/lbm/refinement_stepper.py

@@ -0,0 +1,302 @@
+# CelerisLab/lbm/refinement_stepper.py
+"""Opt-in fixed patch backend with optional stationary fine-interior circles.
+
+The caller owns the CudaContext and must close this backend before that context.
+Only synchronized postcollision_g leaf populations are public. No checkpoint,
+snapshot, forcing, runtime parameter mutation or dynamic refinement is offered.
+"""
+from dataclasses import dataclass, replace
+from numbers import Integral
+import struct
+from threading import RLock
+
+import numpy as np
+import pycuda.driver as cuda
+
+from ..cuda import compiler_v2 as compiler
+from .refinement import ACTIVE, HALO, FixedPatch2D, prepare_fine_boundary
+from .descriptors import SOLID
+
+
+@dataclass
+class _Level:
+    shape: tuple
+    module: object
+    current: object
+    spare: object
+    kernel: object
+    grid: tuple
+
+
+class FixedPatchStepper:
+    """D2Q9/SRT/FP32, one unforced internal 2:1 patch in an x-periodic box.
+
+    ``coarse_boundary`` is ``periodic`` (default) or ``slot_reservoir``:
+    a bottom on-site slot and top fixed-density reservoir from ``cfg``.
+    Periodic mode does not apply the wall/slot boundary fields in ``cfg``.
+    Outer walls remain coarse; fine halo donors must be ordinary fluid.
+    ``fine_objects`` accepts fixed circles in coarse physical coordinates;
+    velocities are in coarse lattice units (dx=dt=1). Independent local copies
+    are projected at fine cell centers. Their complete stencil must be active
+    and disjoint from the actual two-hop transfer footprint.
+    These use the edge/cell-center frame (coarse centers I+.5), whereas ordinary Simulation objects use the integer-node frame.
+
+    ``initialize(coarse_g, fine_g)`` accepts physical Q-major (9,ny,nx) arrays
+    already in postcollision phase. It does not collide or build channel flags.
+    Every ``step()`` publishes only after both levels complete on one stream.
+    On a device/launch failure the instance becomes unusable: close and recreate
+    it, rather than pretending the partially advanced fine buffer rolled back.
+    """
+    def __init__(self, cfg, bounds, *, context, coarse_boundary="periodic", fine_objects=()):
+        if HALO != 4:
+            raise RuntimeError("refinement.cu requires a four-cell fine halo")
+        if coarse_boundary not in ("periodic", "slot_reservoir"):
+            raise ValueError("coarse_boundary must be periodic or slot_reservoir")
+        self.patch = FixedPatch2D.from_config(cfg, bounds)
+        self._wall_aware = coarse_boundary == "slot_reservoir"
+        if self._wall_aware:
+            if (cfg.y_wall_bc != "on_site_velocity" or cfg.top_wall_bc != "reservoir"
+                    or cfg.velocity != 0 or not 2 <= cfg.wall_slot_width < cfg.nx):
+                raise ValueError("slot_reservoir requires on_site_velocity, reservoir, "
+                                 "velocity=0 and 2 <= wall_slot_width < nx")
+            if self.patch.bounds[1] < 4 or self.patch.bounds[3] > cfg.ny - 4:
+                raise ValueError("slot_reservoir requires y0 >= 4 and y1 <= ny-4")
+        self.fine_boundary = prepare_fine_boundary(
+            cfg, self.patch, tuple(fine_objects), wall_aware=self._wall_aware)
+        self._fine_leaf = self.patch.fine_kind == ACTIVE
+        if self.fine_boundary is not None:
+            self._fine_leaf &= (self.fine_boundary.flags.reshape(self.patch.fine_shape) & SOLID) == 0
+        self.ctx = context
+        self._lock = RLock()
+        self._status = "new"
+        self._step_count = 0
+        self._allocations = []
+        self._levels = []
+        self._stream = None
+        self._fine_flags = self._fine_curved = self._curved = None
+        self._block = (cfg.threads_per_block, 1, 1)
+        fine_cfg = replace(cfg, nx=self.patch.fine_shape[1], ny=self.patch.fine_shape[0],
+                           viscosity=2 * cfg.viscosity)
+        if self._wall_aware:
+            fine_cfg = replace(fine_cfg, y_wall_bc="bounce_back", top_wall_bc="velocity",
+                               wall_slot_width=0, wall_normal_velocity=0.0, velocity=0.0)
+        with self.ctx.activate():
+            try:
+                self._stream = cuda.Stream()
+                for config, name, tau in ((replace(cfg), "PatchCoarseStep", self.patch.tau_coarse),
+                                          (fine_cfg, "PatchFineStep", self.patch.tau_fine)):
+                    config.validate()
+                    if not np.isclose(config.omega, 1 / tau, rtol=1e-14, atol=0):
+                        raise RuntimeError("level omega/tau mismatch")
+                    n_objects = (len(self.fine_boundary.objects)
+                                 if name == "PatchFineStep" and self.fine_boundary is not None else 0)
+                    build = compiler.generate_config(config, n_objects=n_objects)
+                    try:
+                        ptx = compiler.compile_kernel(self.ctx.sm_arch, source="refinement.cu",
+                                                      build_dir=build)
+                        module = compiler.load_module(ptx)
+                    finally:
+                        compiler.cleanup_build_dir(build)
+                    ptr, size = module.get_global("d_params")
+                    if size != struct.calcsize("=ffffI"):
+                        raise RuntimeError("LBMParams host/device layout mismatch")
+                    cuda.memcpy_htod(ptr, struct.pack("=ffffI", config.omega, 0, 0, 0, n_objects))
+                    shape = (9, config.ny, config.nx)
+                    size = int(np.prod(shape)) * 4
+                    self._levels.append(_Level(shape, module, self._allocate(size),
+                                              self._allocate(size), module.get_function(name),
+                                              ((config.nx + self._block[0] - 1) // self._block[0],
+                                               config.ny, 1)))
+                p = self.patch
+                order = np.argsort(p.coalesce_coarse)
+                self._ghosts = self._upload(p.ghost_fine)
+                self._donors = self._upload(p.ghost_coarse)
+                self._targets = self._upload(p.coalesce_coarse[order])
+                self._children = self._upload(p.coalesce_fine[order])
+                self._explode = self._levels[1].module.get_function("PatchExplode")
+                if self.fine_boundary is not None:
+                    from .curved_links import CurvedLinkSoA
+                    boundary = self.fine_boundary
+                    self._fine_flags = self._upload(boundary.flags)
+                    # SoA owns its allocations; never also register them in _allocations.
+                    self._fine_curved = CurvedLinkSoA()
+                    self._fine_curved.assign_host(*boundary.curved_host)
+                    self._fine_curved.assign_surface(boundary.objects, fine_cfg, boundary.flags)
+                    self._fine_curved.upload()
+                    self._curved = self._levels[1].module.get_function("CurvedBoundaryKernel")
+            except BaseException:
+                try:
+                    self.close()
+                except Exception:
+                    pass  # Preserve the original construction failure.
+                raise
+
+    def _allocate(self, size):
+        pointer = cuda.mem_alloc(size)
+        self._allocations.append(pointer)
+        return pointer
+
+    def _upload(self, array):
+        array = np.ascontiguousarray(array)
+        pointer = self._allocate(array.nbytes)
+        cuda.memcpy_htod(pointer, array)
+        return pointer
+
+    def _require(self, *states):
+        if self._status not in states:
+            raise RuntimeError(f"Fixed patch is {self._status}; expected {' or '.join(states)}")
+
+    @property
+    def step_count(self):
+        return self._step_count
+
+    def _mask_nonleaves(self, coarse, fine):
+        x0, y0, x1, y1 = self.patch.bounds
+        coarse[:, y0:y1, x0:x1] = np.nan
+        if self._wall_aware:
+            coarse[:, (0, -1), :] = np.nan
+        fine[:, ~self._fine_leaf] = np.nan
+
+    def initialize(self, coarse_g, fine_g):
+        """Upload an explicit G_0; invalid non-leaves are discarded."""
+        with self._lock, self.ctx.activate():
+            self._require("new", "ready")
+            arrays = [np.array(a, dtype=np.float32, order="C", copy=True)
+                      for a in (coarse_g, fine_g)]
+            if any(a.shape != level.shape for a, level in zip(arrays, self._levels)):
+                raise ValueError("Expected Q-major (9,ny,nx) arrays for both levels")
+            x0, y0, x1, y1 = self.patch.bounds
+            leaf = np.ones(self.patch.coarse_shape, dtype=bool)
+            leaf[y0:y1, x0:x1] = False
+            if self._wall_aware:
+                leaf[[0, -1]] = False
+            for array, mask in zip(arrays, (leaf, self._fine_leaf)):
+                if not np.isfinite(array[:, mask]).all() or (array[:, mask].sum(axis=0) <= 0).any():
+                    raise ValueError("Leaf populations must be finite with positive density")
+            self._mask_nonleaves(*arrays)
+            try:
+                for level, array in zip(self._levels, arrays):
+                    cuda.memcpy_htod(level.current, array)
+                    cuda.memcpy_htod(level.spare, array)
+                if self._fine_curved is not None and self._fine_curved.link_exchange_gpu is not None:
+                    cuda.memset_d8(self._fine_curved.link_exchange_gpu, 0, self._fine_curved.count * 12)
+                self._step_count = 0
+                self._status = "ready"
+            except BaseException:
+                self._status = "failed"
+                raise
+
+    def _apply_fine_boundary(self, buffer, **kw):
+        """Stage solid source slots; exchange buffer holds this fine substep only."""
+        cl = self._fine_curved
+        if cl is None:
+            if self.fine_boundary is not None:
+                raise RuntimeError("Fine boundary device buffers are missing")
+            return
+        if self._curved is None or cl.fluid_idx_gpu is None:
+            raise RuntimeError("Fine boundary device buffers are missing")
+        driven = any(obj.surface_actuation is not None for obj in self.fine_boundary.objects)
+        if driven and (cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None):
+            raise RuntimeError("Surface boundary device buffers are missing")
+        self._curved(buffer, cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu,
+                     cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu,
+                     np.uintp(0), np.uintp(0), np.uint32(cl.count),
+                     cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0),
+                     cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0),
+                     grid=((cl.count + self._block[0] - 1) // self._block[0], 1), **kw)
+
+    def step(self, n=1):
+        """Advance G→G: one explosion, two fine collisions, one coarse collision."""
+        if isinstance(n, bool) or not isinstance(n, Integral) or n < 0:
+            raise ValueError("n must be a nonnegative integer")
+        with self._lock, self.ctx.activate():
+            self._require("ready")
+            c, f = self._levels
+            p = self.patch
+            kw = dict(block=self._block, stream=self._stream)
+            try:
+                self._status = "advancing"
+                for _ in range(n):
+                    self._explode(c.current, f.current, self._ghosts, self._donors,
+                                  np.uint32(len(p.ghost_fine)), np.uint64(np.prod(p.coarse_shape)),
+                                  grid=((len(p.ghost_fine) + self._block[0] - 1) // self._block[0], 1),
+                                  **kw)
+                    flags = self._fine_flags if self._fine_flags is not None else np.uintp(0)
+                    if self.fine_boundary is not None and self._fine_flags is None:
+                        raise RuntimeError("Fine boundary flags are missing")
+                    self._apply_fine_boundary(f.current, **kw)
+                    f.kernel(f.current, f.spare, flags, grid=f.grid, **kw)
+                    self._apply_fine_boundary(f.spare, **kw)
+                    f.kernel(f.spare, f.current, flags, grid=f.grid, **kw)
+                    c.kernel(c.current, c.spare, f.current, self._targets, self._children,
+                             np.uint32(len(p.coalesce_coarse)),
+                             *(np.uint32(v) for v in p.bounds), np.uint32(self._wall_aware),
+                             grid=c.grid, **kw)
+                    self._stream.synchronize()
+                    c.current, c.spare = c.spare, c.current
+                    self._step_count += 1
+                self._status = "ready"
+            except BaseException:
+                self._status = "failed"
+                raise
+
+    def download(self):
+        """Return two synchronized G arrays, NaN on all non-leaves; no host cache."""
+        with self._lock, self.ctx.activate():
+            self._require("ready")
+            arrays = [np.empty(level.shape, np.float32) for level in self._levels]
+            try:
+                for level, array in zip(self._levels, arrays):
+                    cuda.memcpy_dtoh(array, level.current)
+            except BaseException:
+                self._status = "failed"
+                raise
+            self._mask_nonleaves(*arrays)
+            return tuple(arrays)
+
+    def close(self):
+        """Release owned resources; retain them for retry if context activation fails."""
+        with self._lock:
+            if self._status == "closed":
+                return
+            self._status = "failed"
+            with self.ctx.activate():
+                self._status = "closed"
+                error = None
+                try:
+                    if self._stream is not None:
+                        self._stream.synchronize()
+                except Exception as exc:
+                    error = exc
+                self._levels.clear()
+                self._explode = None
+                if getattr(self, "_fine_curved", None) is not None:
+                    try:
+                        self._fine_curved.free()
+                    except Exception as exc:
+                        error = error or exc
+                    self._fine_curved = None
+                self._curved = self._fine_flags = None
+                for pointer in self._allocations:
+                    try:
+                        pointer.free()
+                    except Exception as exc:
+                        error = error or exc
+                self._allocations.clear()
+                self._stream = None
+                if error is not None:
+                    raise error
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, *exc):
+        try:
+            self.close()
+        except Exception:
+            if exc[0] is None:
+                raise
+
+
+# Preserve the original entry point and its periodic default.
+PeriodicPatchStepper = FixedPatchStepper

src/CelerisLab/lbm/stepper.py

@@ -24,6 +24,8 @@ class LBMStepper:
         if self._esopull:
             self.step_fn = module.get_function("EsoPullStep")
             self.init_fn = module.get_function("InitEsoPull")
+            self.prepare_boundaries_fn = module.get_function("PrepareEsoPullBoundaries")
+            self.boundary_values_gpu = cuda.mem_alloc((6 * cfg.nx + 2 * cfg.ny * cfg.nq) * 4)
         else:
             self.step_fn = module.get_function("OneStep")
             self.init_fn = module.get_function("InitTubeFlow_v2")
@@ -41,6 +43,7 @@ class LBMStepper:
         else:
             self.grid = ((cfg.nx + tpb - 1) // tpb, cfg.ny, cfg.nz)
 
+        self.boundary_grid = ((max(cfg.nx, cfg.ny) + tpb - 1) // tpb, 2, 1)
         self._step_count = 0
 
     def initialize(self):
@@ -54,7 +57,9 @@ class LBMStepper:
             _gpu_bytes = f.n * f.nq * f.store_bytes
             cuda.memcpy_dtod(f.temp_gpu, f.ddf_gpu, _gpu_bytes)
         cuda.memcpy_dtoh(f.flag, f.flag_gpu)
-        f.download_ddf(step_id=self._step_count, force=True)
+        # Initialization publishes a new t=0 state, including EsoPull parity.
+        self._step_count = 0
+        f.download_ddf(step_id=0, force=True)
 
     def step(
         self,
@@ -71,12 +76,17 @@ class LBMStepper:
             launch_kw["stream"] = stream
 
         for _ in range(n):
+            if self._esopull:
+                self.prepare_boundaries_fn(
+                    f.ddf_gpu, self.boundary_values_gpu, np.uint64(self._step_count),
+                    block=self.block, grid=self.boundary_grid, **launch_kw,
+                )
             self._launch_curved(action_gpu, obs_gpu, **launch_kw)
 
             if self._esopull:
                 self.step_fn(
                     f.ddf_gpu, f.flag_gpu,
-                    action_gpu, obs_gpu,
+                    action_gpu, obs_gpu, self.boundary_values_gpu,
                     np.uint64(self._step_count),
                     block=self.block, grid=self.grid,
                     **launch_kw,
@@ -137,12 +147,17 @@ class LBMStepper:
             )
         else:
             assert cl.fluid_idx_gpu is not None
+            if cl.surface_velocity is not None and (
+                    cl.surface_velocity_gpu is None or cl.link_exchange_gpu is None):
+                raise RuntimeError("Surface boundary device buffers are missing")
             self.curved_fn(
                 f.ddf_gpu,
                 cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu,
                 cl.rx_gpu, cl.ry_gpu, cl.rz_gpu, cl.fallback_class_gpu, cl.body_id_gpu,
                 action_gpu, obs_gpu,
                 np.uint32(f.n_curved),
+                cl.surface_velocity_gpu if cl.surface_velocity_gpu is not None else np.uintp(0),
+                cl.link_exchange_gpu if cl.link_exchange_gpu is not None else np.uintp(0),
                 block=(tpb, 1, 1), grid=grid_c,
                 **launch_kw,
             )

src/CelerisLab/simulation.py

@@ -18,12 +18,13 @@ from __future__ import annotations
 from typing import Dict, Optional, Tuple, Any
 
 import os
+import operator
 
 import numpy as np
 import pycuda.driver as cuda
 
 from .config import LBMConfig, BodyConfig, load_lbm_config, load_body_config
-from .cuda.context import CudaContext
+from .cuda.context import CudaContext, context_bound
 from .cuda import compiler_v2 as compiler
 from .lbm.field import LBMField
 from .lbm.stepper import LBMStepper
@@ -45,39 +46,50 @@ class Simulation:
                  device_id: int = 0):
         # Internal stream for async DTOH coordination
         self._stream: Optional[cuda.Stream] = None
+        self._build_dir = None
+        self._closed = False
 
         # Load configs
         self.lbm_cfg = load_lbm_config(lbm_config_path)
+        self.lbm_cfg.validate_runtime()
         self.body_cfg = load_body_config(body_config_path)
 
         # CUDA context
         self.ctx = CudaContext(device_id)
-        arch = self._resolve_compile_arch()
+        try:
+            with self.ctx.activate():
+                arch = self._resolve_compile_arch()
+                self._build_dir = compiler.generate_config(self.lbm_cfg, n_objects=0)
+                self._ptx_path = compiler.compile_kernel(arch=arch, build_dir=self._build_dir)
+                self._module = compiler.load_module(self._ptx_path)
+                self._compiled_n_objects = 0
+                self._allocation_signature = self.lbm_cfg.allocation_signature()
+                self._compiled_viscosity = self.lbm_cfg.viscosity
 
-        # Compile kernel
+                self.field = LBMField(self.lbm_cfg, self._module)
-        compiler.generate_config(self.lbm_cfg, n_objects=0)
+                self.field.ctx = self.ctx
-        self._compiled_n_objects: int = 0
+                self.stepper = LBMStepper(self.field, self._module, self.lbm_cfg)
-        self._ptx_path = compiler.compile_kernel(arch=arch)
+                self.bodies = ObjectManager(
-        self._module = compiler.load_module(self._ptx_path)
+                    self.lbm_cfg.nx, self.lbm_cfg.ny, self.lbm_cfg.nz,
-
+                    self.lbm_cfg.nq, self.lbm_cfg,
-        # LBM field & stepper
+                )
-        self.field = LBMField(self.lbm_cfg, self._module)
+                self._initialized = False
-        self.stepper = LBMStepper(self.field, self._module, self.lbm_cfg)
+                self._load_objects_from_body_config()
-
+                if self.bodies.count > 0:
-        # Body manager
+                    self.recompile()
-        self.bodies = ObjectManager(
+                self._assert_object_count_contract(expected_count=self.bodies.count)
-            self.lbm_cfg.nx, self.lbm_cfg.ny, self.lbm_cfg.nz,
+        except Exception:
-            self.lbm_cfg.nq, self.lbm_cfg,
+            try:
-        )
+                self.close()
-
+            except Exception:
-        self._initialized = False
+                # Construction must report its original failure. Cleanup is
-        self._load_objects_from_body_config()
+                # best effort here; an explicit close failure remains retryable.
-        if self.bodies.count > 0:
+                pass
-            self.recompile()
+            raise
-        self._assert_object_count_contract(expected_count=self.bodies.count)
 
     # -- Stream access -------------------------------------------------------
     @property
+    @context_bound
     def stream(self) -> cuda.Stream:
         """Internal CUDA stream for async DTOH coordination."""
         if self._stream is None:
@@ -86,7 +98,7 @@ class Simulation:
 
     # -- Object management ---------------------------------------------------
     def add_body(self, type: str = "circle", center=None,
-                 radius: float = 0.0, **kwargs) -> int:
+                 radius: float = 0.0, *, surface_actuation=None, **kwargs) -> int:
         """Add a simulation body. Returns body_id.
 
         Before ``initialize()`` the object is added to the formal registry
@@ -98,12 +110,20 @@ class Simulation:
                 (future: ``"polygon"``, ``"mesh"``).
             center: (x, y) center coordinates (lattice units).
             radius: body radius (lattice units).
+            surface_actuation: Fixed circle aperture dict/SurfaceActuation:
+                angle, half_width (radians), normal_velocity (positive outward).
+                D2Q9/FP32/unshifted/DB only; set before initialization or stage
+                a new body. No runtime motion or amplitude setter.
             **kwargs: passed to the geometry constructor.
 
         Returns:
             int body_id (formal id before initialize; -1 if staged).
         """
+        if self.lbm_cfg.y_wall_bc == "on_site_velocity":
+            raise ValueError("on_site_velocity currently requires an object-free domain")
         type_lower = str(type).strip().lower()
+        if surface_actuation is not None and type_lower != "circle":
+            raise ValueError("Surface actuation is only supported on circles")
 
         if type_lower == "circle":
             from .body.geometry.circle import CircleGeometry
@@ -111,7 +131,7 @@ class Simulation:
             geom = CircleGeometry(cx, cy, radius)
             obj = SimObject(obj_id=-1, geometry=geom,
                             center=center, radius=radius,
-                            is_sensor=False)
+                            is_sensor=False, surface_actuation=surface_actuation)
         elif type_lower == "sensor":
             from .body.geometry.circle import CircleGeometry
             cx, cy = float(center[0]), float(center[1])
@@ -142,7 +162,8 @@ class Simulation:
 
         Before ``initialize()`` the object is removed from the formal registry
         immediately.  After ``initialize()`` the removal is staged and will be
-        committed at the next ``sync_bodies()``.
+        committed at the next ``sync_bodies()``. Removal compacts numeric IDs;
+        obtain the remaining IDs again after immediate removal or runtime sync.
 
         Args:
             id: body_id from ``add_body()``.
@@ -174,6 +195,8 @@ class Simulation:
             center = tuple(spec.get("center", ()))
             radius = float(spec.get("radius", 0.0))
             omega = spec.get("omega", None)
+            if spec.get("surface_actuation") is not None and obj_type != "cylinder":
+                raise ValueError("Surface actuation is only supported on cylinders in body config")
 
             if obj_type == "":
                 raise ValueError(
@@ -185,9 +208,10 @@ class Simulation:
                 )
 
             if obj_type == "cylinder":
-                body_id = self.add_body("circle", center=center, radius=radius)
+                body_id = self.add_body("circle", center=center, radius=radius,
+                                        surface_actuation=spec.get("surface_actuation"))
                 if omega is not None:
-                    self.bodies.get(body_id).state.omega = float(omega)
+                    self.bodies.set_body_state(body_id, omega=float(omega))
             elif obj_type == "sensor":
                 self.add_body("sensor", center=center, radius=radius)
             else:
@@ -210,6 +234,10 @@ class Simulation:
             vy: translation velocity y (future).  ``None`` = no change.
         """
         obj = self.bodies.get(id)
+        proposed = dict(vars(obj.state))
+        proposed.update({k: float(v) for k, v in (("omega", omega), ("vx", vx), ("vy", vy))
+                         if v is not None})
+        obj.validate_surface_actuation(self.lbm_cfg, proposed)
         changed = False
         if omega is not None:
             obj.state.omega = float(omega)
@@ -276,6 +304,19 @@ class Simulation:
             raise RuntimeError("No obs buffer. Call run() first.")
         return self.bodies.read_sensor(id, normalize=normalize)
 
+    @context_bound
+    def read_surface_exchange(self, id: int):
+        """Latest step's wall/aperture (Fx,Fy,Tz) exchange, not total nozzle force.
+
+        Unlike read_body this is not accumulated/time-averaged. It resets after
+        geometry sync or checkpoint load; advance before interpreting it.
+        """
+        if self.bodies.get(id).surface_actuation is None:
+            raise ValueError("Body has no surface actuation")
+        self.stream.synchronize()
+        return self.field.curved.surface_exchange(id)
+
+    @context_bound
     def read_body(self, id: int, *,
                   stream: cuda.Stream | None = None,
                   normalize: bool = True):
@@ -301,6 +342,7 @@ class Simulation:
             stream.synchronize()
         return self.bodies.read_body(id, normalize=normalize)
 
+    @context_bound
     def read_bodies(self, *, stream: cuda.Stream | None = None) -> np.ndarray:
         """Return all bodies' telemetry as a flat float32 array.
 
@@ -320,31 +362,39 @@ class Simulation:
         return self.bodies._obs_array()
 
     # -- Compilation ---------------------------------------------------------
+    @context_bound
     def _recompile(self, expected_n_objects: int | None = None) -> None:
         """Re-generate config headers and recompile the kernel.
 
         If *expected_n_objects* is ``None``, uses ``self.bodies.count``.
         Preserves ``step_count`` and GPU DDF/flag buffers across the
-        module reload.  Does NOT preserve ``action_gpu`` / ``obs_gpu``
+        module reload. Packed action/observation buffers are managed separately
-        (those are recreated by ``sync_to_gpu`` or ``run()``).
+        by the topology synchronization caller.
 
         This is safe to call before or after ``initialize()``.
         """
         if expected_n_objects is None:
             expected_n_objects = self.bodies.count
+        self.lbm_cfg.validate_runtime()
+        self._assert_allocation_contract()
         arch = self._resolve_compile_arch()
-        compiler.generate_config(self.lbm_cfg, n_objects=expected_n_objects)
+        next_dir = compiler.generate_config(self.lbm_cfg, n_objects=expected_n_objects)
+        try:
+            next_ptx = compiler.compile_kernel(arch=arch, build_dir=next_dir)
+            next_module = compiler.load_module(next_ptx)
+        except Exception:
+            compiler.cleanup_build_dir(next_dir)
+            raise
+        old_dir = self._build_dir
+        self._build_dir, self._ptx_path, self._module = next_dir, next_ptx, next_module
         self._compiled_n_objects = int(expected_n_objects)
-        # Remove stale PTX to prevent PyCUDA module handle conflicts.
-        import os as _os
-        if _os.path.exists(self._ptx_path):
-            _os.remove(self._ptx_path)
-        self._ptx_path = compiler.compile_kernel(arch=arch)
-        self._module = compiler.load_module(self._ptx_path)
 
         # Reconnect field and stepper to new module
         self.field.module = self._module
         self.field.invalidate_module_constant_cache()
+        if self.lbm_cfg.viscosity != self._compiled_viscosity:
+            self.field._runtime_params["omega"] = self.lbm_cfg.omega
+        self._compiled_viscosity = self.lbm_cfg.viscosity
         self.field._upload_params()
 
         _prev_step_count = self.stepper._step_count
@@ -354,6 +404,7 @@ class Simulation:
         self.stepper._step_count = _prev_step_count
 
         self._assert_object_count_contract(expected_count=expected_n_objects)
+        compiler.cleanup_build_dir(old_dir)
 
     def recompile(self):
         """Re-generate config headers and recompile kernel.
@@ -369,11 +420,12 @@ class Simulation:
         self._recompile()
 
     # -- Initialization ------------------------------------------------------
+    @context_bound
     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._assert_allocation_contract()
-           self._read_generated_n_objects() != self.bodies.count:
+        if self._read_generated_n_objects() != self.bodies.count:
             self.recompile()
 
         self.field.flag = self.field.build_channel_flags()
@@ -390,6 +442,7 @@ class Simulation:
         self._initialized = True
 
     # -- Runtime body topology sync -------------------------------------------
+    @context_bound
     def sync_bodies(self) -> None:
         """Apply pending body edits (add/remove) to a running simulation.
 
@@ -400,6 +453,7 @@ class Simulation:
         if not self._initialized:
             raise RuntimeError("Call initialize() before sync_bodies()")
 
+        self._assert_allocation_contract()
         # 1. Synchronize GPU
         self.stream.synchronize()
 
@@ -470,6 +524,7 @@ class Simulation:
         self._assert_runtime_contracts()
 
     # -- Stepping ------------------------------------------------------------
+    @context_bound
     def run(self, steps: int, *,
             stream: cuda.Stream | None = None,
             upload_act: bool = True,
@@ -479,18 +534,24 @@ class Simulation:
         """Advance simulation by *steps* time steps.
 
         Args:
-            steps: Number of LBM steps.
+            steps: Nonnegative integer number of LBM steps.
-            stream: CUDA stream for async operations.  ``None`` = internal stream.
+            stream: CUDA stream for async operations. ``None`` = internal stream.
+                External streams must be created under ``with sim.ctx.activate():``.
             upload_act: If True, upload host action array to ``action_gpu``
                 before the step group.
             sync_obs: If True, download ``obs_gpu`` to host pinned buffer
                 after the step group.
             zero_obs: If True (default), zero all obs segments (force, torque,
                 sensor) on GPU and reset the step accumulator before stepping.
-            checkpoint_interval: If >0, save checkpoint every N steps.
+            checkpoint_interval: Nonnegative integer; if >0, save every N steps.
         """
         if not self._initialized:
             raise RuntimeError("Call initialize() first")
+        steps = operator.index(steps)
+        checkpoint_interval = operator.index(checkpoint_interval)
+        if steps < 0 or checkpoint_interval < 0:
+            raise ValueError("steps and checkpoint_interval must be nonnegative")
+        self._assert_runtime_contracts()
         # Discard any uncommitted body edits before stepping.
         if self.bodies.has_pending_edit():
             self.bodies.clear_pending_edits()
@@ -507,7 +568,6 @@ class Simulation:
         if zero_obs:
             self.bodies.zero_obs_async(stream)
 
-        self._assert_runtime_contracts()
         if checkpoint_interval > 0:
             done = 0
             while done < steps:
@@ -518,8 +578,10 @@ class Simulation:
                     obs_gpu=self.bodies.obs_gpu,
                     stream=stream,
                 )
+                self.bodies.increment_obs_steps(batch)
                 done += batch
                 if done < steps or done % checkpoint_interval == 0:
+                    stream.synchronize()
                     self.save_checkpoint()
         else:
             self.stepper.step(
@@ -529,9 +591,9 @@ class Simulation:
                 stream=stream,
             )
 
-        # Accumulate step count for obs time-normalisation
+            # Count each batch before any checkpoint captures observations.
-        if steps > 0:
+            if steps > 0:
-            self.bodies.increment_obs_steps(steps)
+                self.bodies.increment_obs_steps(steps)
 
         # Async download obs
         if sync_obs:
@@ -544,10 +606,12 @@ class Simulation:
         self.run(n)
 
     # -- Data access ---------------------------------------------------------
+    @context_bound
     def get_macroscopic(self) -> Dict[str, np.ndarray]:
         """Download DDF and return rho, ux, uy [, uz]."""
         return self.field.get_macroscopic(step_id=self.stepper.step_count)
 
+    @context_bound
     def get_ddf(self) -> np.ndarray:
         self.field.download_ddf(step_id=self.stepper.step_count)
         return self.field.ddf.copy()
@@ -556,6 +620,7 @@ class Simulation:
         return self.field.flag.copy()
 
     # -- Runtime parameter updates -------------------------------------------
+    @context_bound
     def update_runtime_params(self, **kwargs):
         """Update __constant__ d_params without recompiling.
 
@@ -564,29 +629,54 @@ class Simulation:
         self.field.update_params(**kwargs)
 
     # -- Snapshots -----------------------------------------------------------
+    @context_bound
     def snapshot(self):
-        self.field.snapshot()
+        """Snapshot the flow field; restore keeps the current clock and controls."""
+        self.field.snapshot(step_id=self.stepper.step_count)
 
+    @context_bound
     def restore(self):
-        self.field.restore()
+        """Restore the saved flow field at the current streaming parity."""
+        self.field.restore(step_id=self.stepper.step_count)
 
     # -- Checkpoint / Restore (HDF5) ----------------------------------------
+    @context_bound
     def save_checkpoint(self, path: str = None):
         """Save full simulation state to an HDF5 file."""
         from .common.checkpoint import save_checkpoint
         return save_checkpoint(self.field, self.stepper, self.lbm_cfg,
                                self.bodies, path)
 
+    @context_bound
     def load_checkpoint(self, path: str):
         """Restore simulation state from an HDF5 checkpoint."""
         from .common.checkpoint import load_checkpoint
+        if not self._initialized:
+            self.initialize()
         load_checkpoint(path, self.field, self.stepper, self.lbm_cfg,
                         self.bodies)
         self._initialized = True
 
     # -- Cleanup -------------------------------------------------------------
     def close(self):
-        self.ctx.close()
+        if getattr(self, "_closed", False):
+            return
+        ctx = getattr(self, "ctx", None)
+        if ctx is not None and ctx.is_active:
+            raise RuntimeError("Cannot close Simulation during an active CUDA operation")
+        if ctx is not None and not ctx.is_closed:
+            with ctx.activate():
+                cuda.Context.synchronize()
+                # Release owned GPU resources while their context is current.
+                self.stepper = None
+                self.bodies = None
+                self.field = None
+                self._module = None
+                self._stream = None
+            ctx.close()
+        compiler.cleanup_build_dir(getattr(self, "_build_dir", None))
+        self._build_dir = None
+        self._closed = True
 
     def __enter__(self):
         return self
@@ -658,9 +748,16 @@ class Simulation:
 
     def _assert_runtime_contracts(self) -> None:
         """Validate object-count and obs layout contracts before stepping."""
+        self._assert_allocation_contract()
         self._assert_object_count_contract()
         self._assert_obs_layout_contract()
 
+    def _assert_allocation_contract(self) -> None:
+        if self.lbm_cfg.allocation_signature() != self._allocation_signature:
+            raise ValueError(
+                "Grid, lattice, streaming or storage configuration changed after "
+                "allocation. Construct a new Simulation for structural changes.")
+
 
 # -- Module-level helpers -------------------------------------------------------
 
@@ -668,9 +765,9 @@ def _build_patch_set(added_solid_mask, released_fluid_mask, curved_fluid_idx,
                      nx=None, ny=None):
     """Return sorted uint32 array of cell indices that need DDF patching.
 
-    Includes a 1-cell halo because ``store_physical_node_esopull`` writes
+    Includes a 1-cell halo for the old-fluid donors used by released-cell
-    to neighbor slots (D2Q9 directions 1,3,5,7).  Neighbors must be
+    reconstruction. EsoPull same-time encoding itself has disjoint ownership
-    decoded/encoded together to avoid corrupting the backing layout.
+    and does not require a halo to avoid storage collisions.
     """
     s = set()
     s.update(np.where(added_solid_mask)[0].tolist())
@@ -699,12 +796,14 @@ def _build_patch_set(added_solid_mask, released_fluid_mask, curved_fluid_idx,
 def _scatter_to_field(field, scratch_host, patch_indices):
     """Copy decoded physical DDF from scratch buffer into field.ddf."""
     nq = field.nq
+    populations = field.ddf.reshape(nq, -1)
     for j, idx in enumerate(patch_indices):
-        field.ddf[idx * nq:(idx + 1) * nq] = scratch_host[j * nq:(j + 1) * nq]
+        populations[:, idx] = scratch_host[j * nq:(j + 1) * nq]
 
 
 def _gather_from_field(field, scratch_host, patch_indices):
     """Copy patched physical DDF from field.ddf back to scratch buffer."""
     nq = field.nq
+    populations = field.ddf.reshape(nq, -1)
     for j, idx in enumerate(patch_indices):
-        scratch_host[j * nq:(j + 1) * nq] = field.ddf[idx * nq:(idx + 1) * nq]
+        scratch_host[j * nq:(j + 1) * nq] = populations[:, idx]