kernels: shared block-window to global coordinate helper

BlockStep (slot-reservoir wall rows) and BlockChannelStep computed the
global cell of a block-window cell with the same four lines.
block_global_xy (core/layout.cuh, 2D, __forceinline__) now holds them;
both kernels call it with their own domain width. Same integer
operations: no numerical, ABI or scheduling change. The slot-reservoir
kernel layout is unchanged (no split).

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

src/CelerisLab/lbm/kernels/blocks.cu

@@ -168,21 +168,15 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal
         stream_pull_load(k, f, in + base, j);
         if (cell_kind == KIND_ACTIVE) {
             if (d_block_wall.nx != 0) {
-                const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo;
+                unsigned int gx, gy;
-                const unsigned int nbx = d_block_wall.nx / bx;
+                block_global_xy(x, y, halo, d_block_wall.nx, gx, gy);
-                const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo;
-                const unsigned int gy = (blockIdx.z / nbx) * by + y - halo;
                 if (gy == 1)
                     reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(
                         gx, d_block_wall.nx, d_block_wall.slot_width, d_block_wall.amplitude));
                 else if (gy == d_block_wall.ny - 2) reconstruct_top_reservoir(f, RHO);
             }
 #if COLLISION_MODEL == 2
-            // Uniform MRT sequence (collide_dispatch, zero forcing); SRT path below unchanged.
+            collide_mrt_unforced(f);  // SRT path below unchanged.
-            float rho, ux, uy;
-            compute_rho_u(f, rho, ux, uy);
-            collide_mrt_no_force(f, rho, ux, uy,
-                fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
 #else
             float rho, ux, uy, feq[NQ];
             compute_rho_u(f, rho, ux, uy);

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

@@ -53,6 +53,18 @@ __device__ __forceinline__ void index_from_thread(unsigned int& x,
     k = (unsigned long)y * NX + x;
 }
 
+// Global cell of window cell (x, y) in block blockIdx.z of a row-major,
+// full-domain block layout domain_nx cells wide. NX/NY are the local strides.
+__device__ __forceinline__ void block_global_xy(unsigned int x, unsigned int y,
+                                                unsigned int halo, unsigned int domain_nx,
+                                                unsigned int& gx, unsigned int& gy)
+{
+    const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo;
+    const unsigned int nbx = domain_nx / bx;
+    gx = (blockIdx.z % nbx) * bx + x - halo;
+    gy = (blockIdx.z / nbx) * by + y - halo;
+}
+
 // ---------------------------------------------------------------------------
 // Neighbor list (periodic BCs via modular arithmetic)
 // ---------------------------------------------------------------------------

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

@@ -114,6 +114,16 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
     collide_mrt(g, rho, ux, uy, Fin, omega);
 }
 
+// Unforced MRT with the clamped device omega, as in collide_dispatch with zero
+// forcing. Shared by the blocked, channel-root and dense-patch collisions.
+__device__ __forceinline__ void collide_mrt_unforced(float* __restrict__ g)
+{
+    float rho, ux, uy;
+    compute_rho_u(g, rho, ux, uy);
+    collide_mrt_no_force(g, rho, ux, uy,
+        fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
+}
+
 #elif LATTICE_MODEL == LATTICE_D3Q19
 
 __device__ __forceinline__ void collide_mrt(float* __restrict__ g,

src/CelerisLab/lbm/kernels/refinement.cu

@@ -37,10 +37,7 @@ __device__ unsigned int incoming_counts[NQ * TOTAL_CELLS];
 
 __device__ __forceinline__ void patch_collide(float* f, unsigned long k) {
 #if COLLISION_MODEL == 2
-    float rho, ux, uy;
+    collide_mrt_unforced(f);
-    compute_rho_u(f, rho, ux, uy);
-    collide_mrt_no_force(f, rho, ux, uy,
-        fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
 #else
     float rho, ux, uy, feq[NQ];
     compute_rho_u(f, rho, ux, uy);

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

@@ -15,10 +15,8 @@ extern "C" __global__ void BlockChannelStep(
     unsigned int x, y; unsigned long k;
     index_from_thread(x, y, k);
     if (x < halo || y < halo || x >= NX - halo || y >= NY - halo) return;
-    const unsigned int bx = NX - 2 * halo, by = NY - 2 * halo;
+    unsigned int gx, gy;
-    const unsigned int nbx = d_block_channel.nx / bx;
+    block_global_xy(x, y, halo, d_block_channel.nx, gx, gy);
-    const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo;
-    const unsigned int gy = (blockIdx.z / nbx) * by + y - halo;
     const unsigned long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS;
     const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k;
     float f[NQ];
@@ -55,11 +53,7 @@ extern "C" __global__ void BlockChannelStep(
         }
         if (!wall && !outlet) {
 #if COLLISION_MODEL == 2
-            // Match uniform MRT arithmetic (collide_dispatch, zero forcing).
+            collide_mrt_unforced(f);
-            float rho, ux, uy;
-            compute_rho_u(f, rho, ux, uy);
-            collide_mrt_no_force(f, rho, ux, uy,
-                fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, d_params.omega)));
 #else
             // Match uniform SRT arithmetic, including its shifted equilibrium path.
             float rho, ux, uy, feq[NQ], forcing[NQ] = {};

tests/unit/test_annulus_ledger.py

@@ -1,101 +0,0 @@
-"""Synthetic checks of the annulus mass ledger. Host-only, stdlib; not solver evidence."""
-
-import random
-import sys
-from pathlib import Path
-import unittest
-
-sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "validation"))
-import annulus_ledger as al  # noqa: E402
-
-NX = NY = 24
-# Annulus 4 <= r <= 8 around (12, 12); dyadic values keep every sum exact.
-A = {(x, y) for x in range(NX) for y in range(NY) if 16 <= (x - 12) ** 2 + (y - 12) ** 2 <= 64}
-DELTA = 0.375
-
-
-def owner(seed=1):
-    rng = random.Random(seed)
-    return [[[rng.randint(1, 1024) / 1024 for _ in range(NX)] for _ in range(NY)] for _ in range(9)]
-
-
-def collide(P, source=None):
-    """Mass-conserving stand-in: rotate populations within each cell (+ optional source)."""
-    D = [[[P[(q + 1) % 9][y][x] for x in range(NX)] for y in range(NY)] for q in range(9)]
-    if source:
-        (x, y), s = source
-        D[0][y][x] += s
-    return D
-
-
-def link(interior):
-    """(x, y, q) with x in A and x - c_q inside (seam) or outside (edge) A."""
-    for (x, y) in sorted(A):
-        for q in range(1, 9):
-            src = ((x - al.CX[q]) % NX, (y - al.CY[q]) % NY)
-            if (src in A) == interior:
-                return x, y, q
-    raise AssertionError
-
-
-class Ledger(unittest.TestCase):
-    def setUp(self):
-        self.O = owner()
-        self.P = al.pull(self.O, NX, NY)
-
-    def run_ledger(self, P, D):
-        return al.ledger(self.O, P, D, A, NX, NY)
-
-    def test_ideal_pull_has_no_h_and_closes(self):
-        r = self.run_ledger(self.P, collide(self.P))
-        self.assertEqual((r["H_seam"], r["H_edge"], r["C"], r["residual"]), (0, 0, 0, 0))
-        self.assertEqual(r["dM"], -r["Q_exact"])
-        self.assertNotEqual(r["Q_exact"], 0)
-
-    def test_collision_source_is_c(self):
-        x, y = sorted(A)[3]
-        r = self.run_ledger(self.P, collide(self.P, ((x, y), DELTA)))
-        self.assertEqual((r["C"], r["H"], r["residual"]), (DELTA, 0, 0))
-
-    def test_interior_seam_perturbation(self):
-        x, y, q = link(interior=True)
-        self.P[q][y][x] += DELTA
-        r = self.run_ledger(self.P, collide(self.P))
-        self.assertEqual((r["H_seam"], r["H_edge"], r["residual"]), (DELTA, 0, 0))
-        # Invisible to any flux through the edge of A; only H explains the mass change.
-        self.assertEqual(al.staged_flux(self.O, self.P, A, NX, NY), r["Q_exact"])
-        self.assertEqual(r["dM"] + r["Q_exact"], DELTA)
-
-    def test_edge_perturbation_and_staged_flux_double_count(self):
-        x, y, q = link(interior=False)
-        self.P[q][y][x] += DELTA
-        r = self.run_ledger(self.P, collide(self.P))
-        self.assertEqual((r["H_seam"], r["H_edge"], r["residual"]), (0, DELTA, 0))
-        staged = al.staged_flux(self.O, self.P, A, NX, NY)
-        self.assertEqual(staged, r["Q_exact"] - DELTA)
-        # Staged flux + full H counts the staging change twice.
-        self.assertEqual(r["dM"] + staged - r["C"] - r["H"], -DELTA)
-        # Staged flux alone closes only without the edge part of H.
-        self.assertEqual(r["dM"] + staged - r["C"] - r["H_seam"], 0)
-
-    def test_owner_snapshot_after_staging_hides_the_change(self):
-        # Snapshot taken after staging wrote the same slot into the owner image.
-        x, y, q = link(interior=False)
-        sx, sy = (x - al.CX[q]) % NX, (y - al.CY[q]) % NY
-        self.P[q][y][x] += DELTA
-        late = [[row[:] for row in plane] for plane in self.O]
-        late[q][sy][sx] += DELTA
-        r = al.ledger(late, self.P, collide(self.P), A, NX, NY)
-        self.assertEqual((r["H"], r["residual"]), (0, 0))
-
-    def test_wrong_pull_direction_is_flagged(self):
-        # Indexing error (push instead of pull): H is nonzero with no physical cause.
-        bad = [[[self.O[q][(y + al.CY[q]) % NY][(x + al.CX[q]) % NX] for x in range(NX)]
-                for y in range(NY)] for q in range(9)]
-        r = self.run_ledger(bad, collide(bad))
-        self.assertNotEqual(r["H"], 0)
-        self.assertEqual(r["residual"], 0)
-
-
-if __name__ == "__main__":
-    unittest.main()

tests/unit/test_refinement_mrt.py

@@ -37,8 +37,13 @@ class RefinementMRTHost(unittest.TestCase):
             self.assertIn("STORE_PRECISION != 0", source)
             self.assertIn('#if COLLISION_MODEL == 2\n#include "operators/collision_mrt.cuh"\n#endif', source)
         self.assertIn("USE_DDF_SHIFTING ||", refinement)  # dense/mixed shifted stays rejected
+        # One shared unforced MRT body; each blocked/patch path calls it exactly once.
+        mrt = (KERNELS / "operators/collision_mrt.cuh").read_text()
+        self.assertEqual(mrt.count("void collide_mrt_unforced("), 1)
         for source in (blocks, channel, refinement):
-            self.assertEqual(source.count("collide_mrt_no_force("), 1)
+            self.assertEqual(source.count("collide_mrt_unforced(f);"), 1)
+            self.assertNotIn("collide_mrt_no_force(", source)
+            self.assertNotIn("compute_rho_u(f, rho, ux, uy);\n            collide_mrt", source)
 
 
 if __name__ == "__main__": unittest.main()

tests/validation/annulus_ledger.py

@@ -1,75 +0,0 @@
-"""Per-hop mass ledger of a fluid cell set A (host-only, stdlib, D2Q9 pull scheme).
-
-Diagnostic definition only; no solver code, no GPU. Arrays are indexed [q][y][x] on a
-periodic nx x ny lattice; A is a set of (x, y) ACTIVE cells.
-
-  O  owner populations before any staging (canonical owner slots, start of the hop)
-  P  populations each cell actually pulls: P[q](x) is what x reads for direction q
-     (after halo copy, BFL/coalesce staging, ...)
-  D  populations written by the cell after collision
-
-With M0 = sum_{A,q} O and M1 = sum_{A,q} D (both owner snapshots of A):
-
-  M1 - M0 + Q_exact = C + H                                            (exact)
-
-  C       = sum_{x in A, q} (D_q(x) - P_q(x))              collision source
-  H       = sum_{x in A, q} (P_q(x) - O_q(x - c_q))        actual minus ideal owner input,
-            split into H_seam (x - c_q in A: internal pulls, e.g. tile seams) and
-            H_edge (x - c_q not in A: inputs crossing the edge of A)
-  Q_exact = sum_{y in A, y+c_q not in A} O_q(y) - sum_{y not in A, y+c_q in A} O_q(y)
-            owner-based net outward flux (outgoing minus incoming)
-
-Proof: sum_{A,q} O_q(x - c_q) - M0 = -Q_exact (relabel y = x - c_q per direction), and
-sum_{A,q} D = C + H + sum_{A,q} O_q(x - c_q). Q_exact must use O, not staged values: a
-flux built from staged inputs already contains H_edge, so adding the full H counts it
-twice (see staged_flux). The edge of A must not cross cut links, otherwise O at the solid
-end enters Q_exact.
-"""
-
-CX = (0, 1, -1, 0, 0, 1, -1, 1, -1)
-CY = (0, 0, 0, 1, -1, 1, -1, -1, 1)
-
-
-def _src(x, y, q, nx, ny):
-    return (x - CX[q]) % nx, (y - CY[q]) % ny
-
-
-def pull(O, nx, ny):
-    """Ideal pulled input: P[q](x) = O[q](x - c_q)."""
-    return [[[O[q][_src(x, y, q, nx, ny)[1]][_src(x, y, q, nx, ny)[0]]
-              for x in range(nx)] for y in range(ny)] for q in range(9)]
-
-
-def ledger(O, P, D, A, nx, ny):
-    A = set(A)
-    out = dict(dM=0.0, C=0.0, H_seam=0.0, H_edge=0.0, Q_exact=0.0)
-    for (x, y) in A:
-        for q in range(9):
-            out["dM"] += D[q][y][x] - O[q][y][x]
-            out["C"] += D[q][y][x] - P[q][y][x]
-            sx, sy = _src(x, y, q, nx, ny)
-            key = "H_seam" if (sx, sy) in A else "H_edge"
-            out[key] += P[q][y][x] - O[q][sy][sx]
-            if (sx, sy) not in A:
-                out["Q_exact"] -= O[q][sy][sx]          # incoming owner population
-            tx, ty = (x + CX[q]) % nx, (y + CY[q]) % ny
-            if (tx, ty) not in A:
-                out["Q_exact"] += O[q][y][x]            # outgoing owner population
-    out["H"] = out["H_seam"] + out["H_edge"]
-    out["residual"] = out["dM"] + out["Q_exact"] - out["C"] - out["H"]
-    return out
-
-
-def staged_flux(O, P, A, nx, ny):
-    """Outward flux with incoming terms from the staged inputs P (the wrong combination)."""
-    A = set(A)
-    q_out = 0.0
-    for (x, y) in A:
-        for q in range(9):
-            sx, sy = _src(x, y, q, nx, ny)
-            if (sx, sy) not in A:
-                q_out -= P[q][y][x]
-            tx, ty = (x + CX[q]) % nx, (y + CY[q]) % ny
-            if (tx, ty) not in A:
-                q_out += O[q][y][x]
-    return q_out