D3Q19 SRT

.gitignore

@@ -66,3 +66,6 @@ venv.bak/
 *.tmp
 *.bak
 *.log
+
+# reference:
+ref/

src/CelerisLab/cuda/compiler.py

@@ -51,6 +51,19 @@ def compile_kernel():
         ]
     )
 
+
+def compile_kernel_v2():
+    """Compile the new modular kernel (kernel_v2.cu → kernel_v2.ptx)."""
+    subprocess.run(
+        [
+            "nvcc",
+            "-ptx",
+            kernel_path("kernel_v2.cu"),
+            "-o",
+            kernel_path("kernel_v2.ptx"),
+        ]
+    )
+
 def config_kernal(config_cuda: CudaConfig, config_field: FlowFieldConfig):
     lines = read_lines(kernel_path("macros.h"))
     lines = modify_macro(lines, "MULT_GPU", config_cuda.multi_gpu)
@@ -77,6 +90,31 @@ def config_kernal(config_cuda: CudaConfig, config_field: FlowFieldConfig):
 
     write_lines(kernel_path("macros.h"), lines)
 
+
+def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig,
+                     collision_model: int = 2,
+                     streaming_model: int = 0,
+                     store_precision: int = 0,
+                     use_ddf_shifting: int = 0):
+    """Configure macros.h for the new modular kernel architecture.
+
+    Args:
+        collision_model:  0=SRT, 1=TRT, 2=MRT (default)
+        streaming_model:  0=double-buffer (default), 1=Esoteric-Pull
+        store_precision:  0=FP32 (default), 1=FP16S, 2=FP16C
+        use_ddf_shifting: 0=off (default), 1=on
+    """
+    # First apply legacy config
+    config_kernal(config_cuda, config_field)
+
+    # Then apply new architecture macros
+    lines = read_lines(kernel_path("macros.h"))
+    lines = modify_macro(lines, "COLLISION_MODEL", collision_model)
+    lines = modify_macro(lines, "STREAMING_MODEL", streaming_model)
+    lines = modify_macro(lines, "STORE_PRECISION", store_precision)
+    lines = modify_macro(lines, "USE_DDF_SHIFTING", use_ddf_shifting)
+    write_lines(kernel_path("macros.h"), lines)
+
 def config_object(n_obj: int):
     lines = read_lines(kernel_path("macros.h"))
     lines = modify_macro(lines, "N_OBJS", n_obj)

src/CelerisLab/lbm/configs/config_gym.json

@@ -1,3 +0,0 @@
-{
-    
-}

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

@@ -0,0 +1,131 @@
+// CelerisLab – boundary/bounce_back.cuh
+// Full-way bounce-back for grid-aligned solid walls.
+//
+// For nodes adjacent to a solid wall, the incoming population from the
+// wall direction is reflected:
+//   f[i] = f[opp(i)]   (+ wall velocity correction for moving walls)
+// ============================================================================
+
+#ifndef CELERIS_BOUNDARY_BOUNCE_BACK_CUH
+#define CELERIS_BOUNDARY_BOUNCE_BACK_CUH
+
+// ---------------------------------------------------------------------------
+// Simple bounce-back:  replace f[i] with f[opp(i)] for solid neighbors
+// This version checks the flag of each neighbor.
+// ---------------------------------------------------------------------------
+__device__ inline void apply_bounce_back(float* __restrict__ f,
+                                         const unsigned long* __restrict__ j,
+                                         const uint8_t* __restrict__ flags_arr)
+{
+    // For each direction pair (i, i+1):
+    //   If the neighbor in direction i is SOLID → f[i] should bounce
+    //   from the opposite direction.
+    //
+    // With paired ordering, opp(i) = i+1 and opp(i+1) = i.
+    // In pull streaming, f[i] was loaded from j[opp(i)] = j[i+1].
+    // If j[i+1] is solid, we need bounce-back.
+
+    for (int i = 1; i < NQ; i += 2) {
+        // Direction i: check if source (j[i+1]) is solid
+        if (flags_arr[j[i + 1]] & LEGACY_SOLID) {
+            // Bounce: f[i] = f_post[opp(i)] at current node = pre-collision f[i+1]
+            // But in pull-streaming context, we already loaded from neighbor.
+            // For simple BB on a flat wall, just swap f[i] and f[i+1]:
+            float temp = f[i];
+            f[i] = f[i + 1];
+            f[i + 1] = temp;
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Top/bottom wall bounce-back (specialized for channel flow)
+// For nodes at y=1 (adjacent to y=0 wall) or y=NY-2 (adjacent to y=NY-1 wall)
+// ---------------------------------------------------------------------------
+#if NQ == 9
+__device__ inline void apply_wall_bb_d2q9(unsigned int y,
+                                           float* __restrict__ f)
+{
+    if (y == 1) {
+        float temp;
+        temp = f[3]; f[3] = f[4]; f[4] = temp;   // ±y swap
+        temp = f[5]; f[5] = f[6]; f[6] = temp;   // ±(x+y) swap
+        temp = f[7]; f[7] = f[8]; f[8] = temp;   // ±(x-y) swap
+    }
+    else if (y == NY - 2) {
+        float temp;
+        temp = f[3]; f[3] = f[4]; f[4] = temp;
+        temp = f[5]; f[5] = f[6]; f[6] = temp;
+        temp = f[7]; f[7] = f[8]; f[8] = temp;
+    }
+}
+#endif
+
+// ---------------------------------------------------------------------------
+// D3Q19 wall bounce-back on y=0/NY-1 (channel walls)
+// Pairs with non-zero cy: (3,4)±y, (7,8)±(x+y), (11,12)±(y+z),
+//                          (13,14)±(x-y), (17,18)±(y-z)
+// ---------------------------------------------------------------------------
+// ---------------------------------------------------------------------------
+// D3Q19 y-wall half-way bounce-back for pull double-buffer streaming.
+//
+// Instead of swapping f[i] <-> f[opp(i)] (which corrupts the toward-wall
+// direction with garbage from the wall node), we read the opposite direction
+// at the SAME node from the INPUT buffer.  That value is the post-collision
+// population from the previous time step — exactly the correct half-way BB.
+//
+// At y=1  (wall at y=0):  directions that arrive FROM y=0 have cy_src = -1.
+//   The pull loaded them from the wall (garbage).  Replace with fi_in[k, opp].
+//
+// At y=NY-2 (wall at y=NY-1): directions from y=NY-1, cy_src = +1.
+// ---------------------------------------------------------------------------
+#if NQ == 19
+__device__ inline void apply_wall_bb_d3q19_y_pull(unsigned int y,
+                                                   float* __restrict__ f,
+                                                   const fpxx* __restrict__ fi_in,
+                                                   unsigned long k)
+{
+    if (y == 1) {
+        // Directions whose pull-source is at y=0 (wall):
+        //   f[3]  (+y)    ← source j[4]  = (x, 0, z)
+        //   f[7]  (+x+y)  ← source j[8]  = (x-1, 0, z)
+        //   f[11] (+y+z)  ← source j[12] = (x, 0, z-1)
+        //   f[14] (-x+y)  ← source j[13] = (x+1, 0, z)
+        //   f[17] (+y-z)  ← source j[18] = (x, 0, z+1)
+        f[3]  = load_ddf(fi_in, index_f(k, 4u));
+        f[7]  = load_ddf(fi_in, index_f(k, 8u));
+        f[11] = load_ddf(fi_in, index_f(k, 12u));
+        f[14] = load_ddf(fi_in, index_f(k, 13u));
+        f[17] = load_ddf(fi_in, index_f(k, 18u));
+    }
+    else if (y == (unsigned int)(NY - 2)) {
+        // Directions whose pull-source is at y=NY-1 (wall):
+        //   f[4]  (-y)    ← source j[3]  = (x, NY-1, z)
+        //   f[8]  (-x-y)  ← source j[7]  = (x+1, NY-1, z)
+        //   f[12] (-y-z)  ← source j[11] = (x, NY-1, z+1)
+        //   f[13] (+x-y)  ← source j[14] = (x-1, NY-1, z)
+        //   f[18] (-y+z)  ← source j[17] = (x, NY-1, z-1)
+        f[4]  = load_ddf(fi_in, index_f(k, 3u));
+        f[8]  = load_ddf(fi_in, index_f(k, 7u));
+        f[12] = load_ddf(fi_in, index_f(k, 11u));
+        f[13] = load_ddf(fi_in, index_f(k, 14u));
+        f[18] = load_ddf(fi_in, index_f(k, 17u));
+    }
+}
+
+// Keep the old swap version for backward compatibility
+__device__ inline void apply_wall_bb_d3q19_y(unsigned int y,
+                                              float* __restrict__ f)
+{
+    if (y == 1 || y == NY - 2) {
+        float temp;
+        temp = f[3];  f[3]  = f[4];  f[4]  = temp;   // ±y
+        temp = f[7];  f[7]  = f[8];  f[8]  = temp;   // ±(x+y)
+        temp = f[11]; f[11] = f[12]; f[12] = temp;   // ±(y+z)
+        temp = f[13]; f[13] = f[14]; f[14] = temp;   // ±(x-y)
+        temp = f[17]; f[17] = f[18]; f[18] = temp;   // ±(y-z)
+    }
+}
+#endif
+
+#endif // CELERIS_BOUNDARY_BOUNCE_BACK_CUH

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

@@ -0,0 +1,111 @@
+// CelerisLab – boundary/curved_boundary.cuh
+// Interpolated curved-boundary bounce-back.
+//
+// For curved solid surfaces that don't align with the grid, the
+// reflected population is interpolated using the fractional distance q
+// between the fluid node and the wall intersection.
+//
+// Migrated from existing kernel.cu + ref/main.cu.
+//
+// Delta pool layout per curved-boundary node (11 floats for D2Q9):
+//   delta[offset + 0]  : encoded object id  (bitcast int → float)
+//   delta[offset + 1..NQ-1] : q values for each direction (0 if no wall hit)
+//   delta[offset + NQ]  : normal_y  (wall normal y-component / R)
+//   delta[offset + NQ+1]: normal_x  (wall normal x-component / R)
+// ============================================================================
+
+#ifndef CELERIS_BOUNDARY_CURVED_BOUNDARY_CUH
+#define CELERIS_BOUNDARY_CURVED_BOUNDARY_CUH
+
+#if NQ == 9
+
+// ---------------------------------------------------------------------------
+// apply_curved_boundary:
+//   Interpolated bounce-back with wall velocity correction.
+//   Operates on a SOLID+INTERFACE node's neighbors.
+//
+// Parameters:
+//   n       – current node (solid + interface)
+//   x, y    – coordinates of current node
+//   f_out   – the output DDF buffer (double-buffer scheme) being modified
+//   delta   – parameter pool
+//   id_off  – offset into delta for this node (= indx[n])
+//   Uw, Vw  – wall velocity at this node
+//   obs_fx, obs_fy – accumulated force observation (atomicAdd targets)
+// ---------------------------------------------------------------------------
+__device__ inline void apply_curved_boundary(
+    unsigned long n, unsigned int x, unsigned int y,
+    fpxx* __restrict__ f_out,
+    const float* __restrict__ delta,
+    int id_off,
+    float Uw, float Vw,
+    float* __restrict__ obs_fx,
+    float* __restrict__ obs_fy)
+{
+    // New paired direction ordering:
+    //   cx = {0, 1,-1, 0, 0, 1,-1, 1,-1}
+    //   cy = {0, 0, 0, 1,-1, 1,-1,-1, 1}
+
+    for (int i = 1; i < NQ; i++) {
+        int x_neb = x + d_cx[i];
+        int y_neb = y + d_cy[i];
+
+        // Check bounds (skip if neighbor is outside domain)
+        if (x_neb < 0 || x_neb >= NX || y_neb < 0 || y_neb >= NY) continue;
+
+        unsigned long k_neb = linear_index((unsigned int)x_neb, (unsigned int)y_neb);
+
+        // Only process if neighbor is FLUID
+        // (We read the flag from global memory – this is a boundary kernel,
+        //  called infrequently per node, so the extra read is acceptable.)
+        // The caller should ensure this node IS solid+interface.
+
+        float q = delta[id_off + i];  // fractional distance (0 if no wall hit)
+        if (q <= 0.0f) continue;       // no wall intersection in this direction
+
+        int oi = opp_dir(i);
+        float ci_dot_uw = (float)d_cx[i] * Uw + (float)d_cy[i] * Vw;
+        float wall_term = 6.0f * d_w[i] * ci_dot_uw;
+
+        // Second neighbor for quadratic interpolation
+        int x_neb2 = x + 2 * d_cx[i];
+        int y_neb2 = y + 2 * d_cy[i];
+        // Clamp to domain (simple, could be improved)
+        x_neb2 = max(0, min(NX - 1, x_neb2));
+        y_neb2 = max(0, min(NY - 1, y_neb2));
+        unsigned long k_neb2 = linear_index((unsigned int)x_neb2, (unsigned int)y_neb2);
+
+        // Read current DDF values from output buffer
+        float f_self_opp  = load_ddf(f_out, index_f(n,      (unsigned int)oi));
+        float f_neb_opp   = load_ddf(f_out, index_f(k_neb,  (unsigned int)oi));
+        float f_neb2_fwd  = load_ddf(f_out, index_f(k_neb2, (unsigned int)i));
+
+        // Interpolated bounce-back (Yu, Mei, Shyy, 2003)
+        float f_reflected = (q * f_self_opp + (1.0f - q) * f_neb_opp
+                             + q * f_neb2_fwd + wall_term) / (1.0f + q);
+
+        store_ddf(f_out, index_f(k_neb, (unsigned int)i), f_reflected);
+
+        // Force observation (momentum exchange)
+        float f_neb_fwd = load_ddf(f_out, index_f(k_neb, (unsigned int)i));
+        float f_self_fwd = load_ddf(f_out, index_f(n, (unsigned int)i));
+        float f_sum = f_neb_fwd + f_self_opp;
+
+        int x_back = x - d_cx[i];
+        int y_back = y - d_cy[i];
+        x_back = max(0, min(NX - 1, x_back));
+        y_back = max(0, min(NY - 1, y_back));
+
+        // Accumulate force
+        if (obs_fx != nullptr) {
+            atomicAdd(obs_fx, -f_sum * (float)d_cx[i] + f_self_fwd);
+            atomicAdd(obs_fy, -f_sum * (float)d_cy[i] +
+                      load_ddf(f_out, index_f(linear_index((unsigned int)x_back, (unsigned int)y_back),
+                                              (unsigned int)i)));
+        }
+    }
+}
+
+#endif // NQ == 9
+
+#endif // CELERIS_BOUNDARY_CURVED_BOUNDARY_CUH

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

@@ -0,0 +1,174 @@
+// CelerisLab – boundary/ibm_kernels.cuh
+// Immersed Boundary Method (IBM) kernels: Euler↔Lagrangian coupling.
+// Migrated from ref/main.cu: Eul2Lag, Lag2Eul, Update_Lagrangian.
+//
+// Uses Peskin's 4-point delta function for interpolation/spreading.
+//
+// Data layout (per rigid body):
+//   RigidBodyState2D: {x, y, theta, vx, vy, omega}
+//   Lagrangian points: motionL[NT_LAG * 4]
+//     motionL[4*k + 0] = target vx at point k
+//     motionL[4*k + 1] = target vy at point k
+//     motionL[4*k + 2] = XL (x-position)
+//     motionL[4*k + 3] = YL (y-position)
+//   forceL[NT_LAG * 2] = Lagrangian force (fx, fy) per point
+// ============================================================================
+
+#ifndef CELERIS_BOUNDARY_IBM_KERNELS_CUH
+#define CELERIS_BOUNDARY_IBM_KERNELS_CUH
+
+#include "../core/params.cuh"  // RigidBodyState2D, RigidBodyControl2D
+
+// IBM is currently 2D only. Guard all kernels.
+#if DIM == 2
+
+// ---------------------------------------------------------------------------
+// Peskin 4-point delta function  (1D component)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ float peskin_delta_1d(float r) {
+    float ar = fabsf(r);
+    if (ar < 1.0f)
+        return 0.125f * (3.0f - 2.0f * ar + sqrtf(1.0f + 4.0f * ar - 4.0f * ar * ar));
+    else if (ar < 2.0f)
+        return 0.125f * (5.0f - 2.0f * ar - sqrtf(-7.0f + 12.0f * ar - 4.0f * ar * ar));
+    else
+        return 0.0f;
+}
+
+// ---------------------------------------------------------------------------
+// Peskin 4-point delta function  (2D = product of 1D components)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ float peskin_delta_2d(float dx, float dy) {
+    return peskin_delta_1d(dx) * peskin_delta_1d(dy);
+}
+
+// ---------------------------------------------------------------------------
+// update_lagrangian_control:
+//   Update Lagrangian point positions using RigidBodyState2D.
+//   Supports circular body (points on circumference).
+//
+//   Thread k ∈ [0, NT_LAG)  — one thread per Lagrangian point.
+// ---------------------------------------------------------------------------
+__global__ void update_lagrangian_control(
+    float* __restrict__       motionL,   // [NT_LAG * 4]
+    const RigidBodyState2D*   state,     // single rigid body state
+    float                     radius,    // body radius
+    int                       NT_LAG)    // number of Lagrangian points
+{
+    int k = threadIdx.x + blockIdx.x * blockDim.x;
+    if (k >= NT_LAG) return;
+
+    float angle = state->theta + (float)k * 2.0f * 3.14159265358979323846f / (float)NT_LAG;
+
+    // Target velocity at this Lagrangian point (rigid-body kinematics)
+    float r_sin = radius * sinf(angle);
+    float r_cos = radius * cosf(angle);
+
+    motionL[4 * k + 0] = state->vx - state->omega * r_sin;  // target vx
+    motionL[4 * k + 1] = state->vy + state->omega * r_cos;  // target vy
+    motionL[4 * k + 2] = state->x + r_cos;                  // XL
+    motionL[4 * k + 3] = state->y + r_sin;                  // YL
+}
+
+// ---------------------------------------------------------------------------
+// euler_to_lagrangian (Eul2Lag):
+//   Interpolate Eulerian velocity to Lagrangian points.
+//   Compute Lagrangian force = penalty * (target_vel - interpolated_vel) * dL.
+//
+//   Thread k ∈ [0, NT_LAG)
+// ---------------------------------------------------------------------------
+__global__ void euler_to_lagrangian(
+    float* __restrict__       forceL,     // [NT_LAG * 2] output forces
+    const float* __restrict__ motionL,    // [NT_LAG * 4]
+    const float* __restrict__ u_field,    // [DIM * N] SoA velocity field
+    int                       NT_LAG)
+{
+    int k = threadIdx.x + blockIdx.x * blockDim.x;
+    if (k >= NT_LAG) return;
+
+    float XL = motionL[4 * k + 2];
+    float YL = motionL[4 * k + 3];
+    int ix = (int)XL;
+    int iy = (int)YL;
+
+    // Arc-length element dL (central difference)
+    float dL;
+    if (k > 0 && k < NT_LAG - 1) {
+        float dx1 = motionL[4*k+2] - motionL[4*(k-1)+2];
+        float dy1 = motionL[4*k+3] - motionL[4*(k-1)+3];
+        float dx2 = motionL[4*(k+1)+2] - motionL[4*k+2];
+        float dy2 = motionL[4*(k+1)+3] - motionL[4*k+3];
+        dL = 0.5f * (sqrtf(dx1*dx1 + dy1*dy1) + sqrtf(dx2*dx2 + dy2*dy2));
+    } else if (k == 0) {
+        float dx1 = motionL[4*1+2] - motionL[4*0+2];
+        float dy1 = motionL[4*1+3] - motionL[4*0+3];
+        dL = sqrtf(dx1*dx1 + dy1*dy1);
+    } else {
+        float dx1 = motionL[4*k+2] - motionL[4*(k-1)+2];
+        float dy1 = motionL[4*k+3] - motionL[4*(k-1)+3];
+        dL = sqrtf(dx1*dx1 + dy1*dy1);
+    }
+
+    // Interpolate velocity from Euler grid using Peskin 4-point delta
+    float UL = 0.0f, VL = 0.0f;
+    for (int dy = -2; dy <= 2; dy++) {
+        for (int dx = -2; dx <= 2; dx++) {
+            int ex = ix + dx;
+            int ey = iy + dy;
+            if (ex < 0 || ex >= NX || ey < 0 || ey >= NY) continue;
+
+            float ww = peskin_delta_2d((float)ex - XL, (float)ey - YL);
+            unsigned long en = (unsigned long)ey * NX + (unsigned long)ex;
+            UL += u_field[en] * ww;                       // ux
+            VL += u_field[TOTAL_CELLS + en] * ww;         // uy
+        }
+    }
+
+    // Force = (target_vel - interpolated_vel) * dL
+    forceL[2 * k + 0] = 2.0f * (motionL[4*k+0] - UL) * dL;   // fx
+    forceL[2 * k + 1] = 2.0f * (motionL[4*k+1] - VL) * dL;   // fy
+}
+
+// ---------------------------------------------------------------------------
+// lagrangian_to_euler (Lag2Eul):
+//   Spread Lagrangian forces back to the Eulerian grid.
+//
+//   Thread (x, y) — one thread per Euler node.
+//   Only processes nodes near the body (bounding box check).
+// ---------------------------------------------------------------------------
+__global__ void lagrangian_to_euler(
+    float* __restrict__       forceE,     // [DIM * N] output Euler force field
+    const float* __restrict__ forceL,     // [NT_LAG * 2]
+    const float* __restrict__ motionL,    // [NT_LAG * 4]
+    float  xc, float yc, float L0,       // body center and extent
+    int    NT_LAG)
+{
+    unsigned int x, y;
+    unsigned long k;
+#if DIM == 2
+    index_from_thread(x, y, k);
+#endif
+
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY) return;
+
+    // Bounding box check
+    if ((float)x < xc - L0 || (float)x > xc + L0 ||
+        (float)y < yc - L0 || (float)y > yc + L0) return;
+
+    float forcex = 0.0f, forcey = 0.0f;
+
+    for (int kl = 0; kl < NT_LAG; kl++) {
+        float XL = motionL[4 * kl + 2];
+        float YL = motionL[4 * kl + 3];
+        float ww = peskin_delta_2d((float)x - XL, (float)y - YL);
+        forcex += forceL[2 * kl + 0] * ww;
+        forcey += forceL[2 * kl + 1] * ww;
+    }
+
+    forceE[k]               = forcex;
+    forceE[TOTAL_CELLS + k] = forcey;
+}
+
+#endif // DIM == 2
+
+#endif // CELERIS_BOUNDARY_IBM_KERNELS_CUH

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

@@ -0,0 +1,173 @@
+// CelerisLab – boundary/inlet_outlet.cuh
+// Inlet and outlet boundary conditions  (D2Q9).
+//
+// Parabolic inlet (non-equilibrium extrapolation, Zou-He style):
+//   Left wall (x=0): reconstruct cx>0 populations (i=1,5,7)
+//
+// Pressure outlet (non-equilibrium extrapolation):
+//   Right wall (x=NX-1): reconstruct cx<0 populations (i=2,6,8)
+//
+// New paired D2Q9 ordering:
+//   cx = {0, 1,-1, 0, 0, 1,-1, 1,-1}
+//   cy = {0, 0, 0, 1,-1, 1,-1,-1, 1}
+// ============================================================================
+
+#ifndef CELERIS_BOUNDARY_INLET_OUTLET_CUH
+#define CELERIS_BOUNDARY_INLET_OUTLET_CUH
+
+#if NQ == 9
+
+// ---------------------------------------------------------------------------
+// Parabolic inlet  (x = 0, non-equilibrium extrapolation)
+//
+// f, f_neb are local DDF arrays:
+//   f      = populations at the boundary node (x=0)
+//   f_neb  = populations at the interior neighbor (x=1)
+//   y      = y-coordinate of the boundary node
+//
+// Reconstructs f[1], f[5], f[7]  (cx > 0 directions in new ordering)
+// using: f_bc[i] = f_neb[i] - feq(rho_neb, u_neb)[i] + feq(rho_neb, u_target)[i]
+// ---------------------------------------------------------------------------
+__device__ inline void apply_parabolic_inlet(float* __restrict__ f,
+                                             const float* __restrict__ f_neb,
+                                             float y_coord)
+{
+    // Neighbor macros
+    float p_neb = (f_neb[0]+f_neb[1]+f_neb[2]+f_neb[3]+f_neb[4]
+                  +f_neb[5]+f_neb[6]+f_neb[7]+f_neb[8]) / 3.0f;
+
+    // Target velocity (parabolic profile)
+    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
+    float u_target = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+    float v_target = 0.0f;
+
+    // Neighbor velocity
+    float u_neb = (f_neb[1]-f_neb[2]+f_neb[5]-f_neb[6]+f_neb[7]-f_neb[8]) / RHO;
+    float v_neb = (f_neb[3]-f_neb[4]+f_neb[5]-f_neb[6]-f_neb[7]+f_neb[8]) / RHO;
+
+    // feq for direction i=1 (cx=1, cy=0), w=1/9:
+    //   feq = (2p + RHO*(2u² + 2u - v²)) / 6
+    float feq1_target = (2.0f*p_neb + RHO*(2.0f*u_target*u_target + 2.0f*u_target - v_target*v_target)) / 6.0f;
+    float feq1_neb    = (2.0f*p_neb + RHO*(2.0f*u_neb*u_neb + 2.0f*u_neb - v_neb*v_neb)) / 6.0f;
+
+    // feq for direction i=5 (cx=1, cy=1), w=1/36:
+    //   feq = (p + RHO*(u² + 3uv + u + v² + v)) / 12
+    float feq5_target = (p_neb + RHO*(u_target*u_target + 3.0f*u_target*v_target + u_target + v_target*v_target + v_target)) / 12.0f;
+    float feq5_neb    = (p_neb + RHO*(u_neb*u_neb + 3.0f*u_neb*v_neb + u_neb + v_neb*v_neb + v_neb)) / 12.0f;
+
+    // feq for direction i=7 (cx=1, cy=-1), w=1/36:
+    //   feq = (p + RHO*(u² - 3uv + u + v² - v)) / 12
+    float feq7_target = (p_neb + RHO*(u_target*u_target - 3.0f*u_target*v_target + u_target + v_target*v_target - v_target)) / 12.0f;
+    float feq7_neb    = (p_neb + RHO*(u_neb*u_neb - 3.0f*u_neb*v_neb + u_neb + v_neb*v_neb - v_neb)) / 12.0f;
+
+    // Non-equilibrium extrapolation
+    f[1] = f_neb[1] - feq1_neb + feq1_target;
+    f[5] = f_neb[5] - feq5_neb + feq5_target;
+    f[7] = f_neb[7] - feq7_neb + feq7_target;
+}
+
+// ---------------------------------------------------------------------------
+// Pressure outlet  (x = NX-1, non-equilibrium extrapolation)
+//
+// Reconstructs f[2], f[6], f[8]  (cx < 0 directions in new ordering)
+// p_out = 0 (gauge pressure), uses velocity from neighbor.
+// ---------------------------------------------------------------------------
+__device__ inline void apply_pressure_outlet(float* __restrict__ f,
+                                             const float* __restrict__ f_neb,
+                                             float y_coord)
+{
+    float p_out = 0.0f;
+
+    // Target velocity (parabolic, same as inlet for consistency)
+    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
+    float u_target = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+    float v_target = 0.0f;
+
+    // Neighbor velocity
+    float u_neb = (f_neb[1]-f_neb[2]+f_neb[5]-f_neb[6]+f_neb[7]-f_neb[8]) / RHO;
+    float v_neb = (f_neb[3]-f_neb[4]+f_neb[5]-f_neb[6]-f_neb[7]+f_neb[8]) / RHO;
+
+    // feq for direction i=2 (cx=-1, cy=0), w=1/9:
+    //   feq = (2p - RHO*(-2u² + 2u + v²)) / 6
+    float feq2_target = (2.0f*p_out - RHO*(-2.0f*u_target*u_target + 2.0f*u_target + v_target*v_target)) / 6.0f;
+    float feq2_neb    = (2.0f*p_out - RHO*(-2.0f*u_neb*u_neb + 2.0f*u_neb + v_neb*v_neb)) / 6.0f;
+
+    // feq for direction i=8 (cx=-1, cy=1), w=1/36:
+    //   feq = (p + RHO*(u² - 3uv - u + v² + v)) / 12
+    float feq8_target = (p_out + RHO*(u_target*u_target - 3.0f*u_target*v_target - u_target + v_target*v_target + v_target)) / 12.0f;
+    float feq8_neb    = (p_out + RHO*(u_neb*u_neb - 3.0f*u_neb*v_neb - u_neb + v_neb*v_neb + v_neb)) / 12.0f;
+
+    // feq for direction i=6 (cx=-1, cy=-1), w=1/36:
+    //   feq = (p + RHO*(u² + 3uv - u + v² - v)) / 12
+    float feq6_target = (p_out + RHO*(u_target*u_target + 3.0f*u_target*v_target - u_target + v_target*v_target - v_target)) / 12.0f;
+    float feq6_neb    = (p_out + RHO*(u_neb*u_neb + 3.0f*u_neb*v_neb - u_neb + v_neb*v_neb - v_neb)) / 12.0f;
+
+    f[2] = f_neb[2] - feq2_neb + feq2_target;
+    f[8] = f_neb[8] - feq8_neb + feq8_target;
+    f[6] = f_neb[6] - feq6_neb + feq6_target;
+}
+
+#endif // NQ == 9
+
+// ============================================================================
+// D3Q19 inlet / outlet  (non-equilibrium extrapolation)
+//
+// Parabolic inlet (x=0):  reconstruct cx>0 populations i=1,7,9,13,15
+// Pressure outlet (x=NX-1): reconstruct cx<0 populations i=2,8,10,14,16
+//
+// Uses generic feq computation from macro.cuh to avoid hand-expanded formulas.
+// ============================================================================
+#if NQ == 19
+
+__device__ inline void apply_parabolic_inlet_3d(float* __restrict__ f,
+                                                const float* __restrict__ f_neb,
+                                                float y_coord)
+{
+    // Neighbor macros
+    float rho_neb, un, vn, wn;
+    compute_rho_u(f_neb, rho_neb, un, vn, wn);
+
+    // Target velocity (parabolic in y, uniform in z)
+    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
+    float u_tar = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+
+    // feq arrays
+    float feq_tar[19], feq_neb[19];
+    compute_feq(rho_neb, u_tar, 0.0f, 0.0f, feq_tar);
+    compute_feq(rho_neb, un, vn, wn, feq_neb);
+
+    // Reconstruct cx>0 directions: i = 1, 7, 9, 13, 15
+    f[1]  = f_neb[1]  - feq_neb[1]  + feq_tar[1];
+    f[7]  = f_neb[7]  - feq_neb[7]  + feq_tar[7];
+    f[9]  = f_neb[9]  - feq_neb[9]  + feq_tar[9];
+    f[13] = f_neb[13] - feq_neb[13] + feq_tar[13];
+    f[15] = f_neb[15] - feq_neb[15] + feq_tar[15];
+}
+
+__device__ inline void apply_pressure_outlet_3d(float* __restrict__ f,
+                                                const float* __restrict__ f_neb,
+                                                float y_coord)
+{
+    // Neighbor macros
+    float rho_neb, un, vn, wn;
+    compute_rho_u(f_neb, rho_neb, un, vn, wn);
+
+    // Target: p_out = 0 gauge → rho = 1.0  (using neighbor velocity)
+    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
+    float u_tar = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+
+    float feq_tar[19], feq_neb[19];
+    compute_feq(RHO, u_tar, 0.0f, 0.0f, feq_tar);
+    compute_feq(RHO, un, vn, wn, feq_neb);
+
+    // Reconstruct cx<0 directions: i = 2, 8, 10, 14, 16
+    f[2]  = f_neb[2]  - feq_neb[2]  + feq_tar[2];
+    f[8]  = f_neb[8]  - feq_neb[8]  + feq_tar[8];
+    f[10] = f_neb[10] - feq_neb[10] + feq_tar[10];
+    f[14] = f_neb[14] - feq_neb[14] + feq_tar[14];
+    f[16] = f_neb[16] - feq_neb[16] + feq_tar[16];
+}
+
+#endif // NQ == 19
+
+#endif // CELERIS_BOUNDARY_INLET_OUTLET_CUH

src/CelerisLab/lbm/kernels/core/descriptors.cuh

@@ -0,0 +1,75 @@
+// CelerisLab – core/descriptors.cuh
+// Lattice descriptor constants.
+// Direction ordering: paired for Esoteric-Pull compatibility.
+//   i=0 rest; for i>=1, (i, i^1|1) are opposite pairs:
+//     opp(0) = 0
+//     opp(i) = ((i-1) ^ 1) + 1   for i >= 1
+//
+// D2Q9 :  (0)rest  (1,2)±x  (3,4)±y  (5,6)±(x+y)  (7,8)±(x-y)
+// D3Q19:  (0)rest  (1,2)±x  (3,4)±y  (5,6)±z
+//         (7,8)±(x+y)  (9,10)±(x+z)  (11,12)±(y+z)
+//         (13,14)±(x-y)  (15,16)±(x-z)  (17,18)±(y-z)
+// ============================================================================
+
+#ifndef CELERIS_CORE_DESCRIPTORS_CUH
+#define CELERIS_CORE_DESCRIPTORS_CUH
+
+// NQ and DIM come from macros.h (included before this header)
+
+// ---------------------------------------------------------------------------
+// D2Q9
+// ---------------------------------------------------------------------------
+#if NQ == 9
+
+__constant__ int d_cx[9] = { 0,  1, -1,  0,  0,  1, -1,  1, -1};
+__constant__ int d_cy[9] = { 0,  0,  0,  1, -1,  1, -1, -1,  1};
+
+__constant__ float d_w[9] = {
+    4.0f/9.0f,
+    1.0f/9.0f,  1.0f/9.0f,  1.0f/9.0f,  1.0f/9.0f,
+    1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f
+};
+
+// Named weight macros (avoid __constant__ reads in hot paths)
+#define W0_VAL  (4.0f/9.0f)
+#define WS_VAL  (1.0f/9.0f)
+#define WE_VAL  (1.0f/36.0f)
+
+// ---------------------------------------------------------------------------
+// D3Q19
+// ---------------------------------------------------------------------------
+#elif NQ == 19
+
+__constant__ int d_cx[19] = { 0,  1,-1,  0, 0,  0, 0,  1,-1,  1,-1,  0, 0,  1,-1,  1,-1,  0, 0};
+__constant__ int d_cy[19] = { 0,  0, 0,  1,-1,  0, 0,  1,-1,  0, 0,  1,-1, -1, 1,  0, 0,  1,-1};
+__constant__ int d_cz[19] = { 0,  0, 0,  0, 0,  1,-1,  0, 0,  1,-1,  1,-1,  0, 0, -1, 1, -1, 1};
+
+__constant__ float d_w[19] = {
+    1.0f/3.0f,
+    1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f,
+    1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f,
+    1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f
+};
+
+#define W0_VAL  (1.0f/3.0f)
+#define WS_VAL  (1.0f/18.0f)
+#define WE_VAL  (1.0f/36.0f)
+
+#else
+#error "Unsupported NQ. Use 9 (D2Q9) or 19 (D3Q19)."
+#endif
+
+// ---------------------------------------------------------------------------
+// Lattice sound speed
+// ---------------------------------------------------------------------------
+#define CS2      (1.0f / 3.0f)
+#define CS2_INV  (3.0f)
+
+// ---------------------------------------------------------------------------
+// Opposite direction (works for any paired-ordered lattice)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ int opp_dir(int i) {
+    return (i == 0) ? 0 : (((i - 1) ^ 1) + 1);
+}
+
+#endif // CELERIS_CORE_DESCRIPTORS_CUH

src/CelerisLab/lbm/kernels/core/flags.cuh

@@ -0,0 +1,86 @@
+// CelerisLab – core/flags.cuh
+// Cell flag definitions (uint32_t, layered bytes) and helper functions.
+// ============================================================================
+
+#ifndef CELERIS_CORE_FLAGS_CUH
+#define CELERIS_CORE_FLAGS_CUH
+
+#include <stdint.h>
+
+// ---------------------------------------------------------------------------
+// Byte 0 [0:7]  – basic cell type
+// ---------------------------------------------------------------------------
+#define FLAG_FLUID      0x01u
+#define FLAG_SOLID      0x02u
+#define FLAG_GAS        0x04u
+#define FLAG_INTERFACE  0x08u
+#define FLAG_SENSOR     0x10u
+
+// ---------------------------------------------------------------------------
+// Byte 1 [8:15] – boundary condition type
+// ---------------------------------------------------------------------------
+#define FLAG_BB         (0x01u << 8)   // full-way bounce-back
+#define FLAG_EQ_BC      (0x02u << 8)   // equilibrium (Dirichlet) boundary
+#define FLAG_CURVED     (0x04u << 8)   // curved-boundary interpolation
+#define FLAG_IBM        (0x08u << 8)   // immersed boundary marker
+#define FLAG_MOVING     (0x10u << 8)   // moving-wall velocity
+
+// ---------------------------------------------------------------------------
+// Byte 2 [16:23] – multi-GPU / AMR / ghost
+// ---------------------------------------------------------------------------
+#define FLAG_GHOST      (0x01u << 16)
+#define FLAG_HALO       (0x02u << 16)
+#define FLAG_AMR_FINE   (0x04u << 16)
+#define FLAG_AMR_COARSE (0x08u << 16)
+
+// ---------------------------------------------------------------------------
+// Byte 3 [24:31] – reserved for extensions (phase-field, particles …)
+// ---------------------------------------------------------------------------
+
+// ---------------------------------------------------------------------------
+// Masks
+// ---------------------------------------------------------------------------
+#define MASK_CELL_TYPE  0x000000FFu
+#define MASK_BC_TYPE    0x0000FF00u
+#define MASK_MULTIGPU   0x00FF0000u
+#define MASK_EXTENSION  0xFF000000u
+
+// ---------------------------------------------------------------------------
+// Legacy compatibility (current driver.py uses uint8 with these bits)
+// ---------------------------------------------------------------------------
+#define LEGACY_FLUID     0x01
+#define LEGACY_SOLID     0x02
+#define LEGACY_GAS       0x04
+#define LEGACY_OBSTACLE  0x04   // obstacle / immersed body (triggers BB at adjacent fluid)
+#define LEGACY_INTERFACE 0x08
+#define LEGACY_SENSOR    0x10
+
+// ---------------------------------------------------------------------------
+// Device helper functions
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ bool is_fluid(uint32_t f)     { return (f & FLAG_FLUID) != 0; }
+__device__ __forceinline__ bool is_solid(uint32_t f)     { return (f & FLAG_SOLID) != 0; }
+__device__ __forceinline__ bool is_gas(uint32_t f)       { return (f & FLAG_GAS) != 0; }
+__device__ __forceinline__ bool is_interface(uint32_t f) { return (f & FLAG_INTERFACE) != 0; }
+__device__ __forceinline__ bool is_sensor(uint32_t f)    { return (f & FLAG_SENSOR) != 0; }
+
+__device__ __forceinline__ bool is_bb(uint32_t f)        { return (f & FLAG_BB) != 0; }
+__device__ __forceinline__ bool is_curved(uint32_t f)    { return (f & FLAG_CURVED) != 0; }
+__device__ __forceinline__ bool is_ibm(uint32_t f)       { return (f & FLAG_IBM) != 0; }
+__device__ __forceinline__ bool is_moving(uint32_t f)    { return (f & FLAG_MOVING) != 0; }
+__device__ __forceinline__ bool is_eq_bc(uint32_t f)     { return (f & FLAG_EQ_BC) != 0; }
+
+__device__ __forceinline__ bool is_boundary(uint32_t f)  { return (f & MASK_BC_TYPE) != 0; }
+__device__ __forceinline__ bool is_wall(uint32_t f)      { return (f & (FLAG_BB | FLAG_CURVED | FLAG_MOVING)) != 0; }
+
+// ---------------------------------------------------------------------------
+// Legacy flag promotion  (uint8 flags → uint32 with inferred BC bits)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ uint32_t promote_legacy_flag(uint8_t legacy) {
+    uint32_t f = (uint32_t)legacy & 0x1Fu;              // copy low 5 bits
+    if ((f & FLAG_SOLID) && (f & FLAG_INTERFACE))        // SOLID + INTERFACE
+        f |= FLAG_CURVED;                                //  → mark as curved BC
+    return f;
+}
+
+#endif // CELERIS_CORE_FLAGS_CUH

src/CelerisLab/lbm/kernels/core/layout.cuh

@@ -0,0 +1,155 @@
+// CelerisLab – core/layout.cuh
+// SoA memory layout: index_f, coordinates, and neighbor computation.
+// Depends on: NQ, NX, NY, NZ from macros.h; descriptors.cuh for d_cx/d_cy.
+// ============================================================================
+
+#ifndef CELERIS_CORE_LAYOUT_CUH
+#define CELERIS_CORE_LAYOUT_CUH
+
+// ---------------------------------------------------------------------------
+// Total cell count (compile-time)
+// ---------------------------------------------------------------------------
+#define TOTAL_CELLS ((unsigned long)(NX) * (unsigned long)(NY) * (unsigned long)(NZ))
+
+// ---------------------------------------------------------------------------
+// SoA index:  fi[ direction * N + node ]   ("Q-major")
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ unsigned long index_f(unsigned long n, unsigned int i) {
+    return (unsigned long)i * TOTAL_CELLS + n;
+}
+
+// ---------------------------------------------------------------------------
+// SoA index for DIM-component field:  u[ component * N + node ]
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ unsigned long index_u(unsigned long n, unsigned int d) {
+    return (unsigned long)d * TOTAL_CELLS + n;
+}
+
+// ---------------------------------------------------------------------------
+// 2D linear index → (x, y)  and  (x, y) → linear
+// ---------------------------------------------------------------------------
+#if DIM == 2
+
+__device__ __forceinline__ void coordinates(unsigned long n,
+                                            unsigned int& x,
+                                            unsigned int& y)
+{
+    x = (unsigned int)(n % NX);
+    y = (unsigned int)(n / NX);
+}
+
+__device__ __forceinline__ unsigned long linear_index(unsigned int x,
+                                                      unsigned int y)
+{
+    return (unsigned long)y * NX + x;
+}
+
+__device__ __forceinline__ void index_from_thread(unsigned int& x,
+                                                  unsigned int& y,
+                                                  unsigned long& k)
+{
+    x = threadIdx.x + NT * blockIdx.x;
+    y = blockIdx.y;
+    k = (unsigned long)y * NX + x;
+}
+
+// ---------------------------------------------------------------------------
+// Neighbor list (periodic BCs via modular arithmetic)
+// ---------------------------------------------------------------------------
+__device__ inline void compute_neighbors(unsigned long n, unsigned long* j) {
+    unsigned int x, y;
+    coordinates(n, x, y);
+
+    unsigned int xp = (x + 1u) % NX;
+    unsigned int xm = (x + NX - 1u) % NX;
+    unsigned int yp = (y + 1u) % NY;
+    unsigned int ym = (y + NY - 1u) % NY;
+
+    // i=0: self
+    j[0] = n;
+
+#if NQ == 9
+    // paired: (1,2)±x  (3,4)±y  (5,6)±(x+y)  (7,8)±(x-y)
+    j[1] = linear_index(xp, y );   // +x
+    j[2] = linear_index(xm, y );   // -x
+    j[3] = linear_index(x , yp);   // +y
+    j[4] = linear_index(x , ym);   // -y
+    j[5] = linear_index(xp, yp);   // +x +y
+    j[6] = linear_index(xm, ym);   // -x -y
+    j[7] = linear_index(xp, ym);   // +x -y
+    j[8] = linear_index(xm, yp);   // -x +y
+#elif NQ == 19
+    #error "D3Q19 requires DIM == 3. Set DIM=3 in macros.h."
+#endif
+}
+
+#elif DIM == 3
+
+__device__ __forceinline__ void coordinates(unsigned long n,
+                                            unsigned int& x,
+                                            unsigned int& y,
+                                            unsigned int& z)
+{
+    x = (unsigned int)(n % NX);
+    y = (unsigned int)(n / NX % NY);
+    z = (unsigned int)(n / ((unsigned long)NX * NY));
+}
+
+__device__ __forceinline__ unsigned long linear_index(unsigned int x,
+                                                      unsigned int y,
+                                                      unsigned int z)
+{
+    return (unsigned long)z * NY * NX + (unsigned long)y * NX + x;
+}
+
+__device__ __forceinline__ void index_from_thread(unsigned int& x,
+                                                  unsigned int& y,
+                                                  unsigned int& z,
+                                                  unsigned long& k)
+{
+    x = threadIdx.x + NT * blockIdx.x;
+    y = blockIdx.y;
+    z = blockIdx.z;
+    k = linear_index(x, y, z);
+}
+
+__device__ inline void compute_neighbors(unsigned long n, unsigned long* j) {
+    unsigned int x, y, z;
+    coordinates(n, x, y, z);
+
+    unsigned int xp = (x + 1u) % NX;
+    unsigned int xm = (x + NX - 1u) % NX;
+    unsigned int yp = (y + 1u) % NY;
+    unsigned int ym = (y + NY - 1u) % NY;
+    unsigned int zp = (z + 1u) % NZ;
+    unsigned int zm = (z + NZ - 1u) % NZ;
+
+    j[0]  = n;
+
+#if NQ == 19
+    // ±x, ±y, ±z (straight)
+    j[1]  = linear_index(xp, y,  z );   // +x
+    j[2]  = linear_index(xm, y,  z );   // -x
+    j[3]  = linear_index(x,  yp, z );   // +y
+    j[4]  = linear_index(x,  ym, z );   // -y
+    j[5]  = linear_index(x,  y,  zp);   // +z
+    j[6]  = linear_index(x,  y,  zm);   // -z
+    // diagonal
+    j[7]  = linear_index(xp, yp, z );   // +x+y
+    j[8]  = linear_index(xm, ym, z );   // -x-y
+    j[9]  = linear_index(xp, y,  zp);   // +x+z
+    j[10] = linear_index(xm, y,  zm);   // -x-z
+    j[11] = linear_index(x,  yp, zp);   // +y+z
+    j[12] = linear_index(x,  ym, zm);   // -y-z
+    j[13] = linear_index(xp, ym, z );   // +x-y
+    j[14] = linear_index(xm, yp, z );   // -x+y
+    j[15] = linear_index(xp, y,  zm);   // +x-z
+    j[16] = linear_index(xm, y,  zp);   // -x+z
+    j[17] = linear_index(x,  yp, zm);   // +y-z
+    j[18] = linear_index(x,  ym, zp);   // -y+z
+#endif
+}
+
+#endif // DIM == 3
+
+#endif // CELERIS_CORE_LAYOUT_CUH

src/CelerisLab/lbm/kernels/core/params.cuh

@@ -0,0 +1,62 @@
+// CelerisLab – core/params.cuh
+// Runtime parameter structures transported via __constant__ memory.
+// Includes: LBMParams, RigidBodyState2D, RigidBodyControl2D.
+// ============================================================================
+
+#ifndef CELERIS_CORE_PARAMS_CUH
+#define CELERIS_CORE_PARAMS_CUH
+
+#include <stdint.h>
+
+// ============================================================================
+// LBM runtime parameters  (uploaded once per run or on parameter change)
+// ============================================================================
+struct LBMParams {
+    //--- grid ---
+    unsigned int  Nx, Ny, Nz;
+    unsigned long N;               // Nx * Ny * Nz
+
+    //--- relaxation ---
+    float omega;                   // SRT/TRT主松弛率  w = 1/(3ν + 0.5)
+    float omega_bulk;              // MRT bulk 松弛率   (s_e / s_eps)
+
+    //--- external body force ---
+    float fx, fy, fz;
+
+    //--- reference quantities ---
+    float rho_ref;                 // 参考密度 (通常 1.0)
+    float u_inlet;                 // 入口参考速度
+
+    //--- diagnostics ---
+    unsigned int n_objects;        // 观测对象数量
+};
+
+__constant__ LBMParams d_params;
+
+// ============================================================================
+// Rigid-body state / control  (2D, per architecture §17.3)
+// ============================================================================
+struct RigidBodyState2D {
+    float x, y, theta;            // 位姿 (position + orientation)
+    float vx, vy, omega;          // 速度 (translational + rotational)
+};
+
+struct RigidBodyControl2D {
+    float ax_cmd, ay_cmd;         // 目标加速度（或目标速度/位移, 取决于 mode）
+    float alpha_cmd;              // 目标角加速度 / 角速度
+    int   control_mode;           // 0=力控, 1=速度控, 2=位移控
+};
+
+// ============================================================================
+// Halo plan (multi-GPU / AMR 预留)
+// ============================================================================
+struct HaloPlan {
+    int       face_id;
+    int       peer_gpu;
+    int       level;              // AMR level (0 = base)
+    size_t    count;
+    int*      send_idx;
+    int*      recv_idx;
+};
+
+#endif // CELERIS_CORE_PARAMS_CUH

src/CelerisLab/lbm/kernels/core/precision.cuh

@@ -0,0 +1,91 @@
+// CelerisLab – core/precision.cuh
+// Storage precision abstraction (FP32 / FP16S / FP16C).
+// Compute precision is always float (FP32).
+//
+// Controlled by macros.h:
+//   STORE_PRECISION  0 = FP32 (default)
+//                    1 = FP16S (IEEE-754 half + ×32768 scaling)
+//                    2 = FP16C (custom 1-4-11 format, higher precision)
+// ============================================================================
+
+#ifndef CELERIS_CORE_PRECISION_CUH
+#define CELERIS_CORE_PRECISION_CUH
+
+#ifndef STORE_PRECISION
+#define STORE_PRECISION 0
+#endif
+
+// ============================================================================
+// FP16S – Scaled IEEE-754 Half  (hardware-accelerated conversion)
+// ============================================================================
+#if STORE_PRECISION == 1
+
+#include <cuda_fp16.h>
+using fpxx = __half;
+
+__device__ __forceinline__ float load_ddf(const fpxx* p, unsigned long idx) {
+    return __half2float(p[idx]) * 3.0517578125e-5f;   // ÷ 32768
+}
+__device__ __forceinline__ void store_ddf(fpxx* p, unsigned long idx, float v) {
+    p[idx] = __float2half_rn(v * 32768.0f);
+}
+
+#define FPXX_BYTES 2
+
+// ============================================================================
+// FP16C – Custom 1-4-11 format  (software conversion, ~3.6 dec. digits)
+// ============================================================================
+#elif STORE_PRECISION == 2
+
+using fpxx = unsigned short;
+
+__device__ __forceinline__ float half_custom_to_float(unsigned short x) {
+    unsigned int e = (x & 0x7800) >> 11;
+    unsigned int m = (x & 0x07FF) << 12;
+    unsigned int v = __float_as_uint((float)m) >> 23;
+    return __uint_as_float(
+        (x & 0x8000) << 16 |
+        (e != 0) * ((e + 112) << 23 | m) |
+        ((e == 0) & (m != 0)) * ((v - 37) << 23 | ((m << (150 - v)) & 0x007FF000))
+    );
+}
+
+__device__ __forceinline__ unsigned short float_to_half_custom(float x) {
+    unsigned int b = __float_as_uint(x) + 0x00000800;  // round-to-nearest-even
+    unsigned int e = (b & 0x7F800000) >> 23;
+    unsigned int m = b & 0x007FFFFF;
+    return (unsigned short)(
+        (b & 0x80000000) >> 16 |
+        (e > 112) * ((((e - 112) << 11) & 0x7800) | m >> 12) |
+        ((e < 113) & (e > 100)) * ((((0x007FF800 + m) >> (124 - e)) + 1) >> 1)
+    );
+}
+
+__device__ __forceinline__ float load_ddf(const fpxx* p, unsigned long idx) {
+    return half_custom_to_float(p[idx]);
+}
+__device__ __forceinline__ void store_ddf(fpxx* p, unsigned long idx, float v) {
+    p[idx] = float_to_half_custom(v);
+}
+
+#define FPXX_BYTES 2
+
+// ============================================================================
+// FP32 – No conversion (default)
+// ============================================================================
+#else
+
+using fpxx = float;
+
+__device__ __forceinline__ float load_ddf(const fpxx* p, unsigned long idx) {
+    return p[idx];
+}
+__device__ __forceinline__ void store_ddf(fpxx* p, unsigned long idx, float v) {
+    p[idx] = v;
+}
+
+#define FPXX_BYTES 4
+
+#endif  // STORE_PRECISION
+
+#endif  // CELERIS_CORE_PRECISION_CUH

src/CelerisLab/lbm/kernels/kernel_v2.cu

@@ -0,0 +1,299 @@
+// CelerisLab – kernel_v2.cu
+// ============================================================================
+// Modular compilation entry point  (Stage 1 Architecture)
+//
+// This file includes all modular headers and step kernels,
+// and exports extern "C" functions for PyCUDA.
+//
+// Controlled by macros.h:
+//   NQ              = 9 (D2Q9) or 19 (D3Q19)
+//   COLLISION_MODEL = 0 (SRT) / 1 (TRT) / 2 (MRT)
+//   STREAMING_MODEL = 0 (double-buffer) / 1 (Esoteric-Pull)
+//   STORE_PRECISION = 0 (FP32) / 1 (FP16S) / 2 (FP16C)
+//   USE_DDF_SHIFTING= 0 / 1
+// ============================================================================
+
+#include <stdio.h>
+#include <stdint.h>
+#include <cuda.h>
+
+// ---------------------------------------------------------------------------
+// Layer 0: Configuration (compile-time)
+// ---------------------------------------------------------------------------
+#include "macros.h"
+
+// ---------------------------------------------------------------------------
+// Layer 1: Core primitives
+// ---------------------------------------------------------------------------
+#include "core/precision.cuh"
+#include "core/flags.cuh"
+#include "core/descriptors.cuh"
+#include "core/layout.cuh"
+#include "core/params.cuh"
+
+// ---------------------------------------------------------------------------
+// Layer 2: Operators
+// ---------------------------------------------------------------------------
+#include "operators/macro.cuh"
+#include "operators/collision_srt.cuh"
+#include "operators/collision_trt.cuh"
+#include "operators/collision_mrt.cuh"
+#include "operators/forcing_guo.cuh"
+
+// ---------------------------------------------------------------------------
+// Layer 3: Streaming
+// ---------------------------------------------------------------------------
+#include "streaming/pull_double_buffer.cuh"
+#include "streaming/esopull_single_buffer.cuh"
+
+// ---------------------------------------------------------------------------
+// Layer 4: Boundary conditions
+// ---------------------------------------------------------------------------
+#include "boundary/bounce_back.cuh"
+#include "boundary/curved_boundary.cuh"
+#include "boundary/inlet_outlet.cuh"
+
+// ---------------------------------------------------------------------------
+// Layer 5: Step kernels (selected by STREAMING_MODEL)
+// ---------------------------------------------------------------------------
+#include "step/one_step_double.cu"
+#include "step/one_step_esopull.cu"
+
+// ---------------------------------------------------------------------------
+// Layer 6: IBM kernels (always available, launched separately)
+// ---------------------------------------------------------------------------
+#include "boundary/ibm_kernels.cuh"
+
+// ============================================================================
+// Extern "C" wrappers  (for PyCUDA / ctypes)
+//
+// Naming:
+//   OneStep            – backward-compatible main step (double-buffer)
+//   InitTubeFlow_v2    – channel initialization (new ordering)
+//   StreamCollide      – new API main step (double or esopull)
+//   CurvedBoundary     – post-step curved boundary correction
+//   UpdateFields       – recompute rho/u from DDF
+//   InitEsoPull        – Esoteric-Pull initialization
+//   IBM_UpdateLag      – IBM Lagrangian point update
+//   IBM_Eul2Lag        – Euler → Lagrangian interpolation
+//   IBM_Lag2Eul        – Lagrangian → Euler spreading
+// ============================================================================
+extern "C"
+{
+
+// ----- Backward-compatible init (uses new direction ordering) -----
+__global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi)
+{
+#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;
+
+    float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1))
+                                             * ((float)y - 0.5f * (NY - 1))
+                                / ((float)(NY - 2) * (float)(NY - 2)));
+
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++)
+            store_ddf(fi, index_f(k, (unsigned int)i), 0.0f);
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;
+            float val = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init);
+#if USE_DDF_SHIFTING
+            val -= d_w[i];
+#endif
+            store_ddf(fi, index_f(k, (unsigned int)i), val);
+        }
+    }
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+
+    float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1))
+                                             * ((float)y - 0.5f * (NY - 1))
+                                / ((float)(NY - 2) * (float)(NY - 2)));
+
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++)
+            store_ddf(fi, index_f(k, (unsigned int)i), 0.0f);
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;
+            float val = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init);
+#if USE_DDF_SHIFTING
+            val -= d_w[i];
+#endif
+            store_ddf(fi, index_f(k, (unsigned int)i), val);
+        }
+    }
+#endif
+}
+
+// ----- Main step (double-buffer) -----
+// Signature compatible with driver.py: flag, fi_in, fi_out, indx, delta, action, obs
+__global__ void OneStep(
+    uint8_t*  flag,
+    fpxx*     fi_in,
+    fpxx*     fi_out,
+    int32_t*  indx,
+    float*    delta,
+    float*    action,
+    float*    obs)
+{
+    // Redirect to the modular StreamCollideDouble kernel body
+    // (We inline here to keep a single extern "C" entry point.)
+
+#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;
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+#endif
+
+    uint8_t fl = flag[k];
+
+    unsigned long j[NQ];
+    compute_neighbors(k, j);
+
+    float f[NQ];
+
+#if STREAMING_MODEL == 0
+    stream_pull_load(k, f, fi_in, j);
+#else
+    // Esoteric-Pull for OneStep compat would need timestep — not available
+    // in legacy signature. Use double-buffer mode only for OneStep.
+    stream_pull_load(k, f, fi_in, j);
+#endif
+
+    float rho_n, ux, uy;
+#if NQ == 9
+    compute_rho_u(f, rho_n, ux, uy);
+
+    // Inlet / outlet / wall bounce-back
+    if (fl & LEGACY_SOLID) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x + 1u, y);
+            for (int i = 0; i < NQ; i++)
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            apply_parabolic_inlet(f, f_neb, (float)y);
+        }
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x - 1u, y);
+            for (int i = 0; i < NQ; i++)
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            apply_pressure_outlet(f, f_neb, (float)y);
+        }
+        else {
+            // Wall / corner: full-way bounce-back (swap all pairs)
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+            }
+        }
+    }
+
+    // Collision
+    if (fl & LEGACY_FLUID) {
+        float feq[NQ], Fin[NQ];
+        compute_feq(rho_n, ux, uy, feq);
+        zero_forcing(Fin);
+
+#if COLLISION_MODEL == 0
+        collide_srt(f, feq, Fin, d_params.omega);
+#elif COLLISION_MODEL == 1
+        collide_trt(f, feq, Fin, d_params.omega);
+#elif COLLISION_MODEL == 2
+        collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
+#endif
+    }
+
+    stream_pull_store(k, f, fi_out);
+
+    // Curved boundary processing
+    if ((fl & LEGACY_SOLID) && (fl & LEGACY_INTERFACE)) {
+        int id_off = indx[k];
+        int id_obj = *reinterpret_cast<int*>(&delta[id_off]);
+        float Uw = action[id_obj] * delta[id_off + NQ];
+        float Vw = action[id_obj] * delta[id_off + NQ + 1];
+        float* obs_fx = &obs[DIM * id_obj];
+        float* obs_fy = &obs[DIM * id_obj + 1];
+        apply_curved_boundary(k, x, y, fi_out, delta, id_off, Uw, Vw, obs_fx, obs_fy);
+    }
+
+    // Sensor
+    if (fl & LEGACY_SENSOR) {
+        int id_obj = indx[k];
+        atomicAdd(&obs[DIM * id_obj],     ux);
+        atomicAdd(&obs[DIM * id_obj + 1], uy);
+    }
+
+#elif NQ == 19
+    // ---- Macroscopic (computed for all nodes; only used for fluid collision) ----
+    float uz;
+    compute_rho_u(f, rho_n, ux, uy, uz);
+
+    // ---- Boundary conditions ----
+    if (fl & LEGACY_SOLID) {
+        // Interior of inlet/outlet planes (skip corners at y=0/NY-1)
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
+            // Parabolic velocity inlet
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x + 1u, y, z);
+            for (int i = 0; i < NQ; i++)
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            apply_parabolic_inlet_3d(f, f_neb, (float)y);
+        }
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
+            // Pressure outlet
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x - 1u, y, z);
+            for (int i = 0; i < NQ; i++)
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            apply_pressure_outlet_3d(f, f_neb, (float)y);
+        }
+        else {
+            // Wall nodes and corners: full-way bounce-back.
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+            }
+        }
+    }
+    // Obstacle nodes: full-way bounce-back (same swap)
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+        }
+    }
+
+    // ---- Collision (fluid only) ----
+    if (fl & LEGACY_FLUID) {
+        float feq[NQ], Fin[NQ];
+        compute_feq(rho_n, ux, uy, uz, feq);
+        zero_forcing(Fin);
+        collide_srt(f, feq, Fin, d_params.omega);
+    }
+
+    stream_pull_store(k, f, fi_out);
+#endif
+}
+
+}  // extern "C"

src/CelerisLab/lbm/kernels/macros.h

@@ -3,23 +3,23 @@
 // cuda parameters
 #define MULT_GPU False
 #define NT 128
-#define X_1U 128
+#define X_1U 256
-#define Y_1U 32
+#define Y_1U 128
-#define Z_1U 1
+#define Z_1U 32
 
 // flow parameters
 #define LBtype float
-#define UX 10
+#define UX 1
-#define UY 16
+#define UY 1
 #define UZ 1
-#define NX 1280
+#define NX 256
-#define NY 512
+#define NY 128
-#define NZ 1
+#define NZ 32
-#define DIM 2
+#define DIM 3
-#define NQ 9
+#define NQ 19
-#define VIS 0.004
+#define VIS 0.0096000000
 #define RHO 1.0
-#define U0 0.01
+#define U0 0.04
 
 // constants
 #define PI 3.141592653589793238
@@ -33,5 +33,31 @@
 #define V_TAYLOR    0b00000001
 
 // variables
-#define N_OBJS 7
+#define N_OBJS 0
 // #define N_SENS 2
+
+// ============================================================================
+// New architecture configuration (Stage 1)
+// These defaults are safe for backward compatibility.
+// compiler.py can override any of them via modify_macro().
+// ============================================================================
+
+// Collision model: 0=SRT, 1=TRT, 2=MRT
+#ifndef COLLISION_MODEL
+#define COLLISION_MODEL 0
+#endif
+
+// Streaming model: 0=double-buffer, 1=esoteric-pull
+#ifndef STREAMING_MODEL
+#define STREAMING_MODEL 0
+#endif
+
+// Storage precision: 0=FP32, 1=FP16S, 2=FP16C
+#ifndef STORE_PRECISION
+#define STORE_PRECISION 0
+#endif
+
+// DDF-shifting: 0=off, 1=on
+#ifndef USE_DDF_SHIFTING
+#define USE_DDF_SHIFTING 0
+#endif

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

@@ -0,0 +1,135 @@
+// CelerisLab – operators/collision_mrt.cuh
+// MRT (Multiple-Relaxation-Time) collision operator.
+//
+// D2Q9 Lallemand-Luo (Phys. Rev. E 61, 2000) with NEW paired ordering:
+//   cx = {0, 1,-1, 0, 0, 1,-1, 1,-1}
+//   cy = {0, 0, 0, 1,-1, 1,-1,-1, 1}
+//
+// Moment space:
+//   m[0] = ρ     (conserved,  s₀ = 0 or 1.0 overwrite)
+//   m[1] = e     (energy,     s₁ = s_e)
+//   m[2] = ε     (energy²,   s₂ = s_eps)
+//   m[3] = jx    (conserved,  s₃ = 0 or 1.0 overwrite)
+//   m[4] = qx    (energy flux x, s₄ = s_q)
+//   m[5] = jy    (conserved,  s₅ = 0 or 1.0 overwrite)
+//   m[6] = qy    (energy flux y, s₆ = s_q)
+//   m[7] = pxx   (stress,     s₇ = s_nu = ω = 1/(3ν + 0.5))
+//   m[8] = pxy   (stress,     s₈ = s_nu)
+// ============================================================================
+
+#ifndef CELERIS_OPERATORS_COLLISION_MRT_CUH
+#define CELERIS_OPERATORS_COLLISION_MRT_CUH
+
+#if NQ == 9
+
+// ---------------------------------------------------------------------------
+// D2Q9 MRT  (fully expanded, no loops, no matrix storage)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ void collide_mrt(float* __restrict__ g,
+                                            float rho, float ux, float uy,
+                                            const float* __restrict__ Fin,
+                                            float omega)
+{
+    // ----- Relaxation rates -----
+    //   s_rho = s_jx = s_jy = 1.0  (conserved moments → overwrite to eq)
+    //   s_e   = s_eps = s_q  = 1.2  (kinetic transport)
+    //   s_nu  = omega               (viscosity-related)
+    const float s_rho = 1.0f;   // conserved moment relaxation
+    const float s_e   = 1.2f;
+    const float s_eps = 1.2f;
+    const float s_jx  = 1.0f;
+    const float s_q   = 1.2f;
+    const float s_jy  = 1.0f;
+    const float s_nu  = omega;
+
+    // ----- Pressure (from density) -----
+    const float p = (g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8]) / 3.0f;
+
+    // ----- Forward transform: m = M * g  (new paired ordering) -----
+    float m[9];
+    m[0] =    g[0] +   g[1] +   g[2] +   g[3] +   g[4] +   g[5] +   g[6] +   g[7] +   g[8];
+    m[1] = -4*g[0] -   g[1] -   g[2] -   g[3] -   g[4] + 2*g[5] + 2*g[6] + 2*g[7] + 2*g[8];
+    m[2] =  4*g[0] - 2*g[1] - 2*g[2] - 2*g[3] - 2*g[4] +   g[5] +   g[6] +   g[7] +   g[8];
+    m[3] =            g[1] -   g[2]                      +   g[5] -   g[6] +   g[7] -   g[8];
+    m[4] =         -2*g[1] + 2*g[2]                      +   g[5] -   g[6] +   g[7] -   g[8];
+    m[5] =                            g[3] -   g[4]      +   g[5] -   g[6] -   g[7] +   g[8];
+    m[6] =                         -2*g[3] + 2*g[4]      +   g[5] -   g[6] -   g[7] +   g[8];
+    m[7] =            g[1] +   g[2] -  g[3] -   g[4];
+    m[8] =                                                    g[5] +   g[6] -   g[7] -   g[8];
+
+    // ----- Equilibrium moments -----
+    const float u2 = ux * ux + uy * uy;
+    float meq[9];
+    meq[0] = 3.0f * p;                        // ρ
+    meq[1] = -6.0f * p + 3.0f * RHO * u2;     // e
+    meq[2] =  3.0f * p - 3.0f * RHO * u2;     // ε
+    meq[3] = RHO * ux;                         // jx
+    meq[4] = -RHO * ux;                        // qx
+    meq[5] = RHO * uy;                         // jy
+    meq[6] = -RHO * uy;                        // qy
+    meq[7] = RHO * (ux * ux - uy * uy);        // pxx
+    meq[8] = RHO * ux * uy;                     // pxy
+
+    // ----- Relaxation:  delta_m[i] = s_i * (meq[i] - m[i]) -----
+    float dm[9];
+    dm[0] = s_rho * (meq[0] - m[0]);
+    dm[1] = s_e   * (meq[1] - m[1]);
+    dm[2] = s_eps * (meq[2] - m[2]);
+    dm[3] = s_jx  * (meq[3] - m[3]);
+    dm[4] = s_q   * (meq[4] - m[4]);
+    dm[5] = s_jy  * (meq[5] - m[5]);
+    dm[6] = s_q   * (meq[6] - m[6]);
+    dm[7] = s_nu  * (meq[7] - m[7]);
+    dm[8] = s_nu  * (meq[8] - m[8]);
+
+    // ----- Inverse transform: g += M⁻¹ * dm  (new paired ordering) -----
+    g[0] += (    dm[0] -    dm[1] +    dm[2]                                              ) / 9.0f;
+    g[1] += (4 * dm[0] -    dm[1] - 2* dm[2] + 6* dm[3] - 6* dm[4]               + 9*dm[7]) / 36.0f;
+    g[2] += (4 * dm[0] -    dm[1] - 2* dm[2] - 6* dm[3] + 6* dm[4]               + 9*dm[7]) / 36.0f;
+    g[3] += (4 * dm[0] -    dm[1] - 2* dm[2]                        + 6*dm[5] - 6*dm[6] - 9*dm[7]) / 36.0f;
+    g[4] += (4 * dm[0] -    dm[1] - 2* dm[2]                        - 6*dm[5] + 6*dm[6] - 9*dm[7]) / 36.0f;
+    g[5] += (4 * dm[0] + 2* dm[1] +    dm[2] + 6* dm[3] + 3* dm[4] + 6*dm[5] + 3*dm[6] + 9*dm[8]) / 36.0f;
+    g[6] += (4 * dm[0] + 2* dm[1] +    dm[2] - 6* dm[3] - 3* dm[4] - 6*dm[5] - 3*dm[6] + 9*dm[8]) / 36.0f;
+    g[7] += (4 * dm[0] + 2* dm[1] +    dm[2] + 6* dm[3] + 3* dm[4] - 6*dm[5] - 3*dm[6] - 9*dm[8]) / 36.0f;
+    g[8] += (4 * dm[0] + 2* dm[1] +    dm[2] - 6* dm[3] - 3* dm[4] + 6*dm[5] + 3*dm[6] - 9*dm[8]) / 36.0f;
+
+    // ----- Add forcing (if present) -----
+    #pragma unroll
+    for (int i = 0; i < 9; i++) {
+        g[i] += Fin[i];
+    }
+}
+
+// Convenience wrapper: no external forcing
+__device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
+                                                     float rho, float ux, float uy,
+                                                     float omega)
+{
+    float Fin[9] = {0};
+    collide_mrt(g, rho, ux, uy, Fin, omega);
+}
+
+#elif NQ == 19
+
+// ---------------------------------------------------------------------------
+// D3Q19 MRT  (placeholder – use SRT or TRT until fully validated)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ void collide_mrt(float* __restrict__ g,
+                                            float rho, float ux, float uy, float uz,
+                                            const float* __restrict__ Fin,
+                                            float omega)
+{
+    // TODO: implement D3Q19 MRT with 19-moment Gram–Schmidt basis
+    // Fall back to SRT for now
+    float feq[19];
+    compute_feq(rho, ux, uy, uz, feq);
+    const float c_tau = 1.0f - 0.5f * omega;
+    #pragma unroll
+    for (int i = 0; i < 19; i++) {
+        g[i] = fmaf(1.0f - omega, g[i], fmaf(omega, feq[i], c_tau * Fin[i]));
+    }
+}
+
+#endif  // NQ
+
+#endif  // CELERIS_OPERATORS_COLLISION_MRT_CUH

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

@@ -0,0 +1,35 @@
+// CelerisLab – operators/collision_srt.cuh
+// SRT (BGK) collision operator.
+//   f_out[i] = (1 - ω) f[i] + ω feq[i] + Fin[i]
+// ============================================================================
+
+#ifndef CELERIS_OPERATORS_COLLISION_SRT_CUH
+#define CELERIS_OPERATORS_COLLISION_SRT_CUH
+
+__device__ __forceinline__ void collide_srt(float* __restrict__ f,
+                                            const float* __restrict__ feq,
+                                            const float* __restrict__ Fin,
+                                            float omega)
+{
+    // Pre-compute forcing prefactor (1 - ω/2) for Guo scheme
+    const float c_tau = 1.0f - 0.5f * omega;
+
+    #pragma unroll
+    for (int i = 0; i < NQ; i++) {
+        f[i] = fmaf(1.0f - omega, f[i],
+                     fmaf(omega, feq[i], c_tau * Fin[i]));
+    }
+}
+
+// Variant without forcing (Fin = 0)
+__device__ __forceinline__ void collide_srt_no_force(float* __restrict__ f,
+                                                     const float* __restrict__ feq,
+                                                     float omega)
+{
+    #pragma unroll
+    for (int i = 0; i < NQ; i++) {
+        f[i] = fmaf(1.0f - omega, f[i], omega * feq[i]);
+    }
+}
+
+#endif // CELERIS_OPERATORS_COLLISION_SRT_CUH

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

@@ -0,0 +1,79 @@
+// CelerisLab – operators/collision_trt.cuh
+// TRT (Two-Relaxation-Time) collision operator.
+//
+// Symmetric part   uses ω⁺ (= omega, transport viscosity)
+// Antisymmetric    uses ω⁻ (linked by magic parameter Λ = 3/16)
+//   ω⁻ = 1 / (Λ / (1/ω⁺ - 0.5) + 0.5)
+//
+// Paired direction ordering makes TRT natural:
+//   f_s[i] = 0.5*(f[i] + f[opp(i)])    symmetric
+//   f_a[i] = 0.5*(f[i] - f[opp(i)])    antisymmetric
+// ============================================================================
+
+#ifndef CELERIS_OPERATORS_COLLISION_TRT_CUH
+#define CELERIS_OPERATORS_COLLISION_TRT_CUH
+
+// Magic parameter Λ = 3/16 (optimal for porous-media / bounce-back wall location)
+#define TRT_MAGIC_PARAM (0.1875f)
+
+__device__ __forceinline__ float compute_omega_minus(float omega_plus) {
+    return 1.0f / (TRT_MAGIC_PARAM / (1.0f / omega_plus - 0.5f) + 0.5f);
+}
+
+__device__ __forceinline__ void collide_trt(float* __restrict__ f,
+                                            const float* __restrict__ feq,
+                                            const float* __restrict__ Fin,
+                                            float omega)
+{
+    const float wp = omega;
+    const float wm = compute_omega_minus(wp);
+
+    // Direction 0: rest particle – treat as symmetric-only
+    f[0] = f[0] - wp * (f[0] - feq[0]) + Fin[0];
+
+    // Direction pairs (i, i+1) for i = 1, 3, 5, …
+    #pragma unroll
+    for (int i = 1; i < NQ; i += 2) {
+        const int ib = i + 1;  // opposite (paired layout)
+
+        // Current and opposite
+        const float fi = f[i];
+        const float fb = f[ib];
+
+        // Equilibrium
+        const float feqi = feq[i];
+        const float feqb = feq[ib];
+
+        // Symmetric / antisymmetric non-equilibrium
+        const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb));
+        const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb));
+
+        // Relax + forcing
+        f[i]  = fi  - wp * delta_s - wm * delta_a + Fin[i];
+        f[ib] = fb  - wp * delta_s + wm * delta_a + Fin[ib];
+    }
+}
+
+// Variant without forcing (Fin = 0)
+__device__ __forceinline__ void collide_trt_no_force(float* __restrict__ f,
+                                                     const float* __restrict__ feq,
+                                                     float omega)
+{
+    const float wp = omega;
+    const float wm = compute_omega_minus(wp);
+
+    f[0] = f[0] - wp * (f[0] - feq[0]);
+
+    #pragma unroll
+    for (int i = 1; i < NQ; i += 2) {
+        const int ib = i + 1;
+        const float fi = f[i], fb = f[ib];
+        const float feqi = feq[i], feqb = feq[ib];
+        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;
+        f[ib] = fb  - wp * delta_s + wm * delta_a;
+    }
+}
+
+#endif // CELERIS_OPERATORS_COLLISION_TRT_CUH

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

@@ -0,0 +1,102 @@
+// CelerisLab – operators/forcing_guo.cuh
+// Guo forcing scheme  (Guo, Zheng & Shi, Phys. Rev. E 65, 2002).
+//
+// 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
+//   simplified: Fin[i] = 9·w_i · [ (c_i·F)(c_i·u* + 1/3) - u*·F/3 ]
+// ============================================================================
+
+#ifndef CELERIS_OPERATORS_FORCING_GUO_CUH
+#define CELERIS_OPERATORS_FORCING_GUO_CUH
+
+// ---------------------------------------------------------------------------
+// D2Q9 Guo forcing terms
+// ---------------------------------------------------------------------------
+#if NQ == 9
+
+__device__ __forceinline__ void apply_guo_velocity_correction(
+    float& ux, float& uy,
+    float fx, float fy, float rho)
+{
+    float rho2 = 0.5f / rho;
+    ux = fmaf(fx, rho2, ux);
+    uy = fmaf(fy, rho2, uy);
+}
+
+__device__ __forceinline__ void compute_guo_forcing(
+    float ux, float uy,           // velocity (already corrected: u*)
+    float fx, float fy,           // force density
+    float* __restrict__ Fin)
+{
+    const float uF = -0.33333334f * fmaf(ux, fx, uy * fy);   // -u*·F / 3
+
+    Fin[0] = 9.0f * W0_VAL * uF;
+
+    // Straight directions
+    //   cx = {0, 1,-1, 0, 0, 1,-1, 1,-1}
+    //   cy = {0, 0, 0, 1,-1, 1,-1,-1, 1}
+    Fin[1] = 9.0f * WS_VAL * fmaf( fx,  ux + 0.33333334f, uF);   // +x
+    Fin[2] = 9.0f * WS_VAL * fmaf(-fx, -ux + 0.33333334f, uF);   // -x
+    Fin[3] = 9.0f * WS_VAL * fmaf( fy,  uy + 0.33333334f, uF);   // +y
+    Fin[4] = 9.0f * WS_VAL * fmaf(-fy, -uy + 0.33333334f, uF);   // -y
+
+    // Diagonal directions
+    const float f5x = fx, f5y = fy;   // +x+y
+    const float c5u = ux + uy;
+    const float c5F = f5x + f5y;
+    Fin[5] = 9.0f * WE_VAL * fmaf(c5F, c5u + 0.33333334f, uF);
+
+    const float c6u = -ux - uy;        // -x-y
+    const float c6F = -fx - fy;
+    Fin[6] = 9.0f * WE_VAL * fmaf(c6F, c6u + 0.33333334f, uF);
+
+    const float c7u = ux - uy;         // +x-y
+    const float c7F = fx - fy;
+    Fin[7] = 9.0f * WE_VAL * fmaf(c7F, c7u + 0.33333334f, uF);
+
+    const float c8u = -ux + uy;        // -x+y
+    const float c8F = -fx + fy;
+    Fin[8] = 9.0f * WE_VAL * fmaf(c8F, c8u + 0.33333334f, uF);
+}
+
+#elif NQ == 19
+
+__device__ __forceinline__ void apply_guo_velocity_correction(
+    float& ux, float& uy, float& uz,
+    float fx, float fy, float fz, float rho)
+{
+    float rho2 = 0.5f / rho;
+    ux = fmaf(fx, rho2, ux);
+    uy = fmaf(fy, rho2, uy);
+    uz = fmaf(fz, rho2, uz);
+}
+
+__device__ __forceinline__ void compute_guo_forcing(
+    float ux, float uy, float uz,
+    float fx, float fy, float fz,
+    float* __restrict__ Fin)
+{
+    const float uF = -0.33333334f * fmaf(ux, fx, fmaf(uy, fy, uz * fz));
+
+    Fin[0] = 9.0f * W0_VAL * uF;
+
+    // Generic loop for D3Q19 (compiler will fully unroll)
+    for (int i = 1; i < 19; i++) {
+        float ci_dot_F = d_cx[i]*fx + d_cy[i]*fy + d_cz[i]*fz;
+        float ci_dot_u = d_cx[i]*ux + d_cy[i]*uy + d_cz[i]*uz;
+        Fin[i] = 9.0f * d_w[i] * fmaf(ci_dot_F, ci_dot_u + 0.33333334f, uF);
+    }
+}
+
+#endif // NQ
+
+// ---------------------------------------------------------------------------
+// Zero forcing helper (when no external force)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ void zero_forcing(float* __restrict__ Fin) {
+    #pragma unroll
+    for (int i = 0; i < NQ; i++) Fin[i] = 0.0f;
+}
+
+#endif // CELERIS_OPERATORS_FORCING_GUO_CUH

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

@@ -0,0 +1,168 @@
+// CelerisLab – operators/macro.cuh
+// Macroscopic quantity computation (rho, u) and equilibrium distribution (feq).
+// Supports DDF-shifting (USE_DDF_SHIFTING=1) and standard mode (=0).
+//
+// D2Q9 paired ordering:
+//   cx = {0, 1,-1, 0, 0, 1,-1, 1,-1}
+//   cy = {0, 0, 0, 1,-1, 1,-1,-1, 1}
+// ============================================================================
+
+#ifndef CELERIS_OPERATORS_MACRO_CUH
+#define CELERIS_OPERATORS_MACRO_CUH
+
+#ifndef USE_DDF_SHIFTING
+#define USE_DDF_SHIFTING 0
+#endif
+
+// ---------------------------------------------------------------------------
+// compute_rho_u:  f[NQ] → (rho, ux, uy [, uz])
+// ---------------------------------------------------------------------------
+#if NQ == 9
+
+__device__ __forceinline__ void compute_rho_u(const float* __restrict__ f,
+                                              float& rho,
+                                              float& ux, float& uy)
+{
+    // Density: sum of all populations
+    rho = f[0] + f[1] + f[2] + f[3] + f[4] + f[5] + f[6] + f[7] + f[8];
+#if USE_DDF_SHIFTING
+    rho += 1.0f;   // DDF-shifting: stored f_tilde = f - w, sum(w) = 1
+#endif
+
+    // Momentum (works identically with or without DDF-shifting,
+    // because sum(w_i * cx_i) = 0)
+    float inv_rho = 1.0f / rho;
+    //  cx:  0  1 -1  0  0  1 -1  1 -1
+    ux = (f[1] - f[2] + f[5] - f[6] + f[7] - f[8]) * inv_rho;
+    //  cy:  0  0  0  1 -1  1 -1 -1  1
+    uy = (f[3] - f[4] + f[5] - f[6] - f[7] + f[8]) * inv_rho;
+}
+
+#elif NQ == 19
+
+__device__ __forceinline__ void compute_rho_u(const float* __restrict__ f,
+                                              float& rho,
+                                              float& ux, float& uy, float& uz)
+{
+    rho = f[0];
+    for (int i = 1; i < 19; i++) rho += f[i];
+#if USE_DDF_SHIFTING
+    rho += 1.0f;
+#endif
+
+    float inv_rho = 1.0f / rho;
+    ux = (f[1]-f[2]  + f[7]-f[8]  + f[9]-f[10] + f[13]-f[14] + f[15]-f[16]) * inv_rho;
+    uy = (f[3]-f[4]  + f[7]-f[8]  + f[11]-f[12] + f[14]-f[13] + f[17]-f[18]) * inv_rho;
+    uz = (f[5]-f[6]  + f[9]-f[10] + f[11]-f[12] + f[16]-f[15] + f[18]-f[17]) * inv_rho;
+}
+
+#endif  // NQ
+
+// ---------------------------------------------------------------------------
+// compute_feq:  (rho, u) → feq[NQ]
+// ---------------------------------------------------------------------------
+#if NQ == 9
+
+__device__ __forceinline__ void compute_feq(float rho, float ux, float uy,
+                                            float* __restrict__ feq)
+{
+#if USE_DDF_SHIFTING
+    const float rhom1 = rho - 1.0f;          // perturbation density
+#endif
+    const float ux3 = 3.0f * ux;
+    const float uy3 = 3.0f * uy;
+    const float c3  = -1.5f * (ux * ux + uy * uy);  // = -|u|²/(2cs²)
+
+    // Helper: feq_i = w_i * rho * (1 + 3 ci·u + 4.5 (ci·u)² - 1.5 u²)
+    // With DDF-shifting: feq_tilde_i = feq_i - w_i
+    //   = w_i * (rho * (1 + ...) - 1) = w_i * (rho*(...) + rho - 1)
+
+    #define FEQ_CORE(cu)  (rho * fmaf(0.5f * 3.0f, (cu) * (cu), (cu) + c3))
+
+#if USE_DDF_SHIFTING
+    #define FEQ(w, cu)  ((w) * fmaf(rho, fmaf(0.5f, (cu)*(cu), (cu) + c3), rhom1))
+#else
+    #define FEQ(w, cu)  ((w) * rho * (1.0f + (cu) + 0.5f * (cu) * (cu) + c3))
+#endif
+
+    // i=0: rest
+    feq[0] = FEQ(W0_VAL, 0.0f);
+
+    // (1,2) ±x
+    feq[1] = FEQ(WS_VAL,  ux3);
+    feq[2] = FEQ(WS_VAL, -ux3);
+
+    // (3,4) ±y
+    feq[3] = FEQ(WS_VAL,  uy3);
+    feq[4] = FEQ(WS_VAL, -uy3);
+
+    // Pre-compute diagonal velocity combos
+    const float u0 = ux3 + uy3;    // +x+y
+    const float u1 = ux3 - uy3;    // +x-y
+
+    // (5,6) ±(x+y)
+    feq[5] = FEQ(WE_VAL,  u0);
+    feq[6] = FEQ(WE_VAL, -u0);
+
+    // (7,8) ±(x-y)
+    feq[7] = FEQ(WE_VAL,  u1);
+    feq[8] = FEQ(WE_VAL, -u1);
+
+    #undef FEQ
+    #undef FEQ_CORE
+}
+
+#elif NQ == 19
+
+__device__ __forceinline__ void compute_feq(float rho, float ux, float uy, float uz,
+                                            float* __restrict__ feq)
+{
+#if USE_DDF_SHIFTING
+    const float rhom1 = rho - 1.0f;
+#endif
+    const float ux3 = 3.0f * ux;
+    const float uy3 = 3.0f * uy;
+    const float uz3 = 3.0f * uz;
+    const float c3  = -1.5f * (ux*ux + uy*uy + uz*uz);
+
+#if USE_DDF_SHIFTING
+    #define FEQ(w, cu)  ((w) * fmaf(rho, fmaf(0.5f, (cu)*(cu), (cu) + c3), rhom1))
+#else
+    #define FEQ(w, cu)  ((w) * rho * (1.0f + (cu) + 0.5f * (cu) * (cu) + c3))
+#endif
+
+    feq[0]  = FEQ(W0_VAL, 0.0f);
+
+    // Straight
+    feq[1]  = FEQ(WS_VAL,  ux3);   feq[2]  = FEQ(WS_VAL, -ux3);
+    feq[3]  = FEQ(WS_VAL,  uy3);   feq[4]  = FEQ(WS_VAL, -uy3);
+    feq[5]  = FEQ(WS_VAL,  uz3);   feq[6]  = FEQ(WS_VAL, -uz3);
+
+    // Diagonal
+    const float u0 = ux3+uy3, u1 = ux3+uz3, u2 = uy3+uz3;
+    const float u3 = ux3-uy3, u4 = ux3-uz3, u5 = uy3-uz3;
+
+    feq[7]  = FEQ(WE_VAL,  u0);   feq[8]  = FEQ(WE_VAL, -u0);
+    feq[9]  = FEQ(WE_VAL,  u1);   feq[10] = FEQ(WE_VAL, -u1);
+    feq[11] = FEQ(WE_VAL,  u2);   feq[12] = FEQ(WE_VAL, -u2);
+    feq[13] = FEQ(WE_VAL,  u3);   feq[14] = FEQ(WE_VAL, -u3);
+    feq[15] = FEQ(WE_VAL,  u4);   feq[16] = FEQ(WE_VAL, -u4);
+    feq[17] = FEQ(WE_VAL,  u5);   feq[18] = FEQ(WE_VAL, -u5);
+
+    #undef FEQ
+}
+
+#endif  // NQ
+
+// ---------------------------------------------------------------------------
+// compute_pressure (diagnostic, from state equation p = cs² ρ)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ float compute_pressure(float rho) {
+    return CS2 * rho;
+}
+
+__device__ __forceinline__ float compute_pressure_perturbation(float rho, float rho_ref) {
+    return CS2 * (rho - rho_ref);
+}
+
+#endif // CELERIS_OPERATORS_MACRO_CUH

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

@@ -0,0 +1,332 @@
+// CelerisLab – step/one_step_double.cu
+// Main LBM step kernel using double-buffer PULL streaming.
+//
+// Workflow per node:
+//   1. Pull-stream: gather DDF from neighbors → local f[NQ]
+//   2. Compute macroscopic quantities (ρ, u)
+//   3. Apply boundary conditions (inlet/outlet/wall)
+//   4. Compute equilibrium + forcing
+//   5. Collision (SRT / TRT / MRT)
+//   6. Write to output buffer
+//
+// Kernel signature maintains backward compatibility with driver.py:
+//   flag*   – cell flags (uint8 for legacy, uint32 for new)
+//   fi_in*  – input DDF array  (fpxx)
+//   fi_out* – output DDF array (fpxx)
+//   indx*   – per-cell offset into delta pool
+//   delta*  – curved boundary parameter pool
+//   action* – per-object control input
+//   obs*    – per-object observation output (force / velocity)
+// ============================================================================
+
+#ifndef CELERIS_STEP_ONE_STEP_DOUBLE_CU
+#define CELERIS_STEP_ONE_STEP_DOUBLE_CU
+
+__global__ void StreamCollideDouble(
+    const uint8_t* __restrict__ flag,
+    const fpxx*    __restrict__ fi_in,
+    fpxx*          __restrict__ fi_out,
+    const int32_t* __restrict__ indx,
+    const float*   __restrict__ delta,
+    const float*   __restrict__ action,
+    float*         __restrict__ obs,
+    float*         __restrict__ rho_arr,    // macro output (can be NULL)
+    float*         __restrict__ u_arr)      // macro output (can be NULL)
+{
+    // ----- Thread → node mapping (2D) -----
+#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;
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+#endif
+
+    // ----- Read flag -----
+    uint8_t fl = flag[k];
+
+    // ----- Neighbor indices -----
+    unsigned long j[NQ];
+    compute_neighbors(k, j);
+
+    // ----- Pull-stream: load DDF from neighbors -----
+    float f[NQ];
+    stream_pull_load(k, f, fi_in, j);
+
+    // ----- Compute macroscopic quantities -----
+    float rho_n, ux, uy;
+#if NQ == 9
+    compute_rho_u(f, rho_n, ux, uy);
+#elif NQ == 19
+    float uz;
+    compute_rho_u(f, rho_n, ux, uy, uz);
+#endif
+
+    // ----- Boundary conditions (inlet/outlet/wall) -----
+#if NQ == 9
+    if (fl & LEGACY_SOLID) {
+        if (x == 0) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x + 1u, y);
+            for (int i = 0; i < NQ; i++) {
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            }
+            apply_parabolic_inlet(f, f_neb, (float)y);
+        }
+        else if (x == NX - 1) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x - 1u, y);
+            for (int i = 0; i < NQ; i++) {
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            }
+            apply_pressure_outlet(f, f_neb, (float)y);
+        }
+    }
+#elif NQ == 19
+    if (fl & LEGACY_SOLID) {
+        if (x == 0) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x + 1u, y, z);
+            for (int i = 0; i < NQ; i++) {
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            }
+            apply_parabolic_inlet_3d(f, f_neb, (float)y);
+        }
+        else if (x == NX - 1) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x - 1u, y, z);
+            for (int i = 0; i < NQ; i++) {
+                f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
+            }
+            apply_pressure_outlet_3d(f, f_neb, (float)y);
+        }
+    }
+#endif
+
+    // ----- Wall bounce-back (top/bottom walls) -----
+#if NQ == 9
+    if (y == 1 || y == NY - 2) {
+        apply_wall_bb_d2q9(y, f);
+    }
+#elif NQ == 19
+    if (y == 1 || y == NY - 2) {
+        apply_wall_bb_d3q19_y(y, f);
+    }
+#endif
+
+    // ----- Compute equilibrium -----
+    if (fl & LEGACY_FLUID) {
+        float feq[NQ];
+        float Fin[NQ];
+
+#if NQ == 9
+        compute_feq(rho_n, ux, uy, feq);
+        zero_forcing(Fin);
+
+        // ----- Collision -----
+#if COLLISION_MODEL == 0
+        collide_srt(f, feq, Fin, d_params.omega);
+#elif COLLISION_MODEL == 1
+        collide_trt(f, feq, Fin, d_params.omega);
+#elif COLLISION_MODEL == 2
+        collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
+#endif
+
+#elif NQ == 19
+        compute_feq(rho_n, ux, uy, uz, feq);
+        zero_forcing(Fin);
+        collide_srt(f, feq, Fin, d_params.omega);
+#endif
+    }
+
+    // ----- Write to output buffer -----
+    stream_pull_store(k, f, fi_out);
+
+    // ----- Write macroscopic fields (optional) -----
+    if (rho_arr != nullptr) {
+        rho_arr[k] = rho_n;
+    }
+    if (u_arr != nullptr) {
+#if DIM == 2
+        u_arr[k]               = ux;
+        u_arr[TOTAL_CELLS + k] = uy;
+#elif DIM == 3
+        u_arr[k]                   = ux;
+        u_arr[TOTAL_CELLS + k]     = uy;
+        u_arr[2 * TOTAL_CELLS + k] = uz;
+#endif
+    }
+
+    // ----- Sensor observation -----
+    if (fl & LEGACY_SENSOR) {
+        int id_obj = indx[k];
+        atomicAdd(&obs[DIM * id_obj],     ux);
+        atomicAdd(&obs[DIM * id_obj + 1], uy);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// Curved boundary post-processing kernel (separate launch)
+// Processes SOLID+INTERFACE nodes → modifies fi_out at neighbors.
+// ---------------------------------------------------------------------------
+__global__ void CurvedBoundaryPost(
+    const uint8_t* __restrict__ flag,
+    fpxx*          __restrict__ fi_out,
+    const int32_t* __restrict__ indx,
+    const float*   __restrict__ delta,
+    const float*   __restrict__ action,
+    float*         __restrict__ obs)
+{
+#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;
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+#endif
+
+    uint8_t fl = flag[k];
+
+    if ((fl & LEGACY_SOLID) && (fl & LEGACY_INTERFACE)) {
+#if DIM == 2
+        int id_off = indx[k];
+        int id_obj = *reinterpret_cast<const int*>(&delta[id_off]);
+
+        // Wall velocity from action + normal
+        float Uw = action[id_obj] * delta[id_off + NQ];
+        float Vw = action[id_obj] * delta[id_off + NQ + 1];
+
+        float* obs_fx = (obs != nullptr) ? &obs[DIM * id_obj]     : nullptr;
+        float* obs_fy = (obs != nullptr) ? &obs[DIM * id_obj + 1] : nullptr;
+
+        apply_curved_boundary(k, x, y, fi_out, delta, id_off, Uw, Vw, obs_fx, obs_fy);
+#endif // DIM == 2 curved BC
+    }
+}
+
+// ---------------------------------------------------------------------------
+// InitTubeFlow:  Initialize parabolic tube flow  (backward compatible)
+// ---------------------------------------------------------------------------
+__global__ void InitTubeFlow(uint8_t* flag, fpxx* fi)
+{
+#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;
+
+    float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1))
+                                             * ((float)y - 0.5f * (NY - 1))
+                                / ((float)(NY - 2) * (float)(NY - 2)));
+
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++) {
+            store_ddf(fi, index_f(k, (unsigned int)i), 0.0f);
+        }
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        // feq with rho=RHO, u=(u_init, 0)
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;  // cy=0 for parabolic flow
+            float val = d_w[i] * RHO * (1.0f + 3.0f * cu + 4.5f * cu * cu - 1.5f * u_init * u_init);
+#if USE_DDF_SHIFTING
+            val -= d_w[i];   // store f_tilde = f - w
+#endif
+            store_ddf(fi, index_f(k, (unsigned int)i), val);
+        }
+    }
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+
+    float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1))
+                                             * ((float)y - 0.5f * (NY - 1))
+                                / ((float)(NY - 2) * (float)(NY - 2)));
+
+    // Walls: y=0, y=NY-1, x=0, x=NX-1;  z is periodic
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++) {
+            store_ddf(fi, index_f(k, (unsigned int)i), 0.0f);
+        }
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;
+            float val = d_w[i] * RHO * (1.0f + 3.0f * cu + 4.5f * cu * cu - 1.5f * u_init * u_init);
+#if USE_DDF_SHIFTING
+            val -= d_w[i];
+#endif
+            store_ddf(fi, index_f(k, (unsigned int)i), val);
+        }
+    }
+#endif
+}
+
+// ---------------------------------------------------------------------------
+// UpdateMacro:  Recompute rho/u from DDF (for diagnostics/output)
+// ---------------------------------------------------------------------------
+__global__ void UpdateMacro(
+    const fpxx*    __restrict__ fi,
+    float*         __restrict__ rho_arr,
+    float*         __restrict__ u_arr,
+    const uint8_t* __restrict__ flag)
+{
+#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;
+
+    if (flag[k] & LEGACY_SOLID) return;
+
+    unsigned long j[NQ];
+    compute_neighbors(k, j);
+
+    float f[NQ];
+    // For double-buffer, just read from the current buffer (no streaming)
+    for (int i = 0; i < NQ; i++) {
+        f[i] = load_ddf(fi, index_f(k, (unsigned int)i));
+    }
+
+    float rho_n, ux, uy;
+    compute_rho_u(f, rho_n, ux, uy);
+
+    rho_arr[k]               = rho_n;
+    u_arr[k]                 = ux;
+    u_arr[TOTAL_CELLS + k]   = uy;
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+
+    if (flag[k] & LEGACY_SOLID) return;
+
+    float f[NQ];
+    for (int i = 0; i < NQ; i++) {
+        f[i] = load_ddf(fi, index_f(k, (unsigned int)i));
+    }
+
+    float rho_n, ux, uy, uz;
+    compute_rho_u(f, rho_n, ux, uy, uz);
+
+    rho_arr[k]                   = rho_n;
+    u_arr[k]                     = ux;
+    u_arr[TOTAL_CELLS + k]       = uy;
+    u_arr[2 * TOTAL_CELLS + k]   = uz;
+#endif
+}
+
+#endif // CELERIS_STEP_ONE_STEP_DOUBLE_CU

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

@@ -0,0 +1,198 @@
+// CelerisLab – step/one_step_esopull.cu
+// Main LBM step kernel using Esoteric-Pull single-buffer streaming.
+//
+// Workflow per node:
+//   1. load_f_esopull: read DDF (streaming second half)
+//   2. Compute macroscopic quantities (ρ, u)
+//   3. Boundary conditions
+//   4. Equilibrium + forcing + collision
+//   5. store_f_esopull: write DDF (streaming first half)
+//
+// Compared to double-buffer:
+//   - Uses HALF the memory (single fi array)
+//   - Timestep counter t is required
+//   - fi_alt pointer is unused (can be NULL)
+// ============================================================================
+
+#ifndef CELERIS_STEP_ONE_STEP_ESOPULL_CU
+#define CELERIS_STEP_ONE_STEP_ESOPULL_CU
+
+__global__ void StreamCollideEsoPull(
+    fpxx*          __restrict__ fi,          // single DDF buffer (read/write)
+    const uint8_t* __restrict__ flag,
+    const int32_t* __restrict__ indx,
+    const float*   __restrict__ delta,
+    const float*   __restrict__ action,
+    float*         __restrict__ obs,
+    float*         __restrict__ rho_arr,
+    float*         __restrict__ u_arr,
+    const float*   __restrict__ force_field, // Euler force field [DIM*N], or NULL
+    unsigned long                t)          // current timestep
+{
+    // ----- Thread → node mapping -----
+#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;
+#elif DIM == 3
+    unsigned int x, y, z;
+    unsigned long k;
+    index_from_thread(x, y, z, k);
+    if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
+#endif
+
+    uint8_t fl = flag[k];
+
+    // Skip pure solid / gas
+    if ((fl & LEGACY_SOLID) && !(fl & LEGACY_INTERFACE) && !(fl & LEGACY_SENSOR)) {
+        // For inlet/outlet solid nodes, we still process below
+        if (x != 0 && x != (unsigned int)(NX - 1)) return;
+    }
+
+    // ----- Neighbor indices -----
+    unsigned long j[NQ];
+    compute_neighbors(k, j);
+
+    // ----- Esoteric-Pull: load DDF (streaming second half) -----
+    float f[NQ];
+    load_f_esopull(k, f, fi, j, t);
+
+    // ----- Compute macroscopic quantities -----
+    float rho_n, ux, uy;
+#if NQ == 9
+    compute_rho_u(f, rho_n, ux, uy);
+#elif NQ == 19
+    float uz;
+    compute_rho_u(f, rho_n, ux, uy, uz);
+#endif
+
+    // ----- Boundary conditions -----
+#if NQ == 9
+    if (fl & LEGACY_SOLID) {
+        if (x == 0) {
+            float f_neb[NQ];
+            unsigned long k_neb = linear_index(x + 1u, y);
+            load_f_esopull(k_neb, f_neb, fi, j, t);  // approximate: reuse j
+            // More accurate: compute j_neb separately
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_parabolic_inlet(f, f_neb, (float)y);
+        }
+        else if (x == (unsigned int)(NX - 1)) {
+            unsigned long k_neb = linear_index(x - 1u, y);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_pressure_outlet(f, f_neb, (float)y);
+        }
+    }
+
+    if (y == 1 || y == (unsigned int)(NY - 2)) {
+        apply_wall_bb_d2q9(y, f);
+    }
+#endif
+
+    // ----- Forcing -----
+    float Fin[NQ];
+    float fxn = d_params.fx, fyn = d_params.fy;
+
+    if (force_field != nullptr) {
+        fxn += force_field[k];
+        fyn += force_field[TOTAL_CELLS + k];
+    }
+
+    if (fxn != 0.0f || fyn != 0.0f) {
+#if NQ == 9
+        apply_guo_velocity_correction(ux, uy, fxn, fyn, rho_n);
+        compute_guo_forcing(ux, uy, fxn, fyn, Fin);
+#endif
+    } else {
+        zero_forcing(Fin);
+    }
+
+    // ----- Clamp velocity for stability -----
+    const float CS_LIMIT = 0.57735027f;  // 1/sqrt(3)
+    ux = fminf(fmaxf(ux, -CS_LIMIT), CS_LIMIT);
+    uy = fminf(fmaxf(uy, -CS_LIMIT), CS_LIMIT);
+
+    // ----- Collision -----
+    if (fl & LEGACY_FLUID) {
+        float feq[NQ];
+#if NQ == 9
+        compute_feq(rho_n, ux, uy, feq);
+
+#if COLLISION_MODEL == 0
+        collide_srt(f, feq, Fin, d_params.omega);
+#elif COLLISION_MODEL == 1
+        collide_trt(f, feq, Fin, d_params.omega);
+#elif COLLISION_MODEL == 2
+        collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
+#endif
+
+#elif NQ == 19
+        float fzn = d_params.fz;
+        if (force_field != nullptr) fzn += force_field[2*TOTAL_CELLS + k];
+        apply_guo_velocity_correction(ux, uy, uz, fxn, fyn, fzn, rho_n);
+        compute_feq(rho_n, ux, uy, uz, feq);
+        collide_srt(f, feq, Fin, d_params.omega);
+#endif
+    }
+
+    // ----- Esoteric-Pull: write DDF (streaming first half) -----
+    store_f_esopull(k, f, fi, j, t);
+
+    // ----- Write macroscopic fields -----
+    if (rho_arr != nullptr) rho_arr[k] = rho_n;
+    if (u_arr != nullptr) {
+        u_arr[k]               = ux;
+        u_arr[TOTAL_CELLS + k] = uy;
+    }
+
+    // ----- Sensor observation -----
+    if (fl & LEGACY_SENSOR) {
+        int id_obj = indx[k];
+        atomicAdd(&obs[DIM * id_obj],     ux);
+        atomicAdd(&obs[DIM * id_obj + 1], uy);
+    }
+}
+
+// ---------------------------------------------------------------------------
+// InitializeEsoPull:  Initialize DDF via Esoteric-Pull store (t=1)
+// ---------------------------------------------------------------------------
+__global__ void InitializeEsoPull(uint8_t* flag, fpxx* fi)
+{
+#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(k, j);
+
+    float feq[NQ];
+    if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
+        flag[k] = LEGACY_SOLID;
+        for (int i = 0; i < NQ; i++) feq[i] = 0.0f;
+    } else {
+        flag[k] = (uint8_t)LEGACY_FLUID;
+        float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f*(NY-1))
+                                                  * ((float)y - 0.5f*(NY-1))
+                                    / ((float)(NY-2) * (float)(NY-2)));
+        for (int i = 0; i < NQ; i++) {
+            float cu = (float)d_cx[i] * u_init;
+            feq[i] = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init);
+#if USE_DDF_SHIFTING
+            feq[i] -= d_w[i];
+#endif
+        }
+    }
+    // Write via Esoteric-Pull store at t=1 (odd step) for correct alignment
+    store_f_esopull(k, feq, fi, j, 1ul);
+#endif
+}
+
+#endif // CELERIS_STEP_ONE_STEP_ESOPULL_CU

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

@@ -0,0 +1,79 @@
+// CelerisLab – streaming/esopull_single_buffer.cuh
+// Esoteric-Pull single-buffer streaming  (Geier & Schönherr, 2017).
+//
+// Uses ONE DDF array (half the memory of double-buffer).
+// Relies on paired direction ordering: (i, i+1) are opposites.
+//
+// Key idea: load_f and store_f alternate read/write positions
+// based on timestep parity, so they never conflict in a single step.
+//
+// Even step (t%2 == 0):
+//   Read:   f[i]   ← fi[n,    i+1]   (local, reverse slot)
+//           f[i+1] ← fi[j[i], i  ]   (neighbor, forward slot)
+//   Write:  fi[j[i], i+1] ← f[i]     (neighbor, reverse slot)
+//           fi[n,    i  ] ← f[i+1]   (local, forward slot)
+//
+// Odd step (t%2 == 1):
+//   Read:   f[i]   ← fi[n,    i  ]   (local, forward slot)
+//           f[i+1] ← fi[j[i], i+1]   (neighbor, reverse slot)
+//   Write:  fi[j[i], i+1] ← f[i]     (neighbor, forward slot → stored in reverse)
+//           fi[n,    i  ] ← f[i+1]   (local, reverse slot → stored in forward)
+//
+// Wait, let me use the exact FluidX3D convention:
+//   Odd step read:  f[i] ← fi[n, i],   f[i+1] ← fi[j[i], i+1]
+//   Odd step write: fi[j[i], i+1] ← f[i],  fi[n, i] ← f[i+1]
+// ============================================================================
+
+#ifndef CELERIS_STREAMING_ESOPULL_SINGLE_BUFFER_CUH
+#define CELERIS_STREAMING_ESOPULL_SINGLE_BUFFER_CUH
+
+// ---------------------------------------------------------------------------
+// load_f_esopull:  Esoteric-Pull read phase (streaming second half)
+// ---------------------------------------------------------------------------
+__device__ inline void load_f_esopull(unsigned long n,
+                                      float* __restrict__ f,
+                                      const fpxx* __restrict__ fi,
+                                      const unsigned long* __restrict__ j,
+                                      unsigned long t)
+{
+    // i=0: rest particle always from self
+    f[0] = load_ddf(fi, index_f(n, 0u));
+
+    for (unsigned int i = 1; i < NQ; i += 2) {
+        if (t & 1) {
+            // Odd step
+            f[i]     = load_ddf(fi, index_f(n,    i));
+            f[i + 1] = load_ddf(fi, index_f(j[i], i + 1));
+        } else {
+            // Even step
+            f[i]     = load_ddf(fi, index_f(n,    i + 1));
+            f[i + 1] = load_ddf(fi, index_f(j[i], i));
+        }
+    }
+}
+
+// ---------------------------------------------------------------------------
+// store_f_esopull:  Esoteric-Pull write phase (streaming first half)
+// ---------------------------------------------------------------------------
+__device__ inline void store_f_esopull(unsigned long n,
+                                       const float* __restrict__ f,
+                                       fpxx* __restrict__ fi,
+                                       const unsigned long* __restrict__ j,
+                                       unsigned long t)
+{
+    store_ddf(fi, index_f(n, 0u), f[0]);
+
+    for (unsigned int i = 1; i < NQ; i += 2) {
+        if (t & 1) {
+            // Odd step
+            store_ddf(fi, index_f(j[i], i + 1), f[i]);
+            store_ddf(fi, index_f(n,    i),     f[i + 1]);
+        } else {
+            // Even step
+            store_ddf(fi, index_f(j[i], i),     f[i]);
+            store_ddf(fi, index_f(n,    i + 1), f[i + 1]);
+        }
+    }
+}
+
+#endif // CELERIS_STREAMING_ESOPULL_SINGLE_BUFFER_CUH

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

@@ -0,0 +1,56 @@
+// CelerisLab – streaming/pull_double_buffer.cuh
+// Double-buffer PULL streaming.
+//   f_local[i] = f_in[ neighbor(n, -c_i), i ]
+//   f_out[n, i] = f_collided[i]   (after collision, written to second buffer)
+//
+// This replaces the old PUSH scheme where each thread writes to neighbors.
+// Pull streaming provides better memory coalescing on GPU.
+// ============================================================================
+
+#ifndef CELERIS_STREAMING_PULL_DOUBLE_BUFFER_CUH
+#define CELERIS_STREAMING_PULL_DOUBLE_BUFFER_CUH
+
+// ---------------------------------------------------------------------------
+// stream_pull_load:  gather DDF from neighboring nodes into local f[NQ]
+//   n    – current node
+//   f    – output: local distribution (compute precision, float)
+//   fi   – input:  source DDF array (storage precision, fpxx)
+//   j    – neighbor indices (j[i] = neighbor in direction i)
+//
+// Pull semantic:  f[i] at node n = fi[ j[opp(i)], i ]
+//   i.e. the population traveling in direction i at node n
+//   arrived from the node at n - c_i  = neighbor in direction opp(i)
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ void stream_pull_load(unsigned long n,
+                                                 float* __restrict__ f,
+                                                 const fpxx* __restrict__ fi,
+                                                 const unsigned long* __restrict__ j)
+{
+    // i=0 (rest): always from self
+    f[0] = load_ddf(fi, index_f(n, 0u));
+
+    // For paired ordering, opp(i) swaps within each pair:
+    //   opp(1)=2, opp(2)=1, opp(3)=4, opp(4)=3, ...
+    #pragma unroll
+    for (int i = 1; i < NQ; i += 2) {
+        // population i arrives from direction opp(i) = i+1
+        f[i]     = load_ddf(fi, index_f(j[i+1], (unsigned int)i));
+        // population i+1 arrives from direction opp(i+1) = i
+        f[i + 1] = load_ddf(fi, index_f(j[i],   (unsigned int)(i + 1)));
+    }
+}
+
+// ---------------------------------------------------------------------------
+// stream_pull_store:  write collided f[NQ] to the OUTPUT buffer at node n
+// ---------------------------------------------------------------------------
+__device__ __forceinline__ void stream_pull_store(unsigned long n,
+                                                  const float* __restrict__ f,
+                                                  fpxx* __restrict__ fi_out)
+{
+    #pragma unroll
+    for (int i = 0; i < NQ; i++) {
+        store_ddf(fi_out, index_f(n, (unsigned int)i), f[i]);
+    }
+}
+
+#endif // CELERIS_STREAMING_PULL_DOUBLE_BUFFER_CUH

tests/test_d2q9_visual.py

@@ -0,0 +1,362 @@
+#!/usr/bin/env python3
+"""
+D2Q9 Regression Test — Poiseuille Channel + Cylinder Flow
+==========================================================
+Uses the ORIGINAL kernel.cu with same grid / BCs.
+Produces matplotlib figures for visual validation.
+
+Usage:
+    python tests/test_d2q9_visual.py --device 2
+    python tests/test_d2q9_visual.py --device 2 --cylinder
+"""
+
+import sys, os, argparse, time
+
+# Ensure CelerisLab package is importable
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src'))
+
+import numpy as np
+import pycuda.driver as cuda
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+from matplotlib.colors import Normalize
+
+from CelerisLab.cuda import compiler
+from CelerisLab.common import preprocess as preproc
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Configuration  ━━━━━━━━━━━━━━━━━━━━━━━
+NX, NY = 1280, 512
+NQ = 9
+NT = 128
+DIM = 2
+VIS = 0.002
+U0 = 0.01
+RHO = 1.0
+TOTAL = NX * NY
+
+# Original direction vectors (const.h ordering)
+E = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]], dtype=np.int32)
+OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6], dtype=np.int32)
+
+FLUID_FLAG = 0b00000001
+SOLID_FLAG = 0b00000010
+INTERFACE_FLAG = 0b00001000
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Helpers  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
+
+def configure_macros(n_objs=0):
+    """Write macros.h to match original kernel settings."""
+    lines = compiler.read_lines(compiler.kernel_path("macros.h"))
+    defs = {
+        'MULT_GPU': 'False', 'NT': NT,
+        'X_1U': 128, 'Y_1U': 32, 'Z_1U': 1,
+        'LBtype': 'float',
+        'UX': 10, 'UY': 16, 'UZ': 1,
+        'NX': NX, 'NY': NY, 'NZ': 1,
+        'DIM': DIM, 'NQ': NQ,
+        'VIS': VIS, 'RHO': f'{RHO}', 'U0': U0,
+        'N_OBJS': n_objs,
+    }
+    for name, val in defs.items():
+        lines = compiler.modify_macro(lines, name, val)
+    compiler.write_lines(compiler.kernel_path("macros.h"), lines)
+
+
+def extract_fields(ddf_host):
+    """Compute rho, u, v from host DDF (original const.h ordering).
+
+    The original kernel uses DDF-shifting: stores f_shifted = f - w_i*RHO.
+    So sum(f_shifted) = rho - RHO  (~0 for incompressible flow),
+    and momentum = sum(e_x * f_shifted) works because sum(w_i * e_x) = 0.
+    """
+    f = ddf_host.reshape(NQ, NY, NX)
+    rho = np.sum(f, axis=0) + RHO          # un-shift: rho = sum(f_shifted) + RHO
+    u = (f[1] + f[5] + f[8] - f[3] - f[6] - f[7]) / RHO
+    v = (f[2] + f[5] + f[6] - f[4] - f[7] - f[8]) / RHO
+    return rho, u, v
+
+
+def analytical_poiseuille(y_arr):
+    """Analytical parabolic profile matching InitTubeFlow."""
+    yy = (y_arr - 0.5 * (NY - 1)) / (NY - 2.0)
+    return U0 * 1.5 * (1 - 4 * yy**2)
+
+
+def build_cylinder_data(cx, cy, radius):
+    """Replicate driver.py add_cylinder logic for flag / delta / indx."""
+    flag = np.ones(TOTAL, dtype=np.uint8)   # init all FLUID
+    indx = np.zeros(TOTAL, dtype=np.int32)
+    delta_list = []
+    index_offset = 0
+
+    # Build Poiseuille flag first (walls + solid borders)
+    for y in range(NY):
+        for x in range(NX):
+            k = x + y * NX
+            if y == 0 or y == NY - 1 or x == 0 or x == NX - 1:
+                flag[k] = SOLID_FLAG
+
+    # Add cylinder
+    for x in range(int(cx - radius) - 1, int(cx + radius) + 1):
+        for y in range(int(cy - radius) - 1, int(cy + radius) + 1):
+            if (x - cx)**2 + (y - cy)**2 < radius**2:
+                k = x + y * NX
+                flag[k] = SOLID_FLAG
+                dt = np.zeros(11, dtype=np.float32)
+                dt[0] = np.int32(0).view(np.float32)  # id_object = 0
+                has_interface = False
+                for i in range(NQ):
+                    xn = x + E[i][0]
+                    yn = y + E[i][1]
+                    if (xn - cx)**2 + (yn - cy)**2 >= radius**2:
+                        has_interface = True
+                        xi, yi = preproc.find_circle_intersection(
+                            x, y, xn, yn, cx, cy, radius)
+                        d_neb = np.sqrt((xi - xn)**2 + (yi - yn)**2)
+                        e_len = np.sqrt(E[i][0]**2 + E[i][1]**2)
+                        if e_len > 0:
+                            dt[i] = d_neb / e_len
+                if has_interface:
+                    flag[k] |= INTERFACE_FLAG
+                    dt[9] = (cy - y) / radius
+                    dt[10] = (x - cx) / radius
+                    indx[k] = index_offset
+                    delta_list.append(dt)
+                    index_offset += 11
+
+    delta = np.concatenate(delta_list) if delta_list else np.zeros(1, dtype=np.float32)
+    return flag, indx, delta
+
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Simulation  ━━━━━━━━━━━━━━━━━━━━━━━━━━
+
+def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host):
+    """Compile kernel, run LBM, return DDF on host."""
+    cuda.init()
+    dev = cuda.Device(device_id)
+    ctx = dev.make_context()
+    print(f"[GPU {device_id}] {dev.name()}")
+
+    try:
+        configure_macros(n_objs)
+        compiler.compile_kernel()
+        ptx_path = compiler.kernel_path("kernel.ptx")
+        mod = cuda.module_from_file(ptx_path)
+        step_fn = mod.get_function("OneStep")
+        init_fn = mod.get_function("InitTubeFlow")
+
+        nbytes_ddf = TOTAL * NQ * 4
+        ddf_gpu   = cuda.mem_alloc(nbytes_ddf)
+        temp_gpu  = cuda.mem_alloc(nbytes_ddf)
+        flag_gpu  = cuda.mem_alloc(flag_host.nbytes)
+        indx_gpu  = cuda.mem_alloc(indx_host.nbytes)
+        delta_gpu = cuda.mem_alloc(max(delta_host.nbytes, 4))
+
+        action_host = np.zeros(max(n_objs, 1), dtype=np.float32)
+        obs_host    = np.zeros(max(n_objs * DIM, 1), dtype=np.float32)
+        action_gpu = cuda.mem_alloc(action_host.nbytes)
+        obs_gpu    = cuda.mem_alloc(obs_host.nbytes)
+
+        cuda.memcpy_htod(action_gpu, action_host)
+        cuda.memcpy_htod(obs_gpu, obs_host)
+
+        block = (NT, 1, 1)
+        grid  = (NX // NT, NY, 1)
+
+        # Init Poiseuille
+        init_fn(flag_gpu, ddf_gpu, block=block, grid=grid)
+        ctx.synchronize()
+
+        # Overwrite flag / indx / delta for cylinder case
+        cuda.memcpy_htod(flag_gpu, flag_host)
+        cuda.memcpy_htod(indx_gpu, indx_host)
+        cuda.memcpy_htod(delta_gpu, delta_host)
+
+        # Step loop
+        t0 = time.time()
+        for i in range(n_steps):
+            step_fn(flag_gpu, ddf_gpu, temp_gpu, indx_gpu, delta_gpu,
+                    action_gpu, obs_gpu, block=block, grid=grid)
+            ddf_gpu, temp_gpu = temp_gpu, ddf_gpu
+        ctx.synchronize()
+        dt = time.time() - t0
+        mlups = TOTAL * n_steps / dt / 1e6
+        print(f"  {n_steps} steps in {dt:.2f}s  ({mlups:.1f} MLUPS)")
+
+        # Copy back
+        ddf = np.zeros(TOTAL * NQ, dtype=np.float32)
+        cuda.memcpy_dtoh(ddf, ddf_gpu)
+        flag_out = np.zeros(TOTAL, dtype=np.uint8)
+        cuda.memcpy_dtoh(flag_out, flag_gpu)
+
+        return ddf, flag_out
+    finally:
+        ctx.pop()
+
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Visualization  ━━━━━━━━━━━━━━━━━━━━━━━
+
+def plot_poiseuille(ddf, flag, out_path):
+    """3-panel figure: velocity mag, u(y) profile, pressure along centerline."""
+    rho, u, v = extract_fields(ddf)
+    vel_mag = np.sqrt(u**2 + v**2)
+
+    # Mask solid cells for display
+    mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool)
+    vel_mag_masked = np.ma.array(vel_mag, mask=mask)
+
+    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
+
+    # (a) Velocity magnitude heatmap
+    ax = axes[0]
+    im = ax.imshow(vel_mag_masked, origin='lower', aspect='auto',
+                   cmap='jet', extent=[0, NX, 0, NY])
+    plt.colorbar(im, ax=ax, label='|u|')
+    ax.set_title('Velocity magnitude')
+    ax.set_xlabel('x'); ax.set_ylabel('y')
+
+    # (b) u(y) profile at x = NX/2  vs. analytical
+    ax = axes[1]
+    x_mid = NX // 2
+    y_arr = np.arange(NY, dtype=float)
+    ax.plot(u[:, x_mid], y_arr, 'b-', lw=2, label='LBM')
+    ax.plot(analytical_poiseuille(y_arr), y_arr, 'r--', lw=1.5, label='Analytical')
+    ax.set_xlabel('u_x'); ax.set_ylabel('y')
+    ax.set_title(f'u(y) at x={x_mid}')
+    ax.legend()
+    ax.grid(True, alpha=0.3)
+
+    # (c) Pressure along centerline y = NY/2
+    ax = axes[2]
+    y_mid = NY // 2
+    p = rho / 3.0
+    ax.plot(np.arange(NX), p[y_mid, :], 'g-', lw=1.5)
+    ax.set_xlabel('x'); ax.set_ylabel('p = ρ/3')
+    ax.set_title(f'Pressure along centerline (y={y_mid})')
+    ax.grid(True, alpha=0.3)
+
+    fig.suptitle(f'D2Q9 Poiseuille – NX={NX}, NY={NY}, VIS={VIS}, U0={U0}', fontsize=13)
+    fig.tight_layout()
+    fig.savefig(out_path, dpi=150)
+    print(f"  Saved: {out_path}")
+    plt.close(fig)
+
+
+def plot_cylinder(ddf, flag, cx, cy, radius, out_path):
+    """3-panel figure: velocity magnitude (zoom), vorticity, streamlines."""
+    rho, u, v = extract_fields(ddf)
+    vel_mag = np.sqrt(u**2 + v**2)
+    mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool)
+
+    # Zoom window around cylinder
+    pad = int(radius * 8)
+    x0 = max(int(cx - pad), 0)
+    x1 = min(int(cx + pad * 2), NX)
+    y0 = max(int(cy - pad), 0)
+    y1 = min(int(cy + pad), NY)
+
+    fig, axes = plt.subplots(1, 3, figsize=(20, 6))
+
+    # (a) Velocity magnitude (zoomed)
+    ax = axes[0]
+    vm_z = np.ma.array(vel_mag[y0:y1, x0:x1], mask=mask[y0:y1, x0:x1])
+    im = ax.imshow(vm_z, origin='lower', aspect='equal',
+                   cmap='jet', extent=[x0, x1, y0, y1])
+    circ = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7)
+    ax.add_patch(circ)
+    plt.colorbar(im, ax=ax, label='|u|')
+    ax.set_title('Velocity magnitude')
+
+    # (b) Vorticity
+    ax = axes[1]
+    dvdx = np.gradient(v, axis=1)
+    dudy = np.gradient(u, axis=0)
+    omega = dvdx - dudy
+    om_z = np.ma.array(omega[y0:y1, x0:x1], mask=mask[y0:y1, x0:x1])
+    vmax = np.percentile(np.abs(omega[~mask]), 99)
+    im = ax.imshow(om_z, origin='lower', aspect='equal',
+                   cmap='RdBu_r', extent=[x0, x1, y0, y1],
+                   vmin=-vmax, vmax=vmax)
+    circ2 = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7)
+    ax.add_patch(circ2)
+    plt.colorbar(im, ax=ax, label='ω')
+    ax.set_title('Vorticity')
+
+    # (c) Streamlines
+    ax = axes[2]
+    X, Y = np.meshgrid(np.arange(x0, x1), np.arange(y0, y1))
+    u_z = u[y0:y1, x0:x1].copy()
+    v_z = v[y0:y1, x0:x1].copy()
+    u_z[mask[y0:y1, x0:x1]] = 0
+    v_z[mask[y0:y1, x0:x1]] = 0
+    speed = np.sqrt(u_z**2 + v_z**2)
+    ax.streamplot(X, Y, u_z, v_z, color=speed, cmap='jet',
+                  density=2.0, linewidth=0.8)
+    circ3 = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7)
+    ax.add_patch(circ3)
+    ax.set_xlim(x0, x1); ax.set_ylim(y0, y1)
+    ax.set_aspect('equal')
+    ax.set_title('Streamlines')
+
+    fig.suptitle(f'D2Q9 Cylinder Flow – Re_D={U0*1.5*2*radius/VIS:.0f}, D={2*radius}', fontsize=13)
+    fig.tight_layout()
+    fig.savefig(out_path, dpi=150)
+    print(f"  Saved: {out_path}")
+    plt.close(fig)
+
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Main  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
+
+def main():
+    parser = argparse.ArgumentParser(description='D2Q9 Regression Test')
+    parser.add_argument('--device', type=int, default=2,
+                        help='CUDA device ID (default: 2)')
+    parser.add_argument('--cylinder', action='store_true',
+                        help='Also run cylinder flow test')
+    parser.add_argument('--steps-pois', type=int, default=5000,
+                        help='Steps for Poiseuille (default: 5000)')
+    parser.add_argument('--steps-cyl', type=int, default=30000,
+                        help='Steps for cylinder (default: 30000)')
+    args = parser.parse_args()
+
+    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output')
+    os.makedirs(out_dir, exist_ok=True)
+
+    # ---- Test 1: Poiseuille ----
+    print("\n===== Test 1: Poiseuille Channel Flow =====")
+    flag_pois = np.ones(TOTAL, dtype=np.uint8)
+    indx_pois = np.zeros(TOTAL, dtype=np.int32)
+    delta_pois = np.zeros(1, dtype=np.float32)
+
+    ddf, flag = run_simulation(args.device, args.steps_pois, 0,
+                               flag_pois, indx_pois, delta_pois)
+    plot_poiseuille(ddf, flag, os.path.join(out_dir, 'poiseuille_d2q9.png'))
+
+    # Error metric
+    rho, u, v = extract_fields(ddf)
+    y_arr = np.arange(NY, dtype=float)
+    u_ana = analytical_poiseuille(y_arr)
+    x_mid = NX // 2
+    u_num = u[:, x_mid]
+    # Interior cells only (skip walls)
+    err = np.max(np.abs(u_num[2:-2] - u_ana[2:-2])) / np.max(np.abs(u_ana[2:-2]))
+    print(f"  L∞ relative error at x={x_mid}: {err:.2e}")
+
+    # ---- Test 2: Cylinder ----
+    if args.cylinder:
+        print("\n===== Test 2: Flow Around Cylinder =====")
+        cyl_cx, cyl_cy, cyl_r = 256.0, 256.0, 32.0
+        flag_cyl, indx_cyl, delta_cyl = build_cylinder_data(cyl_cx, cyl_cy, cyl_r)
+
+        ddf2, flag2 = run_simulation(args.device, args.steps_cyl, 1,
+                                     flag_cyl, indx_cyl, delta_cyl)
+        plot_cylinder(ddf2, flag2, cyl_cx, cyl_cy, cyl_r,
+                      os.path.join(out_dir, 'cylinder_d2q9.png'))
+
+    print("\nDone.")
+
+
+if __name__ == '__main__':
+    main()

tests/test_d3q19_cylinder.py

@@ -0,0 +1,327 @@
+#!/usr/bin/env python3
+"""
+D3Q19 SRT — Cylinder Wake Flow (Periodic Z)
+=============================================
+Tests the v2 modular kernel (kernel_v2.cu) in 3D.
+Cylinder axis along z, parabolic inlet, pressure outlet, no-slip y-walls.
+Produces cross-section visualizations at z=NZ/2.
+
+Usage:
+    python tests/test_d3q19_cylinder.py --device 3
+    python tests/test_d3q19_cylinder.py --device 3 --re 200 --steps 50000
+"""
+
+import sys, os, argparse, time, struct
+
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src'))
+
+import numpy as np
+import pycuda.driver as cuda
+
+import matplotlib
+matplotlib.use('Agg')
+import matplotlib.pyplot as plt
+
+from CelerisLab.cuda import compiler
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Configuration  ━━━━━━━━━━━━━━━━━━━━━━━
+# Reasonable 3D grid — fits in < 500 MB GPU memory
+NX, NY, NZ = 256, 128, 32
+NQ  = 19
+NT  = 128
+DIM = 3
+RHO = 1.0
+CYL_CX, CYL_CY = 64.0, 64.0   # Cylinder center (x,y)
+CYL_R = 12.0                   # Cylinder radius
+TOTAL = NX * NY * NZ
+
+# D3Q19 paired direction ordering (from descriptors.cuh)
+#   0:rest  (1,2)±x  (3,4)±y  (5,6)±z
+#   (7,8)±(x+y) (9,10)±(x+z) (11,12)±(y+z)
+#   (13,14)±(x-y) (15,16)±(x-z) (17,18)±(y-z)
+CX = np.array([0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0], dtype=np.int32)
+CY = np.array([0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1], dtype=np.int32)
+CZ = np.array([0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1], dtype=np.int32)
+W  = np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float32)
+
+FLUID_FLAG = 0x01
+SOLID_FLAG = 0x02
+OBSTACLE_FLAG = 0x04   # triggers half-way BB at adjacent fluid nodes
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Helpers  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
+
+def compute_vis_omega(Re, D, U0):
+    """Compute viscosity and omega from target Reynolds number."""
+    vis = U0 * D / Re
+    omega = 1.0 / (3.0 * vis + 0.5)
+    return vis, omega
+
+
+def configure_macros_3d(vis, u0, n_objs=0):
+    """Write macros.h for D3Q19 SRT."""
+    lines = compiler.read_lines(compiler.kernel_path("macros.h"))
+    defs = {
+        'MULT_GPU': 'False', 'NT': NT,
+        'X_1U': NX, 'Y_1U': NY, 'Z_1U': NZ,
+        'LBtype': 'float',
+        'UX': 1, 'UY': 1, 'UZ': 1,
+        'NX': NX, 'NY': NY, 'NZ': NZ,
+        'DIM': DIM, 'NQ': NQ,
+        'VIS': f'{vis:.10f}', 'RHO': f'{RHO}', 'U0': u0,
+        'N_OBJS': n_objs,
+        'COLLISION_MODEL': 0,     # SRT
+        'STREAMING_MODEL': 0,     # double-buffer
+        'STORE_PRECISION': 0,     # FP32
+        'USE_DDF_SHIFTING': 0,
+    }
+    for name, val in defs.items():
+        lines = compiler.modify_macro(lines, name, val)
+    compiler.write_lines(compiler.kernel_path("macros.h"), lines)
+
+
+def build_cylinder_3d(cx, cy, radius):
+    """Build flag array for 3D cylinder (axis along z, periodic z)."""
+    flag = np.ones(TOTAL, dtype=np.uint8) * FLUID_FLAG
+    for z in range(NZ):
+        for y in range(NY):
+            for x in range(NX):
+                k = z * NY * NX + y * NX + x
+                # Channel walls
+                if y == 0 or y == NY - 1 or x == 0 or x == NX - 1:
+                    flag[k] = SOLID_FLAG
+                # Cylinder body (obstacle — triggers BB at fluid neighbors)
+                elif (x - cx)**2 + (y - cy)**2 < radius**2:
+                    flag[k] = OBSTACLE_FLAG
+    return flag
+
+
+def extract_fields_3d(ddf_host, z_slice):
+    """Extract rho, u, v, w at a given z-slice from D3Q19 DDF (v2 ordering)."""
+    # DDF layout: f[i * TOTAL + k]  where k = z*NY*NX + y*NX + x
+    f = ddf_host.reshape(NQ, NZ, NY, NX)
+    fz = f[:, z_slice, :, :]  # shape (NQ, NY, NX)
+    rho = np.sum(fz, axis=0)
+    ux = np.zeros_like(rho)
+    uy = np.zeros_like(rho)
+    uz_field = np.zeros_like(rho)
+    for i in range(NQ):
+        ux += CX[i] * fz[i]
+        uy += CY[i] * fz[i]
+        uz_field += CZ[i] * fz[i]
+    ux /= rho
+    uy /= rho
+    uz_field /= rho
+    return rho, ux, uy, uz_field
+
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Simulation  ━━━━━━━━━━━━━━━━━━━━━━━━━━
+
+def run_d3q19(device_id, n_steps, vis, u0, flag_host):
+    """Compile v2 kernel, run D3Q19 SRT, return DDF."""
+    omega = 1.0 / (3.0 * vis + 0.5)
+    cuda.init()
+    dev = cuda.Device(device_id)
+    ctx = dev.make_context()
+    print(f"[GPU {device_id}] {dev.name()}")
+
+    try:
+        configure_macros_3d(vis, u0)
+        compiler.compile_kernel_v2()
+        ptx_path = compiler.kernel_path("kernel_v2.ptx")
+        mod = cuda.module_from_file(ptx_path)
+
+        # Get kernels (extern "C" entries from kernel_v2.cu)
+        init_fn = mod.get_function("InitTubeFlow_v2")
+        step_fn = mod.get_function("OneStep")
+
+        # Set d_params.omega via __constant__ memory
+        params_ptr, params_size = mod.get_global("d_params")
+        # LBMParams struct layout (see params.cuh):
+        # Nx(4) Ny(4) Nz(4) N(8) omega(4) omega_bulk(4) fx(4) fy(4) fz(4)
+        # rho_ref(4) u_inlet(4) n_objects(4)
+        # Pack: unsigned int Nx, Ny, Nz; unsigned long N; float omega, omega_bulk, fx, fy, fz, rho_ref, u_inlet; unsigned int n_objects
+        params_data = struct.pack('IIIQfffffffI',
+                                  NX, NY, NZ,
+                                  TOTAL,
+                                  omega, 0.0,      # omega, omega_bulk
+                                  0.0, 0.0, 0.0,   # fx, fy, fz
+                                  RHO, u0,          # rho_ref, u_inlet
+                                  0)                 # n_objects
+        # Pad to match struct size
+        if len(params_data) < params_size:
+            params_data += b'\x00' * (params_size - len(params_data))
+        cuda.memcpy_htod(params_ptr, params_data)
+
+        # Allocate
+        nbytes_ddf = TOTAL * NQ * 4
+        ddf_gpu  = cuda.mem_alloc(nbytes_ddf)
+        temp_gpu = cuda.mem_alloc(nbytes_ddf)
+        flag_gpu = cuda.mem_alloc(flag_host.nbytes)
+        indx_gpu = cuda.mem_alloc(TOTAL * 4)
+        delta_gpu = cuda.mem_alloc(4)
+        action_gpu = cuda.mem_alloc(4)
+        obs_gpu = cuda.mem_alloc(4)
+
+        # Dummy arrays
+        cuda.memset_d32(indx_gpu, 0, TOTAL)
+        cuda.memset_d32(delta_gpu, 0, 1)
+        cuda.memset_d32(action_gpu, 0, 1)
+        cuda.memset_d32(obs_gpu, 0, 1)
+
+        block = (NT, 1, 1)
+        grid  = (NX // NT, NY, NZ)
+
+        # Initialize parabolic flow
+        init_fn(flag_gpu, ddf_gpu, block=block, grid=grid)
+        ctx.synchronize()
+
+        # Overwrite flags with cylinder geometry
+        cuda.memcpy_htod(flag_gpu, flag_host)
+
+        # Step loop
+        print(f"  Running {n_steps} steps (NX={NX}, NY={NY}, NZ={NZ}, omega={omega:.4f})...")
+        t0 = time.time()
+        for i in range(n_steps):
+            step_fn(flag_gpu, ddf_gpu, temp_gpu, indx_gpu, delta_gpu,
+                    action_gpu, obs_gpu,
+                    block=block, grid=grid)
+            ddf_gpu, temp_gpu = temp_gpu, ddf_gpu
+            if (i + 1) % 5000 == 0:
+                ctx.synchronize()
+                elapsed = time.time() - t0
+                mlups = TOTAL * (i + 1) / elapsed / 1e6
+                print(f"    step {i+1}/{n_steps}  ({mlups:.1f} MLUPS)")
+        ctx.synchronize()
+        dt = time.time() - t0
+        mlups = TOTAL * n_steps / dt / 1e6
+        print(f"  Done: {dt:.1f}s, {mlups:.1f} MLUPS")
+
+        # Copy back
+        ddf = np.zeros(TOTAL * NQ, dtype=np.float32)
+        cuda.memcpy_dtoh(ddf, ddf_gpu)
+        return ddf
+
+    finally:
+        ctx.pop()
+
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Visualization  ━━━━━━━━━━━━━━━━━━━━━━━
+
+def plot_d3q19_cylinder(ddf, flag, Re, u0, out_path):
+    """4-panel figure at z=NZ/2: vel-mag, vorticity, streamlines, u(y) profile."""
+    z_mid = NZ // 2
+    rho, ux, uy, uz = extract_fields_3d(ddf, z_mid)
+    vel_mag = np.sqrt(ux**2 + uy**2 + uz**2)
+
+    mask2d = ((flag.reshape(NZ, NY, NX)[z_mid] & (SOLID_FLAG | OBSTACLE_FLAG)) != 0)
+    vel_masked = np.ma.array(vel_mag, mask=mask2d)
+
+    fig, axes = plt.subplots(2, 2, figsize=(16, 10))
+
+    # (a) Velocity magnitude
+    ax = axes[0, 0]
+    im = ax.imshow(vel_masked, origin='lower', aspect='auto',
+                   cmap='jet', extent=[0, NX, 0, NY])
+    circ = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7)
+    ax.add_patch(circ)
+    plt.colorbar(im, ax=ax, label='|u|')
+    ax.set_title(f'Velocity magnitude (z={z_mid})')
+    ax.set_xlabel('x'); ax.set_ylabel('y')
+
+    # (b) Vorticity  ω_z = ∂v/∂x − ∂u/∂y
+    ax = axes[0, 1]
+    dvdx = np.gradient(uy, axis=1)
+    dudy = np.gradient(ux, axis=0)
+    omega_z = dvdx - dudy
+    om_masked = np.ma.array(omega_z, mask=mask2d)
+    vmax = np.percentile(np.abs(omega_z[~mask2d]), 99) if np.any(~mask2d) else 1e-3
+    im = ax.imshow(om_masked, origin='lower', aspect='auto',
+                   cmap='RdBu_r', extent=[0, NX, 0, NY],
+                   vmin=-vmax, vmax=vmax)
+    circ2 = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7)
+    ax.add_patch(circ2)
+    plt.colorbar(im, ax=ax, label='ω_z')
+    ax.set_title('Vorticity ω_z')
+    ax.set_xlabel('x'); ax.set_ylabel('y')
+
+    # (c) Streamlines
+    ax = axes[1, 0]
+    X, Y = np.meshgrid(np.arange(NX), np.arange(NY))
+    ux_c = ux.copy(); ux_c[mask2d] = 0
+    uy_c = uy.copy(); uy_c[mask2d] = 0
+    speed = np.sqrt(ux_c**2 + uy_c**2)
+    ax.streamplot(X, Y, ux_c, uy_c, color=speed, cmap='jet',
+                  density=2.5, linewidth=0.7)
+    circ3 = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7)
+    ax.add_patch(circ3)
+    ax.set_xlim(0, NX); ax.set_ylim(0, NY)
+    ax.set_aspect('auto')
+    ax.set_title('Streamlines')
+    ax.set_xlabel('x'); ax.set_ylabel('y')
+
+    # (d) u_x(y) profiles at different x stations
+    ax = axes[1, 1]
+    y_arr = np.arange(NY)
+    # Analytical parabolic inlet
+    yy = (y_arr - 0.5 * (NY - 1)) / (NY - 2.0)
+    u_ana = u0 * 1.5 * (1 - 4 * yy**2)
+
+    x_stations = [NX // 8, NX // 4, NX // 2, 3 * NX // 4]
+    for xs in x_stations:
+        ax.plot(ux[:, xs], y_arr, label=f'x={xs}')
+    ax.plot(u_ana, y_arr, 'k--', lw=1.5, label='Analytical inlet')
+    ax.set_xlabel('u_x'); ax.set_ylabel('y')
+    ax.set_title('u_x(y) profiles')
+    ax.legend(fontsize=8)
+    ax.grid(True, alpha=0.3)
+
+    fig.suptitle(f'D3Q19 SRT Cylinder — Re={Re:.0f}, D={2*CYL_R:.0f}, '
+                 f'Grid={NX}×{NY}×{NZ}', fontsize=13)
+    fig.tight_layout()
+    fig.savefig(out_path, dpi=150)
+    print(f"  Saved: {out_path}")
+    plt.close(fig)
+
+
+# ━━━━━━━━━━━━━━━━━━━━━━━  Main  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
+
+def main():
+    parser = argparse.ArgumentParser(description='D3Q19 SRT Cylinder Flow')
+    parser.add_argument('--device', type=int, default=0,
+                        help='CUDA device ID (default: 0)')
+    parser.add_argument('--re', type=float, default=100.0,
+                        help='Reynolds number based on diameter (default: 100)')
+    parser.add_argument('--u0', type=float, default=0.04,
+                        help='Inlet characteristic velocity (default: 0.04)')
+    parser.add_argument('--steps', type=int, default=30000,
+                        help='Number of LBM steps (default: 30000)')
+    args = parser.parse_args()
+
+    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output')
+    os.makedirs(out_dir, exist_ok=True)
+
+    D = 2 * CYL_R
+    vis, omega = compute_vis_omega(args.re, D, args.u0)
+    print(f"\n===== D3Q19 SRT Cylinder Flow =====")
+    print(f"  Re = {args.re:.0f},  D = {D:.0f},  U0 = {args.u0}")
+    print(f"  ν  = {vis:.6f},  ω = {omega:.4f}")
+    if omega > 1.95:
+        print(f"  WARNING: omega={omega:.4f} is close to 2.0, stability may be poor.")
+        print(f"  Consider reducing U0 or Re.")
+
+    flag = build_cylinder_3d(CYL_CX, CYL_CY, CYL_R)
+    n_solid = np.sum(flag == SOLID_FLAG)
+    n_fluid = np.sum(flag == FLUID_FLAG)
+    print(f"  Grid: {NX}×{NY}×{NZ} = {TOTAL} cells  (fluid: {n_fluid}, solid: {n_solid})")
+    print(f"  Memory: ~{2 * TOTAL * NQ * 4 / 1e6:.0f} MB for double-buffer DDF")
+
+    ddf = run_d3q19(args.device, args.steps, vis, args.u0, flag)
+
+    plot_d3q19_cylinder(ddf, flag, args.re, args.u0,
+                        os.path.join(out_dir, f'cylinder_d3q19_re{int(args.re)}.png'))
+
+    print("\nDone.")
+
+
+if __name__ == '__main__':
+    main()