From bf5088f84c01c2994066504290be936b85e9658b Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 22 Oct 2026 01:46:49 +0000 Subject: [PATCH] kernels (K9): blocks.cu becomes an include-only entry; kernels under step/ blocks.cu keeps its header, includes and #error guards and then includes the kernels, as kernel_v2.cu does: - step/block_transfer.cu: BlockHaloCopy, BlockExplode, BlockExplodeProjected, BlockCoalesce - step/block_boundary.cu: BlockCurvedBoundary, BlockAccumulateExchange - step/block_step.cu: BlockStep (rejects BLOCK_ROOT_BC == 2) - step/block_channel.cu: unchanged (mode 2 only) Kernel text moved verbatim (checked against the previous blocks.cu) and in the same definition order, so the translation unit is the same sequence of kernels. The two source-reading tests follow the move: test_projected_explosion reads step/block_transfer.cu, test_refinement_mrt reads blocks.cu + step/block_step.cu; their assertions are unchanged. Co-Authored-By: Claude Opus 5.5 --- src/CelerisLab/lbm/kernels/blocks.cu | 153 +----------------- .../lbm/kernels/step/block_boundary.cu | 42 +++++ src/CelerisLab/lbm/kernels/step/block_step.cu | 35 ++++ .../lbm/kernels/step/block_transfer.cu | 78 +++++++++ tests/unit/test_projected_explosion.py | 2 +- tests/unit/test_refinement_mrt.py | 2 +- 6 files changed, 163 insertions(+), 149 deletions(-) create mode 100644 src/CelerisLab/lbm/kernels/step/block_boundary.cu create mode 100644 src/CelerisLab/lbm/kernels/step/block_step.cu create mode 100644 src/CelerisLab/lbm/kernels/step/block_transfer.cu diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 44fc8c2..65d1320 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/src/CelerisLab/lbm/kernels/blocks.cu @@ -38,153 +38,12 @@ #include "operators/block_collision.cuh" #endif -// 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) -{ - 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; - for (unsigned int q = 0; q < NQ; ++q) - to[index_f(d % TOTAL_CELLS, q)] = from[index_f(s % TOTAL_CELLS, q)]; -} - -// Fixed-patch explosion into blocked fine storage: target = owner block cell -// id of a patch-ghost cell; donors are parent block cell ids. Dense parent -// storage is the one-window case with coarse_cells = full domain size. -extern "C" __global__ void BlockExplode( - const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor, - unsigned int count, unsigned long coarse_cells) -{ - 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; - const uint64_t d = donor[row]; - const fpxx* from = coarse + (d / coarse_cells) * NQ * coarse_cells; - for (unsigned int q = 0; q < NQ; ++q) - to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells]; -} - -// Experimental Schornbaum--Ruede (arXiv:1508.07982, Algorithm 2) -// and Chen et al. (doi:10.1016/j.physa.2005.09.036) population interpolation. -// Independent implementation of the mathematical stencil; no upstream source copied. -// Stencil cells are owner cores; a disabled axis repeats its center twice. -extern "C" __global__ void BlockExplodeProjected( - const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor, - unsigned int count, unsigned long coarse_cells, - const uint64_t* neighbors, const float* offsets) -{ - 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; - 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]; - } - float gx = .5f * (v[1] - v[0]), gy = .5f * (v[3] - v[2]); - const float ex = d_cx[q], ey = d_cy[q]; - if (q != 0) { - const float along = (ex * gx + ey * gy) / (ex * ex + ey * ey); - gx -= ex * along; - gy -= ey * along; - } - to[index_f(t % TOTAL_CELLS, q)] = - coarse[(d / coarse_cells) * NQ * coarse_cells + q * coarse_cells + d % coarse_cells] - + offsets[2ul * row] * gx + offsets[2ul * row + 1] * gy; - } -} - -// Targets are receiver-local SOURCE slots for the immediately following pull. -// Keep child summation order identical to PatchCoarseStep; no collision here. -extern "C" __global__ void BlockCoalesce( - const fpxx* fine, fpxx* coarse, const uint64_t* target, const uint64_t* children, - unsigned int count) -{ - unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; - if (row >= count) return; - const uint64_t* s = children + 4ul * row; - coarse[target[row]] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); -} - -// Static cut links, receiver-local staging: cell is the receiver's local core -// cell, so the existing BFL writes the solid-source window image that this -// block's next pull reads (no halo copy between this kernel and BlockStep). -// surface_velocity (nullable, 2 floats per link): prescribed fluid velocity on -// fixed geometry, same convention as CurvedBoundaryKernel. Optional raw rows -// (jx,jy,torque,dm) are unique physical links, never body-observation slots. -extern "C" __global__ void BlockCurvedBoundary( - fpxx* f, const uint32_t* block, const uint32_t* cell, const unsigned char* dir, - const float* q, const unsigned char* fallback, unsigned int count, - const float* surface_velocity, const float* rx, const float* ry, float* raw) -{ - unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; - if (row >= count) return; - float Uw = 0.0f, Vw = 0.0f; - if (surface_velocity != nullptr) { - const float u = surface_velocity[2 * row], v = surface_velocity[2 * row + 1]; - if (u != 0.0f || v != 0.0f) { Uw = u; Vw = v; } - } - apply_bouzidi_link(dir[row], q[row], cell[row], f + (unsigned long)block[row] * NQ * TOTAL_CELLS, - Uw, Vw, rx ? rx[row] : 0.0f, ry ? ry[row] : 0.0f, - fallback[row], nullptr, -1, raw ? raw + 4ul * row : nullptr, - raw ? raw + 4ul * row + 3 : nullptr); -} - -// Same-stream BFL -> this row-owned accumulation -> pull step. No atomics, -// no cross-row reduction and no history allocation proportional to substeps. -extern "C" __global__ void BlockAccumulateExchange( - const float* raw, double* integral, unsigned int count, double h, double dt) -{ - unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; - if (row >= count) return; - const unsigned long k = 4ul * row; - const double impulse = h * dt; - integral[k] += impulse * (double)raw[k]; - integral[k + 1] += impulse * (double)raw[k + 1]; - integral[k + 2] += h * impulse * (double)raw[k + 2]; - integral[k + 3] += h * h * (double)raw[k + 3]; -} - +// Kernels, in definition order: transfers, boundary staging, then the one full +// step selected for this module (BlockStep, or BlockChannelStep for a channel root). +#include "step/block_transfer.cu" +#include "step/block_boundary.cu" #if BLOCK_ROOT_BC != 2 -// Pull/collide/store core cells only; halo cells of out are left untouched. -// flags (block cell ids, nullable): solid cores store NaN and never collide. -// kind (block cell ids, nullable = all ACTIVE): OUTER skips the pull and -// stores NaN, INNER pulls without collision, ACTIVE pulls and collides. -extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int halo, - const uint16_t* flags, const unsigned char* kind) { - 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 long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS; - const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k; - 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]))) { - fill_invalid(f); - } else { - unsigned long j[NQ]; - compute_neighbors(x, y, j); // Core neighbors stay inside the window: no wrap. - stream_pull_load(k, f, in + base, j); - if (cell_kind == KIND_ACTIVE) { -#if BLOCK_ROOT_BC == 1 - block_wall(f, x, y, halo); -#endif - block_collide(f); - } - } - stream_pull_store(k, f, out + base); -} - -#endif - -// Root channel is a replacement full step, selected once by the host. -#if BLOCK_ROOT_BC == 2 +#include "step/block_step.cu" +#else #include "step/block_channel.cu" #endif diff --git a/src/CelerisLab/lbm/kernels/step/block_boundary.cu b/src/CelerisLab/lbm/kernels/step/block_boundary.cu new file mode 100644 index 0000000..9b89451 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/block_boundary.cu @@ -0,0 +1,42 @@ +// CelerisLab/lbm/kernels/step/block_boundary.cu +// Receiver-local curved-boundary staging and its raw-exchange accumulation. +// Part of the blocks.cu compilation entry; relies on its includes and guards. + +// Static cut links, receiver-local staging: cell is the receiver's local core +// cell, so the existing BFL writes the solid-source window image that this +// block's next pull reads (no halo copy between this kernel and BlockStep). +// surface_velocity (nullable, 2 floats per link): prescribed fluid velocity on +// fixed geometry, same convention as CurvedBoundaryKernel. Optional raw rows +// (jx,jy,torque,dm) are unique physical links, never body-observation slots. +extern "C" __global__ void BlockCurvedBoundary( + fpxx* f, const uint32_t* block, const uint32_t* cell, const unsigned char* dir, + const float* q, const unsigned char* fallback, unsigned int count, + const float* surface_velocity, const float* rx, const float* ry, float* raw) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + float Uw = 0.0f, Vw = 0.0f; + if (surface_velocity != nullptr) { + const float u = surface_velocity[2 * row], v = surface_velocity[2 * row + 1]; + if (u != 0.0f || v != 0.0f) { Uw = u; Vw = v; } + } + apply_bouzidi_link(dir[row], q[row], cell[row], f + (unsigned long)block[row] * NQ * TOTAL_CELLS, + Uw, Vw, rx ? rx[row] : 0.0f, ry ? ry[row] : 0.0f, + fallback[row], nullptr, -1, raw ? raw + 4ul * row : nullptr, + raw ? raw + 4ul * row + 3 : nullptr); +} + +// Same-stream BFL -> this row-owned accumulation -> pull step. No atomics, +// no cross-row reduction and no history allocation proportional to substeps. +extern "C" __global__ void BlockAccumulateExchange( + const float* raw, double* integral, unsigned int count, double h, double dt) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + const unsigned long k = 4ul * row; + const double impulse = h * dt; + integral[k] += impulse * (double)raw[k]; + integral[k + 1] += impulse * (double)raw[k + 1]; + integral[k + 2] += h * impulse * (double)raw[k + 2]; + integral[k + 3] += h * h * (double)raw[k + 3]; +} diff --git a/src/CelerisLab/lbm/kernels/step/block_step.cu b/src/CelerisLab/lbm/kernels/step/block_step.cu new file mode 100644 index 0000000..dc1a38c --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/block_step.cu @@ -0,0 +1,35 @@ +// CelerisLab/lbm/kernels/step/block_step.cu +// Generic blocked pull/collide/store for periodic/fine and slot/reservoir levels. +// Part of the blocks.cu compilation entry; relies on its includes and guards. +#if BLOCK_ROOT_BC == 2 +#error "BlockStep is not built for the channel root (BLOCK_ROOT_BC == 2)" +#endif + +// Pull/collide/store core cells only; halo cells of out are left untouched. +// flags (block cell ids, nullable): solid cores store NaN and never collide. +// kind (block cell ids, nullable = all ACTIVE): OUTER skips the pull and +// stores NaN, INNER pulls without collision, ACTIVE pulls and collides. +extern "C" __global__ void BlockStep(const fpxx* in, fpxx* out, unsigned int halo, + const uint16_t* flags, const unsigned char* kind) { + 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 long base = (unsigned long)blockIdx.z * NQ * TOTAL_CELLS; + const unsigned long cell = (unsigned long)blockIdx.z * TOTAL_CELLS + k; + 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]))) { + fill_invalid(f); + } else { + unsigned long j[NQ]; + compute_neighbors(x, y, j); // Core neighbors stay inside the window: no wrap. + stream_pull_load(k, f, in + base, j); + if (cell_kind == KIND_ACTIVE) { +#if BLOCK_ROOT_BC == 1 + block_wall(f, x, y, halo); +#endif + block_collide(f); + } + } + stream_pull_store(k, f, out + base); +} diff --git a/src/CelerisLab/lbm/kernels/step/block_transfer.cu b/src/CelerisLab/lbm/kernels/step/block_transfer.cu new file mode 100644 index 0000000..b2b09cb --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/block_transfer.cu @@ -0,0 +1,78 @@ +// CelerisLab/lbm/kernels/step/block_transfer.cu +// Same-level halo refresh and 2:1 explosion/coalescence between blocked levels. +// Part of the blocks.cu compilation entry; relies on its includes and guards. + +// 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) +{ + 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; + for (unsigned int q = 0; q < NQ; ++q) + to[index_f(d % TOTAL_CELLS, q)] = from[index_f(s % TOTAL_CELLS, q)]; +} + +// Fixed-patch explosion into blocked fine storage: target = owner block cell +// id of a patch-ghost cell; donors are parent block cell ids. Dense parent +// storage is the one-window case with coarse_cells = full domain size. +extern "C" __global__ void BlockExplode( + const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor, + unsigned int count, unsigned long coarse_cells) +{ + 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; + const uint64_t d = donor[row]; + const fpxx* from = coarse + (d / coarse_cells) * NQ * coarse_cells; + for (unsigned int q = 0; q < NQ; ++q) + to[index_f(t % TOTAL_CELLS, q)] = from[q * coarse_cells + d % coarse_cells]; +} + +// Experimental Schornbaum--Ruede (arXiv:1508.07982, Algorithm 2) +// and Chen et al. (doi:10.1016/j.physa.2005.09.036) population interpolation. +// Independent implementation of the mathematical stencil; no upstream source copied. +// Stencil cells are owner cores; a disabled axis repeats its center twice. +extern "C" __global__ void BlockExplodeProjected( + const fpxx* coarse, fpxx* fine, const uint64_t* target, const uint64_t* donor, + unsigned int count, unsigned long coarse_cells, + const uint64_t* neighbors, const float* offsets) +{ + 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; + 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]; + } + float gx = .5f * (v[1] - v[0]), gy = .5f * (v[3] - v[2]); + const float ex = d_cx[q], ey = d_cy[q]; + if (q != 0) { + const float along = (ex * gx + ey * gy) / (ex * ex + ey * ey); + gx -= ex * along; + gy -= ey * along; + } + to[index_f(t % TOTAL_CELLS, q)] = + coarse[(d / coarse_cells) * NQ * coarse_cells + q * coarse_cells + d % coarse_cells] + + offsets[2ul * row] * gx + offsets[2ul * row + 1] * gy; + } +} + +// Targets are receiver-local SOURCE slots for the immediately following pull. +// Keep child summation order identical to PatchCoarseStep; no collision here. +extern "C" __global__ void BlockCoalesce( + const fpxx* fine, fpxx* coarse, const uint64_t* target, const uint64_t* children, + unsigned int count) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + const uint64_t* s = children + 4ul * row; + coarse[target[row]] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); +} diff --git a/tests/unit/test_projected_explosion.py b/tests/unit/test_projected_explosion.py index 3c16634..bf6f59b 100644 --- a/tests/unit/test_projected_explosion.py +++ b/tests/unit/test_projected_explosion.py @@ -162,7 +162,7 @@ class ProjectedKernelTests(unittest.TestCase): targets = child.owner.ravel()[p.ghost_fine] donors = parent.owner.ravel()[p.ghost_coarse] source = (Path(__file__).resolve().parents[2] / - "src/CelerisLab/lbm/kernels/blocks.cu").read_text() + "src/CelerisLab/lbm/kernels/step/block_transfer.cu").read_text() 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): diff --git a/tests/unit/test_refinement_mrt.py b/tests/unit/test_refinement_mrt.py index 9302b72..c150728 100644 --- a/tests/unit/test_refinement_mrt.py +++ b/tests/unit/test_refinement_mrt.py @@ -22,7 +22,7 @@ class RefinementMRTHost(unittest.TestCase): def test_srt_statements_and_guards_preserved(self): """MRT is an additive compile-time branch; SRT statements and other guards remain.""" - blocks = (KERNELS / "blocks.cu").read_text() + blocks = (KERNELS / "blocks.cu").read_text() + (KERNELS / "step/block_step.cu").read_text() channel = (KERNELS / "step/block_channel.cu").read_text() block_collision = (KERNELS / "operators/block_collision.cuh").read_text() self.assertIn('#include "operators/block_collision.cuh"', blocks)