diff --git a/src/CelerisLab/lbm/kernels/refinement.cu b/src/CelerisLab/lbm/kernels/refinement.cu index 3efd0fd..5d6ac9d 100644 --- a/src/CelerisLab/lbm/kernels/refinement.cu +++ b/src/CelerisLab/lbm/kernels/refinement.cu @@ -36,93 +36,5 @@ __device__ unsigned int incoming_counts[NQ * TOTAL_CELLS]; } #endif -__device__ __forceinline__ void patch_collide(float* f, unsigned long k) { -#if COLLISION_MODEL == 2 - collide_mrt_unforced(f); -#else - float rho, ux, uy, feq[NQ]; - compute_rho_u(f, rho, ux, uy); - compute_feq(rho, ux, uy, feq); - collide_srt_no_force(f, feq, d_params.omega); -#endif -#ifdef CELERIS_REFINEMENT_AUDIT - ++collision_counts[k]; -#endif -} - -extern "C" __global__ void PatchExplode( - const float* coarse, float* fine, const uint32_t* targets, - const uint32_t* donors, unsigned int count, unsigned long coarse_cells) -{ - unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; - if (row >= count) return; - for (unsigned int q = 0; q < NQ; ++q) - fine[index_f(targets[row], q)] = coarse[q * coarse_cells + donors[row]]; -#ifdef CELERIS_REFINEMENT_AUDIT - ++explosion_counts[targets[row]]; -#endif -} - -extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint16_t* flags) { - unsigned int x, y; unsigned long k; - index_from_thread(x, y, k); - if (x >= NX || y >= NY) return; - float f[NQ]; - // Outer halo is invalid after EACH hop, including when reusing the buffer - // 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]))) { - fill_invalid(f); - } else { - unsigned long j[NQ]; - compute_neighbors(x, y, j); // These interior targets never wrap. - stream_pull_load(k, f, in, j); - if (x >= 4 && y >= 4 && x < NX - 4 && y < NY - 4) - patch_collide(f, k); - } - stream_pull_store(k, f, out); -} - -extern "C" __global__ void PatchCoarseStep( - const float* in, float* out, const float* fine, - const uint64_t* targets, const uint64_t* children, unsigned int count, - unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, - unsigned int wall_aware) -{ - unsigned int x, y; unsigned long k; - index_from_thread(x, y, k); - if (x >= NX || y >= NY) return; - // Both buffers were poisoned at initialization; covered cells stay invalid. - if (x >= x0 && x < x1 && y >= y0 && y < y1) return; - if (wall_aware && (y == 0 || y == NY - 1)) return; - unsigned long j[NQ]; - compute_neighbors(x, y, j); - float f[NQ]; - stream_pull_load(k, f, in, j); - if (x >= x0 - 1 && x <= x1 && y >= y0 - 1 && y <= y1) { - for (unsigned int q = 1; q < NQ; ++q) { - uint64_t slot = index_f(k, q); - unsigned int lo = 0, hi = count; - while (lo < hi) { - unsigned int mid = lo + (hi - lo) / 2; - if (targets[mid] < slot) lo = mid + 1; - else hi = mid; - } - if (lo < count && targets[lo] == slot) { - const uint64_t* s = children + 4ul * lo; - f[q] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); -#ifdef CELERIS_REFINEMENT_AUDIT - ++incoming_counts[slot]; -#endif - } - } - } - // Includes the eight grazing-corner links with UNCOVERED coarse donors. - // All replacement arrivals must be in place before computing rho/u. - if (wall_aware) { - if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x)); - else if (y == NY - 2) reconstruct_top_reservoir(f, RHO); - } - patch_collide(f, k); - stream_pull_store(k, f, out); -} +// Kernels: patch collision helper, PatchExplode, PatchFineStep, PatchCoarseStep. +#include "step/patch_step.cu" diff --git a/src/CelerisLab/lbm/kernels/step/patch_step.cu b/src/CelerisLab/lbm/kernels/step/patch_step.cu new file mode 100644 index 0000000..1f64ea7 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/patch_step.cu @@ -0,0 +1,95 @@ +// CelerisLab/lbm/kernels/step/patch_step.cu +// Dense fixed-patch collision, explosion and the fine/coarse steps. +// Part of the refinement.cu compilation entry; relies on its includes, guards +// and the optional CELERIS_REFINEMENT_AUDIT counters declared there. + +__device__ __forceinline__ void patch_collide(float* f, unsigned long k) { +#if COLLISION_MODEL == 2 + collide_mrt_unforced(f); +#else + float rho, ux, uy, feq[NQ]; + compute_rho_u(f, rho, ux, uy); + compute_feq(rho, ux, uy, feq); + collide_srt_no_force(f, feq, d_params.omega); +#endif +#ifdef CELERIS_REFINEMENT_AUDIT + ++collision_counts[k]; +#endif +} + +extern "C" __global__ void PatchExplode( + const float* coarse, float* fine, const uint32_t* targets, + const uint32_t* donors, unsigned int count, unsigned long coarse_cells) +{ + unsigned int row = blockIdx.x * blockDim.x + threadIdx.x; + if (row >= count) return; + for (unsigned int q = 0; q < NQ; ++q) + fine[index_f(targets[row], q)] = coarse[q * coarse_cells + donors[row]]; +#ifdef CELERIS_REFINEMENT_AUDIT + ++explosion_counts[targets[row]]; +#endif +} + +extern "C" __global__ void PatchFineStep(const float* in, float* out, const uint16_t* flags) { + unsigned int x, y; unsigned long k; + index_from_thread(x, y, k); + if (x >= NX || y >= NY) return; + float f[NQ]; + // Outer halo is invalid after EACH hop, including when reusing the buffer + // 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]))) { + fill_invalid(f); + } else { + unsigned long j[NQ]; + compute_neighbors(x, y, j); // These interior targets never wrap. + stream_pull_load(k, f, in, j); + if (x >= 4 && y >= 4 && x < NX - 4 && y < NY - 4) + patch_collide(f, k); + } + stream_pull_store(k, f, out); +} + +extern "C" __global__ void PatchCoarseStep( + const float* in, float* out, const float* fine, + const uint64_t* targets, const uint64_t* children, unsigned int count, + unsigned int x0, unsigned int y0, unsigned int x1, unsigned int y1, + unsigned int wall_aware) +{ + unsigned int x, y; unsigned long k; + index_from_thread(x, y, k); + if (x >= NX || y >= NY) return; + // Both buffers were poisoned at initialization; covered cells stay invalid. + if (x >= x0 && x < x1 && y >= y0 && y < y1) return; + if (wall_aware && (y == 0 || y == NY - 1)) return; + unsigned long j[NQ]; + compute_neighbors(x, y, j); + float f[NQ]; + stream_pull_load(k, f, in, j); + if (x >= x0 - 1 && x <= x1 && y >= y0 - 1 && y <= y1) { + for (unsigned int q = 1; q < NQ; ++q) { + uint64_t slot = index_f(k, q); + unsigned int lo = 0, hi = count; + while (lo < hi) { + unsigned int mid = lo + (hi - lo) / 2; + if (targets[mid] < slot) lo = mid + 1; + else hi = mid; + } + if (lo < count && targets[lo] == slot) { + const uint64_t* s = children + 4ul * lo; + f[q] = 0.25f * (fine[s[0]] + fine[s[1]] + fine[s[2]] + fine[s[3]]); +#ifdef CELERIS_REFINEMENT_AUDIT + ++incoming_counts[slot]; +#endif + } + } + } + // Includes the eight grazing-corner links with UNCOVERED coarse donors. + // All replacement arrivals must be in place before computing rho/u. + if (wall_aware) { + if (y == 1) reconstruct_bottom_velocity(f, 0.0f, bottom_wall_normal_velocity(x)); + else if (y == NY - 2) reconstruct_top_reservoir(f, RHO); + } + patch_collide(f, k); + stream_pull_store(k, f, out); +} diff --git a/tests/unit/test_refinement_mrt.py b/tests/unit/test_refinement_mrt.py index c150728..dd2f93a 100644 --- a/tests/unit/test_refinement_mrt.py +++ b/tests/unit/test_refinement_mrt.py @@ -27,7 +27,7 @@ class RefinementMRTHost(unittest.TestCase): block_collision = (KERNELS / "operators/block_collision.cuh").read_text() self.assertIn('#include "operators/block_collision.cuh"', blocks) self.assertEqual(blocks.count("block_collide(f);"), 1) - refinement = (KERNELS / "refinement.cu").read_text() + refinement = (KERNELS / "refinement.cu").read_text() + (KERNELS / "step/patch_step.cu").read_text() for text in ("compute_feq_from_momentum(drho, jx, jy, ux, uy, feq);", "compute_feq(rho, ux, uy, feq);", "collide_srt_no_force(f, feq, d_params.omega);"):