From 781333407ed855bcbbcc5aacd460b02dae4e517d Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 01:15:47 +0000 Subject: [PATCH 1/2] kernels: one shared unforced MRT collision for blocked, channel and patch paths BlockStep, BlockChannelStep and refinement.cu patch_collide each carried the same compute_rho_u -> omega clamp -> collide_mrt_no_force sequence. collide_mrt_unforced (collision_mrt.cuh, __forceinline__) now holds it once; the three sites call it. Same operations in the same order: no numerical, ABI or SRT change. The host source contract now requires one helper definition and exactly one call per file, and forbids the old inline sequence. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/kernels/blocks.cu | 6 +----- src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh | 10 ++++++++++ src/CelerisLab/lbm/kernels/refinement.cu | 5 +---- src/CelerisLab/lbm/kernels/step/block_channel.cu | 6 +----- tests/unit/test_refinement_mrt.py | 7 ++++++- 5 files changed, 19 insertions(+), 15 deletions(-) diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 8e3764f..141d23a 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -178,11 +178,7 @@ extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int hal 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. - 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))); + collide_mrt_unforced(f); // SRT path below unchanged. #else float rho, ux, uy, feq[NQ]; compute_rho_u(f, rho, ux, uy); diff --git a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh index efc59c1..80d4637 100644 --- a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh +++ b/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, diff --git a/src/CelerisLab/lbm/kernels/refinement.cu b/src/CelerisLab/lbm/kernels/refinement.cu index 93c79b9..f496b10 100644 --- a/src/CelerisLab/lbm/kernels/refinement.cu +++ b/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; - 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))); + collide_mrt_unforced(f); #else float rho, ux, uy, feq[NQ]; compute_rho_u(f, rho, ux, uy); diff --git a/src/CelerisLab/lbm/kernels/step/block_channel.cu b/src/CelerisLab/lbm/kernels/step/block_channel.cu index 17e90fb..67a613d 100644 --- a/src/CelerisLab/lbm/kernels/step/block_channel.cu +++ b/src/CelerisLab/lbm/kernels/step/block_channel.cu @@ -55,11 +55,7 @@ extern "C" __global__ void BlockChannelStep( } if (!wall && !outlet) { #if COLLISION_MODEL == 2 - // Match uniform MRT arithmetic (collide_dispatch, zero forcing). - 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))); + collide_mrt_unforced(f); #else // Match uniform SRT arithmetic, including its shifted equilibrium path. float rho, ux, uy, feq[NQ], forcing[NQ] = {}; diff --git a/tests/unit/test_refinement_mrt.py b/tests/unit/test_refinement_mrt.py index 5f33782..be4d6f5 100644 --- a/tests/unit/test_refinement_mrt.py +++ b/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() From 1fa498c7d7ba8857c72c9740a0d6fed086e3f528 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 01:24:27 +0000 Subject: [PATCH 2/2] 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 --- src/CelerisLab/lbm/kernels/blocks.cu | 6 ++---- src/CelerisLab/lbm/kernels/core/layout.cuh | 12 ++++++++++++ src/CelerisLab/lbm/kernels/step/block_channel.cu | 6 ++---- 3 files changed, 16 insertions(+), 8 deletions(-) diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 141d23a..b8c7ef4 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -168,10 +168,8 @@ 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; - const unsigned int nbx = d_block_wall.nx / bx; - const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo; - const unsigned int gy = (blockIdx.z / nbx) * by + y - halo; + unsigned int gx, gy; + block_global_xy(x, y, halo, d_block_wall.nx, gx, gy); 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)); diff --git a/src/CelerisLab/lbm/kernels/core/layout.cuh b/src/CelerisLab/lbm/kernels/core/layout.cuh index 56a6d61..3b0f76a 100644 --- a/src/CelerisLab/lbm/kernels/core/layout.cuh +++ b/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) // --------------------------------------------------------------------------- diff --git a/src/CelerisLab/lbm/kernels/step/block_channel.cu b/src/CelerisLab/lbm/kernels/step/block_channel.cu index 67a613d..d31bd0d 100644 --- a/src/CelerisLab/lbm/kernels/step/block_channel.cu +++ b/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; - const unsigned int nbx = d_block_channel.nx / bx; - const unsigned int gx = (blockIdx.z % nbx) * bx + x - halo; - const unsigned int gy = (blockIdx.z / nbx) * by + y - halo; + unsigned int gx, gy; + block_global_xy(x, y, halo, d_block_channel.nx, gx, gy); 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];