Merge K4 test harness fix into K5

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

src/CelerisLab/lbm/kernels/blocks.cu

@@ -13,6 +13,7 @@
 #include "core/params.cuh"
 #include "core/flags.cuh"
 #include "core/obs.cuh"
+#include "core/cell_kind.cuh"
 #include "operators/macro.cuh"
 #include "operators/collision_srt.cuh"
 #if COLLISION_MODEL == 2
@@ -20,6 +21,10 @@
 #endif
 #include "streaming/pull_double_buffer.cuh"
 #include "boundary/curved_boundary.cuh"
+
+#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0
+#error "Blocks require FP32 D2Q9 SRT or MRT double buffers without LES"
+#endif
 #ifndef BLOCK_ROOT_BC
 #error "Blocks require a generated BLOCK_ROOT_BC mode"
 #elif BLOCK_ROOT_BC < 0 || BLOCK_ROOT_BC > 2
@@ -33,13 +38,6 @@
 #include "operators/block_collision.cuh"
 #endif
 
-#if DIM != 2 || NQ != 9 || (COLLISION_MODEL != 0 && COLLISION_MODEL != 2) || STORE_PRECISION != 0 || USE_LES || STREAMING_MODEL != 0
-#error "Blocks require FP32 D2Q9 SRT or MRT double buffers without LES"
-#endif
-
-// Fixed-patch cell kinds, matching refinement.py ACTIVE/INNER_GHOST/OUTER_GHOST.
-constexpr unsigned char KIND_ACTIVE = 0u, KIND_INNER = 1u, KIND_OUTER = 2u;
-
 // Halo refresh: block cell ids b*TOTAL_CELLS+c, sources are always core cells.
 extern "C" __global__ void BlockHaloCopy(
     fpxx* f, const uint64_t* dst, const uint64_t* src, unsigned int count)
@@ -47,8 +45,8 @@ extern "C" __global__ void BlockHaloCopy(
     unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
     if (row >= count) return;
     const uint64_t d = dst[row], s = src[row];
-    fpxx* to = f + (d / TOTAL_CELLS) * NQ * TOTAL_CELLS;
-    const fpxx* from = f + (s / TOTAL_CELLS) * NQ * TOTAL_CELLS;
+    fpxx* to = block_ddf(f, d);
+    const fpxx* from = block_ddf(f, s);
     for (unsigned int q = 0; q < NQ; ++q)
         to[index_f(d % TOTAL_CELLS, q)] = from[index_f(s % TOTAL_CELLS, q)];
 }
@@ -63,9 +61,9 @@ extern "C" __global__ void BlockExplode(
     unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
     if (row >= count) return;
     const uint64_t t = target[row];
-    fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS;
+    fpxx* to = block_ddf(fine, t);
     const uint64_t d = donor[row];
-    const fpxx* from = coarse + (d / coarse_cells) * NQ * coarse_cells;
+    const fpxx* from = block_ddf(coarse, d, coarse_cells);
     for (unsigned int q = 0; q < NQ; ++q)
         to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells];
 }
@@ -82,13 +80,12 @@ extern "C" __global__ void BlockExplodeProjected(
     unsigned int row = blockIdx.x * blockDim.x + threadIdx.x;
     if (row >= count) return;
     const uint64_t t = target[row], d = donor[row];
-    fpxx* to = fine + (t / TOTAL_CELLS) * NQ * TOTAL_CELLS;
+    fpxx* to = block_ddf(fine, t);
     for (unsigned int q = 0; q < NQ; ++q) {
         float v[4];
         for (unsigned int a = 0; a < 4; ++a) {
             const uint64_t n = neighbors[4ul * row + a];
-            v[a] = coarse[(n / coarse_cells) * NQ * coarse_cells
-                          + q * coarse_cells + n % coarse_cells];
+            v[a] = coarse[block_slot(n, coarse_cells, q)];
         }
         float gx = .5f * (v[1] - v[0]), gy = .5f * (v[3] - v[2]);
         const float ex = d_cx[q], ey = d_cy[q];
@@ -98,7 +95,7 @@ extern "C" __global__ void BlockExplodeProjected(
             gy -= ey * along;
         }
         to[index_f(t % TOTAL_CELLS, q)] =
-            coarse[(d / coarse_cells) * NQ * coarse_cells + q * coarse_cells + d % coarse_cells]
+            coarse[block_slot(d, coarse_cells, q)]
             + offsets[2ul * row] * gx + offsets[2ul * row + 1] * gy;
     }
 }
@@ -169,7 +166,7 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
     const unsigned char cell_kind = kind != nullptr ? kind[cell] : KIND_ACTIVE;
     float f[NQ];
     if (cell_kind == KIND_OUTER || (flags != nullptr && is_solid(flags[cell]))) {
-        for (int q = 0; q < NQ; ++q) f[q] = NAN;
+        fill_invalid(f);
     } else {
         unsigned long j[NQ];
         compute_neighbors(x, y, j);  // Core neighbors stay inside the window: no wrap.

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

@@ -15,10 +15,13 @@
 // D2Q9 y-wall half-way bounce-back for pull double-buffer streaming.
 // ---------------------------------------------------------------------------
 #if LATTICE_MODEL == LATTICE_D2Q9
-__device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y,
-                                                  float* __restrict__ f,
-                                                  const fpxx* __restrict__ fi_in,
-                                                  unsigned long k)
+// Rows 1 and top_row (= physical ny - 2) next to walls at rows 0 and ny - 1.
+// y is the physical row; k addresses the node in fi_in (a block window may
+// have NX/NY strides different from the physical extent).
+__device__ inline void apply_wall_bb_d2q9_rows(unsigned int y, unsigned int top_row,
+                                               float* __restrict__ f,
+                                               const fpxx* __restrict__ fi_in,
+                                               unsigned long k)
 {
     if (y == 1) {
         // Directions sourced from y=0 wall in pull step: +y, +x+y, -x+y.
@@ -26,14 +29,22 @@ __device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y,
         f[5] = load_ddf(fi_in, index_f(k, 6u));
         f[8] = load_ddf(fi_in, index_f(k, 7u));
     }
-    else if (y == (unsigned int)(NY - 2)) {
-        // Directions sourced from y=NY-1 wall in pull step: -y, -x-y, +x-y.
+    else if (y == top_row) {
+        // Directions sourced from the y = ny - 1 wall in pull step: -y, -x-y, +x-y.
         f[4] = load_ddf(fi_in, index_f(k, 3u));
         f[6] = load_ddf(fi_in, index_f(k, 5u));
         f[7] = load_ddf(fi_in, index_f(k, 8u));
     }
 }
 
+__device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y,
+                                                  float* __restrict__ f,
+                                                  const fpxx* __restrict__ fi_in,
+                                                  unsigned long k)
+{
+    apply_wall_bb_d2q9_rows(y, (unsigned int)(NY - 2), f, fi_in, k);
+}
+
 // ---------------------------------------------------------------------------
 // D2Q9 y-wall free-slip (specular reflection) for pull double-buffer streaming.
 //

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

@@ -0,0 +1,16 @@
+// CelerisLab – core/cell_kind.cuh
+// Fixed-patch cell kinds stored per block cell, matching lbm/refinement.py
+// ACTIVE/INNER_GHOST/OUTER_GHOST: OUTER is invalid (stores NaN, no pull),
+// INNER streams without collision, ACTIVE pulls and collides.
+#ifndef CELERIS_CORE_CELL_KIND_CUH
+#define CELERIS_CORE_CELL_KIND_CUH
+
+constexpr unsigned char KIND_ACTIVE = 0u, KIND_INNER = 1u, KIND_OUTER = 2u;
+
+// NaN populations mark a cell without valid state (OUTER ghost or solid core).
+__device__ __forceinline__ void fill_invalid(float* f)
+{
+    for (int q = 0; q < NQ; ++q) f[q] = NAN;
+}
+
+#endif

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

@@ -65,6 +65,21 @@ __device__ __forceinline__ void block_global_xy(unsigned int x, unsigned int y,
     gy = (blockIdx.z / nbx) * by + y - halo;
 }
 
+// Block-contiguous storage: block b holds NQ*cells values (Q-major inside) and
+// block cell id b*cells + c names cell c of block b. cells is TOTAL_CELLS for a
+// block window, or the whole domain for a dense (one-window) parent.
+template <typename T>
+__device__ __forceinline__ T* block_ddf(T* f, uint64_t block_cell, unsigned long cells = TOTAL_CELLS)
+{
+    return f + (block_cell / cells) * NQ * cells;
+}
+
+// Offset of population q of block cell block_cell from the start of storage.
+__device__ __forceinline__ uint64_t block_slot(uint64_t block_cell, unsigned long cells, unsigned int q)
+{
+    return (block_cell / cells) * NQ * cells + q * cells + block_cell % cells;
+}
+
 // ---------------------------------------------------------------------------
 // Neighbor list (periodic BCs via modular arithmetic)
 // ---------------------------------------------------------------------------

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

@@ -2,6 +2,13 @@
 #ifndef CELERIS_OPERATORS_BLOCK_COLLISION_CUH
 #define CELERIS_OPERATORS_BLOCK_COLLISION_CUH
 
+#include "../core/params.cuh"
+#include "macro.cuh"
+#include "collision_srt.cuh"
+#if COLLISION_MODEL == 2
+#include "collision_mrt.cuh"
+#endif
+
 __device__ __forceinline__ void block_collide(float* f) {
 #if COLLISION_MODEL == 2
     collide_mrt_unforced(f);  // SRT path below unchanged.

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

@@ -0,0 +1,14 @@
+// CelerisLab – operators/bounce_back_swap.cuh
+// Full-way bounce-back at a wall node: swap every opposite pair (paired
+// order q, q+1). Needs only NQ; shared by the uniform and blocked steps.
+#ifndef CELERIS_OPERATORS_BOUNCE_BACK_SWAP_CUH
+#define CELERIS_OPERATORS_BOUNCE_BACK_SWAP_CUH
+
+__device__ __forceinline__ void bounce_back_swap(float* f) {
+    #pragma unroll
+    for (int i = 1; i < NQ; i += 2) {
+        float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+    }
+}
+
+#endif

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

@@ -5,12 +5,7 @@
 #ifndef CELERIS_OPERATORS_HELPERS_CUH
 #define CELERIS_OPERATORS_HELPERS_CUH
 
-__device__ __forceinline__ void bounce_back_swap(float* f) {
-    #pragma unroll
-    for (int i = 1; i < NQ; i += 2) {
-        float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
-    }
-}
+#include "bounce_back_swap.cuh"
 
 __device__ __forceinline__ void collide_dispatch(
     float* f, float rho_n,

src/CelerisLab/lbm/kernels/refinement.cu

@@ -9,6 +9,7 @@
 #include "core/params.cuh"
 #include "core/flags.cuh"
 #include "core/obs.cuh"
+#include "core/cell_kind.cuh"
 #include "operators/macro.cuh"
 #include "operators/collision_srt.cuh"
 #if COLLISION_MODEL == 2
@@ -71,7 +72,7 @@ extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint
     // that held explosion data. Inner halo only streams; it never collides.
     if (x < 2 || y < 2 || x >= NX - 2 || y >= NY - 2
             || (flags != nullptr && is_solid(flags[k]))) {
-        for (int q = 0; q < NQ; ++q) f[q] = NAN;
+        fill_invalid(f);
     } else {
         unsigned long j[NQ];
         compute_neighbors(x, y, j);  // These interior targets never wrap.

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

@@ -4,6 +4,9 @@
 
 // Root-only blocked counterpart of OneStep's D2Q9/SRT channel closure.
 // NX/NY are LOCAL strides; all physical boundary decisions use global coordinates.
+#include "../core/cell_kind.cuh"
+#include "../operators/bounce_back_swap.cuh"
+#include "../boundary/bounce_back.cuh"
 #include "../boundary/inlet/common.cuh"
 #include "../boundary/inlet/zou_he_local.cuh"
 #include "../boundary/outlet/pressure_neq.cuh"
@@ -25,7 +28,7 @@ extern "C" __global__ void BlockChannelStep(
     const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k;
     float f[NQ];
     if (kind != nullptr && kind[cell] == KIND_OUTER) {
-        for (int q = 0; q < NQ; ++q) f[q] = NAN;
+        fill_invalid(f);
     } else {
         unsigned long j[NQ];
         compute_neighbors(x, y, j);
@@ -33,27 +36,17 @@ extern "C" __global__ void BlockChannelStep(
         const bool wall = gy == 0 || gy == d_block_channel.ny - 1;
         const bool outlet = gx == d_block_channel.nx - 1;
         if (wall) {
-            // Wall owns all four corners. Ghosts pull and swap, never collide.
-            for (int q = 1; q < NQ; q += 2) {
-                const float t = f[q]; f[q] = f[q + 1]; f[q + 1] = t;
-            }
+            bounce_back_swap(f);  // Wall owns all four corners; ghosts never collide.
         } else if (gx == 0) {
             apply_zou_he_left_velocity_inlet_d2q9(
                 f, inlet_target_u((float)gy, d_block_channel.ny, d_block_channel.velocity), 0.0f);
         } else if (outlet) {
             // Immutable input donor, receiver-local (also valid across a block seam).
             float donor[NQ];
-            for (int q = 0; q < NQ; ++q)
-                donor[q] = load_ddf(in + base, index_f(k - 1, q));
+            load_node(in + base, k - 1, donor);
             apply_pressure_outlet_d2q9(f, donor);
-        } else if (gy == 1) {
-            f[3] = load_ddf(in + base, index_f(k, 4));
-            f[5] = load_ddf(in + base, index_f(k, 6));
-            f[8] = load_ddf(in + base, index_f(k, 7));
-        } else if (gy == d_block_channel.ny - 2) {
-            f[4] = load_ddf(in + base, index_f(k, 3));
-            f[6] = load_ddf(in + base, index_f(k, 5));
-            f[7] = load_ddf(in + base, index_f(k, 8));
+        } else {
+            apply_wall_bb_d2q9_rows(gy, d_block_channel.ny - 2, f, in + base, k);
         }
         if (!wall && !outlet) {
 #if COLLISION_MODEL == 2

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

@@ -21,6 +21,14 @@
 //   i.e. the population traveling in direction i at node n
 //   arrived from the node at n - c_i  = neighbor in direction opp(i)
 // ---------------------------------------------------------------------------
+// load_node: all NQ populations of node n as stored (no streaming).
+__device__ __forceinline__ void load_node(const fpxx* __restrict__ fi, unsigned long n,
+                                          float* __restrict__ f)
+{
+    for (int q = 0; q < NQ; ++q)
+        f[q] = load_ddf(fi, index_f(n, q));
+}
+
 __device__ __forceinline__ void stream_pull_load(unsigned long n,
                                                  float* __restrict__ f,
                                                  const fpxx* __restrict__ fi,

tests/unit/cpp/channel_root.cpp

@@ -41,7 +41,7 @@ inline float __fdiv_rn(float a,float b){return a/b;}
 #include "boundary/inlet_outlet.cuh"
 #include "step/one_step_double.cu"
 #else
-constexpr unsigned char KIND_OUTER=2;
+#include "core/cell_kind.cuh"
 #include "step/block_channel.cu"
 #endif
 int main(){

tests/unit/test_projected_explosion.py

@@ -163,6 +163,12 @@ class ProjectedKernelTests(unittest.TestCase):
         donors = parent.owner.ravel()[p.ghost_coarse]
         source = (Path(__file__).resolve().parents[2] /
                   "src/CelerisLab/lbm/kernels/blocks.cu").read_text()
+        # The kernel addresses storage through the real core/layout.cuh helpers.
+        layout = (Path(__file__).resolve().parents[2] /
+                  "src/CelerisLab/lbm/kernels/core/layout.cuh").read_text()
+        start = layout.index("template <typename T>\n__device__ __forceinline__ T* block_ddf(")
+        end = layout.index("}\n", layout.index("__device__ __forceinline__ uint64_t block_slot(")) + 2
+        helpers = layout[start:end].replace("__device__ __forceinline__ ", "inline ")
         kernel = source.split('extern "C" __global__ void BlockExplodeProjected(',1)[1]
         kernel = 'void BlockExplodeProjected(' + kernel.split('// Targets are receiver-local SOURCE slots',1)[0]
         def array(name, values, kind):
@@ -183,7 +189,7 @@ class ProjectedKernelTests(unittest.TestCase):
             code += f'constexpr unsigned NQ=9, TOTAL_CELLS={int(np.prod(child.local_shape))};\n'
             code += 'struct Id {unsigned x;}; Id blockIdx{0},blockDim{1},threadIdx{0};\n'
             code += 'unsigned index_f(unsigned k,unsigned q){return q*TOTAL_CELLS+k;}\n'
-            code += 'int d_cx[]={0,1,-1,0,0,1,-1,1,-1},d_cy[]={0,0,0,1,-1,1,-1,-1,1};\n'+kernel
+            code += 'int d_cx[]={0,1,-1,0,0,1,-1,1,-1},d_cy[]={0,0,0,1,-1,1,-1,-1,1};\n'+helpers+kernel
             code += data+array('target',targets,'uint64_t')+array('donor',donors,'uint64_t')
             code += array('neighbor',ids,'uint64_t')+array('offset',offsets,'float')
             code += f'float fine[{child.nblocks*9*int(np.prod(child.local_shape))}]={{}};\n'