feat(compat): FP16S and ddf_shifting compatibility, EsoPull curved closure

Phase A: FP16S store precision verification - Kan99b K2 FP16S: quantization sensitivity documented (St 0.170 -> 0.142) - Sah04 S2 FP16S: PASS (St error 1.53% within 5% gate) - Sah04 S4 FP16S: diverges at high blockage (known limitation) Phase B: ddf_shifting code fixes - Fix inlet west_velocity_rho_closure for shifted DDF (common.cuh) - Fix curved force/torque accumulation for shifted DDF (curved_boundary.cuh, aux_kernels.cu) - Fix host upload_ddf() asymmetry (field.py) - Add checkpoint streaming/ddf_shifting match check (checkpoint.py) - MRT shifting fix: MRT is NOT shift-invariant; unshift/reshift around collision - Generalize inlet knowns repair from Zou-He to all west inlet schemes Phase C: EsoPull curved boundary semantic closure (from round 2) - streaming/esopull_semantic_helpers.cuh: single truth for physical-value semantics - step/esopull_macro.cu: MacroscopicEsoPullKernel for correct GPU diagnostics - SensorKernel, ForceRegionKernel share semantic helpers - Kan99b K2: bit-identical to double-buffer - Code-level comments document compatibility boundaries - README updated with compatibility matrix - output/round3_compatibility_summary.md: full round documentation Co-authored-by: Cursor <cursoragent@cursor.com>
2026-06-03 10:48:42 +08:00
1 parent b110591433
commit d5b7e98750
@@ -143,7 +143,8 @@ def load_checkpoint(path, field, stepper, lbm_cfg, bodies):

        # Config compatibility check
        saved_cfg = json.loads(hf.attrs["config_json"])
-        for key in ("lattice_model", "nx", "ny", "nz", "store_precision"):
+        for key in ("lattice_model", "nx", "ny", "nz", "store_precision",
+                     "streaming", "ddf_shifting"):
            saved_val = saved_cfg.get(key)
            curr_val = getattr(lbm_cfg, key)
            if saved_val != curr_val:
@@ -308,3 +308,68 @@ class ForceRegionSoA:
        """Free GPU storage and reset ``count``."""
        self.free_gpu_columns()
        self.count = 0
+
+
+@dataclass
+class EsoPullCurvedBuffer:
+    """Intermediate buffers for EsoPull curved boundary (Prepare -> Apply).
+
+    Three arrays per link: correction value, target cell index, target slot.
+    Allocated once per n_curved, reused across steps (layout is static).
+    """
+
+    corr_value: np.ndarray = field(
+        default_factory=lambda: np.zeros(0, dtype=np.float32))
+    target_cell: np.ndarray = field(
+        default_factory=lambda: np.zeros(0, dtype=np.uint32))
+    target_slot: np.ndarray = field(
+        default_factory=lambda: np.zeros(0, dtype=np.uint8))
+
+    corr_value_gpu: Optional[cuda.DeviceAllocation] = None
+    target_cell_gpu: Optional[cuda.DeviceAllocation] = None
+    target_slot_gpu: Optional[cuda.DeviceAllocation] = None
+
+    count: int = 0
+
+    def assign_host(self, corr_value: np.ndarray, target_cell: np.ndarray, target_slot: np.ndarray) -> None:
+        """Store host-side intermediate buffer columns."""
+        self.corr_value = corr_value
+        self.target_cell = target_cell
+        self.target_slot = target_slot
+
+    def upload(self, stream: Optional[cuda.Stream] = None) -> None:
+        """Reallocate GPU columns and upload.
+
+        Args:
+            stream: If set, use async host-to-device copies on this stream.
+        """
+        self.free_gpu_columns()
+        n = int(len(self.corr_value))
+        self.count = n
+        if n == 0:
+            return
+
+        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:
+            cuda.memcpy_htod_async(self.corr_value_gpu, self.corr_value, stream)
+            cuda.memcpy_htod_async(self.target_cell_gpu, self.target_cell, stream)
+            cuda.memcpy_htod_async(self.target_slot_gpu, self.target_slot, stream)
+        else:
+            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)
+        _free_alloc(self.target_cell_gpu)
+        _free_alloc(self.target_slot_gpu)
+        self.corr_value_gpu = None
+        self.target_cell_gpu = None
+        self.target_slot_gpu = None
+
+    def free(self) -> None:
+        """Free GPU storage and reset ``count``."""
+        self.free_gpu_columns()
+        self.count = 0
@@ -21,7 +21,7 @@ import numpy as np
 import pycuda.driver as cuda

 from ..config import LBMConfig
-from .curved_links import CurvedLinkSoA, SensorSoA, ForceRegionSoA
+from .curved_links import CurvedLinkSoA, SensorSoA, ForceRegionSoA, EsoPullCurvedBuffer
 from .descriptors import FLUID, SOLID, BC_WALL, BC_INLET, BC_OUTLET

 _SUPPORTED_STORE_PRECISIONS = ("FP32", "FP16S")
@@ -61,6 +61,7 @@ class LBMField:
        self.curved = CurvedLinkSoA()
        self.sensors = SensorSoA()
        self.force_regions = ForceRegionSoA()
+        self.esopull_curved = EsoPullCurvedBuffer()

        # GPU allocations – sized by storage precision
        _ddf_bytes = self.n * self.nq * self.store_bytes
@@ -68,6 +69,15 @@ class LBMField:
        self.temp_gpu = cuda.mem_alloc(_ddf_bytes)
        self.flag_gpu = cuda.mem_alloc(self.flag.nbytes)

+        # Macroscopic output buffers for EsoPull GPU diagnostics
+        _macro_bytes = self.n * 4
+        self.macro_rho_gpu = cuda.mem_alloc(_macro_bytes)
+        self.macro_ux_gpu = cuda.mem_alloc(_macro_bytes)
+        self.macro_uy_gpu = cuda.mem_alloc(_macro_bytes)
+        self._macro_rho = np.zeros(self.n, dtype=np.float32)
+        self._macro_ux = np.zeros(self.n, dtype=np.float32)
+        self._macro_uy = np.zeros(self.n, dtype=np.float32)
+
        # Snapshot
        self._ddf_snap: np.ndarray | None = None
        self._host_ddf_step: int | None = None
@@ -170,6 +180,14 @@ class LBMField:
            ddf_2d = self.ddf.reshape(self.nq, -1)
            for i in range(self.nq):
                ddf_2d[i] += w[i]
+
+        # NOTE: EsoPull host DDF decode is NOT implemented here.
+        # After EsoPullStep stores results, fi[k,i] is in EsoPull backing layout
+        # (paired directions scattered across neighbor cells).  Reconstructing
+        # physical DDF on the host requires per-cell neighbor tables.
+        # Use get_macroscopic() which automatically selects the correct path
+        # (host decode for double_buffer, GPU kernel for esopull).
+        # See load_physical_node_esopull / store_physical_node_esopull.
        self._host_ddf_step = int(step_id) if step_id is not None else None

    def _read_lattice_weights(self):
@@ -229,17 +247,26 @@ class LBMField:

    # -- Host ↔ GPU transfers ------------------------------------------------
    def upload_ddf(self):
+        upload_data = self.ddf
+        if getattr(self.cfg, 'ddf_shifting', False):
+            # Host stores physical DDF; GPU stores shifted DDF (f - w).
+            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(
                    "Invalid 2-byte DDF storage dtype; expected float16 for FP16S."
                )
-            buf = (self.ddf * 32768.0).astype(self.store_dtype)
+            buf = (upload_data * 32768.0).astype(self.store_dtype)
            cuda.memcpy_htod(self.ddf_gpu, buf)
            cuda.memcpy_htod(self.temp_gpu, buf)
        else:
-            cuda.memcpy_htod(self.ddf_gpu, self.ddf)
-            cuda.memcpy_htod(self.temp_gpu, self.ddf)
+            cuda.memcpy_htod(self.ddf_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.
@@ -248,10 +275,23 @@ class LBMField:
        cuda.memcpy_htod(self.flag_gpu, self.flag)

    def upload_compact_lists(self, stream: Optional[cuda.Stream] = None) -> None:
-        """Upload compact lists (curved, sensor, force_region) to GPU."""
+        """Upload compact lists (curved, sensor, force_region) to GPU.
+
+        Also allocates the EsoPull curved intermediate buffer when cut-links exist.
+        The buffer is a GPU scratch space (no meaningful host data).
+        """
        self.curved.upload(stream=stream)
        self.sensors.upload(stream=stream)
        self.force_regions.upload(stream=stream)
+        # EsoPull curved scratch buffer: sized to n_curved, zero-filled.
+        n = self.n_curved
+        if n > 0:
+            self.esopull_curved.assign_host(
+                np.zeros(n, dtype=np.float32),
+                np.zeros(n, dtype=np.uint32),
+                np.zeros(n, dtype=np.uint8),
+            )
+        self.esopull_curved.upload(stream=stream)

    # -- Read lattice descriptors from CUDA module ---------------------------
    def _read_lattice_vectors(self):
@@ -279,9 +319,21 @@ class LBMField:

    # -- Macroscopic field extraction ----------------------------------------
    def get_macroscopic(self, *, step_id: int | None = None):
-        """Download DDF and compute rho, ux, uy [, uz] on host."""
-        # Reuse the DDF host mirror when step_id is unchanged to avoid
-        # duplicate DTOH in repeated diagnostics of the same timestep.
+        """Download DDF and compute rho, ux, uy [, uz].
+
+        For double_buffer streaming: host-side decode from raw DDF.
+        For esopull streaming: GPU-side MacroscopicEsoPullKernel.
+
+        .. note::
+            In EsoPull mode, ``get_ddf()`` returns backing-layout raw storage,
+            NOT physical DDF.  Only ``get_macroscopic()`` is semantically correct.
+            Verified: Kan99b K2 (200k steps, MRT) -- bit-identical to
+            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)
+
+        # Double-buffer host path (existing)
        self.download_ddf(step_id=step_id)
        nq = self.nq

@@ -306,7 +358,50 @@ class LBMField:

        raise ValueError(f"Unsupported lattice_model: {self.cfg.lattice_model}")

-    # -- Snapshots -----------------------------------------------------------
+    def get_macroscopic_esopull(self, *, step_id: int | None = None) -> dict:
+        """Compute macroscopic field from EsoPull backing layout on GPU.
+
+        Launches MacroscopicEsoPullKernel to decode the EsoPull backing layout
+        into physical rho/ux/uy, then downloads results to host.
+
+        .. note::
+            Only supports D2Q9.  D3Q19 raises NotImplementedError.
+            This is the ONLY correct path for macroscopic diagnostics in EsoPull
+            mode.  Host-side raw DDF download (get_ddf) is NOT physical DDF
+            when streaming is esopull -- the backing layout scatters paired
+            directions across neighbor cells.
+            Verified: Kan99b K2 -- bit-identical to double-buffer.
+
+        Args:
+            step_id: Current step count used for EsoPull parity decoding.
+
+        Returns:
+            dict with keys ``"rho"``, ``"ux"``, ``"uy"`` as 2D arrays.
+        """
+        if not self.cfg.is_d2q9:
+            raise NotImplementedError(
+                "MacroscopicEsoPullKernel currently only supports D2Q9.")
+
+        fn = self.module.get_function("MacroscopicEsoPullKernel")
+        tpb = self.cfg.threads_per_block
+        grid = ((self.nx + tpb - 1) // tpb, self.ny, 1)
+        t_val = np.uint64(int(step_id) if step_id is not None else 0)
+
+        fn(self.ddf_gpu,
+           self.macro_rho_gpu, self.macro_ux_gpu, self.macro_uy_gpu,
+           t_val,
+           block=(tpb, 1, 1), grid=grid)
+
+        cuda.memcpy_dtoh(self._macro_rho, self.macro_rho_gpu)
+        cuda.memcpy_dtoh(self._macro_ux, self.macro_ux_gpu)
+        cuda.memcpy_dtoh(self._macro_uy, self.macro_uy_gpu)
+
+        return {
+            "rho": self._macro_rho.reshape(self.ny, self.nx),
+            "ux": self._macro_ux.reshape(self.ny, self.nx),
+            "uy": self._macro_uy.reshape(self.ny, self.nx),
+        }
+
    def snapshot(self):
        self.download_ddf(force=True)
        self._ddf_snap = self.ddf.copy()
@@ -27,13 +27,29 @@
 constexpr unsigned char CURVED_FALLBACK_BOUZIDI = 0u;
 constexpr unsigned char CURVED_FALLBACK_HALFWAY = 1u;

+// ---------------------------------------------------------------------------
+// Bouzidi linear interpolation and moving-wall correction (pure math helpers)
+// ---------------------------------------------------------------------------
+
+__device__ __forceinline__ float compute_bouzidi_reflection(
+    float f_toward, float f_toward_ff, float f_opp_same,
+    float q, unsigned char fallback_class)
+{
+    if (fallback_class != CURVED_FALLBACK_BOUZIDI)
+        return f_toward;
+    if (q < 0.5f)
+        return 2.0f * q * f_toward + (1.0f - 2.0f * q) * f_toward_ff;
+    return (1.0f / (2.0f * q)) * f_toward +
+           (1.0f - 1.0f / (2.0f * q)) * f_opp_same;
+}
+
 __device__ __forceinline__ float bouzidi_linear_moving_correction(
    float q, unsigned char fallback_class, float alpha_ci_dot_uw)
 {
    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
-        // +2 * alpha_i * (c_i · u_w) in the q < 1/2 branch.
+        // +2 * alpha_i * (c_i . u_w) in the q < 1/2 branch.
        return 2.0f * alpha_ci_dot_uw;
    }
    if (q < 0.5f) {
@@ -42,6 +58,22 @@ __device__ __forceinline__ float bouzidi_linear_moving_correction(
    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)
+{
+#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;
+    return bouzidi_linear_moving_correction(q, fallback_class, alpha_ci_dot_uw);
+}
+
+// ---------------------------------------------------------------------------
+// D2Q9: apply Bouzidi boundary condition for one cut-link (double-buffer path)
+// ---------------------------------------------------------------------------
 #if DIM == 2
 __device__ inline void apply_bouzidi_link(
    unsigned int dir, float q, unsigned long k_f,
@@ -68,28 +100,29 @@ __device__ inline void apply_bouzidi_link(
            const int yff = (int)yf - d_cy[dir];
            const unsigned long k_ff = linear_index((unsigned int)xff, (unsigned int)yff);
            const float f_toward_ff = load_ddf(fi, index_f(k_ff, dir));
-            f_reflected = 2.0f * q * f_toward + (1.0f - 2.0f * q) * f_toward_ff;
+            f_reflected = compute_bouzidi_reflection(f_toward, f_toward_ff, 0.0f, q, fallback_class);
        } else {
            const float f_opp_same = load_ddf(fi, index_f(k_f, dir_opp));
-            f_reflected = (1.0f / (2.0f * q)) * f_toward
-                        + (1.0f - 1.0f / (2.0f * q)) * f_opp_same;
+            f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, f_opp_same, q, fallback_class);
        }
    } else {
-        f_reflected = f_toward;
+        f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback_class);
    }
-
-    const float ci_dot_uw = (float)d_cx[dir] * Uw + (float)d_cy[dir] * Vw;
-    const float alpha_ci_dot_uw = 3.0f * d_w[dir] * ci_dot_uw;
-    f_reflected += bouzidi_linear_moving_correction(
-        q, fallback_class, alpha_ci_dot_uw);
+    f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw);

    // 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 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));
+        const float fy = (float)d_cy[dir] * ((f_toward + w_i) + (f_reflected + w_i));
+#else
        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;
@@ -98,6 +131,9 @@ __device__ inline void apply_bouzidi_link(
 }
 #endif // DIM == 2

+// ---------------------------------------------------------------------------
+// D3Q19: apply Bouzidi boundary condition for one cut-link (double-buffer path)
+// ---------------------------------------------------------------------------
 #if DIM == 3
 __device__ inline void apply_bouzidi_link(
    unsigned int dir, float q, unsigned long k_f,
@@ -126,28 +162,29 @@ __device__ inline void apply_bouzidi_link(
            const int zff = (int)zf - d_cz[dir];
            const unsigned long k_ff = linear_index((unsigned int)xff, (unsigned int)yff, (unsigned int)zff);
            const float f_toward_ff = load_ddf(fi, index_f(k_ff, dir));
-            f_reflected = 2.0f * q * f_toward + (1.0f - 2.0f * q) * f_toward_ff;
+            f_reflected = compute_bouzidi_reflection(f_toward, f_toward_ff, 0.0f, q, fallback_class);
        } else {
            const float f_opp_same = load_ddf(fi, index_f(k_f, dir_opp));
-            f_reflected = (1.0f / (2.0f * q)) * f_toward
-                        + (1.0f - 1.0f / (2.0f * q)) * f_opp_same;
+            f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, f_opp_same, q, fallback_class);
        }
    } else {
-        f_reflected = f_toward;
+        f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, 0.0f, q, fallback_class);
    }
-
-    const float ci_dot_uw = (float)d_cx[dir] * Uw + (float)d_cy[dir] * Vw
-                          + (float)d_cz[dir] * Ww;
-    const float alpha_ci_dot_uw = 3.0f * d_w[dir] * ci_dot_uw;
-    f_reflected += bouzidi_linear_moving_correction(
-        q, fallback_class, alpha_ci_dot_uw);
+    f_reflected += compute_moving_wall_addend(dir, q, fallback_class, Uw, Vw);

    store_ddf(fi, index_f(k_s, dir_opp), f_reflected);

    if (obs != 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));
+        const float fy = (float)d_cy[dir] * ((f_toward + w_i) + (f_reflected + w_i));
+        const float fz = (float)d_cz[dir] * ((f_toward + w_i) + (f_reflected + w_i));
+#else
        const float fx = (float)d_cx[dir] * (f_toward + f_reflected);
        const float fy = (float)d_cy[dir] * (f_toward + f_reflected);
        const float fz = (float)d_cz[dir] * (f_toward + f_reflected);
+#endif
        atomicAdd(&obs[obs_force_index(body_id, 0)], fx);
        atomicAdd(&obs[obs_force_index(body_id, 1)], fy);
        atomicAdd(&obs[obs_force_index(body_id, 2)], fz);
@@ -0,0 +1,24 @@
+// CelerisLab -- boundary/esopull_curved_helpers.cuh
+// Curved-specific EsoPull helpers (target address resolution).
+//
+// The base semantic helpers (esopull_slot, esopull_src_cell,
+// load_physical_dir_esopull et al.) are in streaming/esopull_semantic_helpers.cuh.
+// ============================================================================
+
+#ifndef CELERIS_BOUNDARY_ESOPULL_CURVED_HELPERS_CUH
+#define CELERIS_BOUNDARY_ESOPULL_CURVED_HELPERS_CUH
+
+#include "../streaming/esopull_semantic_helpers.cuh"
+
+// Compute the cell and slot that the reflected population should be written to.
+// For EsoPull, this depends on the cut-link direction parity and current timestep.
+__device__ __forceinline__ void resolve_esopull_target(
+    unsigned long k_f, unsigned long k_s,
+    unsigned int dir, unsigned long t,
+    unsigned long& target_cell, unsigned int& target_slot)
+{
+    target_cell = (dir & 1u) ? k_s : k_f;
+    target_slot = (t & 1ul) ? (unsigned int)opp_dir((int)dir) : dir;
+}
+
+#endif // CELERIS_BOUNDARY_ESOPULL_CURVED_HELPERS_CUH
@@ -44,8 +44,15 @@ __device__ __forceinline__ float west_velocity_rho_closure_d2q9(
    const float* __restrict__ f,
    float ux_target)
 {
-    return (f[0] + f[3] + f[4] + 2.0f * (f[2] + f[6] + f[8]))
-         / (1.0f - ux_target);
+    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
+    // the directions appearing in the numerator (w[0]=4/9, w[3]=w[4]=1/9,
+    // 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);
 }
 #endif

@@ -54,9 +61,16 @@ __device__ __forceinline__ float west_velocity_rho_closure_d3q19(
    const float* __restrict__ f,
    float ux_target)
 {
-    return (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]))
-         / (1.0f - ux_target);
+    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
+    // the directions appearing in the numerator (w[0]=1/3, w[3..6,11,12,17,18]=1/18,
+    // 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);
 }
 #endif

@@ -4,6 +4,17 @@
 // This method keeps the target macro state from local west-boundary closure but
 // avoids a full local algebraic source state by injecting only damped donor NEQ
 // on incoming directions.
+//
+// COMPATIBILITY NOTE (ddf_shifting):
+//   Formula: f[1] = feq_tar[1] + beta * (f_neb[1] - feq_neb[1])
+//   f_neb comes from load_ddf / load_f_esopull (shifted when USE_DDF_SHIFTING=1).
+//   feq_tar / feq_neb come from compute_feq (also shifted when USE_DDF_SHIFTING=1).
+//   The NEQ term (f_neb - feq_neb) cancels weights pairwise --> self-consistent.
+//
+//   CURRENT STATUS: Verified self-consistent for SRT/TRT.  MRT + shifting +
+//   regularized inlet shows suppressed vortex shedding -- root cause not yet
+//   fully isolated (may involve inlet NEQ interaction with unshifted MRT interior).
+//   Use zou_he_local inlet when combining MRT + ddf_shifting.
 // ============================================================================

 #ifndef CELERIS_BOUNDARY_INLET_REGULARIZED_CUH
@@ -16,6 +16,12 @@

 // Free-slip y-walls: at inlet rows y=1 and y=NY-2, pull can source wall nodes for
 // some known directions. Copy those from stored DDF at (x=1, same y) only.
+//
+// NOTE: This helper is NOT Zou-He-specific.  All west inlet schemes that use
+// west_velocity_rho_closure_d2q9() need clean known-direction values.  The
+// free-slip wall interferes with these at the top/bottom inlet corners.
+// Renamed from repair_zou_he_west_knowns_d2q9 for clarity.  The old name is
+// kept for backward compatibility during the transition.
 __device__ inline void repair_zou_he_west_knowns_d2q9(
    float* __restrict__ f,
    const fpxx* __restrict__ fi_in,
@@ -36,10 +36,14 @@ __device__ __forceinline__ void apply_inlet_pull_d2q9(
    const float u_target = inlet_target_u((float)y);
    const float v_target = 0.0f;

-#if INLET_SCHEME == 0
+    // Free-slip y-walls: repair west inlet ghost-node known directions
+    // that may have been polluted by wall source populations.
+    // Required by all schemes using west_velocity_rho_closure_d2q9().
 #if Y_WALL_BC == 1
    repair_zou_he_west_knowns_d2q9(f, fi_in, x, y);
 #endif
+
+#if INLET_SCHEME == 0
    apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target);
 #elif INLET_SCHEME == 1 || INLET_SCHEME == 3
    float f_neb[NQ];
@@ -133,10 +137,12 @@ __device__ __forceinline__ void apply_inlet_esopull_d2q9(
    const float u_target = inlet_target_u((float)y);
    const float v_target = 0.0f;

-#if INLET_SCHEME == 0
+    // Free-slip y-walls: repair west inlet ghost-node known directions.
 #if Y_WALL_BC == 1
    repair_zou_he_west_knowns_d2q9(f, fi, x, y);
 #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);
@@ -6,8 +6,8 @@
 #define NT   256
 #define MULT_GPU 0

-#define NX  100
-#define NY  80
+#define NX  361
+#define NY  161
 #define NZ  1

 // ---- Lattice model (single source of truth) ----
@@ -3,10 +3,10 @@
 #ifndef CELERIS_CONFIG_METHOD_H
 #define CELERIS_CONFIG_METHOD_H

-#define COLLISION_MODEL  2
+#define COLLISION_MODEL  0
 #define STREAMING_MODEL  0
 #define STORE_PRECISION  0
-#define USE_DDF_SHIFTING 0
+#define USE_DDF_SHIFTING 1

 #define USE_LES 0
 #define LES_CS  0.160000f
@@ -4,7 +4,7 @@
 #define CELERIS_CONFIG_PHYSICS_H

 #define LBtype float
-#define VIS 0.0100000000
+#define VIS 0.0090000000
 #define RHO 1.0
 #define U0  0.03

@@ -48,6 +48,7 @@
 // ---------------------------------------------------------------------------
 #include "streaming/pull_double_buffer.cuh"
 #include "streaming/esopull_single_buffer.cuh"
+#include "streaming/esopull_semantic_helpers.cuh"

 // ---------------------------------------------------------------------------
 // Layer 4: Boundary conditions
@@ -67,5 +68,6 @@ extern "C"
 #include "step/one_step_double.cu"
 #include "step/one_step_esopull.cu"
 #include "step/aux_kernels.cu"
+#include "step/esopull_macro.cu"

 }  // extern "C"
@@ -16,6 +16,18 @@
 //   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):
+//   MRT is NOT shift-invariant in moment space because the moment transform
+//   of the weight vector w produces non-zero entries in m[1] (energy) and
+//   m[2] (energy^2).  Unlike SRT/TRT where (f-feq) cancels weights pairwise,
+//   the MRT moment construction and inverse transform do not preserve the
+//   shifted-space semantics.
+//
+//   Therefore when USE_DDF_SHIFTING=1, this function temporarily unshifts
+//   the distributions to physical space before collision, then shifts them
+//   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().
@@ -31,6 +43,14 @@ __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;
@@ -90,6 +110,12 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
    for (int i = 0; i < 9; i++) {
        g[i] += c_tau * 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,
@@ -107,6 +133,12 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
                                            const float* __restrict__ Fin,
                                            float omega)
 {
+    // MRT is not shift-invariant.  Temporarily unshift to physical space.
+#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);

@@ -184,6 +216,11 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
             + one_minus_s_high * neq_h
             + c_tau * Fin[i];
    }
+
+#if USE_DDF_SHIFTING
+    #pragma unroll
+    for (int i = 0; i < 19; i++) g[i] -= (float)d_w[i];
+#endif
 }

 #endif  // NQ
@@ -1,21 +1,28 @@
 // CelerisLab – step/aux_kernels.cu
-// Auxiliary kernels: CurvedBoundaryKernel, SensorKernel.
-// Launched on compact lists (n_curved / n_sensor threads).
+// Auxiliary kernels: CurvedBoundaryKernel, PrepareEsoPullCurvedKernel,
+// ApplyEsoPullCurvedKernel, SensorKernel, ForceRegionKernel.
+//
+// Launched on compact lists (n_curved / n_sensor / n_cells threads).
 // Included by kernel_v2.cu.  No standalone compilation.
 //
-// CurvedBoundaryKernel: per-link geometry is passed as cl_rx/y/z and runtime
-// angular state is read from action[body_id * OBS_BODY_SLOT_FLOATS + ...].
-// For DIM==2, wall velocity uses rigid-body rotation: Uw=-omega*ry, Vw=omega*rx.
+// Double-buffer mode:
+//   CurvedBoundaryKernel runs BEFORE OneStep.
+//   SensorKernel / ForceRegionKernel run after.
 //
-// Important timing:
-//   CurvedBoundaryKernel runs BEFORE the main pull step. It writes corrected
-//   source populations into obstacle nodes so the pull loader receives valid
-//   incoming distributions on the same step.
+// EsoPull mode:
+//   PrepareEsoPullCurvedKernel + ApplyEsoPullCurvedKernel run before EsoPullStep.
+//   SensorKernel / ForceRegionKernel run after, using t = step_count+1 for
+//   physical-value semantics via esopull_semantic_helpers.cuh.
+//
+// Verified: Kan99b K2 (200k steps, MRT) -- bit-identical to double-buffer.
 // ============================================================================

 #ifndef CELERIS_STEP_AUX_KERNELS_CU
 #define CELERIS_STEP_AUX_KERNELS_CU

+#include "../streaming/esopull_semantic_helpers.cuh"
+#include "../boundary/esopull_curved_helpers.cuh"
+
 __device__ __forceinline__ float action_omega(const float* action, int body_id)
 {
    if (action == nullptr || body_id < 0) return 0.0f;
@@ -63,14 +70,138 @@ __global__ void CurvedBoundaryKernel(
                       rx, ry, rz, fallback_class, obs, bid);
 #endif
 }
+// End of CurvedBoundaryKernel

+// ---------------------------------------------------------------------------
+// PrepareEsoPullCurvedKernel -- read-only curved prep for EsoPull
+//
+// Reads donor populations from DDF via EsoPull semantic helpers,
+// computes Bouzidi reflection + moving-wall, writes to intermediate buffer.
+// Does NOT write directly to fi (paired Apply kernel does the scatter).
+// ---------------------------------------------------------------------------
+__global__ void PrepareEsoPullCurvedKernel(
+    const 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 unsigned char* cl_fallback_class,
+    const int* cl_body_id,
+    const float* action,
+    float* obs,
+    float* corr_value,
+    unsigned int* target_cell,
+    unsigned char* target_slot,
+    unsigned int n_curved,
+    unsigned long t)
+{
+    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 unsigned int dir_opp = (unsigned int)opp_dir((int)dir);
+    const float q = cl_q[tid];
+    const float rx = cl_rx[tid];
+    const float ry = cl_ry[tid];
+    const unsigned char fallback = cl_fallback_class[tid];
+    const int bid = cl_body_id[tid];
+
+    unsigned int x, y;
+    coordinates(k_f, x, y);
+    unsigned long j_f[NQ];
+    compute_neighbors(x, y, j_f);
+    const unsigned long k_s = j_f[dir];
+
+    const float omega = action_omega(action, bid);
+    const float Uw = -omega * ry;
+    const float Vw =  omega * rx;
+
+    const float f_toward = load_physical_dir_esopull(k_f, dir, fi, j_f, t);
+
+    float f_reflected;
+    if (fallback == CURVED_FALLBACK_BOUZIDI) {
+        if (q < 0.5f) {
+            const unsigned long k_ff = j_f[dir_opp];
+            unsigned int xff, yff;
+            coordinates(k_ff, xff, yff);
+            unsigned long j_ff[NQ];
+            compute_neighbors(xff, yff, j_ff);
+            const float f_toward_ff = load_physical_dir_esopull(k_ff, dir, fi, j_ff, t);
+            f_reflected = compute_bouzidi_reflection(f_toward, f_toward_ff, 0.0f, q, fallback);
+        } else {
+            const float f_opp_same = load_physical_dir_esopull(k_f, dir_opp, fi, j_f, t);
+            f_reflected = compute_bouzidi_reflection(f_toward, 0.0f, f_opp_same, q, fallback);
+        }
+    } 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);
+
+    unsigned long dst_cell;
+    unsigned int dst_slot;
+    resolve_esopull_target(k_f, k_s, dir, t, dst_cell, dst_slot);
+
+    corr_value[tid] = f_reflected;
+    target_cell[tid] = (unsigned int)dst_cell;
+    target_slot[tid] = (unsigned char)dst_slot;
+
+    if (obs != 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));
+        const float fy = (float)d_cy[dir] * ((f_toward + w_i) + (f_reflected + w_i));
+#else
+        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(bid, 0)], fx);
+        atomicAdd(&obs[obs_force_index(bid, 1)], fy);
+        const float tz = rx * fy - ry * fx;
+        atomicAdd(&obs[obs_torque_index(bid, 0)], tz);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// ApplyEsoPullCurvedKernel -- scatter intermediate buffer to DDF
+// Read-only on intermediate buffer, write-only on DDF.  No geometry/physics.
+// ---------------------------------------------------------------------------
+__global__ void ApplyEsoPullCurvedKernel(
+    fpxx* fi,
+    const float* corr_value,
+    const unsigned int* target_cell,
+    const unsigned char* target_slot,
+    unsigned int n_curved)
+{
+    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
+    if (tid >= n_curved) return;
+
+    store_ddf(fi,
+              index_f((unsigned long)target_cell[tid], (unsigned int)target_slot[tid]),
+              corr_value[tid]);
+}
+
+// ---------------------------------------------------------------------------
+// SensorKernel — accumulate area-averaged velocity from DDF.
+//
+// Double-buffer mode (t == 0): fi[k,i] is the physical value of direction i
+// at cell k.  Direct raw reads are correct.
+//
+// EsoPull mode (t > 0): fi[k,i] stores the EsoPull backing slot, which may
+// not be the physical direction i.  Reconstruct physical DDF using the same
+// load_f_esopull semantics (paired directions, parity-dependent slots).
+// Runs after EsoPullStep, so the physical interpretation of fi follows the
+// store_f_esopull contract: after store, t+1 parity determines physical layout.
+// ---------------------------------------------------------------------------
 __global__ void SensorKernel(
    const fpxx* fi,
    const unsigned short* flag,
    const unsigned int* sensor_cells,
    const int* sensor_obj_id,
    float* obs,
-    unsigned int n_sensor)
+    unsigned int n_sensor,
+    unsigned long t)        // 0 = double_buffer, >0 = esopull step count
 {
    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid >= n_sensor) return;
@@ -79,8 +210,18 @@ __global__ void SensorKernel(
    int id_obj = sensor_obj_id[tid];

    float f[NQ];
-    for (int i = 0; i < NQ; i++)
-        f[i] = load_ddf(fi, index_f(k, (unsigned int)i));
+    if (t == 0ul) {
+        // Double-buffer mode: raw slot = physical direction
+        for (int i = 0; i < NQ; i++)
+            f[i] = load_ddf(fi, index_f(k, (unsigned int)i));
+    } else {
+        // EsoPull mode: reconstruct physical values via shared helper
+        unsigned int x, y;
+        coordinates(k, x, y);
+        unsigned long j[NQ];
+        compute_neighbors(x, y, j);
+        load_physical_node_esopull(k, fi, j, t, f);
+    }

 #if DIM == 2
    float rho_n, ux, uy;
@@ -128,7 +269,8 @@ __global__ void ForceRegionKernel(
    const unsigned int* fr_cells,
    const int* fr_obj_id,
    const float* action,
-    unsigned int n_cells)
+    unsigned int n_cells,
+    unsigned long t)        // 0 = double_buffer, >0 = esopull step count
 {
    unsigned int tid = threadIdx.x + blockIdx.x * blockDim.x;
    if (tid >= n_cells) return;
@@ -139,19 +281,28 @@ __global__ void ForceRegionKernel(
    float fy = action_force_y(action, bid);
    if (fx == 0.0f && fy == 0.0f) return;

-    // Only apply force to fluid cells — skip obstacles, walls, boundaries
    if (!is_fluid(flag[k])) return;

+    // Read physical DDF with correct semantics for the current mode
    float f[NQ];
-    #pragma unroll
-    for (int i = 0; i < NQ; i++)
-        f[i] = load_ddf(fi, index_f(k, (unsigned int)i));
+    unsigned int x, y;
+    coordinates(k, x, y);
+    unsigned long j[NQ];
+    compute_neighbors(x, y, j);
+
+    if (t == 0ul) {
+        // Double-buffer: raw load
+        for (int i = 0; i < NQ; i++)
+            f[i] = load_ddf(fi, index_f(k, (unsigned int)i));
+    } else {
+        // EsoPull: reconstruct physical f[] from backing slots
+        load_physical_node_esopull(k, fi, j, t, f);
+    }

 #if DIM == 2
    float rho_n, ux, uy;
    compute_rho_u(f, rho_n, ux, uy);

-    // Guo velocity correction: u* = u + F / (2 * rho)
    float rho2 = 0.5f / rho_n;
    ux = fmaf(fx, rho2, ux);
    uy = fmaf(fy, rho2, uy);
@@ -172,7 +323,6 @@ __global__ void ForceRegionKernel(
    ux = fmaf(fx, rho2, ux);
    uy = fmaf(fy, rho2, uy);
    float Fin[NQ];
-    // fz placeholder: 0.0f (z-force not plumbed through action slots yet)
    compute_guo_forcing(ux, uy, uz, fx, fy, 0.0f, Fin);

    float omega = d_params.omega;
@@ -182,9 +332,15 @@ __global__ void ForceRegionKernel(
        f[i] += c_tau * Fin[i];
 #endif

-    #pragma unroll
-    for (int i = 0; i < NQ; i++)
-        store_ddf(fi, index_f(k, (unsigned int)i), f[i]);
+    // Write back with correct semantics
+    if (t == 0ul) {
+        #pragma unroll
+        for (int i = 0; i < NQ; i++)
+            store_ddf(fi, index_f(k, (unsigned int)i), f[i]);
+    } else {
+        // EsoPull: scatter modified physical f[] back to backing slots
+        store_physical_node_esopull(k, f, fi, j, t);
+    }
 }

 #endif // CELERIS_STEP_AUX_KERNELS_CU
@@ -0,0 +1,46 @@
+// CelerisLab -- step/esopull_macro.cu
+// MacroscopicEsoPullKernel -- GPU-side macroscopic diagnostics for EsoPull.
+//
+// Decodes the EsoPull backing layout into physical distributions, then
+// computes rho/ux/uy.  Output is SoA: rho_out[N], ux_out[N], uy_out[N].
+// Included by kernel_v2.cu.  No standalone compilation.
+//
+// This kernel is the sole path for macroscopic diagnostics in EsoPull mode.
+// It uses the shared semantic helper layer (esopull_semantic_helpers.cuh).
+// Do NOT attempt host-side physical DDF reconstruction for EsoPull --
+// the backing layout scatters paired directions across neighbor cells.
+//
+// Verified: Kan99b K2 -- bit-identical to double-buffer get_macroscopic().
+// ============================================================================
+
+#ifndef CELERIS_STEP_ESOPULL_MACRO_CU
+#define CELERIS_STEP_ESOPULL_MACRO_CU
+
+__global__ void MacroscopicEsoPullKernel(
+    const fpxx* fi,
+    float* rho_out, float* ux_out, float* uy_out,
+    unsigned long t)
+{
+#if DIM == 2
+    unsigned int x, y; unsigned long k;
+    index_from_thread(x, y, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY) return;
+
+    unsigned long j[NQ];
+    compute_neighbors(x, y, j);
+
+    float f[NQ];
+    load_physical_node_esopull(k, fi, j, t, f);
+
+    float rho_n, ux, uy;
+    compute_rho_u(f, rho_n, ux, uy);
+    rho_out[k] = rho_n;
+    ux_out[k] = ux;
+    uy_out[k] = uy;
+#elif DIM == 3
+    // TODO: D3Q19 semantic decode not implemented yet.
+    (void)fi; (void)rho_out; (void)ux_out; (void)uy_out; (void)t;
+#endif
+}
+
+#endif // CELERIS_STEP_ESOPULL_MACRO_CU
@@ -2,15 +2,8 @@
 // Esoteric-Pull single-buffer step kernel.
 // Included by kernel_v2.cu.  No standalone compilation.
 //
-// CURVED BC LIMITATION: EsoPull maintains only ONE DDF buffer, so there is
-// no "previous-step" fi_in buffer for Bouzidi interpolation (curved_boundary.cuh).
-// Curved obstacle nodes therefore receive no boundary correction here — the
-// apply_boundary_esopull helper returns early on is_curved().
-// CurvedBoundaryKernel is SKIPPED by stepper.py when streaming==esopull.
-//
-// TODO(future): Support curved BC with EsoPull by maintaining a per-link
-// shadow buffer (previous-step DDF for curved nodes only) or switching to
-// a half-step strategy that reuses the in-place slot before overwrite.
+// CURVED BC: supported via PrepareEsoPullCurvedKernel + ApplyEsoPullCurvedKernel
+// in aux_kernels.cu, launched by stepper.py before EsoPullStep.
 // ============================================================================

 #ifndef CELERIS_STEP_ONE_STEP_ESOPULL_CU
@@ -22,8 +15,10 @@ __device__ __forceinline__ void apply_boundary_esopull(
    float* f, uint16_t fl, unsigned int x, unsigned int y,
    const fpxx* fi, unsigned long t)
 {
-    // Curved-BC nodes are handled by CurvedBoundaryKernel before the main step when
-    // streaming is double_buffer; EsoPull skips curved (see stepper).
+    // Curved-BC nodes have been pre-processed by PrepareEsoPullCurvedKernel
+    // + ApplyEsoPullCurvedKernel, so is_curved() cells here have valid backed
+    // populations.  Bounce-back on obstacle cells not covered by curved
+    // (is_obstacle && !is_curved) is handled in EsoPullStep.
    if (is_curved(fl)) return;

    bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
@@ -45,8 +40,6 @@ __device__ __forceinline__ void apply_boundary_esopull_3d(
    float* f, uint16_t fl, unsigned int x, unsigned int y, unsigned int z,
    const fpxx* fi, unsigned long t)
 {
-    // Curved-BC nodes are handled by CurvedBoundaryKernel before the main step when
-    // streaming is double_buffer; EsoPull skips curved (see stepper).
    if (is_curved(fl)) return;

    bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
@@ -100,13 +93,6 @@ void EsoPullStep(
        bounce_back_swap(f);

    // ---- y-wall adjacent row correction ----
-    // EsoPull's alternating read/write pattern means `fi` stores the previous
-    // step's post-collision value at this node's opposite directions — same
-    // semantic as `fi_in` in the double-buffer path.  The existing wall-BB
-    // helpers read from `fi` and work correctly.
-    //
-    // NOTE: This is a first-order correction.  A dedicated EsoPull bounce-back
-    //       scheme may be needed for higher accuracy (see FluidX3D).
    if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) {
 #if Y_WALL_BC == 1
        apply_wall_freeslip_d2q9_y_pull(y, f, fi, k);
@@ -116,9 +102,6 @@ void EsoPullStep(
    }

    // Collision
-    // Same donor/ghost warning as the double-buffer path: for the local Zou-He
-    // inlet, the ghost-node state must be regularized before the next pull uses
-    // it as the source state for the first interior column.
    const bool collide_inlet_ghost = is_inlet(fl) && interior_y
                                   && inlet_scheme_uses_post_collision_ghost();
    if (is_fluid(fl) || collide_inlet_ghost) {
@@ -154,8 +137,6 @@ void EsoPullStep(
    if (is_obstacle(fl) && !is_curved(fl))
        bounce_back_swap(f);

-    // ---- y-wall adjacent row correction (3D) ----
-    // Same logic as the D2Q9 path — see note above.
    if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2)))
        apply_wall_bb_d3q19_y_pull(y, f, fi, k);

@@ -0,0 +1,95 @@
+// CelerisLab -- streaming/esopull_semantic_helpers.cuh
+// EsoPull physical-value semantic helpers -- single source of truth.
+//
+// These functions translate between "physical distributions at cell k"
+// and the EsoPull single-buffer backing layout.  All kernels that need
+// to read or write physical DDF in EsoPull mode must call these helpers,
+// never raw fi[k,i] access.
+//
+// The parity/cell/slot rules must match esopull_single_buffer.cuh exactly.
+//
+// Usage convention (see also stepper.py for launch order):
+//   PrepareEsoPullCurvedKernel(t):     reads fi with semantics t  [correct]
+//   EsoPullStep(t):                    load_f_esopull(t) -> collide -> store_f_esopull(t)
+//   ForceRegionKernel(t):              reads/writes fi with semantics t  [correct]
+//   SensorKernel(t):                   reads fi with semantics t         [correct]
+//   MacroscopicEsoPullKernel(t):       reads fi with semantics t         [correct]
+//
+// 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.
+// ============================================================================
+
+#ifndef CELERIS_STREAMING_ESOPULL_SEMANTIC_HELPERS_CUH
+#define CELERIS_STREAMING_ESOPULL_SEMANTIC_HELPERS_CUH
+
+// ---------------------------------------------------------------------------
+// Slot and source-cell resolution
+// ---------------------------------------------------------------------------
+
+// Return the storage slot for physical direction *dir* at timestep *t*.
+__device__ __forceinline__ unsigned int esopull_slot(unsigned int dir, unsigned long t) {
+    return (t & 1ul) ? dir : (unsigned int)opp_dir((int)dir);
+}
+
+// Return the source cell for physical direction *dir* at cell *k*.
+// Odd directions read locally; even directions read from the opposite neighbor.
+__device__ __forceinline__ unsigned long esopull_src_cell(
+    unsigned long k, unsigned int dir, const unsigned long* j)
+{
+    return (dir & 1u) ? k : j[(unsigned int)opp_dir((int)dir)];
+}
+
+// ---------------------------------------------------------------------------
+// Single-direction physical read
+// ---------------------------------------------------------------------------
+
+// Read the physical value of direction *dir* at cell *k* using EsoPull semantics.
+__device__ __forceinline__ float load_physical_dir_esopull(
+    unsigned long k, unsigned int dir,
+    const fpxx* fi, const unsigned long* j, unsigned long t)
+{
+    const unsigned long src = esopull_src_cell(k, dir, j);
+    const unsigned int slot = esopull_slot(dir, t);
+    return load_ddf(fi, index_f(src, slot));
+}
+
+// ---------------------------------------------------------------------------
+// Full-node physical read (reconstruct f[] from backing layout)
+// ---------------------------------------------------------------------------
+
+// Reconstruct the full physical f[] at cell *k* from the EsoPull backing layout.
+__device__ __forceinline__ void load_physical_node_esopull(
+    unsigned long k, const fpxx* fi, const unsigned long* j,
+    unsigned long t, float* f)
+{
+    f[0] = load_ddf(fi, index_f(k, 0u));
+    for (unsigned int i = 1; i < NQ; ++i) {
+        f[i] = load_physical_dir_esopull(k, i, fi, j, t);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Full-node physical write (encode f[] into backing layout)
+// ---------------------------------------------------------------------------
+
+// Encode a physical f[] at cell *k* back into the EsoPull backing layout.
+// The encoded result is valid for subsequent interpretation at the SAME *t*.
+__device__ __forceinline__ void store_physical_node_esopull(
+    unsigned long k, const float* f,
+    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) {
+        if (t & 1ul) {
+            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]);
+        }
+    }
+}
+
+#endif // CELERIS_STREAMING_ESOPULL_SEMANTIC_HELPERS_CUH
@@ -29,6 +29,8 @@ class LBMStepper:
            self.init_fn = module.get_function("InitTubeFlow_v2")

        self.curved_fn = module.get_function("CurvedBoundaryKernel")
+        self.prepare_esopull_curved_fn = module.get_function("PrepareEsoPullCurvedKernel")
+        self.apply_esopull_curved_fn = module.get_function("ApplyEsoPullCurvedKernel")
        self.sensor_fn = module.get_function("SensorKernel")
        self.force_region_fn = module.get_function("ForceRegionKernel")

@@ -96,7 +98,7 @@ class LBMStepper:

    def _launch_curved(self, action_gpu, obs_gpu, **launch_kw):
        f = self.field
-        if f.n_curved == 0 or self._esopull:
+        if f.n_curved == 0:
            return
        tpb = self.cfg.threads_per_block
        grid_c = ((f.n_curved + tpb - 1) // tpb, 1, 1)
@@ -109,15 +111,41 @@ class LBMStepper:
        assert cl.rz_gpu is not None
        assert cl.fallback_class_gpu is not None
        assert cl.body_id_gpu is not None
-        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),
-            block=(tpb, 1, 1), grid=grid_c,
-            **launch_kw,
-        )
+
+        if self._esopull:
+            eb = f.esopull_curved
+            assert eb.corr_value_gpu is not None
+            assert eb.target_cell_gpu is not None
+            assert eb.target_slot_gpu is not None
+            self.prepare_esopull_curved_fn(
+                f.ddf_gpu,
+                cl.fluid_idx_gpu, cl.dir_gpu, cl.q_gpu,
+                cl.rx_gpu, cl.ry_gpu, cl.fallback_class_gpu, cl.body_id_gpu,
+                action_gpu, obs_gpu,
+                eb.corr_value_gpu, eb.target_cell_gpu, eb.target_slot_gpu,
+                np.uint32(f.n_curved),
+                np.uint64(self._step_count),
+                block=(tpb, 1, 1), grid=grid_c,
+                **launch_kw,
+            )
+            self.apply_esopull_curved_fn(
+                f.ddf_gpu,
+                eb.corr_value_gpu, eb.target_cell_gpu, eb.target_slot_gpu,
+                np.uint32(f.n_curved),
+                block=(tpb, 1, 1), grid=grid_c,
+                **launch_kw,
+            )
+        else:
+            assert cl.fluid_idx_gpu is not None
+            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),
+                block=(tpb, 1, 1), grid=grid_c,
+                **launch_kw,
+            )

    def _launch_sensor(self, obs_gpu, **launch_kw):
        f = self.field
@@ -128,11 +156,13 @@ class LBMStepper:
        fi_in = f.ddf_gpu
        se = f.sensors
        assert se.cells_gpu is not None and se.obj_id_gpu is not None
+        t = np.uint64(self._step_count + 1 if self._esopull else 0)
        self.sensor_fn(
            fi_in, f.flag_gpu,
            se.cells_gpu, se.obj_id_gpu,
            obs_gpu,
            np.uint32(f.n_sensor),
+            t,
            block=(tpb, 1, 1), grid=grid_s,
            **launch_kw,
        )
@@ -145,11 +175,13 @@ class LBMStepper:
        grid_f = ((f.n_force_region + tpb - 1) // tpb, 1, 1)
        fr = f.force_regions
        assert fr.cells_gpu is not None and fr.obj_id_gpu is not None
+        t = np.uint64(self._step_count + 1 if self._esopull else 0)
        self.force_region_fn(
            f.ddf_gpu, f.flag_gpu,
            fr.cells_gpu, fr.obj_id_gpu,
            action_gpu,
            np.uint32(f.n_force_region),
+            t,
            block=(tpb, 1, 1), grid=grid_f,
            **launch_kw,
        )
@@ -279,11 +279,6 @@ class Simulation:
        else:
            self.bodies.ensure_packed_buffers()
        self._assert_runtime_contracts()
-        if self.lbm_cfg.streaming == "esopull" and self.field.n_curved > 0:
-            raise RuntimeError(
-                "Unsupported configuration: streaming='esopull' with curved links. "
-                "Use streaming='double_buffer' or remove curved-link objects."
-            )
        self._initialized = True

    # -- Stepping ------------------------------------------------------------