kernels (K4): blocked storage address helpers and one invalid-cell fill

- core/layout.cuh: block_ddf(f, block_cell[, cells]) returns the start of a
  block's populations and block_slot(block_cell, cells, q) the offset of one
  population; BlockHaloCopy, BlockExplode and BlockExplodeProjected use them
  instead of repeating (c / cells) * NQ * cells [+ q * cells + c % cells].
  BlockExplode keeps its loop-hoisted parent base; types stay uint64_t /
  unsigned long; the expressions are the same ones, now named.
- core/cell_kind.cuh: fill_invalid(f) writes the NaN marker; BlockStep,
  BlockChannelStep and PatchFineStep call it instead of three copies of the
  loop.
No change to operation order, types, launches or arithmetic.

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

src/CelerisLab/lbm/kernels/blocks.cu

@@ -45,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)];
 }
@@ -61,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];
 }
@@ -80,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];
@@ -96,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;
     }
 }
@@ -167,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/core/cell_kind.cuh

@@ -7,4 +7,10 @@
 
 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/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

@@ -26,7 +26,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);