diff --git a/src/CelerisLab/lbm/kernels/blocks.cu b/src/CelerisLab/lbm/kernels/blocks.cu index 2ff9f89..21f61b3 100644 --- a/src/CelerisLab/lbm/kernels/blocks.cu +++ b/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. diff --git a/src/CelerisLab/lbm/kernels/core/cell_kind.cuh b/src/CelerisLab/lbm/kernels/core/cell_kind.cuh index 16a0c61..c90a147 100644 --- a/src/CelerisLab/lbm/kernels/core/cell_kind.cuh +++ b/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 diff --git a/src/CelerisLab/lbm/kernels/core/layout.cuh b/src/CelerisLab/lbm/kernels/core/layout.cuh index 3b0f76a..1fce767 100644 --- a/src/CelerisLab/lbm/kernels/core/layout.cuh +++ b/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 +__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) // --------------------------------------------------------------------------- diff --git a/src/CelerisLab/lbm/kernels/refinement.cu b/src/CelerisLab/lbm/kernels/refinement.cu index f496b10..3efd0fd 100644 --- a/src/CelerisLab/lbm/kernels/refinement.cu +++ b/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. diff --git a/src/CelerisLab/lbm/kernels/step/block_channel.cu b/src/CelerisLab/lbm/kernels/step/block_channel.cu index e570a8c..5791a50 100644 --- a/src/CelerisLab/lbm/kernels/step/block_channel.cu +++ b/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);