diff --git a/.gitignore b/.gitignore index f2ebeeb..ee1fdd0 100644 --- a/.gitignore +++ b/.gitignore @@ -66,3 +66,6 @@ venv.bak/ *.tmp *.bak *.log + +# reference: +ref/ \ No newline at end of file diff --git a/GiteaSyncTest b/GiteaSyncTest deleted file mode 100644 index e69de29..0000000 diff --git a/output/cylinder_d2q9.png b/output/cylinder_d2q9.png new file mode 100644 index 0000000..c75dfb4 Binary files /dev/null and b/output/cylinder_d2q9.png differ diff --git a/output/cylinder_d3q19_re100.png b/output/cylinder_d3q19_re100.png new file mode 100644 index 0000000..4617ce0 Binary files /dev/null and b/output/cylinder_d3q19_re100.png differ diff --git a/output/cylinder_d3q19_re20.png b/output/cylinder_d3q19_re20.png new file mode 100644 index 0000000..87355fd Binary files /dev/null and b/output/cylinder_d3q19_re20.png differ diff --git a/output/poiseuille_d2q9.png b/output/poiseuille_d2q9.png new file mode 100644 index 0000000..752f3a1 Binary files /dev/null and b/output/poiseuille_d2q9.png differ diff --git a/src/CelerisLab/cuda/compiler.py b/src/CelerisLab/cuda/compiler.py index 3734d75..2164b71 100644 --- a/src/CelerisLab/cuda/compiler.py +++ b/src/CelerisLab/cuda/compiler.py @@ -51,6 +51,19 @@ def compile_kernel(): ] ) + +def compile_kernel_v2(): + """Compile the new modular kernel (kernel_v2.cu → kernel_v2.ptx).""" + subprocess.run( + [ + "nvcc", + "-ptx", + kernel_path("kernel_v2.cu"), + "-o", + kernel_path("kernel_v2.ptx"), + ] + ) + def config_kernal(config_cuda: CudaConfig, config_field: FlowFieldConfig): lines = read_lines(kernel_path("macros.h")) lines = modify_macro(lines, "MULT_GPU", config_cuda.multi_gpu) @@ -77,6 +90,31 @@ def config_kernal(config_cuda: CudaConfig, config_field: FlowFieldConfig): write_lines(kernel_path("macros.h"), lines) + +def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig, + collision_model: int = 2, + streaming_model: int = 0, + store_precision: int = 0, + use_ddf_shifting: int = 0): + """Configure macros.h for the new modular kernel architecture. + + Args: + collision_model: 0=SRT, 1=TRT, 2=MRT (default) + streaming_model: 0=double-buffer (default), 1=Esoteric-Pull + store_precision: 0=FP32 (default), 1=FP16S, 2=FP16C + use_ddf_shifting: 0=off (default), 1=on + """ + # First apply legacy config + config_kernal(config_cuda, config_field) + + # Then apply new architecture macros + lines = read_lines(kernel_path("macros.h")) + lines = modify_macro(lines, "COLLISION_MODEL", collision_model) + lines = modify_macro(lines, "STREAMING_MODEL", streaming_model) + lines = modify_macro(lines, "STORE_PRECISION", store_precision) + lines = modify_macro(lines, "USE_DDF_SHIFTING", use_ddf_shifting) + write_lines(kernel_path("macros.h"), lines) + def config_object(n_obj: int): lines = read_lines(kernel_path("macros.h")) lines = modify_macro(lines, "N_OBJS", n_obj) diff --git a/src/CelerisLab/lbm/configs/config_gym.json b/src/CelerisLab/lbm/configs/config_gym.json deleted file mode 100644 index 544b7b4..0000000 --- a/src/CelerisLab/lbm/configs/config_gym.json +++ /dev/null @@ -1,3 +0,0 @@ -{ - -} \ No newline at end of file diff --git a/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh b/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh new file mode 100644 index 0000000..eb9e21e --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh @@ -0,0 +1,131 @@ +// CelerisLab – boundary/bounce_back.cuh +// Full-way bounce-back for grid-aligned solid walls. +// +// For nodes adjacent to a solid wall, the incoming population from the +// wall direction is reflected: +// f[i] = f[opp(i)] (+ wall velocity correction for moving walls) +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_BOUNCE_BACK_CUH +#define CELERIS_BOUNDARY_BOUNCE_BACK_CUH + +// --------------------------------------------------------------------------- +// Simple bounce-back: replace f[i] with f[opp(i)] for solid neighbors +// This version checks the flag of each neighbor. +// --------------------------------------------------------------------------- +__device__ inline void apply_bounce_back(float* __restrict__ f, + const unsigned long* __restrict__ j, + const uint8_t* __restrict__ flags_arr) +{ + // For each direction pair (i, i+1): + // If the neighbor in direction i is SOLID → f[i] should bounce + // from the opposite direction. + // + // With paired ordering, opp(i) = i+1 and opp(i+1) = i. + // In pull streaming, f[i] was loaded from j[opp(i)] = j[i+1]. + // If j[i+1] is solid, we need bounce-back. + + for (int i = 1; i < NQ; i += 2) { + // Direction i: check if source (j[i+1]) is solid + if (flags_arr[j[i + 1]] & LEGACY_SOLID) { + // Bounce: f[i] = f_post[opp(i)] at current node = pre-collision f[i+1] + // But in pull-streaming context, we already loaded from neighbor. + // For simple BB on a flat wall, just swap f[i] and f[i+1]: + float temp = f[i]; + f[i] = f[i + 1]; + f[i + 1] = temp; + } + } +} + +// --------------------------------------------------------------------------- +// Top/bottom wall bounce-back (specialized for channel flow) +// For nodes at y=1 (adjacent to y=0 wall) or y=NY-2 (adjacent to y=NY-1 wall) +// --------------------------------------------------------------------------- +#if NQ == 9 +__device__ inline void apply_wall_bb_d2q9(unsigned int y, + float* __restrict__ f) +{ + if (y == 1) { + float temp; + temp = f[3]; f[3] = f[4]; f[4] = temp; // ±y swap + temp = f[5]; f[5] = f[6]; f[6] = temp; // ±(x+y) swap + temp = f[7]; f[7] = f[8]; f[8] = temp; // ±(x-y) swap + } + else if (y == NY - 2) { + float temp; + temp = f[3]; f[3] = f[4]; f[4] = temp; + temp = f[5]; f[5] = f[6]; f[6] = temp; + temp = f[7]; f[7] = f[8]; f[8] = temp; + } +} +#endif + +// --------------------------------------------------------------------------- +// D3Q19 wall bounce-back on y=0/NY-1 (channel walls) +// Pairs with non-zero cy: (3,4)±y, (7,8)±(x+y), (11,12)±(y+z), +// (13,14)±(x-y), (17,18)±(y-z) +// --------------------------------------------------------------------------- +// --------------------------------------------------------------------------- +// D3Q19 y-wall half-way bounce-back for pull double-buffer streaming. +// +// Instead of swapping f[i] <-> f[opp(i)] (which corrupts the toward-wall +// direction with garbage from the wall node), we read the opposite direction +// at the SAME node from the INPUT buffer. That value is the post-collision +// population from the previous time step — exactly the correct half-way BB. +// +// At y=1 (wall at y=0): directions that arrive FROM y=0 have cy_src = -1. +// The pull loaded them from the wall (garbage). Replace with fi_in[k, opp]. +// +// At y=NY-2 (wall at y=NY-1): directions from y=NY-1, cy_src = +1. +// --------------------------------------------------------------------------- +#if NQ == 19 +__device__ inline void apply_wall_bb_d3q19_y_pull(unsigned int y, + float* __restrict__ f, + const fpxx* __restrict__ fi_in, + unsigned long k) +{ + if (y == 1) { + // Directions whose pull-source is at y=0 (wall): + // f[3] (+y) ← source j[4] = (x, 0, z) + // f[7] (+x+y) ← source j[8] = (x-1, 0, z) + // f[11] (+y+z) ← source j[12] = (x, 0, z-1) + // f[14] (-x+y) ← source j[13] = (x+1, 0, z) + // f[17] (+y-z) ← source j[18] = (x, 0, z+1) + f[3] = load_ddf(fi_in, index_f(k, 4u)); + f[7] = load_ddf(fi_in, index_f(k, 8u)); + f[11] = load_ddf(fi_in, index_f(k, 12u)); + f[14] = load_ddf(fi_in, index_f(k, 13u)); + f[17] = load_ddf(fi_in, index_f(k, 18u)); + } + else if (y == (unsigned int)(NY - 2)) { + // Directions whose pull-source is at y=NY-1 (wall): + // f[4] (-y) ← source j[3] = (x, NY-1, z) + // f[8] (-x-y) ← source j[7] = (x+1, NY-1, z) + // f[12] (-y-z) ← source j[11] = (x, NY-1, z+1) + // f[13] (+x-y) ← source j[14] = (x-1, NY-1, z) + // f[18] (-y+z) ← source j[17] = (x, NY-1, z-1) + f[4] = load_ddf(fi_in, index_f(k, 3u)); + f[8] = load_ddf(fi_in, index_f(k, 7u)); + f[12] = load_ddf(fi_in, index_f(k, 11u)); + f[13] = load_ddf(fi_in, index_f(k, 14u)); + f[18] = load_ddf(fi_in, index_f(k, 17u)); + } +} + +// Keep the old swap version for backward compatibility +__device__ inline void apply_wall_bb_d3q19_y(unsigned int y, + float* __restrict__ f) +{ + if (y == 1 || y == NY - 2) { + float temp; + temp = f[3]; f[3] = f[4]; f[4] = temp; // ±y + temp = f[7]; f[7] = f[8]; f[8] = temp; // ±(x+y) + temp = f[11]; f[11] = f[12]; f[12] = temp; // ±(y+z) + temp = f[13]; f[13] = f[14]; f[14] = temp; // ±(x-y) + temp = f[17]; f[17] = f[18]; f[18] = temp; // ±(y-z) + } +} +#endif + +#endif // CELERIS_BOUNDARY_BOUNCE_BACK_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh new file mode 100644 index 0000000..fffa618 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/curved_boundary.cuh @@ -0,0 +1,111 @@ +// CelerisLab – boundary/curved_boundary.cuh +// Interpolated curved-boundary bounce-back. +// +// For curved solid surfaces that don't align with the grid, the +// reflected population is interpolated using the fractional distance q +// between the fluid node and the wall intersection. +// +// Migrated from existing kernel.cu + ref/main.cu. +// +// Delta pool layout per curved-boundary node (11 floats for D2Q9): +// delta[offset + 0] : encoded object id (bitcast int → float) +// delta[offset + 1..NQ-1] : q values for each direction (0 if no wall hit) +// delta[offset + NQ] : normal_y (wall normal y-component / R) +// delta[offset + NQ+1]: normal_x (wall normal x-component / R) +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_CURVED_BOUNDARY_CUH +#define CELERIS_BOUNDARY_CURVED_BOUNDARY_CUH + +#if NQ == 9 + +// --------------------------------------------------------------------------- +// apply_curved_boundary: +// Interpolated bounce-back with wall velocity correction. +// Operates on a SOLID+INTERFACE node's neighbors. +// +// Parameters: +// n – current node (solid + interface) +// x, y – coordinates of current node +// f_out – the output DDF buffer (double-buffer scheme) being modified +// delta – parameter pool +// id_off – offset into delta for this node (= indx[n]) +// Uw, Vw – wall velocity at this node +// obs_fx, obs_fy – accumulated force observation (atomicAdd targets) +// --------------------------------------------------------------------------- +__device__ inline void apply_curved_boundary( + unsigned long n, unsigned int x, unsigned int y, + fpxx* __restrict__ f_out, + const float* __restrict__ delta, + int id_off, + float Uw, float Vw, + float* __restrict__ obs_fx, + float* __restrict__ obs_fy) +{ + // New paired direction ordering: + // cx = {0, 1,-1, 0, 0, 1,-1, 1,-1} + // cy = {0, 0, 0, 1,-1, 1,-1,-1, 1} + + for (int i = 1; i < NQ; i++) { + int x_neb = x + d_cx[i]; + int y_neb = y + d_cy[i]; + + // Check bounds (skip if neighbor is outside domain) + if (x_neb < 0 || x_neb >= NX || y_neb < 0 || y_neb >= NY) continue; + + unsigned long k_neb = linear_index((unsigned int)x_neb, (unsigned int)y_neb); + + // Only process if neighbor is FLUID + // (We read the flag from global memory – this is a boundary kernel, + // called infrequently per node, so the extra read is acceptable.) + // The caller should ensure this node IS solid+interface. + + float q = delta[id_off + i]; // fractional distance (0 if no wall hit) + if (q <= 0.0f) continue; // no wall intersection in this direction + + int oi = opp_dir(i); + float ci_dot_uw = (float)d_cx[i] * Uw + (float)d_cy[i] * Vw; + float wall_term = 6.0f * d_w[i] * ci_dot_uw; + + // Second neighbor for quadratic interpolation + int x_neb2 = x + 2 * d_cx[i]; + int y_neb2 = y + 2 * d_cy[i]; + // Clamp to domain (simple, could be improved) + x_neb2 = max(0, min(NX - 1, x_neb2)); + y_neb2 = max(0, min(NY - 1, y_neb2)); + unsigned long k_neb2 = linear_index((unsigned int)x_neb2, (unsigned int)y_neb2); + + // Read current DDF values from output buffer + float f_self_opp = load_ddf(f_out, index_f(n, (unsigned int)oi)); + float f_neb_opp = load_ddf(f_out, index_f(k_neb, (unsigned int)oi)); + float f_neb2_fwd = load_ddf(f_out, index_f(k_neb2, (unsigned int)i)); + + // Interpolated bounce-back (Yu, Mei, Shyy, 2003) + float f_reflected = (q * f_self_opp + (1.0f - q) * f_neb_opp + + q * f_neb2_fwd + wall_term) / (1.0f + q); + + store_ddf(f_out, index_f(k_neb, (unsigned int)i), f_reflected); + + // Force observation (momentum exchange) + float f_neb_fwd = load_ddf(f_out, index_f(k_neb, (unsigned int)i)); + float f_self_fwd = load_ddf(f_out, index_f(n, (unsigned int)i)); + float f_sum = f_neb_fwd + f_self_opp; + + int x_back = x - d_cx[i]; + int y_back = y - d_cy[i]; + x_back = max(0, min(NX - 1, x_back)); + y_back = max(0, min(NY - 1, y_back)); + + // Accumulate force + if (obs_fx != nullptr) { + atomicAdd(obs_fx, -f_sum * (float)d_cx[i] + f_self_fwd); + atomicAdd(obs_fy, -f_sum * (float)d_cy[i] + + load_ddf(f_out, index_f(linear_index((unsigned int)x_back, (unsigned int)y_back), + (unsigned int)i))); + } + } +} + +#endif // NQ == 9 + +#endif // CELERIS_BOUNDARY_CURVED_BOUNDARY_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/ibm_kernels.cuh b/src/CelerisLab/lbm/kernels/boundary/ibm_kernels.cuh new file mode 100644 index 0000000..f139fb3 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/ibm_kernels.cuh @@ -0,0 +1,174 @@ +// CelerisLab – boundary/ibm_kernels.cuh +// Immersed Boundary Method (IBM) kernels: Euler↔Lagrangian coupling. +// Migrated from ref/main.cu: Eul2Lag, Lag2Eul, Update_Lagrangian. +// +// Uses Peskin's 4-point delta function for interpolation/spreading. +// +// Data layout (per rigid body): +// RigidBodyState2D: {x, y, theta, vx, vy, omega} +// Lagrangian points: motionL[NT_LAG * 4] +// motionL[4*k + 0] = target vx at point k +// motionL[4*k + 1] = target vy at point k +// motionL[4*k + 2] = XL (x-position) +// motionL[4*k + 3] = YL (y-position) +// forceL[NT_LAG * 2] = Lagrangian force (fx, fy) per point +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_IBM_KERNELS_CUH +#define CELERIS_BOUNDARY_IBM_KERNELS_CUH + +#include "../core/params.cuh" // RigidBodyState2D, RigidBodyControl2D + +// IBM is currently 2D only. Guard all kernels. +#if DIM == 2 + +// --------------------------------------------------------------------------- +// Peskin 4-point delta function (1D component) +// --------------------------------------------------------------------------- +__device__ __forceinline__ float peskin_delta_1d(float r) { + float ar = fabsf(r); + if (ar < 1.0f) + return 0.125f * (3.0f - 2.0f * ar + sqrtf(1.0f + 4.0f * ar - 4.0f * ar * ar)); + else if (ar < 2.0f) + return 0.125f * (5.0f - 2.0f * ar - sqrtf(-7.0f + 12.0f * ar - 4.0f * ar * ar)); + else + return 0.0f; +} + +// --------------------------------------------------------------------------- +// Peskin 4-point delta function (2D = product of 1D components) +// --------------------------------------------------------------------------- +__device__ __forceinline__ float peskin_delta_2d(float dx, float dy) { + return peskin_delta_1d(dx) * peskin_delta_1d(dy); +} + +// --------------------------------------------------------------------------- +// update_lagrangian_control: +// Update Lagrangian point positions using RigidBodyState2D. +// Supports circular body (points on circumference). +// +// Thread k ∈ [0, NT_LAG) — one thread per Lagrangian point. +// --------------------------------------------------------------------------- +__global__ void update_lagrangian_control( + float* __restrict__ motionL, // [NT_LAG * 4] + const RigidBodyState2D* state, // single rigid body state + float radius, // body radius + int NT_LAG) // number of Lagrangian points +{ + int k = threadIdx.x + blockIdx.x * blockDim.x; + if (k >= NT_LAG) return; + + float angle = state->theta + (float)k * 2.0f * 3.14159265358979323846f / (float)NT_LAG; + + // Target velocity at this Lagrangian point (rigid-body kinematics) + float r_sin = radius * sinf(angle); + float r_cos = radius * cosf(angle); + + motionL[4 * k + 0] = state->vx - state->omega * r_sin; // target vx + motionL[4 * k + 1] = state->vy + state->omega * r_cos; // target vy + motionL[4 * k + 2] = state->x + r_cos; // XL + motionL[4 * k + 3] = state->y + r_sin; // YL +} + +// --------------------------------------------------------------------------- +// euler_to_lagrangian (Eul2Lag): +// Interpolate Eulerian velocity to Lagrangian points. +// Compute Lagrangian force = penalty * (target_vel - interpolated_vel) * dL. +// +// Thread k ∈ [0, NT_LAG) +// --------------------------------------------------------------------------- +__global__ void euler_to_lagrangian( + float* __restrict__ forceL, // [NT_LAG * 2] output forces + const float* __restrict__ motionL, // [NT_LAG * 4] + const float* __restrict__ u_field, // [DIM * N] SoA velocity field + int NT_LAG) +{ + int k = threadIdx.x + blockIdx.x * blockDim.x; + if (k >= NT_LAG) return; + + float XL = motionL[4 * k + 2]; + float YL = motionL[4 * k + 3]; + int ix = (int)XL; + int iy = (int)YL; + + // Arc-length element dL (central difference) + float dL; + if (k > 0 && k < NT_LAG - 1) { + float dx1 = motionL[4*k+2] - motionL[4*(k-1)+2]; + float dy1 = motionL[4*k+3] - motionL[4*(k-1)+3]; + float dx2 = motionL[4*(k+1)+2] - motionL[4*k+2]; + float dy2 = motionL[4*(k+1)+3] - motionL[4*k+3]; + dL = 0.5f * (sqrtf(dx1*dx1 + dy1*dy1) + sqrtf(dx2*dx2 + dy2*dy2)); + } else if (k == 0) { + float dx1 = motionL[4*1+2] - motionL[4*0+2]; + float dy1 = motionL[4*1+3] - motionL[4*0+3]; + dL = sqrtf(dx1*dx1 + dy1*dy1); + } else { + float dx1 = motionL[4*k+2] - motionL[4*(k-1)+2]; + float dy1 = motionL[4*k+3] - motionL[4*(k-1)+3]; + dL = sqrtf(dx1*dx1 + dy1*dy1); + } + + // Interpolate velocity from Euler grid using Peskin 4-point delta + float UL = 0.0f, VL = 0.0f; + for (int dy = -2; dy <= 2; dy++) { + for (int dx = -2; dx <= 2; dx++) { + int ex = ix + dx; + int ey = iy + dy; + if (ex < 0 || ex >= NX || ey < 0 || ey >= NY) continue; + + float ww = peskin_delta_2d((float)ex - XL, (float)ey - YL); + unsigned long en = (unsigned long)ey * NX + (unsigned long)ex; + UL += u_field[en] * ww; // ux + VL += u_field[TOTAL_CELLS + en] * ww; // uy + } + } + + // Force = (target_vel - interpolated_vel) * dL + forceL[2 * k + 0] = 2.0f * (motionL[4*k+0] - UL) * dL; // fx + forceL[2 * k + 1] = 2.0f * (motionL[4*k+1] - VL) * dL; // fy +} + +// --------------------------------------------------------------------------- +// lagrangian_to_euler (Lag2Eul): +// Spread Lagrangian forces back to the Eulerian grid. +// +// Thread (x, y) — one thread per Euler node. +// Only processes nodes near the body (bounding box check). +// --------------------------------------------------------------------------- +__global__ void lagrangian_to_euler( + float* __restrict__ forceE, // [DIM * N] output Euler force field + const float* __restrict__ forceL, // [NT_LAG * 2] + const float* __restrict__ motionL, // [NT_LAG * 4] + float xc, float yc, float L0, // body center and extent + int NT_LAG) +{ + unsigned int x, y; + unsigned long k; +#if DIM == 2 + index_from_thread(x, y, k); +#endif + + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; + + // Bounding box check + if ((float)x < xc - L0 || (float)x > xc + L0 || + (float)y < yc - L0 || (float)y > yc + L0) return; + + float forcex = 0.0f, forcey = 0.0f; + + for (int kl = 0; kl < NT_LAG; kl++) { + float XL = motionL[4 * kl + 2]; + float YL = motionL[4 * kl + 3]; + float ww = peskin_delta_2d((float)x - XL, (float)y - YL); + forcex += forceL[2 * kl + 0] * ww; + forcey += forceL[2 * kl + 1] * ww; + } + + forceE[k] = forcex; + forceE[TOTAL_CELLS + k] = forcey; +} + +#endif // DIM == 2 + +#endif // CELERIS_BOUNDARY_IBM_KERNELS_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh new file mode 100644 index 0000000..3a87a5a --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh @@ -0,0 +1,173 @@ +// CelerisLab – boundary/inlet_outlet.cuh +// Inlet and outlet boundary conditions (D2Q9). +// +// Parabolic inlet (non-equilibrium extrapolation, Zou-He style): +// Left wall (x=0): reconstruct cx>0 populations (i=1,5,7) +// +// Pressure outlet (non-equilibrium extrapolation): +// Right wall (x=NX-1): reconstruct cx<0 populations (i=2,6,8) +// +// New paired D2Q9 ordering: +// cx = {0, 1,-1, 0, 0, 1,-1, 1,-1} +// cy = {0, 0, 0, 1,-1, 1,-1,-1, 1} +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_INLET_OUTLET_CUH +#define CELERIS_BOUNDARY_INLET_OUTLET_CUH + +#if NQ == 9 + +// --------------------------------------------------------------------------- +// Parabolic inlet (x = 0, non-equilibrium extrapolation) +// +// f, f_neb are local DDF arrays: +// f = populations at the boundary node (x=0) +// f_neb = populations at the interior neighbor (x=1) +// y = y-coordinate of the boundary node +// +// Reconstructs f[1], f[5], f[7] (cx > 0 directions in new ordering) +// using: f_bc[i] = f_neb[i] - feq(rho_neb, u_neb)[i] + feq(rho_neb, u_target)[i] +// --------------------------------------------------------------------------- +__device__ inline void apply_parabolic_inlet(float* __restrict__ f, + const float* __restrict__ f_neb, + float y_coord) +{ + // Neighbor macros + float p_neb = (f_neb[0]+f_neb[1]+f_neb[2]+f_neb[3]+f_neb[4] + +f_neb[5]+f_neb[6]+f_neb[7]+f_neb[8]) / 3.0f; + + // Target velocity (parabolic profile) + float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f); + float u_target = U0 * 1.5f * (1.0f - 4.0f * yy * yy); + float v_target = 0.0f; + + // Neighbor velocity + float u_neb = (f_neb[1]-f_neb[2]+f_neb[5]-f_neb[6]+f_neb[7]-f_neb[8]) / RHO; + float v_neb = (f_neb[3]-f_neb[4]+f_neb[5]-f_neb[6]-f_neb[7]+f_neb[8]) / RHO; + + // feq for direction i=1 (cx=1, cy=0), w=1/9: + // feq = (2p + RHO*(2u² + 2u - v²)) / 6 + float feq1_target = (2.0f*p_neb + RHO*(2.0f*u_target*u_target + 2.0f*u_target - v_target*v_target)) / 6.0f; + float feq1_neb = (2.0f*p_neb + RHO*(2.0f*u_neb*u_neb + 2.0f*u_neb - v_neb*v_neb)) / 6.0f; + + // feq for direction i=5 (cx=1, cy=1), w=1/36: + // feq = (p + RHO*(u² + 3uv + u + v² + v)) / 12 + float feq5_target = (p_neb + RHO*(u_target*u_target + 3.0f*u_target*v_target + u_target + v_target*v_target + v_target)) / 12.0f; + float feq5_neb = (p_neb + RHO*(u_neb*u_neb + 3.0f*u_neb*v_neb + u_neb + v_neb*v_neb + v_neb)) / 12.0f; + + // feq for direction i=7 (cx=1, cy=-1), w=1/36: + // feq = (p + RHO*(u² - 3uv + u + v² - v)) / 12 + float feq7_target = (p_neb + RHO*(u_target*u_target - 3.0f*u_target*v_target + u_target + v_target*v_target - v_target)) / 12.0f; + float feq7_neb = (p_neb + RHO*(u_neb*u_neb - 3.0f*u_neb*v_neb + u_neb + v_neb*v_neb - v_neb)) / 12.0f; + + // Non-equilibrium extrapolation + f[1] = f_neb[1] - feq1_neb + feq1_target; + f[5] = f_neb[5] - feq5_neb + feq5_target; + f[7] = f_neb[7] - feq7_neb + feq7_target; +} + +// --------------------------------------------------------------------------- +// Pressure outlet (x = NX-1, non-equilibrium extrapolation) +// +// Reconstructs f[2], f[6], f[8] (cx < 0 directions in new ordering) +// p_out = 0 (gauge pressure), uses velocity from neighbor. +// --------------------------------------------------------------------------- +__device__ inline void apply_pressure_outlet(float* __restrict__ f, + const float* __restrict__ f_neb, + float y_coord) +{ + float p_out = 0.0f; + + // Target velocity (parabolic, same as inlet for consistency) + float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f); + float u_target = U0 * 1.5f * (1.0f - 4.0f * yy * yy); + float v_target = 0.0f; + + // Neighbor velocity + float u_neb = (f_neb[1]-f_neb[2]+f_neb[5]-f_neb[6]+f_neb[7]-f_neb[8]) / RHO; + float v_neb = (f_neb[3]-f_neb[4]+f_neb[5]-f_neb[6]-f_neb[7]+f_neb[8]) / RHO; + + // feq for direction i=2 (cx=-1, cy=0), w=1/9: + // feq = (2p - RHO*(-2u² + 2u + v²)) / 6 + float feq2_target = (2.0f*p_out - RHO*(-2.0f*u_target*u_target + 2.0f*u_target + v_target*v_target)) / 6.0f; + float feq2_neb = (2.0f*p_out - RHO*(-2.0f*u_neb*u_neb + 2.0f*u_neb + v_neb*v_neb)) / 6.0f; + + // feq for direction i=8 (cx=-1, cy=1), w=1/36: + // feq = (p + RHO*(u² - 3uv - u + v² + v)) / 12 + float feq8_target = (p_out + RHO*(u_target*u_target - 3.0f*u_target*v_target - u_target + v_target*v_target + v_target)) / 12.0f; + float feq8_neb = (p_out + RHO*(u_neb*u_neb - 3.0f*u_neb*v_neb - u_neb + v_neb*v_neb + v_neb)) / 12.0f; + + // feq for direction i=6 (cx=-1, cy=-1), w=1/36: + // feq = (p + RHO*(u² + 3uv - u + v² - v)) / 12 + float feq6_target = (p_out + RHO*(u_target*u_target + 3.0f*u_target*v_target - u_target + v_target*v_target - v_target)) / 12.0f; + float feq6_neb = (p_out + RHO*(u_neb*u_neb + 3.0f*u_neb*v_neb - u_neb + v_neb*v_neb - v_neb)) / 12.0f; + + f[2] = f_neb[2] - feq2_neb + feq2_target; + f[8] = f_neb[8] - feq8_neb + feq8_target; + f[6] = f_neb[6] - feq6_neb + feq6_target; +} + +#endif // NQ == 9 + +// ============================================================================ +// D3Q19 inlet / outlet (non-equilibrium extrapolation) +// +// Parabolic inlet (x=0): reconstruct cx>0 populations i=1,7,9,13,15 +// Pressure outlet (x=NX-1): reconstruct cx<0 populations i=2,8,10,14,16 +// +// Uses generic feq computation from macro.cuh to avoid hand-expanded formulas. +// ============================================================================ +#if NQ == 19 + +__device__ inline void apply_parabolic_inlet_3d(float* __restrict__ f, + const float* __restrict__ f_neb, + float y_coord) +{ + // Neighbor macros + float rho_neb, un, vn, wn; + compute_rho_u(f_neb, rho_neb, un, vn, wn); + + // Target velocity (parabolic in y, uniform in z) + float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f); + float u_tar = U0 * 1.5f * (1.0f - 4.0f * yy * yy); + + // feq arrays + float feq_tar[19], feq_neb[19]; + compute_feq(rho_neb, u_tar, 0.0f, 0.0f, feq_tar); + compute_feq(rho_neb, un, vn, wn, feq_neb); + + // Reconstruct cx>0 directions: i = 1, 7, 9, 13, 15 + f[1] = f_neb[1] - feq_neb[1] + feq_tar[1]; + f[7] = f_neb[7] - feq_neb[7] + feq_tar[7]; + f[9] = f_neb[9] - feq_neb[9] + feq_tar[9]; + f[13] = f_neb[13] - feq_neb[13] + feq_tar[13]; + f[15] = f_neb[15] - feq_neb[15] + feq_tar[15]; +} + +__device__ inline void apply_pressure_outlet_3d(float* __restrict__ f, + const float* __restrict__ f_neb, + float y_coord) +{ + // Neighbor macros + float rho_neb, un, vn, wn; + compute_rho_u(f_neb, rho_neb, un, vn, wn); + + // Target: p_out = 0 gauge → rho = 1.0 (using neighbor velocity) + float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f); + float u_tar = U0 * 1.5f * (1.0f - 4.0f * yy * yy); + + float feq_tar[19], feq_neb[19]; + compute_feq(RHO, u_tar, 0.0f, 0.0f, feq_tar); + compute_feq(RHO, un, vn, wn, feq_neb); + + // Reconstruct cx<0 directions: i = 2, 8, 10, 14, 16 + f[2] = f_neb[2] - feq_neb[2] + feq_tar[2]; + f[8] = f_neb[8] - feq_neb[8] + feq_tar[8]; + f[10] = f_neb[10] - feq_neb[10] + feq_tar[10]; + f[14] = f_neb[14] - feq_neb[14] + feq_tar[14]; + f[16] = f_neb[16] - feq_neb[16] + feq_tar[16]; +} + +#endif // NQ == 19 + +#endif // CELERIS_BOUNDARY_INLET_OUTLET_CUH diff --git a/src/CelerisLab/lbm/kernels/core/descriptors.cuh b/src/CelerisLab/lbm/kernels/core/descriptors.cuh new file mode 100644 index 0000000..47398b1 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/core/descriptors.cuh @@ -0,0 +1,75 @@ +// CelerisLab – core/descriptors.cuh +// Lattice descriptor constants. +// Direction ordering: paired for Esoteric-Pull compatibility. +// i=0 rest; for i>=1, (i, i^1|1) are opposite pairs: +// opp(0) = 0 +// opp(i) = ((i-1) ^ 1) + 1 for i >= 1 +// +// D2Q9 : (0)rest (1,2)±x (3,4)±y (5,6)±(x+y) (7,8)±(x-y) +// D3Q19: (0)rest (1,2)±x (3,4)±y (5,6)±z +// (7,8)±(x+y) (9,10)±(x+z) (11,12)±(y+z) +// (13,14)±(x-y) (15,16)±(x-z) (17,18)±(y-z) +// ============================================================================ + +#ifndef CELERIS_CORE_DESCRIPTORS_CUH +#define CELERIS_CORE_DESCRIPTORS_CUH + +// NQ and DIM come from macros.h (included before this header) + +// --------------------------------------------------------------------------- +// D2Q9 +// --------------------------------------------------------------------------- +#if NQ == 9 + +__constant__ int d_cx[9] = { 0, 1, -1, 0, 0, 1, -1, 1, -1}; +__constant__ int d_cy[9] = { 0, 0, 0, 1, -1, 1, -1, -1, 1}; + +__constant__ float d_w[9] = { + 4.0f/9.0f, + 1.0f/9.0f, 1.0f/9.0f, 1.0f/9.0f, 1.0f/9.0f, + 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f +}; + +// Named weight macros (avoid __constant__ reads in hot paths) +#define W0_VAL (4.0f/9.0f) +#define WS_VAL (1.0f/9.0f) +#define WE_VAL (1.0f/36.0f) + +// --------------------------------------------------------------------------- +// D3Q19 +// --------------------------------------------------------------------------- +#elif NQ == 19 + +__constant__ int d_cx[19] = { 0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0}; +__constant__ int d_cy[19] = { 0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1, -1, 1, 0, 0, 1,-1}; +__constant__ int d_cz[19] = { 0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0, -1, 1, -1, 1}; + +__constant__ float d_w[19] = { + 1.0f/3.0f, + 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, 1.0f/18.0f, + 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, + 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f, 1.0f/36.0f +}; + +#define W0_VAL (1.0f/3.0f) +#define WS_VAL (1.0f/18.0f) +#define WE_VAL (1.0f/36.0f) + +#else +#error "Unsupported NQ. Use 9 (D2Q9) or 19 (D3Q19)." +#endif + +// --------------------------------------------------------------------------- +// Lattice sound speed +// --------------------------------------------------------------------------- +#define CS2 (1.0f / 3.0f) +#define CS2_INV (3.0f) + +// --------------------------------------------------------------------------- +// Opposite direction (works for any paired-ordered lattice) +// --------------------------------------------------------------------------- +__device__ __forceinline__ int opp_dir(int i) { + return (i == 0) ? 0 : (((i - 1) ^ 1) + 1); +} + +#endif // CELERIS_CORE_DESCRIPTORS_CUH diff --git a/src/CelerisLab/lbm/kernels/core/flags.cuh b/src/CelerisLab/lbm/kernels/core/flags.cuh new file mode 100644 index 0000000..7524651 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/core/flags.cuh @@ -0,0 +1,86 @@ +// CelerisLab – core/flags.cuh +// Cell flag definitions (uint32_t, layered bytes) and helper functions. +// ============================================================================ + +#ifndef CELERIS_CORE_FLAGS_CUH +#define CELERIS_CORE_FLAGS_CUH + +#include + +// --------------------------------------------------------------------------- +// Byte 0 [0:7] – basic cell type +// --------------------------------------------------------------------------- +#define FLAG_FLUID 0x01u +#define FLAG_SOLID 0x02u +#define FLAG_GAS 0x04u +#define FLAG_INTERFACE 0x08u +#define FLAG_SENSOR 0x10u + +// --------------------------------------------------------------------------- +// Byte 1 [8:15] – boundary condition type +// --------------------------------------------------------------------------- +#define FLAG_BB (0x01u << 8) // full-way bounce-back +#define FLAG_EQ_BC (0x02u << 8) // equilibrium (Dirichlet) boundary +#define FLAG_CURVED (0x04u << 8) // curved-boundary interpolation +#define FLAG_IBM (0x08u << 8) // immersed boundary marker +#define FLAG_MOVING (0x10u << 8) // moving-wall velocity + +// --------------------------------------------------------------------------- +// Byte 2 [16:23] – multi-GPU / AMR / ghost +// --------------------------------------------------------------------------- +#define FLAG_GHOST (0x01u << 16) +#define FLAG_HALO (0x02u << 16) +#define FLAG_AMR_FINE (0x04u << 16) +#define FLAG_AMR_COARSE (0x08u << 16) + +// --------------------------------------------------------------------------- +// Byte 3 [24:31] – reserved for extensions (phase-field, particles …) +// --------------------------------------------------------------------------- + +// --------------------------------------------------------------------------- +// Masks +// --------------------------------------------------------------------------- +#define MASK_CELL_TYPE 0x000000FFu +#define MASK_BC_TYPE 0x0000FF00u +#define MASK_MULTIGPU 0x00FF0000u +#define MASK_EXTENSION 0xFF000000u + +// --------------------------------------------------------------------------- +// Legacy compatibility (current driver.py uses uint8 with these bits) +// --------------------------------------------------------------------------- +#define LEGACY_FLUID 0x01 +#define LEGACY_SOLID 0x02 +#define LEGACY_GAS 0x04 +#define LEGACY_OBSTACLE 0x04 // obstacle / immersed body (triggers BB at adjacent fluid) +#define LEGACY_INTERFACE 0x08 +#define LEGACY_SENSOR 0x10 + +// --------------------------------------------------------------------------- +// Device helper functions +// --------------------------------------------------------------------------- +__device__ __forceinline__ bool is_fluid(uint32_t f) { return (f & FLAG_FLUID) != 0; } +__device__ __forceinline__ bool is_solid(uint32_t f) { return (f & FLAG_SOLID) != 0; } +__device__ __forceinline__ bool is_gas(uint32_t f) { return (f & FLAG_GAS) != 0; } +__device__ __forceinline__ bool is_interface(uint32_t f) { return (f & FLAG_INTERFACE) != 0; } +__device__ __forceinline__ bool is_sensor(uint32_t f) { return (f & FLAG_SENSOR) != 0; } + +__device__ __forceinline__ bool is_bb(uint32_t f) { return (f & FLAG_BB) != 0; } +__device__ __forceinline__ bool is_curved(uint32_t f) { return (f & FLAG_CURVED) != 0; } +__device__ __forceinline__ bool is_ibm(uint32_t f) { return (f & FLAG_IBM) != 0; } +__device__ __forceinline__ bool is_moving(uint32_t f) { return (f & FLAG_MOVING) != 0; } +__device__ __forceinline__ bool is_eq_bc(uint32_t f) { return (f & FLAG_EQ_BC) != 0; } + +__device__ __forceinline__ bool is_boundary(uint32_t f) { return (f & MASK_BC_TYPE) != 0; } +__device__ __forceinline__ bool is_wall(uint32_t f) { return (f & (FLAG_BB | FLAG_CURVED | FLAG_MOVING)) != 0; } + +// --------------------------------------------------------------------------- +// Legacy flag promotion (uint8 flags → uint32 with inferred BC bits) +// --------------------------------------------------------------------------- +__device__ __forceinline__ uint32_t promote_legacy_flag(uint8_t legacy) { + uint32_t f = (uint32_t)legacy & 0x1Fu; // copy low 5 bits + if ((f & FLAG_SOLID) && (f & FLAG_INTERFACE)) // SOLID + INTERFACE + f |= FLAG_CURVED; // → mark as curved BC + return f; +} + +#endif // CELERIS_CORE_FLAGS_CUH diff --git a/src/CelerisLab/lbm/kernels/core/layout.cuh b/src/CelerisLab/lbm/kernels/core/layout.cuh new file mode 100644 index 0000000..23db86b --- /dev/null +++ b/src/CelerisLab/lbm/kernels/core/layout.cuh @@ -0,0 +1,155 @@ +// CelerisLab – core/layout.cuh +// SoA memory layout: index_f, coordinates, and neighbor computation. +// Depends on: NQ, NX, NY, NZ from macros.h; descriptors.cuh for d_cx/d_cy. +// ============================================================================ + +#ifndef CELERIS_CORE_LAYOUT_CUH +#define CELERIS_CORE_LAYOUT_CUH + +// --------------------------------------------------------------------------- +// Total cell count (compile-time) +// --------------------------------------------------------------------------- +#define TOTAL_CELLS ((unsigned long)(NX) * (unsigned long)(NY) * (unsigned long)(NZ)) + +// --------------------------------------------------------------------------- +// SoA index: fi[ direction * N + node ] ("Q-major") +// --------------------------------------------------------------------------- +__device__ __forceinline__ unsigned long index_f(unsigned long n, unsigned int i) { + return (unsigned long)i * TOTAL_CELLS + n; +} + +// --------------------------------------------------------------------------- +// SoA index for DIM-component field: u[ component * N + node ] +// --------------------------------------------------------------------------- +__device__ __forceinline__ unsigned long index_u(unsigned long n, unsigned int d) { + return (unsigned long)d * TOTAL_CELLS + n; +} + +// --------------------------------------------------------------------------- +// 2D linear index → (x, y) and (x, y) → linear +// --------------------------------------------------------------------------- +#if DIM == 2 + +__device__ __forceinline__ void coordinates(unsigned long n, + unsigned int& x, + unsigned int& y) +{ + x = (unsigned int)(n % NX); + y = (unsigned int)(n / NX); +} + +__device__ __forceinline__ unsigned long linear_index(unsigned int x, + unsigned int y) +{ + return (unsigned long)y * NX + x; +} + +__device__ __forceinline__ void index_from_thread(unsigned int& x, + unsigned int& y, + unsigned long& k) +{ + x = threadIdx.x + NT * blockIdx.x; + y = blockIdx.y; + k = (unsigned long)y * NX + x; +} + +// --------------------------------------------------------------------------- +// Neighbor list (periodic BCs via modular arithmetic) +// --------------------------------------------------------------------------- +__device__ inline void compute_neighbors(unsigned long n, unsigned long* j) { + unsigned int x, y; + coordinates(n, x, y); + + unsigned int xp = (x + 1u) % NX; + unsigned int xm = (x + NX - 1u) % NX; + unsigned int yp = (y + 1u) % NY; + unsigned int ym = (y + NY - 1u) % NY; + + // i=0: self + j[0] = n; + +#if NQ == 9 + // paired: (1,2)±x (3,4)±y (5,6)±(x+y) (7,8)±(x-y) + j[1] = linear_index(xp, y ); // +x + j[2] = linear_index(xm, y ); // -x + j[3] = linear_index(x , yp); // +y + j[4] = linear_index(x , ym); // -y + j[5] = linear_index(xp, yp); // +x +y + j[6] = linear_index(xm, ym); // -x -y + j[7] = linear_index(xp, ym); // +x -y + j[8] = linear_index(xm, yp); // -x +y +#elif NQ == 19 + #error "D3Q19 requires DIM == 3. Set DIM=3 in macros.h." +#endif +} + +#elif DIM == 3 + +__device__ __forceinline__ void coordinates(unsigned long n, + unsigned int& x, + unsigned int& y, + unsigned int& z) +{ + x = (unsigned int)(n % NX); + y = (unsigned int)(n / NX % NY); + z = (unsigned int)(n / ((unsigned long)NX * NY)); +} + +__device__ __forceinline__ unsigned long linear_index(unsigned int x, + unsigned int y, + unsigned int z) +{ + return (unsigned long)z * NY * NX + (unsigned long)y * NX + x; +} + +__device__ __forceinline__ void index_from_thread(unsigned int& x, + unsigned int& y, + unsigned int& z, + unsigned long& k) +{ + x = threadIdx.x + NT * blockIdx.x; + y = blockIdx.y; + z = blockIdx.z; + k = linear_index(x, y, z); +} + +__device__ inline void compute_neighbors(unsigned long n, unsigned long* j) { + unsigned int x, y, z; + coordinates(n, x, y, z); + + unsigned int xp = (x + 1u) % NX; + unsigned int xm = (x + NX - 1u) % NX; + unsigned int yp = (y + 1u) % NY; + unsigned int ym = (y + NY - 1u) % NY; + unsigned int zp = (z + 1u) % NZ; + unsigned int zm = (z + NZ - 1u) % NZ; + + j[0] = n; + +#if NQ == 19 + // ±x, ±y, ±z (straight) + j[1] = linear_index(xp, y, z ); // +x + j[2] = linear_index(xm, y, z ); // -x + j[3] = linear_index(x, yp, z ); // +y + j[4] = linear_index(x, ym, z ); // -y + j[5] = linear_index(x, y, zp); // +z + j[6] = linear_index(x, y, zm); // -z + // diagonal + j[7] = linear_index(xp, yp, z ); // +x+y + j[8] = linear_index(xm, ym, z ); // -x-y + j[9] = linear_index(xp, y, zp); // +x+z + j[10] = linear_index(xm, y, zm); // -x-z + j[11] = linear_index(x, yp, zp); // +y+z + j[12] = linear_index(x, ym, zm); // -y-z + j[13] = linear_index(xp, ym, z ); // +x-y + j[14] = linear_index(xm, yp, z ); // -x+y + j[15] = linear_index(xp, y, zm); // +x-z + j[16] = linear_index(xm, y, zp); // -x+z + j[17] = linear_index(x, yp, zm); // +y-z + j[18] = linear_index(x, ym, zp); // -y+z +#endif +} + +#endif // DIM == 3 + +#endif // CELERIS_CORE_LAYOUT_CUH diff --git a/src/CelerisLab/lbm/kernels/core/params.cuh b/src/CelerisLab/lbm/kernels/core/params.cuh new file mode 100644 index 0000000..fbe3d30 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/core/params.cuh @@ -0,0 +1,62 @@ +// CelerisLab – core/params.cuh +// Runtime parameter structures transported via __constant__ memory. +// Includes: LBMParams, RigidBodyState2D, RigidBodyControl2D. +// ============================================================================ + +#ifndef CELERIS_CORE_PARAMS_CUH +#define CELERIS_CORE_PARAMS_CUH + +#include + +// ============================================================================ +// LBM runtime parameters (uploaded once per run or on parameter change) +// ============================================================================ +struct LBMParams { + //--- grid --- + unsigned int Nx, Ny, Nz; + unsigned long N; // Nx * Ny * Nz + + //--- relaxation --- + float omega; // SRT/TRT主松弛率 w = 1/(3ν + 0.5) + float omega_bulk; // MRT bulk 松弛率 (s_e / s_eps) + + //--- external body force --- + float fx, fy, fz; + + //--- reference quantities --- + float rho_ref; // 参考密度 (通常 1.0) + float u_inlet; // 入口参考速度 + + //--- diagnostics --- + unsigned int n_objects; // 观测对象数量 +}; + +__constant__ LBMParams d_params; + +// ============================================================================ +// Rigid-body state / control (2D, per architecture §17.3) +// ============================================================================ +struct RigidBodyState2D { + float x, y, theta; // 位姿 (position + orientation) + float vx, vy, omega; // 速度 (translational + rotational) +}; + +struct RigidBodyControl2D { + float ax_cmd, ay_cmd; // 目标加速度（或目标速度/位移, 取决于 mode） + float alpha_cmd; // 目标角加速度 / 角速度 + int control_mode; // 0=力控, 1=速度控, 2=位移控 +}; + +// ============================================================================ +// Halo plan (multi-GPU / AMR 预留) +// ============================================================================ +struct HaloPlan { + int face_id; + int peer_gpu; + int level; // AMR level (0 = base) + size_t count; + int* send_idx; + int* recv_idx; +}; + +#endif // CELERIS_CORE_PARAMS_CUH diff --git a/src/CelerisLab/lbm/kernels/core/precision.cuh b/src/CelerisLab/lbm/kernels/core/precision.cuh new file mode 100644 index 0000000..dfd9d17 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/core/precision.cuh @@ -0,0 +1,91 @@ +// CelerisLab – core/precision.cuh +// Storage precision abstraction (FP32 / FP16S / FP16C). +// Compute precision is always float (FP32). +// +// Controlled by macros.h: +// STORE_PRECISION 0 = FP32 (default) +// 1 = FP16S (IEEE-754 half + ×32768 scaling) +// 2 = FP16C (custom 1-4-11 format, higher precision) +// ============================================================================ + +#ifndef CELERIS_CORE_PRECISION_CUH +#define CELERIS_CORE_PRECISION_CUH + +#ifndef STORE_PRECISION +#define STORE_PRECISION 0 +#endif + +// ============================================================================ +// FP16S – Scaled IEEE-754 Half (hardware-accelerated conversion) +// ============================================================================ +#if STORE_PRECISION == 1 + +#include +using fpxx = __half; + +__device__ __forceinline__ float load_ddf(const fpxx* p, unsigned long idx) { + return __half2float(p[idx]) * 3.0517578125e-5f; // ÷ 32768 +} +__device__ __forceinline__ void store_ddf(fpxx* p, unsigned long idx, float v) { + p[idx] = __float2half_rn(v * 32768.0f); +} + +#define FPXX_BYTES 2 + +// ============================================================================ +// FP16C – Custom 1-4-11 format (software conversion, ~3.6 dec. digits) +// ============================================================================ +#elif STORE_PRECISION == 2 + +using fpxx = unsigned short; + +__device__ __forceinline__ float half_custom_to_float(unsigned short x) { + unsigned int e = (x & 0x7800) >> 11; + unsigned int m = (x & 0x07FF) << 12; + unsigned int v = __float_as_uint((float)m) >> 23; + return __uint_as_float( + (x & 0x8000) << 16 | + (e != 0) * ((e + 112) << 23 | m) | + ((e == 0) & (m != 0)) * ((v - 37) << 23 | ((m << (150 - v)) & 0x007FF000)) + ); +} + +__device__ __forceinline__ unsigned short float_to_half_custom(float x) { + unsigned int b = __float_as_uint(x) + 0x00000800; // round-to-nearest-even + unsigned int e = (b & 0x7F800000) >> 23; + unsigned int m = b & 0x007FFFFF; + return (unsigned short)( + (b & 0x80000000) >> 16 | + (e > 112) * ((((e - 112) << 11) & 0x7800) | m >> 12) | + ((e < 113) & (e > 100)) * ((((0x007FF800 + m) >> (124 - e)) + 1) >> 1) + ); +} + +__device__ __forceinline__ float load_ddf(const fpxx* p, unsigned long idx) { + return half_custom_to_float(p[idx]); +} +__device__ __forceinline__ void store_ddf(fpxx* p, unsigned long idx, float v) { + p[idx] = float_to_half_custom(v); +} + +#define FPXX_BYTES 2 + +// ============================================================================ +// FP32 – No conversion (default) +// ============================================================================ +#else + +using fpxx = float; + +__device__ __forceinline__ float load_ddf(const fpxx* p, unsigned long idx) { + return p[idx]; +} +__device__ __forceinline__ void store_ddf(fpxx* p, unsigned long idx, float v) { + p[idx] = v; +} + +#define FPXX_BYTES 4 + +#endif // STORE_PRECISION + +#endif // CELERIS_CORE_PRECISION_CUH diff --git a/src/CelerisLab/lbm/kernels/kernel_v2.cu b/src/CelerisLab/lbm/kernels/kernel_v2.cu new file mode 100644 index 0000000..199eefb --- /dev/null +++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu @@ -0,0 +1,299 @@ +// CelerisLab – kernel_v2.cu +// ============================================================================ +// Modular compilation entry point (Stage 1 Architecture) +// +// This file includes all modular headers and step kernels, +// and exports extern "C" functions for PyCUDA. +// +// Controlled by macros.h: +// NQ = 9 (D2Q9) or 19 (D3Q19) +// COLLISION_MODEL = 0 (SRT) / 1 (TRT) / 2 (MRT) +// STREAMING_MODEL = 0 (double-buffer) / 1 (Esoteric-Pull) +// STORE_PRECISION = 0 (FP32) / 1 (FP16S) / 2 (FP16C) +// USE_DDF_SHIFTING= 0 / 1 +// ============================================================================ + +#include +#include +#include + +// --------------------------------------------------------------------------- +// Layer 0: Configuration (compile-time) +// --------------------------------------------------------------------------- +#include "macros.h" + +// --------------------------------------------------------------------------- +// Layer 1: Core primitives +// --------------------------------------------------------------------------- +#include "core/precision.cuh" +#include "core/flags.cuh" +#include "core/descriptors.cuh" +#include "core/layout.cuh" +#include "core/params.cuh" + +// --------------------------------------------------------------------------- +// Layer 2: Operators +// --------------------------------------------------------------------------- +#include "operators/macro.cuh" +#include "operators/collision_srt.cuh" +#include "operators/collision_trt.cuh" +#include "operators/collision_mrt.cuh" +#include "operators/forcing_guo.cuh" + +// --------------------------------------------------------------------------- +// Layer 3: Streaming +// --------------------------------------------------------------------------- +#include "streaming/pull_double_buffer.cuh" +#include "streaming/esopull_single_buffer.cuh" + +// --------------------------------------------------------------------------- +// Layer 4: Boundary conditions +// --------------------------------------------------------------------------- +#include "boundary/bounce_back.cuh" +#include "boundary/curved_boundary.cuh" +#include "boundary/inlet_outlet.cuh" + +// --------------------------------------------------------------------------- +// Layer 5: Step kernels (selected by STREAMING_MODEL) +// --------------------------------------------------------------------------- +#include "step/one_step_double.cu" +#include "step/one_step_esopull.cu" + +// --------------------------------------------------------------------------- +// Layer 6: IBM kernels (always available, launched separately) +// --------------------------------------------------------------------------- +#include "boundary/ibm_kernels.cuh" + +// ============================================================================ +// Extern "C" wrappers (for PyCUDA / ctypes) +// +// Naming: +// OneStep – backward-compatible main step (double-buffer) +// InitTubeFlow_v2 – channel initialization (new ordering) +// StreamCollide – new API main step (double or esopull) +// CurvedBoundary – post-step curved boundary correction +// UpdateFields – recompute rho/u from DDF +// InitEsoPull – Esoteric-Pull initialization +// IBM_UpdateLag – IBM Lagrangian point update +// IBM_Eul2Lag – Euler → Lagrangian interpolation +// IBM_Lag2Eul – Lagrangian → Euler spreading +// ============================================================================ +extern "C" +{ + +// ----- Backward-compatible init (uses new direction ordering) ----- +__global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi) +{ +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; + + float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1)) + * ((float)y - 0.5f * (NY - 1)) + / ((float)(NY - 2) * (float)(NY - 2))); + + if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) { + flag[k] = LEGACY_SOLID; + for (int i = 0; i < NQ; i++) + store_ddf(fi, index_f(k, (unsigned int)i), 0.0f); + } else { + flag[k] = (uint8_t)LEGACY_FLUID; + for (int i = 0; i < NQ; i++) { + float cu = (float)d_cx[i] * u_init; + float val = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init); +#if USE_DDF_SHIFTING + val -= d_w[i]; +#endif + store_ddf(fi, index_f(k, (unsigned int)i), val); + } + } +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; + + float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1)) + * ((float)y - 0.5f * (NY - 1)) + / ((float)(NY - 2) * (float)(NY - 2))); + + if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) { + flag[k] = LEGACY_SOLID; + for (int i = 0; i < NQ; i++) + store_ddf(fi, index_f(k, (unsigned int)i), 0.0f); + } else { + flag[k] = (uint8_t)LEGACY_FLUID; + for (int i = 0; i < NQ; i++) { + float cu = (float)d_cx[i] * u_init; + float val = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init); +#if USE_DDF_SHIFTING + val -= d_w[i]; +#endif + store_ddf(fi, index_f(k, (unsigned int)i), val); + } + } +#endif +} + +// ----- Main step (double-buffer) ----- +// Signature compatible with driver.py: flag, fi_in, fi_out, indx, delta, action, obs +__global__ void OneStep( + uint8_t* flag, + fpxx* fi_in, + fpxx* fi_out, + int32_t* indx, + float* delta, + float* action, + float* obs) +{ + // Redirect to the modular StreamCollideDouble kernel body + // (We inline here to keep a single extern "C" entry point.) + +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; +#endif + + uint8_t fl = flag[k]; + + unsigned long j[NQ]; + compute_neighbors(k, j); + + float f[NQ]; + +#if STREAMING_MODEL == 0 + stream_pull_load(k, f, fi_in, j); +#else + // Esoteric-Pull for OneStep compat would need timestep — not available + // in legacy signature. Use double-buffer mode only for OneStep. + stream_pull_load(k, f, fi_in, j); +#endif + + float rho_n, ux, uy; +#if NQ == 9 + compute_rho_u(f, rho_n, ux, uy); + + // Inlet / outlet / wall bounce-back + if (fl & LEGACY_SOLID) { + bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); + if (x == 0 && interior_y) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x + 1u, y); + for (int i = 0; i < NQ; i++) + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + apply_parabolic_inlet(f, f_neb, (float)y); + } + else if (x == (unsigned int)(NX - 1) && interior_y) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x - 1u, y); + for (int i = 0; i < NQ; i++) + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + apply_pressure_outlet(f, f_neb, (float)y); + } + else { + // Wall / corner: full-way bounce-back (swap all pairs) + #pragma unroll + for (int i = 1; i < NQ; i += 2) { + float t = f[i]; f[i] = f[i+1]; f[i+1] = t; + } + } + } + + // Collision + if (fl & LEGACY_FLUID) { + float feq[NQ], Fin[NQ]; + compute_feq(rho_n, ux, uy, feq); + zero_forcing(Fin); + +#if COLLISION_MODEL == 0 + collide_srt(f, feq, Fin, d_params.omega); +#elif COLLISION_MODEL == 1 + collide_trt(f, feq, Fin, d_params.omega); +#elif COLLISION_MODEL == 2 + collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega); +#endif + } + + stream_pull_store(k, f, fi_out); + + // Curved boundary processing + if ((fl & LEGACY_SOLID) && (fl & LEGACY_INTERFACE)) { + int id_off = indx[k]; + int id_obj = *reinterpret_cast(&delta[id_off]); + float Uw = action[id_obj] * delta[id_off + NQ]; + float Vw = action[id_obj] * delta[id_off + NQ + 1]; + float* obs_fx = &obs[DIM * id_obj]; + float* obs_fy = &obs[DIM * id_obj + 1]; + apply_curved_boundary(k, x, y, fi_out, delta, id_off, Uw, Vw, obs_fx, obs_fy); + } + + // Sensor + if (fl & LEGACY_SENSOR) { + int id_obj = indx[k]; + atomicAdd(&obs[DIM * id_obj], ux); + atomicAdd(&obs[DIM * id_obj + 1], uy); + } + +#elif NQ == 19 + // ---- Macroscopic (computed for all nodes; only used for fluid collision) ---- + float uz; + compute_rho_u(f, rho_n, ux, uy, uz); + + // ---- Boundary conditions ---- + if (fl & LEGACY_SOLID) { + // Interior of inlet/outlet planes (skip corners at y=0/NY-1) + bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); + if (x == 0 && interior_y) { + // Parabolic velocity inlet + float f_neb[NQ]; + unsigned long k_neb = linear_index(x + 1u, y, z); + for (int i = 0; i < NQ; i++) + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + apply_parabolic_inlet_3d(f, f_neb, (float)y); + } + else if (x == (unsigned int)(NX - 1) && interior_y) { + // Pressure outlet + float f_neb[NQ]; + unsigned long k_neb = linear_index(x - 1u, y, z); + for (int i = 0; i < NQ; i++) + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + apply_pressure_outlet_3d(f, f_neb, (float)y); + } + else { + // Wall nodes and corners: full-way bounce-back. + #pragma unroll + for (int i = 1; i < NQ; i += 2) { + float t = f[i]; f[i] = f[i+1]; f[i+1] = t; + } + } + } + // Obstacle nodes: full-way bounce-back (same swap) + if (fl & LEGACY_OBSTACLE) { + #pragma unroll + for (int i = 1; i < NQ; i += 2) { + float t = f[i]; f[i] = f[i+1]; f[i+1] = t; + } + } + + // ---- Collision (fluid only) ---- + if (fl & LEGACY_FLUID) { + float feq[NQ], Fin[NQ]; + compute_feq(rho_n, ux, uy, uz, feq); + zero_forcing(Fin); + collide_srt(f, feq, Fin, d_params.omega); + } + + stream_pull_store(k, f, fi_out); +#endif +} + +} // extern "C" diff --git a/src/CelerisLab/lbm/kernels/macros.h b/src/CelerisLab/lbm/kernels/macros.h index 7a2466d..dda58fc 100644 --- a/src/CelerisLab/lbm/kernels/macros.h +++ b/src/CelerisLab/lbm/kernels/macros.h @@ -3,23 +3,23 @@ // cuda parameters #define MULT_GPU False #define NT 128 -#define X_1U 128 -#define Y_1U 32 -#define Z_1U 1 +#define X_1U 256 +#define Y_1U 128 +#define Z_1U 32 // flow parameters #define LBtype float -#define UX 10 -#define UY 16 +#define UX 1 +#define UY 1 #define UZ 1 -#define NX 1280 -#define NY 512 -#define NZ 1 -#define DIM 2 -#define NQ 9 -#define VIS 0.004 +#define NX 256 +#define NY 128 +#define NZ 32 +#define DIM 3 +#define NQ 19 +#define VIS 0.0096000000 #define RHO 1.0 -#define U0 0.01 +#define U0 0.04 // constants #define PI 3.141592653589793238 @@ -33,5 +33,31 @@ #define V_TAYLOR 0b00000001 // variables -#define N_OBJS 7 -// #define N_SENS 2 \ No newline at end of file +#define N_OBJS 0 +// #define N_SENS 2 + +// ============================================================================ +// New architecture configuration (Stage 1) +// These defaults are safe for backward compatibility. +// compiler.py can override any of them via modify_macro(). +// ============================================================================ + +// Collision model: 0=SRT, 1=TRT, 2=MRT +#ifndef COLLISION_MODEL +#define COLLISION_MODEL 0 +#endif + +// Streaming model: 0=double-buffer, 1=esoteric-pull +#ifndef STREAMING_MODEL +#define STREAMING_MODEL 0 +#endif + +// Storage precision: 0=FP32, 1=FP16S, 2=FP16C +#ifndef STORE_PRECISION +#define STORE_PRECISION 0 +#endif + +// DDF-shifting: 0=off, 1=on +#ifndef USE_DDF_SHIFTING +#define USE_DDF_SHIFTING 0 +#endif \ No newline at end of file diff --git a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh new file mode 100644 index 0000000..5aab699 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh @@ -0,0 +1,135 @@ +// CelerisLab – operators/collision_mrt.cuh +// MRT (Multiple-Relaxation-Time) collision operator. +// +// D2Q9 Lallemand-Luo (Phys. Rev. E 61, 2000) with NEW paired ordering: +// cx = {0, 1,-1, 0, 0, 1,-1, 1,-1} +// cy = {0, 0, 0, 1,-1, 1,-1,-1, 1} +// +// Moment space: +// m[0] = ρ (conserved, s₀ = 0 or 1.0 overwrite) +// m[1] = e (energy, s₁ = s_e) +// m[2] = ε (energy², s₂ = s_eps) +// m[3] = jx (conserved, s₃ = 0 or 1.0 overwrite) +// m[4] = qx (energy flux x, s₄ = s_q) +// m[5] = jy (conserved, s₅ = 0 or 1.0 overwrite) +// m[6] = qy (energy flux y, s₆ = s_q) +// m[7] = pxx (stress, s₇ = s_nu = ω = 1/(3ν + 0.5)) +// m[8] = pxy (stress, s₈ = s_nu) +// ============================================================================ + +#ifndef CELERIS_OPERATORS_COLLISION_MRT_CUH +#define CELERIS_OPERATORS_COLLISION_MRT_CUH + +#if NQ == 9 + +// --------------------------------------------------------------------------- +// D2Q9 MRT (fully expanded, no loops, no matrix storage) +// --------------------------------------------------------------------------- +__device__ __forceinline__ void collide_mrt(float* __restrict__ g, + float rho, float ux, float uy, + const float* __restrict__ Fin, + float omega) +{ + // ----- Relaxation rates ----- + // s_rho = s_jx = s_jy = 1.0 (conserved moments → overwrite to eq) + // s_e = s_eps = s_q = 1.2 (kinetic transport) + // s_nu = omega (viscosity-related) + const float s_rho = 1.0f; // conserved moment relaxation + const float s_e = 1.2f; + const float s_eps = 1.2f; + const float s_jx = 1.0f; + const float s_q = 1.2f; + const float s_jy = 1.0f; + const float s_nu = omega; + + // ----- Pressure (from density) ----- + const float p = (g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8]) / 3.0f; + + // ----- Forward transform: m = M * g (new paired ordering) ----- + float m[9]; + m[0] = g[0] + g[1] + g[2] + g[3] + g[4] + g[5] + g[6] + g[7] + g[8]; + m[1] = -4*g[0] - g[1] - g[2] - g[3] - g[4] + 2*g[5] + 2*g[6] + 2*g[7] + 2*g[8]; + m[2] = 4*g[0] - 2*g[1] - 2*g[2] - 2*g[3] - 2*g[4] + g[5] + g[6] + g[7] + g[8]; + m[3] = g[1] - g[2] + g[5] - g[6] + g[7] - g[8]; + m[4] = -2*g[1] + 2*g[2] + g[5] - g[6] + g[7] - g[8]; + m[5] = g[3] - g[4] + g[5] - g[6] - g[7] + g[8]; + m[6] = -2*g[3] + 2*g[4] + g[5] - g[6] - g[7] + g[8]; + m[7] = g[1] + g[2] - g[3] - g[4]; + m[8] = g[5] + g[6] - g[7] - g[8]; + + // ----- Equilibrium moments ----- + const float u2 = ux * ux + uy * uy; + float meq[9]; + meq[0] = 3.0f * p; // ρ + meq[1] = -6.0f * p + 3.0f * RHO * u2; // e + meq[2] = 3.0f * p - 3.0f * RHO * u2; // ε + meq[3] = RHO * ux; // jx + meq[4] = -RHO * ux; // qx + meq[5] = RHO * uy; // jy + meq[6] = -RHO * uy; // qy + meq[7] = RHO * (ux * ux - uy * uy); // pxx + meq[8] = RHO * ux * uy; // pxy + + // ----- Relaxation: delta_m[i] = s_i * (meq[i] - m[i]) ----- + float dm[9]; + dm[0] = s_rho * (meq[0] - m[0]); + dm[1] = s_e * (meq[1] - m[1]); + dm[2] = s_eps * (meq[2] - m[2]); + dm[3] = s_jx * (meq[3] - m[3]); + dm[4] = s_q * (meq[4] - m[4]); + dm[5] = s_jy * (meq[5] - m[5]); + dm[6] = s_q * (meq[6] - m[6]); + dm[7] = s_nu * (meq[7] - m[7]); + dm[8] = s_nu * (meq[8] - m[8]); + + // ----- Inverse transform: g += M⁻¹ * dm (new paired ordering) ----- + g[0] += ( dm[0] - dm[1] + dm[2] ) / 9.0f; + g[1] += (4 * dm[0] - dm[1] - 2* dm[2] + 6* dm[3] - 6* dm[4] + 9*dm[7]) / 36.0f; + g[2] += (4 * dm[0] - dm[1] - 2* dm[2] - 6* dm[3] + 6* dm[4] + 9*dm[7]) / 36.0f; + g[3] += (4 * dm[0] - dm[1] - 2* dm[2] + 6*dm[5] - 6*dm[6] - 9*dm[7]) / 36.0f; + g[4] += (4 * dm[0] - dm[1] - 2* dm[2] - 6*dm[5] + 6*dm[6] - 9*dm[7]) / 36.0f; + g[5] += (4 * dm[0] + 2* dm[1] + dm[2] + 6* dm[3] + 3* dm[4] + 6*dm[5] + 3*dm[6] + 9*dm[8]) / 36.0f; + g[6] += (4 * dm[0] + 2* dm[1] + dm[2] - 6* dm[3] - 3* dm[4] - 6*dm[5] - 3*dm[6] + 9*dm[8]) / 36.0f; + g[7] += (4 * dm[0] + 2* dm[1] + dm[2] + 6* dm[3] + 3* dm[4] - 6*dm[5] - 3*dm[6] - 9*dm[8]) / 36.0f; + g[8] += (4 * dm[0] + 2* dm[1] + dm[2] - 6* dm[3] - 3* dm[4] + 6*dm[5] + 3*dm[6] - 9*dm[8]) / 36.0f; + + // ----- Add forcing (if present) ----- + #pragma unroll + for (int i = 0; i < 9; i++) { + g[i] += Fin[i]; + } +} + +// Convenience wrapper: no external forcing +__device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g, + float rho, float ux, float uy, + float omega) +{ + float Fin[9] = {0}; + collide_mrt(g, rho, ux, uy, Fin, omega); +} + +#elif NQ == 19 + +// --------------------------------------------------------------------------- +// D3Q19 MRT (placeholder – use SRT or TRT until fully validated) +// --------------------------------------------------------------------------- +__device__ __forceinline__ void collide_mrt(float* __restrict__ g, + float rho, float ux, float uy, float uz, + const float* __restrict__ Fin, + float omega) +{ + // TODO: implement D3Q19 MRT with 19-moment Gram–Schmidt basis + // Fall back to SRT for now + float feq[19]; + compute_feq(rho, ux, uy, uz, feq); + const float c_tau = 1.0f - 0.5f * omega; + #pragma unroll + for (int i = 0; i < 19; i++) { + g[i] = fmaf(1.0f - omega, g[i], fmaf(omega, feq[i], c_tau * Fin[i])); + } +} + +#endif // NQ + +#endif // CELERIS_OPERATORS_COLLISION_MRT_CUH diff --git a/src/CelerisLab/lbm/kernels/operators/collision_srt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_srt.cuh new file mode 100644 index 0000000..1c4b8aa --- /dev/null +++ b/src/CelerisLab/lbm/kernels/operators/collision_srt.cuh @@ -0,0 +1,35 @@ +// CelerisLab – operators/collision_srt.cuh +// SRT (BGK) collision operator. +// f_out[i] = (1 - ω) f[i] + ω feq[i] + Fin[i] +// ============================================================================ + +#ifndef CELERIS_OPERATORS_COLLISION_SRT_CUH +#define CELERIS_OPERATORS_COLLISION_SRT_CUH + +__device__ __forceinline__ void collide_srt(float* __restrict__ f, + const float* __restrict__ feq, + const float* __restrict__ Fin, + float omega) +{ + // Pre-compute forcing prefactor (1 - ω/2) for Guo scheme + const float c_tau = 1.0f - 0.5f * omega; + + #pragma unroll + for (int i = 0; i < NQ; i++) { + f[i] = fmaf(1.0f - omega, f[i], + fmaf(omega, feq[i], c_tau * Fin[i])); + } +} + +// Variant without forcing (Fin = 0) +__device__ __forceinline__ void collide_srt_no_force(float* __restrict__ f, + const float* __restrict__ feq, + float omega) +{ + #pragma unroll + for (int i = 0; i < NQ; i++) { + f[i] = fmaf(1.0f - omega, f[i], omega * feq[i]); + } +} + +#endif // CELERIS_OPERATORS_COLLISION_SRT_CUH diff --git a/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh new file mode 100644 index 0000000..1d0ad1e --- /dev/null +++ b/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh @@ -0,0 +1,79 @@ +// CelerisLab – operators/collision_trt.cuh +// TRT (Two-Relaxation-Time) collision operator. +// +// Symmetric part uses ω⁺ (= omega, transport viscosity) +// Antisymmetric uses ω⁻ (linked by magic parameter Λ = 3/16) +// ω⁻ = 1 / (Λ / (1/ω⁺ - 0.5) + 0.5) +// +// Paired direction ordering makes TRT natural: +// f_s[i] = 0.5*(f[i] + f[opp(i)]) symmetric +// f_a[i] = 0.5*(f[i] - f[opp(i)]) antisymmetric +// ============================================================================ + +#ifndef CELERIS_OPERATORS_COLLISION_TRT_CUH +#define CELERIS_OPERATORS_COLLISION_TRT_CUH + +// Magic parameter Λ = 3/16 (optimal for porous-media / bounce-back wall location) +#define TRT_MAGIC_PARAM (0.1875f) + +__device__ __forceinline__ float compute_omega_minus(float omega_plus) { + return 1.0f / (TRT_MAGIC_PARAM / (1.0f / omega_plus - 0.5f) + 0.5f); +} + +__device__ __forceinline__ void collide_trt(float* __restrict__ f, + const float* __restrict__ feq, + const float* __restrict__ Fin, + float omega) +{ + const float wp = omega; + const float wm = compute_omega_minus(wp); + + // Direction 0: rest particle – treat as symmetric-only + f[0] = f[0] - wp * (f[0] - feq[0]) + Fin[0]; + + // Direction pairs (i, i+1) for i = 1, 3, 5, … + #pragma unroll + for (int i = 1; i < NQ; i += 2) { + const int ib = i + 1; // opposite (paired layout) + + // Current and opposite + const float fi = f[i]; + const float fb = f[ib]; + + // Equilibrium + const float feqi = feq[i]; + const float feqb = feq[ib]; + + // Symmetric / antisymmetric non-equilibrium + const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb)); + const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb)); + + // Relax + forcing + f[i] = fi - wp * delta_s - wm * delta_a + Fin[i]; + f[ib] = fb - wp * delta_s + wm * delta_a + Fin[ib]; + } +} + +// Variant without forcing (Fin = 0) +__device__ __forceinline__ void collide_trt_no_force(float* __restrict__ f, + const float* __restrict__ feq, + float omega) +{ + const float wp = omega; + const float wm = compute_omega_minus(wp); + + f[0] = f[0] - wp * (f[0] - feq[0]); + + #pragma unroll + for (int i = 1; i < NQ; i += 2) { + const int ib = i + 1; + const float fi = f[i], fb = f[ib]; + const float feqi = feq[i], feqb = feq[ib]; + const float delta_s = 0.5f * ((fi - feqi) + (fb - feqb)); + const float delta_a = 0.5f * ((fi - feqi) - (fb - feqb)); + f[i] = fi - wp * delta_s - wm * delta_a; + f[ib] = fb - wp * delta_s + wm * delta_a; + } +} + +#endif // CELERIS_OPERATORS_COLLISION_TRT_CUH diff --git a/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh b/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh new file mode 100644 index 0000000..f337425 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/operators/forcing_guo.cuh @@ -0,0 +1,102 @@ +// CelerisLab – operators/forcing_guo.cuh +// Guo forcing scheme (Guo, Zheng & Shi, Phys. Rev. E 65, 2002). +// +// Velocity correction: u* = u + F·dt / (2ρ) +// Forcing term: +// Fin[i] = (1 - ω/2) · w_i · [ (c_i - u*)/cs² + (c_i · u*)/cs⁴ · c_i ] · F +// simplified: Fin[i] = 9·w_i · [ (c_i·F)(c_i·u* + 1/3) - u*·F/3 ] +// ============================================================================ + +#ifndef CELERIS_OPERATORS_FORCING_GUO_CUH +#define CELERIS_OPERATORS_FORCING_GUO_CUH + +// --------------------------------------------------------------------------- +// D2Q9 Guo forcing terms +// --------------------------------------------------------------------------- +#if NQ == 9 + +__device__ __forceinline__ void apply_guo_velocity_correction( + float& ux, float& uy, + float fx, float fy, float rho) +{ + float rho2 = 0.5f / rho; + ux = fmaf(fx, rho2, ux); + uy = fmaf(fy, rho2, uy); +} + +__device__ __forceinline__ void compute_guo_forcing( + float ux, float uy, // velocity (already corrected: u*) + float fx, float fy, // force density + float* __restrict__ Fin) +{ + const float uF = -0.33333334f * fmaf(ux, fx, uy * fy); // -u*·F / 3 + + Fin[0] = 9.0f * W0_VAL * uF; + + // Straight directions + // cx = {0, 1,-1, 0, 0, 1,-1, 1,-1} + // cy = {0, 0, 0, 1,-1, 1,-1,-1, 1} + Fin[1] = 9.0f * WS_VAL * fmaf( fx, ux + 0.33333334f, uF); // +x + Fin[2] = 9.0f * WS_VAL * fmaf(-fx, -ux + 0.33333334f, uF); // -x + Fin[3] = 9.0f * WS_VAL * fmaf( fy, uy + 0.33333334f, uF); // +y + Fin[4] = 9.0f * WS_VAL * fmaf(-fy, -uy + 0.33333334f, uF); // -y + + // Diagonal directions + const float f5x = fx, f5y = fy; // +x+y + const float c5u = ux + uy; + const float c5F = f5x + f5y; + Fin[5] = 9.0f * WE_VAL * fmaf(c5F, c5u + 0.33333334f, uF); + + const float c6u = -ux - uy; // -x-y + const float c6F = -fx - fy; + Fin[6] = 9.0f * WE_VAL * fmaf(c6F, c6u + 0.33333334f, uF); + + const float c7u = ux - uy; // +x-y + const float c7F = fx - fy; + Fin[7] = 9.0f * WE_VAL * fmaf(c7F, c7u + 0.33333334f, uF); + + const float c8u = -ux + uy; // -x+y + const float c8F = -fx + fy; + Fin[8] = 9.0f * WE_VAL * fmaf(c8F, c8u + 0.33333334f, uF); +} + +#elif NQ == 19 + +__device__ __forceinline__ void apply_guo_velocity_correction( + float& ux, float& uy, float& uz, + float fx, float fy, float fz, float rho) +{ + float rho2 = 0.5f / rho; + ux = fmaf(fx, rho2, ux); + uy = fmaf(fy, rho2, uy); + uz = fmaf(fz, rho2, uz); +} + +__device__ __forceinline__ void compute_guo_forcing( + float ux, float uy, float uz, + float fx, float fy, float fz, + float* __restrict__ Fin) +{ + const float uF = -0.33333334f * fmaf(ux, fx, fmaf(uy, fy, uz * fz)); + + Fin[0] = 9.0f * W0_VAL * uF; + + // Generic loop for D3Q19 (compiler will fully unroll) + for (int i = 1; i < 19; i++) { + float ci_dot_F = d_cx[i]*fx + d_cy[i]*fy + d_cz[i]*fz; + float ci_dot_u = d_cx[i]*ux + d_cy[i]*uy + d_cz[i]*uz; + Fin[i] = 9.0f * d_w[i] * fmaf(ci_dot_F, ci_dot_u + 0.33333334f, uF); + } +} + +#endif // NQ + +// --------------------------------------------------------------------------- +// Zero forcing helper (when no external force) +// --------------------------------------------------------------------------- +__device__ __forceinline__ void zero_forcing(float* __restrict__ Fin) { + #pragma unroll + for (int i = 0; i < NQ; i++) Fin[i] = 0.0f; +} + +#endif // CELERIS_OPERATORS_FORCING_GUO_CUH diff --git a/src/CelerisLab/lbm/kernels/operators/macro.cuh b/src/CelerisLab/lbm/kernels/operators/macro.cuh new file mode 100644 index 0000000..b820bbe --- /dev/null +++ b/src/CelerisLab/lbm/kernels/operators/macro.cuh @@ -0,0 +1,168 @@ +// CelerisLab – operators/macro.cuh +// Macroscopic quantity computation (rho, u) and equilibrium distribution (feq). +// Supports DDF-shifting (USE_DDF_SHIFTING=1) and standard mode (=0). +// +// D2Q9 paired ordering: +// cx = {0, 1,-1, 0, 0, 1,-1, 1,-1} +// cy = {0, 0, 0, 1,-1, 1,-1,-1, 1} +// ============================================================================ + +#ifndef CELERIS_OPERATORS_MACRO_CUH +#define CELERIS_OPERATORS_MACRO_CUH + +#ifndef USE_DDF_SHIFTING +#define USE_DDF_SHIFTING 0 +#endif + +// --------------------------------------------------------------------------- +// compute_rho_u: f[NQ] → (rho, ux, uy [, uz]) +// --------------------------------------------------------------------------- +#if NQ == 9 + +__device__ __forceinline__ void compute_rho_u(const float* __restrict__ f, + float& rho, + float& ux, float& uy) +{ + // Density: sum of all populations + rho = f[0] + f[1] + f[2] + f[3] + f[4] + f[5] + f[6] + f[7] + f[8]; +#if USE_DDF_SHIFTING + rho += 1.0f; // DDF-shifting: stored f_tilde = f - w, sum(w) = 1 +#endif + + // Momentum (works identically with or without DDF-shifting, + // because sum(w_i * cx_i) = 0) + float inv_rho = 1.0f / rho; + // cx: 0 1 -1 0 0 1 -1 1 -1 + ux = (f[1] - f[2] + f[5] - f[6] + f[7] - f[8]) * inv_rho; + // cy: 0 0 0 1 -1 1 -1 -1 1 + uy = (f[3] - f[4] + f[5] - f[6] - f[7] + f[8]) * inv_rho; +} + +#elif NQ == 19 + +__device__ __forceinline__ void compute_rho_u(const float* __restrict__ f, + float& rho, + float& ux, float& uy, float& uz) +{ + rho = f[0]; + for (int i = 1; i < 19; i++) rho += f[i]; +#if USE_DDF_SHIFTING + rho += 1.0f; +#endif + + float inv_rho = 1.0f / rho; + ux = (f[1]-f[2] + f[7]-f[8] + f[9]-f[10] + f[13]-f[14] + f[15]-f[16]) * inv_rho; + uy = (f[3]-f[4] + f[7]-f[8] + f[11]-f[12] + f[14]-f[13] + f[17]-f[18]) * inv_rho; + uz = (f[5]-f[6] + f[9]-f[10] + f[11]-f[12] + f[16]-f[15] + f[18]-f[17]) * inv_rho; +} + +#endif // NQ + +// --------------------------------------------------------------------------- +// compute_feq: (rho, u) → feq[NQ] +// --------------------------------------------------------------------------- +#if NQ == 9 + +__device__ __forceinline__ void compute_feq(float rho, float ux, float uy, + float* __restrict__ feq) +{ +#if USE_DDF_SHIFTING + const float rhom1 = rho - 1.0f; // perturbation density +#endif + const float ux3 = 3.0f * ux; + const float uy3 = 3.0f * uy; + const float c3 = -1.5f * (ux * ux + uy * uy); // = -|u|²/(2cs²) + + // Helper: feq_i = w_i * rho * (1 + 3 ci·u + 4.5 (ci·u)² - 1.5 u²) + // With DDF-shifting: feq_tilde_i = feq_i - w_i + // = w_i * (rho * (1 + ...) - 1) = w_i * (rho*(...) + rho - 1) + + #define FEQ_CORE(cu) (rho * fmaf(0.5f * 3.0f, (cu) * (cu), (cu) + c3)) + +#if USE_DDF_SHIFTING + #define FEQ(w, cu) ((w) * fmaf(rho, fmaf(0.5f, (cu)*(cu), (cu) + c3), rhom1)) +#else + #define FEQ(w, cu) ((w) * rho * (1.0f + (cu) + 0.5f * (cu) * (cu) + c3)) +#endif + + // i=0: rest + feq[0] = FEQ(W0_VAL, 0.0f); + + // (1,2) ±x + feq[1] = FEQ(WS_VAL, ux3); + feq[2] = FEQ(WS_VAL, -ux3); + + // (3,4) ±y + feq[3] = FEQ(WS_VAL, uy3); + feq[4] = FEQ(WS_VAL, -uy3); + + // Pre-compute diagonal velocity combos + const float u0 = ux3 + uy3; // +x+y + const float u1 = ux3 - uy3; // +x-y + + // (5,6) ±(x+y) + feq[5] = FEQ(WE_VAL, u0); + feq[6] = FEQ(WE_VAL, -u0); + + // (7,8) ±(x-y) + feq[7] = FEQ(WE_VAL, u1); + feq[8] = FEQ(WE_VAL, -u1); + + #undef FEQ + #undef FEQ_CORE +} + +#elif NQ == 19 + +__device__ __forceinline__ void compute_feq(float rho, float ux, float uy, float uz, + float* __restrict__ feq) +{ +#if USE_DDF_SHIFTING + const float rhom1 = rho - 1.0f; +#endif + const float ux3 = 3.0f * ux; + const float uy3 = 3.0f * uy; + const float uz3 = 3.0f * uz; + const float c3 = -1.5f * (ux*ux + uy*uy + uz*uz); + +#if USE_DDF_SHIFTING + #define FEQ(w, cu) ((w) * fmaf(rho, fmaf(0.5f, (cu)*(cu), (cu) + c3), rhom1)) +#else + #define FEQ(w, cu) ((w) * rho * (1.0f + (cu) + 0.5f * (cu) * (cu) + c3)) +#endif + + feq[0] = FEQ(W0_VAL, 0.0f); + + // Straight + feq[1] = FEQ(WS_VAL, ux3); feq[2] = FEQ(WS_VAL, -ux3); + feq[3] = FEQ(WS_VAL, uy3); feq[4] = FEQ(WS_VAL, -uy3); + feq[5] = FEQ(WS_VAL, uz3); feq[6] = FEQ(WS_VAL, -uz3); + + // Diagonal + const float u0 = ux3+uy3, u1 = ux3+uz3, u2 = uy3+uz3; + const float u3 = ux3-uy3, u4 = ux3-uz3, u5 = uy3-uz3; + + feq[7] = FEQ(WE_VAL, u0); feq[8] = FEQ(WE_VAL, -u0); + feq[9] = FEQ(WE_VAL, u1); feq[10] = FEQ(WE_VAL, -u1); + feq[11] = FEQ(WE_VAL, u2); feq[12] = FEQ(WE_VAL, -u2); + feq[13] = FEQ(WE_VAL, u3); feq[14] = FEQ(WE_VAL, -u3); + feq[15] = FEQ(WE_VAL, u4); feq[16] = FEQ(WE_VAL, -u4); + feq[17] = FEQ(WE_VAL, u5); feq[18] = FEQ(WE_VAL, -u5); + + #undef FEQ +} + +#endif // NQ + +// --------------------------------------------------------------------------- +// compute_pressure (diagnostic, from state equation p = cs² ρ) +// --------------------------------------------------------------------------- +__device__ __forceinline__ float compute_pressure(float rho) { + return CS2 * rho; +} + +__device__ __forceinline__ float compute_pressure_perturbation(float rho, float rho_ref) { + return CS2 * (rho - rho_ref); +} + +#endif // CELERIS_OPERATORS_MACRO_CUH diff --git a/src/CelerisLab/lbm/kernels/step/one_step_double.cu b/src/CelerisLab/lbm/kernels/step/one_step_double.cu new file mode 100644 index 0000000..960c0f8 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/one_step_double.cu @@ -0,0 +1,332 @@ +// CelerisLab – step/one_step_double.cu +// Main LBM step kernel using double-buffer PULL streaming. +// +// Workflow per node: +// 1. Pull-stream: gather DDF from neighbors → local f[NQ] +// 2. Compute macroscopic quantities (ρ, u) +// 3. Apply boundary conditions (inlet/outlet/wall) +// 4. Compute equilibrium + forcing +// 5. Collision (SRT / TRT / MRT) +// 6. Write to output buffer +// +// Kernel signature maintains backward compatibility with driver.py: +// flag* – cell flags (uint8 for legacy, uint32 for new) +// fi_in* – input DDF array (fpxx) +// fi_out* – output DDF array (fpxx) +// indx* – per-cell offset into delta pool +// delta* – curved boundary parameter pool +// action* – per-object control input +// obs* – per-object observation output (force / velocity) +// ============================================================================ + +#ifndef CELERIS_STEP_ONE_STEP_DOUBLE_CU +#define CELERIS_STEP_ONE_STEP_DOUBLE_CU + +__global__ void StreamCollideDouble( + const uint8_t* __restrict__ flag, + const fpxx* __restrict__ fi_in, + fpxx* __restrict__ fi_out, + const int32_t* __restrict__ indx, + const float* __restrict__ delta, + const float* __restrict__ action, + float* __restrict__ obs, + float* __restrict__ rho_arr, // macro output (can be NULL) + float* __restrict__ u_arr) // macro output (can be NULL) +{ + // ----- Thread → node mapping (2D) ----- +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; +#endif + + // ----- Read flag ----- + uint8_t fl = flag[k]; + + // ----- Neighbor indices ----- + unsigned long j[NQ]; + compute_neighbors(k, j); + + // ----- Pull-stream: load DDF from neighbors ----- + float f[NQ]; + stream_pull_load(k, f, fi_in, j); + + // ----- Compute macroscopic quantities ----- + float rho_n, ux, uy; +#if NQ == 9 + compute_rho_u(f, rho_n, ux, uy); +#elif NQ == 19 + float uz; + compute_rho_u(f, rho_n, ux, uy, uz); +#endif + + // ----- Boundary conditions (inlet/outlet/wall) ----- +#if NQ == 9 + if (fl & LEGACY_SOLID) { + if (x == 0) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x + 1u, y); + for (int i = 0; i < NQ; i++) { + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + } + apply_parabolic_inlet(f, f_neb, (float)y); + } + else if (x == NX - 1) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x - 1u, y); + for (int i = 0; i < NQ; i++) { + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + } + apply_pressure_outlet(f, f_neb, (float)y); + } + } +#elif NQ == 19 + if (fl & LEGACY_SOLID) { + if (x == 0) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x + 1u, y, z); + for (int i = 0; i < NQ; i++) { + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + } + apply_parabolic_inlet_3d(f, f_neb, (float)y); + } + else if (x == NX - 1) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x - 1u, y, z); + for (int i = 0; i < NQ; i++) { + f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); + } + apply_pressure_outlet_3d(f, f_neb, (float)y); + } + } +#endif + + // ----- Wall bounce-back (top/bottom walls) ----- +#if NQ == 9 + if (y == 1 || y == NY - 2) { + apply_wall_bb_d2q9(y, f); + } +#elif NQ == 19 + if (y == 1 || y == NY - 2) { + apply_wall_bb_d3q19_y(y, f); + } +#endif + + // ----- Compute equilibrium ----- + if (fl & LEGACY_FLUID) { + float feq[NQ]; + float Fin[NQ]; + +#if NQ == 9 + compute_feq(rho_n, ux, uy, feq); + zero_forcing(Fin); + + // ----- Collision ----- +#if COLLISION_MODEL == 0 + collide_srt(f, feq, Fin, d_params.omega); +#elif COLLISION_MODEL == 1 + collide_trt(f, feq, Fin, d_params.omega); +#elif COLLISION_MODEL == 2 + collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega); +#endif + +#elif NQ == 19 + compute_feq(rho_n, ux, uy, uz, feq); + zero_forcing(Fin); + collide_srt(f, feq, Fin, d_params.omega); +#endif + } + + // ----- Write to output buffer ----- + stream_pull_store(k, f, fi_out); + + // ----- Write macroscopic fields (optional) ----- + if (rho_arr != nullptr) { + rho_arr[k] = rho_n; + } + if (u_arr != nullptr) { +#if DIM == 2 + u_arr[k] = ux; + u_arr[TOTAL_CELLS + k] = uy; +#elif DIM == 3 + u_arr[k] = ux; + u_arr[TOTAL_CELLS + k] = uy; + u_arr[2 * TOTAL_CELLS + k] = uz; +#endif + } + + // ----- Sensor observation ----- + if (fl & LEGACY_SENSOR) { + int id_obj = indx[k]; + atomicAdd(&obs[DIM * id_obj], ux); + atomicAdd(&obs[DIM * id_obj + 1], uy); + } +} + +// --------------------------------------------------------------------------- +// Curved boundary post-processing kernel (separate launch) +// Processes SOLID+INTERFACE nodes → modifies fi_out at neighbors. +// --------------------------------------------------------------------------- +__global__ void CurvedBoundaryPost( + const uint8_t* __restrict__ flag, + fpxx* __restrict__ fi_out, + const int32_t* __restrict__ indx, + const float* __restrict__ delta, + const float* __restrict__ action, + float* __restrict__ obs) +{ +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; +#endif + + uint8_t fl = flag[k]; + + if ((fl & LEGACY_SOLID) && (fl & LEGACY_INTERFACE)) { +#if DIM == 2 + int id_off = indx[k]; + int id_obj = *reinterpret_cast(&delta[id_off]); + + // Wall velocity from action + normal + float Uw = action[id_obj] * delta[id_off + NQ]; + float Vw = action[id_obj] * delta[id_off + NQ + 1]; + + float* obs_fx = (obs != nullptr) ? &obs[DIM * id_obj] : nullptr; + float* obs_fy = (obs != nullptr) ? &obs[DIM * id_obj + 1] : nullptr; + + apply_curved_boundary(k, x, y, fi_out, delta, id_off, Uw, Vw, obs_fx, obs_fy); +#endif // DIM == 2 curved BC + } +} + +// --------------------------------------------------------------------------- +// InitTubeFlow: Initialize parabolic tube flow (backward compatible) +// --------------------------------------------------------------------------- +__global__ void InitTubeFlow(uint8_t* flag, fpxx* fi) +{ +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; + + float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1)) + * ((float)y - 0.5f * (NY - 1)) + / ((float)(NY - 2) * (float)(NY - 2))); + + if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) { + flag[k] = LEGACY_SOLID; + for (int i = 0; i < NQ; i++) { + store_ddf(fi, index_f(k, (unsigned int)i), 0.0f); + } + } else { + flag[k] = (uint8_t)LEGACY_FLUID; + // feq with rho=RHO, u=(u_init, 0) + for (int i = 0; i < NQ; i++) { + float cu = (float)d_cx[i] * u_init; // cy=0 for parabolic flow + float val = d_w[i] * RHO * (1.0f + 3.0f * cu + 4.5f * cu * cu - 1.5f * u_init * u_init); +#if USE_DDF_SHIFTING + val -= d_w[i]; // store f_tilde = f - w +#endif + store_ddf(fi, index_f(k, (unsigned int)i), val); + } + } +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; + + float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1)) + * ((float)y - 0.5f * (NY - 1)) + / ((float)(NY - 2) * (float)(NY - 2))); + + // Walls: y=0, y=NY-1, x=0, x=NX-1; z is periodic + if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) { + flag[k] = LEGACY_SOLID; + for (int i = 0; i < NQ; i++) { + store_ddf(fi, index_f(k, (unsigned int)i), 0.0f); + } + } else { + flag[k] = (uint8_t)LEGACY_FLUID; + for (int i = 0; i < NQ; i++) { + float cu = (float)d_cx[i] * u_init; + float val = d_w[i] * RHO * (1.0f + 3.0f * cu + 4.5f * cu * cu - 1.5f * u_init * u_init); +#if USE_DDF_SHIFTING + val -= d_w[i]; +#endif + store_ddf(fi, index_f(k, (unsigned int)i), val); + } + } +#endif +} + +// --------------------------------------------------------------------------- +// UpdateMacro: Recompute rho/u from DDF (for diagnostics/output) +// --------------------------------------------------------------------------- +__global__ void UpdateMacro( + const fpxx* __restrict__ fi, + float* __restrict__ rho_arr, + float* __restrict__ u_arr, + const uint8_t* __restrict__ flag) +{ +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; + + if (flag[k] & LEGACY_SOLID) return; + + unsigned long j[NQ]; + compute_neighbors(k, j); + + float f[NQ]; + // For double-buffer, just read from the current buffer (no streaming) + for (int i = 0; i < NQ; i++) { + f[i] = load_ddf(fi, index_f(k, (unsigned int)i)); + } + + float rho_n, ux, uy; + compute_rho_u(f, rho_n, ux, uy); + + rho_arr[k] = rho_n; + u_arr[k] = ux; + u_arr[TOTAL_CELLS + k] = uy; +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; + + if (flag[k] & LEGACY_SOLID) return; + + float f[NQ]; + for (int i = 0; i < NQ; i++) { + f[i] = load_ddf(fi, index_f(k, (unsigned int)i)); + } + + float rho_n, ux, uy, uz; + compute_rho_u(f, rho_n, ux, uy, uz); + + rho_arr[k] = rho_n; + u_arr[k] = ux; + u_arr[TOTAL_CELLS + k] = uy; + u_arr[2 * TOTAL_CELLS + k] = uz; +#endif +} + +#endif // CELERIS_STEP_ONE_STEP_DOUBLE_CU diff --git a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu new file mode 100644 index 0000000..907fd63 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu @@ -0,0 +1,198 @@ +// CelerisLab – step/one_step_esopull.cu +// Main LBM step kernel using Esoteric-Pull single-buffer streaming. +// +// Workflow per node: +// 1. load_f_esopull: read DDF (streaming second half) +// 2. Compute macroscopic quantities (ρ, u) +// 3. Boundary conditions +// 4. Equilibrium + forcing + collision +// 5. store_f_esopull: write DDF (streaming first half) +// +// Compared to double-buffer: +// - Uses HALF the memory (single fi array) +// - Timestep counter t is required +// - fi_alt pointer is unused (can be NULL) +// ============================================================================ + +#ifndef CELERIS_STEP_ONE_STEP_ESOPULL_CU +#define CELERIS_STEP_ONE_STEP_ESOPULL_CU + +__global__ void StreamCollideEsoPull( + fpxx* __restrict__ fi, // single DDF buffer (read/write) + const uint8_t* __restrict__ flag, + const int32_t* __restrict__ indx, + const float* __restrict__ delta, + const float* __restrict__ action, + float* __restrict__ obs, + float* __restrict__ rho_arr, + float* __restrict__ u_arr, + const float* __restrict__ force_field, // Euler force field [DIM*N], or NULL + unsigned long t) // current timestep +{ + // ----- Thread → node mapping ----- +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; +#elif DIM == 3 + unsigned int x, y, z; + unsigned long k; + index_from_thread(x, y, z, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return; +#endif + + uint8_t fl = flag[k]; + + // Skip pure solid / gas + if ((fl & LEGACY_SOLID) && !(fl & LEGACY_INTERFACE) && !(fl & LEGACY_SENSOR)) { + // For inlet/outlet solid nodes, we still process below + if (x != 0 && x != (unsigned int)(NX - 1)) return; + } + + // ----- Neighbor indices ----- + unsigned long j[NQ]; + compute_neighbors(k, j); + + // ----- Esoteric-Pull: load DDF (streaming second half) ----- + float f[NQ]; + load_f_esopull(k, f, fi, j, t); + + // ----- Compute macroscopic quantities ----- + float rho_n, ux, uy; +#if NQ == 9 + compute_rho_u(f, rho_n, ux, uy); +#elif NQ == 19 + float uz; + compute_rho_u(f, rho_n, ux, uy, uz); +#endif + + // ----- Boundary conditions ----- +#if NQ == 9 + if (fl & LEGACY_SOLID) { + if (x == 0) { + float f_neb[NQ]; + unsigned long k_neb = linear_index(x + 1u, y); + load_f_esopull(k_neb, f_neb, fi, j, t); // approximate: reuse j + // More accurate: compute j_neb separately + unsigned long j_neb[NQ]; + compute_neighbors(k_neb, j_neb); + load_f_esopull(k_neb, f_neb, fi, j_neb, t); + apply_parabolic_inlet(f, f_neb, (float)y); + } + else if (x == (unsigned int)(NX - 1)) { + unsigned long k_neb = linear_index(x - 1u, y); + unsigned long j_neb[NQ]; + compute_neighbors(k_neb, j_neb); + float f_neb[NQ]; + load_f_esopull(k_neb, f_neb, fi, j_neb, t); + apply_pressure_outlet(f, f_neb, (float)y); + } + } + + if (y == 1 || y == (unsigned int)(NY - 2)) { + apply_wall_bb_d2q9(y, f); + } +#endif + + // ----- Forcing ----- + float Fin[NQ]; + float fxn = d_params.fx, fyn = d_params.fy; + + if (force_field != nullptr) { + fxn += force_field[k]; + fyn += force_field[TOTAL_CELLS + k]; + } + + if (fxn != 0.0f || fyn != 0.0f) { +#if NQ == 9 + apply_guo_velocity_correction(ux, uy, fxn, fyn, rho_n); + compute_guo_forcing(ux, uy, fxn, fyn, Fin); +#endif + } else { + zero_forcing(Fin); + } + + // ----- Clamp velocity for stability ----- + const float CS_LIMIT = 0.57735027f; // 1/sqrt(3) + ux = fminf(fmaxf(ux, -CS_LIMIT), CS_LIMIT); + uy = fminf(fmaxf(uy, -CS_LIMIT), CS_LIMIT); + + // ----- Collision ----- + if (fl & LEGACY_FLUID) { + float feq[NQ]; +#if NQ == 9 + compute_feq(rho_n, ux, uy, feq); + +#if COLLISION_MODEL == 0 + collide_srt(f, feq, Fin, d_params.omega); +#elif COLLISION_MODEL == 1 + collide_trt(f, feq, Fin, d_params.omega); +#elif COLLISION_MODEL == 2 + collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega); +#endif + +#elif NQ == 19 + float fzn = d_params.fz; + if (force_field != nullptr) fzn += force_field[2*TOTAL_CELLS + k]; + apply_guo_velocity_correction(ux, uy, uz, fxn, fyn, fzn, rho_n); + compute_feq(rho_n, ux, uy, uz, feq); + collide_srt(f, feq, Fin, d_params.omega); +#endif + } + + // ----- Esoteric-Pull: write DDF (streaming first half) ----- + store_f_esopull(k, f, fi, j, t); + + // ----- Write macroscopic fields ----- + if (rho_arr != nullptr) rho_arr[k] = rho_n; + if (u_arr != nullptr) { + u_arr[k] = ux; + u_arr[TOTAL_CELLS + k] = uy; + } + + // ----- Sensor observation ----- + if (fl & LEGACY_SENSOR) { + int id_obj = indx[k]; + atomicAdd(&obs[DIM * id_obj], ux); + atomicAdd(&obs[DIM * id_obj + 1], uy); + } +} + +// --------------------------------------------------------------------------- +// InitializeEsoPull: Initialize DDF via Esoteric-Pull store (t=1) +// --------------------------------------------------------------------------- +__global__ void InitializeEsoPull(uint8_t* flag, fpxx* fi) +{ +#if DIM == 2 + unsigned int x, y; + unsigned long k; + index_from_thread(x, y, k); + if (x >= (unsigned int)NX || y >= (unsigned int)NY) return; + + unsigned long j[NQ]; + compute_neighbors(k, j); + + float feq[NQ]; + if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) { + flag[k] = LEGACY_SOLID; + for (int i = 0; i < NQ; i++) feq[i] = 0.0f; + } else { + flag[k] = (uint8_t)LEGACY_FLUID; + float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f*(NY-1)) + * ((float)y - 0.5f*(NY-1)) + / ((float)(NY-2) * (float)(NY-2))); + for (int i = 0; i < NQ; i++) { + float cu = (float)d_cx[i] * u_init; + feq[i] = d_w[i] * RHO * (1.0f + 3.0f*cu + 4.5f*cu*cu - 1.5f*u_init*u_init); +#if USE_DDF_SHIFTING + feq[i] -= d_w[i]; +#endif + } + } + // Write via Esoteric-Pull store at t=1 (odd step) for correct alignment + store_f_esopull(k, feq, fi, j, 1ul); +#endif +} + +#endif // CELERIS_STEP_ONE_STEP_ESOPULL_CU diff --git a/src/CelerisLab/lbm/kernels/streaming/esopull_single_buffer.cuh b/src/CelerisLab/lbm/kernels/streaming/esopull_single_buffer.cuh new file mode 100644 index 0000000..fc3da4b --- /dev/null +++ b/src/CelerisLab/lbm/kernels/streaming/esopull_single_buffer.cuh @@ -0,0 +1,79 @@ +// CelerisLab – streaming/esopull_single_buffer.cuh +// Esoteric-Pull single-buffer streaming (Geier & Schönherr, 2017). +// +// Uses ONE DDF array (half the memory of double-buffer). +// Relies on paired direction ordering: (i, i+1) are opposites. +// +// Key idea: load_f and store_f alternate read/write positions +// based on timestep parity, so they never conflict in a single step. +// +// Even step (t%2 == 0): +// Read: f[i] ← fi[n, i+1] (local, reverse slot) +// f[i+1] ← fi[j[i], i ] (neighbor, forward slot) +// Write: fi[j[i], i+1] ← f[i] (neighbor, reverse slot) +// fi[n, i ] ← f[i+1] (local, forward slot) +// +// Odd step (t%2 == 1): +// Read: f[i] ← fi[n, i ] (local, forward slot) +// f[i+1] ← fi[j[i], i+1] (neighbor, reverse slot) +// Write: fi[j[i], i+1] ← f[i] (neighbor, forward slot → stored in reverse) +// fi[n, i ] ← f[i+1] (local, reverse slot → stored in forward) +// +// Wait, let me use the exact FluidX3D convention: +// Odd step read: f[i] ← fi[n, i], f[i+1] ← fi[j[i], i+1] +// Odd step write: fi[j[i], i+1] ← f[i], fi[n, i] ← f[i+1] +// ============================================================================ + +#ifndef CELERIS_STREAMING_ESOPULL_SINGLE_BUFFER_CUH +#define CELERIS_STREAMING_ESOPULL_SINGLE_BUFFER_CUH + +// --------------------------------------------------------------------------- +// load_f_esopull: Esoteric-Pull read phase (streaming second half) +// --------------------------------------------------------------------------- +__device__ inline void load_f_esopull(unsigned long n, + float* __restrict__ f, + const fpxx* __restrict__ fi, + const unsigned long* __restrict__ j, + unsigned long t) +{ + // i=0: rest particle always from self + f[0] = load_ddf(fi, index_f(n, 0u)); + + for (unsigned int i = 1; i < NQ; i += 2) { + if (t & 1) { + // Odd step + f[i] = load_ddf(fi, index_f(n, i)); + f[i + 1] = load_ddf(fi, index_f(j[i], i + 1)); + } else { + // Even step + f[i] = load_ddf(fi, index_f(n, i + 1)); + f[i + 1] = load_ddf(fi, index_f(j[i], i)); + } + } +} + +// --------------------------------------------------------------------------- +// store_f_esopull: Esoteric-Pull write phase (streaming first half) +// --------------------------------------------------------------------------- +__device__ inline void store_f_esopull(unsigned long n, + const float* __restrict__ f, + fpxx* __restrict__ fi, + const unsigned long* __restrict__ j, + unsigned long t) +{ + store_ddf(fi, index_f(n, 0u), f[0]); + + for (unsigned int i = 1; i < NQ; i += 2) { + if (t & 1) { + // Odd step + store_ddf(fi, index_f(j[i], i + 1), f[i]); + store_ddf(fi, index_f(n, i), f[i + 1]); + } else { + // Even step + store_ddf(fi, index_f(j[i], i), f[i]); + store_ddf(fi, index_f(n, i + 1), f[i + 1]); + } + } +} + +#endif // CELERIS_STREAMING_ESOPULL_SINGLE_BUFFER_CUH diff --git a/src/CelerisLab/lbm/kernels/streaming/pull_double_buffer.cuh b/src/CelerisLab/lbm/kernels/streaming/pull_double_buffer.cuh new file mode 100644 index 0000000..413ae3c --- /dev/null +++ b/src/CelerisLab/lbm/kernels/streaming/pull_double_buffer.cuh @@ -0,0 +1,56 @@ +// CelerisLab – streaming/pull_double_buffer.cuh +// Double-buffer PULL streaming. +// f_local[i] = f_in[ neighbor(n, -c_i), i ] +// f_out[n, i] = f_collided[i] (after collision, written to second buffer) +// +// This replaces the old PUSH scheme where each thread writes to neighbors. +// Pull streaming provides better memory coalescing on GPU. +// ============================================================================ + +#ifndef CELERIS_STREAMING_PULL_DOUBLE_BUFFER_CUH +#define CELERIS_STREAMING_PULL_DOUBLE_BUFFER_CUH + +// --------------------------------------------------------------------------- +// stream_pull_load: gather DDF from neighboring nodes into local f[NQ] +// n – current node +// f – output: local distribution (compute precision, float) +// fi – input: source DDF array (storage precision, fpxx) +// j – neighbor indices (j[i] = neighbor in direction i) +// +// Pull semantic: f[i] at node n = fi[ j[opp(i)], i ] +// i.e. the population traveling in direction i at node n +// arrived from the node at n - c_i = neighbor in direction opp(i) +// --------------------------------------------------------------------------- +__device__ __forceinline__ void stream_pull_load(unsigned long n, + float* __restrict__ f, + const fpxx* __restrict__ fi, + const unsigned long* __restrict__ j) +{ + // i=0 (rest): always from self + f[0] = load_ddf(fi, index_f(n, 0u)); + + // For paired ordering, opp(i) swaps within each pair: + // opp(1)=2, opp(2)=1, opp(3)=4, opp(4)=3, ... + #pragma unroll + for (int i = 1; i < NQ; i += 2) { + // population i arrives from direction opp(i) = i+1 + f[i] = load_ddf(fi, index_f(j[i+1], (unsigned int)i)); + // population i+1 arrives from direction opp(i+1) = i + f[i + 1] = load_ddf(fi, index_f(j[i], (unsigned int)(i + 1))); + } +} + +// --------------------------------------------------------------------------- +// stream_pull_store: write collided f[NQ] to the OUTPUT buffer at node n +// --------------------------------------------------------------------------- +__device__ __forceinline__ void stream_pull_store(unsigned long n, + const float* __restrict__ f, + fpxx* __restrict__ fi_out) +{ + #pragma unroll + for (int i = 0; i < NQ; i++) { + store_ddf(fi_out, index_f(n, (unsigned int)i), f[i]); + } +} + +#endif // CELERIS_STREAMING_PULL_DOUBLE_BUFFER_CUH diff --git a/tests/test_d2q9_visual.py b/tests/test_d2q9_visual.py new file mode 100644 index 0000000..8c96c84 --- /dev/null +++ b/tests/test_d2q9_visual.py @@ -0,0 +1,362 @@ +#!/usr/bin/env python3 +""" +D2Q9 Regression Test — Poiseuille Channel + Cylinder Flow +========================================================== +Uses the ORIGINAL kernel.cu with same grid / BCs. +Produces matplotlib figures for visual validation. + +Usage: + python tests/test_d2q9_visual.py --device 2 + python tests/test_d2q9_visual.py --device 2 --cylinder +""" + +import sys, os, argparse, time + +# Ensure CelerisLab package is importable +sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src')) + +import numpy as np +import pycuda.driver as cuda + +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt +from matplotlib.colors import Normalize + +from CelerisLab.cuda import compiler +from CelerisLab.common import preprocess as preproc + +# ━━━━━━━━━━━━━━━━━━━━━━━ Configuration ━━━━━━━━━━━━━━━━━━━━━━━ +NX, NY = 1280, 512 +NQ = 9 +NT = 128 +DIM = 2 +VIS = 0.002 +U0 = 0.01 +RHO = 1.0 +TOTAL = NX * NY + +# Original direction vectors (const.h ordering) +E = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]], dtype=np.int32) +OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6], dtype=np.int32) + +FLUID_FLAG = 0b00000001 +SOLID_FLAG = 0b00000010 +INTERFACE_FLAG = 0b00001000 + +# ━━━━━━━━━━━━━━━━━━━━━━━ Helpers ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ + +def configure_macros(n_objs=0): + """Write macros.h to match original kernel settings.""" + lines = compiler.read_lines(compiler.kernel_path("macros.h")) + defs = { + 'MULT_GPU': 'False', 'NT': NT, + 'X_1U': 128, 'Y_1U': 32, 'Z_1U': 1, + 'LBtype': 'float', + 'UX': 10, 'UY': 16, 'UZ': 1, + 'NX': NX, 'NY': NY, 'NZ': 1, + 'DIM': DIM, 'NQ': NQ, + 'VIS': VIS, 'RHO': f'{RHO}', 'U0': U0, + 'N_OBJS': n_objs, + } + for name, val in defs.items(): + lines = compiler.modify_macro(lines, name, val) + compiler.write_lines(compiler.kernel_path("macros.h"), lines) + + +def extract_fields(ddf_host): + """Compute rho, u, v from host DDF (original const.h ordering). + + The original kernel uses DDF-shifting: stores f_shifted = f - w_i*RHO. + So sum(f_shifted) = rho - RHO (~0 for incompressible flow), + and momentum = sum(e_x * f_shifted) works because sum(w_i * e_x) = 0. + """ + f = ddf_host.reshape(NQ, NY, NX) + rho = np.sum(f, axis=0) + RHO # un-shift: rho = sum(f_shifted) + RHO + u = (f[1] + f[5] + f[8] - f[3] - f[6] - f[7]) / RHO + v = (f[2] + f[5] + f[6] - f[4] - f[7] - f[8]) / RHO + return rho, u, v + + +def analytical_poiseuille(y_arr): + """Analytical parabolic profile matching InitTubeFlow.""" + yy = (y_arr - 0.5 * (NY - 1)) / (NY - 2.0) + return U0 * 1.5 * (1 - 4 * yy**2) + + +def build_cylinder_data(cx, cy, radius): + """Replicate driver.py add_cylinder logic for flag / delta / indx.""" + flag = np.ones(TOTAL, dtype=np.uint8) # init all FLUID + indx = np.zeros(TOTAL, dtype=np.int32) + delta_list = [] + index_offset = 0 + + # Build Poiseuille flag first (walls + solid borders) + for y in range(NY): + for x in range(NX): + k = x + y * NX + if y == 0 or y == NY - 1 or x == 0 or x == NX - 1: + flag[k] = SOLID_FLAG + + # Add cylinder + for x in range(int(cx - radius) - 1, int(cx + radius) + 1): + for y in range(int(cy - radius) - 1, int(cy + radius) + 1): + if (x - cx)**2 + (y - cy)**2 < radius**2: + k = x + y * NX + flag[k] = SOLID_FLAG + dt = np.zeros(11, dtype=np.float32) + dt[0] = np.int32(0).view(np.float32) # id_object = 0 + has_interface = False + for i in range(NQ): + xn = x + E[i][0] + yn = y + E[i][1] + if (xn - cx)**2 + (yn - cy)**2 >= radius**2: + has_interface = True + xi, yi = preproc.find_circle_intersection( + x, y, xn, yn, cx, cy, radius) + d_neb = np.sqrt((xi - xn)**2 + (yi - yn)**2) + e_len = np.sqrt(E[i][0]**2 + E[i][1]**2) + if e_len > 0: + dt[i] = d_neb / e_len + if has_interface: + flag[k] |= INTERFACE_FLAG + dt[9] = (cy - y) / radius + dt[10] = (x - cx) / radius + indx[k] = index_offset + delta_list.append(dt) + index_offset += 11 + + delta = np.concatenate(delta_list) if delta_list else np.zeros(1, dtype=np.float32) + return flag, indx, delta + + +# ━━━━━━━━━━━━━━━━━━━━━━━ Simulation ━━━━━━━━━━━━━━━━━━━━━━━━━━ + +def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host): + """Compile kernel, run LBM, return DDF on host.""" + cuda.init() + dev = cuda.Device(device_id) + ctx = dev.make_context() + print(f"[GPU {device_id}] {dev.name()}") + + try: + configure_macros(n_objs) + compiler.compile_kernel() + ptx_path = compiler.kernel_path("kernel.ptx") + mod = cuda.module_from_file(ptx_path) + step_fn = mod.get_function("OneStep") + init_fn = mod.get_function("InitTubeFlow") + + nbytes_ddf = TOTAL * NQ * 4 + ddf_gpu = cuda.mem_alloc(nbytes_ddf) + temp_gpu = cuda.mem_alloc(nbytes_ddf) + flag_gpu = cuda.mem_alloc(flag_host.nbytes) + indx_gpu = cuda.mem_alloc(indx_host.nbytes) + delta_gpu = cuda.mem_alloc(max(delta_host.nbytes, 4)) + + action_host = np.zeros(max(n_objs, 1), dtype=np.float32) + obs_host = np.zeros(max(n_objs * DIM, 1), dtype=np.float32) + action_gpu = cuda.mem_alloc(action_host.nbytes) + obs_gpu = cuda.mem_alloc(obs_host.nbytes) + + cuda.memcpy_htod(action_gpu, action_host) + cuda.memcpy_htod(obs_gpu, obs_host) + + block = (NT, 1, 1) + grid = (NX // NT, NY, 1) + + # Init Poiseuille + init_fn(flag_gpu, ddf_gpu, block=block, grid=grid) + ctx.synchronize() + + # Overwrite flag / indx / delta for cylinder case + cuda.memcpy_htod(flag_gpu, flag_host) + cuda.memcpy_htod(indx_gpu, indx_host) + cuda.memcpy_htod(delta_gpu, delta_host) + + # Step loop + t0 = time.time() + for i in range(n_steps): + step_fn(flag_gpu, ddf_gpu, temp_gpu, indx_gpu, delta_gpu, + action_gpu, obs_gpu, block=block, grid=grid) + ddf_gpu, temp_gpu = temp_gpu, ddf_gpu + ctx.synchronize() + dt = time.time() - t0 + mlups = TOTAL * n_steps / dt / 1e6 + print(f" {n_steps} steps in {dt:.2f}s ({mlups:.1f} MLUPS)") + + # Copy back + ddf = np.zeros(TOTAL * NQ, dtype=np.float32) + cuda.memcpy_dtoh(ddf, ddf_gpu) + flag_out = np.zeros(TOTAL, dtype=np.uint8) + cuda.memcpy_dtoh(flag_out, flag_gpu) + + return ddf, flag_out + finally: + ctx.pop() + + +# ━━━━━━━━━━━━━━━━━━━━━━━ Visualization ━━━━━━━━━━━━━━━━━━━━━━━ + +def plot_poiseuille(ddf, flag, out_path): + """3-panel figure: velocity mag, u(y) profile, pressure along centerline.""" + rho, u, v = extract_fields(ddf) + vel_mag = np.sqrt(u**2 + v**2) + + # Mask solid cells for display + mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool) + vel_mag_masked = np.ma.array(vel_mag, mask=mask) + + fig, axes = plt.subplots(1, 3, figsize=(18, 5)) + + # (a) Velocity magnitude heatmap + ax = axes[0] + im = ax.imshow(vel_mag_masked, origin='lower', aspect='auto', + cmap='jet', extent=[0, NX, 0, NY]) + plt.colorbar(im, ax=ax, label='|u|') + ax.set_title('Velocity magnitude') + ax.set_xlabel('x'); ax.set_ylabel('y') + + # (b) u(y) profile at x = NX/2 vs. analytical + ax = axes[1] + x_mid = NX // 2 + y_arr = np.arange(NY, dtype=float) + ax.plot(u[:, x_mid], y_arr, 'b-', lw=2, label='LBM') + ax.plot(analytical_poiseuille(y_arr), y_arr, 'r--', lw=1.5, label='Analytical') + ax.set_xlabel('u_x'); ax.set_ylabel('y') + ax.set_title(f'u(y) at x={x_mid}') + ax.legend() + ax.grid(True, alpha=0.3) + + # (c) Pressure along centerline y = NY/2 + ax = axes[2] + y_mid = NY // 2 + p = rho / 3.0 + ax.plot(np.arange(NX), p[y_mid, :], 'g-', lw=1.5) + ax.set_xlabel('x'); ax.set_ylabel('p = ρ/3') + ax.set_title(f'Pressure along centerline (y={y_mid})') + ax.grid(True, alpha=0.3) + + fig.suptitle(f'D2Q9 Poiseuille – NX={NX}, NY={NY}, VIS={VIS}, U0={U0}', fontsize=13) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + print(f" Saved: {out_path}") + plt.close(fig) + + +def plot_cylinder(ddf, flag, cx, cy, radius, out_path): + """3-panel figure: velocity magnitude (zoom), vorticity, streamlines.""" + rho, u, v = extract_fields(ddf) + vel_mag = np.sqrt(u**2 + v**2) + mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool) + + # Zoom window around cylinder + pad = int(radius * 8) + x0 = max(int(cx - pad), 0) + x1 = min(int(cx + pad * 2), NX) + y0 = max(int(cy - pad), 0) + y1 = min(int(cy + pad), NY) + + fig, axes = plt.subplots(1, 3, figsize=(20, 6)) + + # (a) Velocity magnitude (zoomed) + ax = axes[0] + vm_z = np.ma.array(vel_mag[y0:y1, x0:x1], mask=mask[y0:y1, x0:x1]) + im = ax.imshow(vm_z, origin='lower', aspect='equal', + cmap='jet', extent=[x0, x1, y0, y1]) + circ = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7) + ax.add_patch(circ) + plt.colorbar(im, ax=ax, label='|u|') + ax.set_title('Velocity magnitude') + + # (b) Vorticity + ax = axes[1] + dvdx = np.gradient(v, axis=1) + dudy = np.gradient(u, axis=0) + omega = dvdx - dudy + om_z = np.ma.array(omega[y0:y1, x0:x1], mask=mask[y0:y1, x0:x1]) + vmax = np.percentile(np.abs(omega[~mask]), 99) + im = ax.imshow(om_z, origin='lower', aspect='equal', + cmap='RdBu_r', extent=[x0, x1, y0, y1], + vmin=-vmax, vmax=vmax) + circ2 = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7) + ax.add_patch(circ2) + plt.colorbar(im, ax=ax, label='ω') + ax.set_title('Vorticity') + + # (c) Streamlines + ax = axes[2] + X, Y = np.meshgrid(np.arange(x0, x1), np.arange(y0, y1)) + u_z = u[y0:y1, x0:x1].copy() + v_z = v[y0:y1, x0:x1].copy() + u_z[mask[y0:y1, x0:x1]] = 0 + v_z[mask[y0:y1, x0:x1]] = 0 + speed = np.sqrt(u_z**2 + v_z**2) + ax.streamplot(X, Y, u_z, v_z, color=speed, cmap='jet', + density=2.0, linewidth=0.8) + circ3 = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7) + ax.add_patch(circ3) + ax.set_xlim(x0, x1); ax.set_ylim(y0, y1) + ax.set_aspect('equal') + ax.set_title('Streamlines') + + fig.suptitle(f'D2Q9 Cylinder Flow – Re_D={U0*1.5*2*radius/VIS:.0f}, D={2*radius}', fontsize=13) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + print(f" Saved: {out_path}") + plt.close(fig) + + +# ━━━━━━━━━━━━━━━━━━━━━━━ Main ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ + +def main(): + parser = argparse.ArgumentParser(description='D2Q9 Regression Test') + parser.add_argument('--device', type=int, default=2, + help='CUDA device ID (default: 2)') + parser.add_argument('--cylinder', action='store_true', + help='Also run cylinder flow test') + parser.add_argument('--steps-pois', type=int, default=5000, + help='Steps for Poiseuille (default: 5000)') + parser.add_argument('--steps-cyl', type=int, default=30000, + help='Steps for cylinder (default: 30000)') + args = parser.parse_args() + + out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output') + os.makedirs(out_dir, exist_ok=True) + + # ---- Test 1: Poiseuille ---- + print("\n===== Test 1: Poiseuille Channel Flow =====") + flag_pois = np.ones(TOTAL, dtype=np.uint8) + indx_pois = np.zeros(TOTAL, dtype=np.int32) + delta_pois = np.zeros(1, dtype=np.float32) + + ddf, flag = run_simulation(args.device, args.steps_pois, 0, + flag_pois, indx_pois, delta_pois) + plot_poiseuille(ddf, flag, os.path.join(out_dir, 'poiseuille_d2q9.png')) + + # Error metric + rho, u, v = extract_fields(ddf) + y_arr = np.arange(NY, dtype=float) + u_ana = analytical_poiseuille(y_arr) + x_mid = NX // 2 + u_num = u[:, x_mid] + # Interior cells only (skip walls) + err = np.max(np.abs(u_num[2:-2] - u_ana[2:-2])) / np.max(np.abs(u_ana[2:-2])) + print(f" L∞ relative error at x={x_mid}: {err:.2e}") + + # ---- Test 2: Cylinder ---- + if args.cylinder: + print("\n===== Test 2: Flow Around Cylinder =====") + cyl_cx, cyl_cy, cyl_r = 256.0, 256.0, 32.0 + flag_cyl, indx_cyl, delta_cyl = build_cylinder_data(cyl_cx, cyl_cy, cyl_r) + + ddf2, flag2 = run_simulation(args.device, args.steps_cyl, 1, + flag_cyl, indx_cyl, delta_cyl) + plot_cylinder(ddf2, flag2, cyl_cx, cyl_cy, cyl_r, + os.path.join(out_dir, 'cylinder_d2q9.png')) + + print("\nDone.") + + +if __name__ == '__main__': + main() diff --git a/tests/test_d3q19_cylinder.py b/tests/test_d3q19_cylinder.py new file mode 100644 index 0000000..bc69c9a --- /dev/null +++ b/tests/test_d3q19_cylinder.py @@ -0,0 +1,327 @@ +#!/usr/bin/env python3 +""" +D3Q19 SRT — Cylinder Wake Flow (Periodic Z) +============================================= +Tests the v2 modular kernel (kernel_v2.cu) in 3D. +Cylinder axis along z, parabolic inlet, pressure outlet, no-slip y-walls. +Produces cross-section visualizations at z=NZ/2. + +Usage: + python tests/test_d3q19_cylinder.py --device 3 + python tests/test_d3q19_cylinder.py --device 3 --re 200 --steps 50000 +""" + +import sys, os, argparse, time, struct + +sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src')) + +import numpy as np +import pycuda.driver as cuda + +import matplotlib +matplotlib.use('Agg') +import matplotlib.pyplot as plt + +from CelerisLab.cuda import compiler + +# ━━━━━━━━━━━━━━━━━━━━━━━ Configuration ━━━━━━━━━━━━━━━━━━━━━━━ +# Reasonable 3D grid — fits in < 500 MB GPU memory +NX, NY, NZ = 256, 128, 32 +NQ = 19 +NT = 128 +DIM = 3 +RHO = 1.0 +CYL_CX, CYL_CY = 64.0, 64.0 # Cylinder center (x,y) +CYL_R = 12.0 # Cylinder radius +TOTAL = NX * NY * NZ + +# D3Q19 paired direction ordering (from descriptors.cuh) +# 0:rest (1,2)±x (3,4)±y (5,6)±z +# (7,8)±(x+y) (9,10)±(x+z) (11,12)±(y+z) +# (13,14)±(x-y) (15,16)±(x-z) (17,18)±(y-z) +CX = np.array([0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0], dtype=np.int32) +CY = np.array([0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1], dtype=np.int32) +CZ = np.array([0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1], dtype=np.int32) +W = np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float32) + +FLUID_FLAG = 0x01 +SOLID_FLAG = 0x02 +OBSTACLE_FLAG = 0x04 # triggers half-way BB at adjacent fluid nodes + +# ━━━━━━━━━━━━━━━━━━━━━━━ Helpers ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ + +def compute_vis_omega(Re, D, U0): + """Compute viscosity and omega from target Reynolds number.""" + vis = U0 * D / Re + omega = 1.0 / (3.0 * vis + 0.5) + return vis, omega + + +def configure_macros_3d(vis, u0, n_objs=0): + """Write macros.h for D3Q19 SRT.""" + lines = compiler.read_lines(compiler.kernel_path("macros.h")) + defs = { + 'MULT_GPU': 'False', 'NT': NT, + 'X_1U': NX, 'Y_1U': NY, 'Z_1U': NZ, + 'LBtype': 'float', + 'UX': 1, 'UY': 1, 'UZ': 1, + 'NX': NX, 'NY': NY, 'NZ': NZ, + 'DIM': DIM, 'NQ': NQ, + 'VIS': f'{vis:.10f}', 'RHO': f'{RHO}', 'U0': u0, + 'N_OBJS': n_objs, + 'COLLISION_MODEL': 0, # SRT + 'STREAMING_MODEL': 0, # double-buffer + 'STORE_PRECISION': 0, # FP32 + 'USE_DDF_SHIFTING': 0, + } + for name, val in defs.items(): + lines = compiler.modify_macro(lines, name, val) + compiler.write_lines(compiler.kernel_path("macros.h"), lines) + + +def build_cylinder_3d(cx, cy, radius): + """Build flag array for 3D cylinder (axis along z, periodic z).""" + flag = np.ones(TOTAL, dtype=np.uint8) * FLUID_FLAG + for z in range(NZ): + for y in range(NY): + for x in range(NX): + k = z * NY * NX + y * NX + x + # Channel walls + if y == 0 or y == NY - 1 or x == 0 or x == NX - 1: + flag[k] = SOLID_FLAG + # Cylinder body (obstacle — triggers BB at fluid neighbors) + elif (x - cx)**2 + (y - cy)**2 < radius**2: + flag[k] = OBSTACLE_FLAG + return flag + + +def extract_fields_3d(ddf_host, z_slice): + """Extract rho, u, v, w at a given z-slice from D3Q19 DDF (v2 ordering).""" + # DDF layout: f[i * TOTAL + k] where k = z*NY*NX + y*NX + x + f = ddf_host.reshape(NQ, NZ, NY, NX) + fz = f[:, z_slice, :, :] # shape (NQ, NY, NX) + rho = np.sum(fz, axis=0) + ux = np.zeros_like(rho) + uy = np.zeros_like(rho) + uz_field = np.zeros_like(rho) + for i in range(NQ): + ux += CX[i] * fz[i] + uy += CY[i] * fz[i] + uz_field += CZ[i] * fz[i] + ux /= rho + uy /= rho + uz_field /= rho + return rho, ux, uy, uz_field + + +# ━━━━━━━━━━━━━━━━━━━━━━━ Simulation ━━━━━━━━━━━━━━━━━━━━━━━━━━ + +def run_d3q19(device_id, n_steps, vis, u0, flag_host): + """Compile v2 kernel, run D3Q19 SRT, return DDF.""" + omega = 1.0 / (3.0 * vis + 0.5) + cuda.init() + dev = cuda.Device(device_id) + ctx = dev.make_context() + print(f"[GPU {device_id}] {dev.name()}") + + try: + configure_macros_3d(vis, u0) + compiler.compile_kernel_v2() + ptx_path = compiler.kernel_path("kernel_v2.ptx") + mod = cuda.module_from_file(ptx_path) + + # Get kernels (extern "C" entries from kernel_v2.cu) + init_fn = mod.get_function("InitTubeFlow_v2") + step_fn = mod.get_function("OneStep") + + # Set d_params.omega via __constant__ memory + params_ptr, params_size = mod.get_global("d_params") + # LBMParams struct layout (see params.cuh): + # Nx(4) Ny(4) Nz(4) N(8) omega(4) omega_bulk(4) fx(4) fy(4) fz(4) + # rho_ref(4) u_inlet(4) n_objects(4) + # Pack: unsigned int Nx, Ny, Nz; unsigned long N; float omega, omega_bulk, fx, fy, fz, rho_ref, u_inlet; unsigned int n_objects + params_data = struct.pack('IIIQfffffffI', + NX, NY, NZ, + TOTAL, + omega, 0.0, # omega, omega_bulk + 0.0, 0.0, 0.0, # fx, fy, fz + RHO, u0, # rho_ref, u_inlet + 0) # n_objects + # Pad to match struct size + if len(params_data) < params_size: + params_data += b'\x00' * (params_size - len(params_data)) + cuda.memcpy_htod(params_ptr, params_data) + + # Allocate + nbytes_ddf = TOTAL * NQ * 4 + ddf_gpu = cuda.mem_alloc(nbytes_ddf) + temp_gpu = cuda.mem_alloc(nbytes_ddf) + flag_gpu = cuda.mem_alloc(flag_host.nbytes) + indx_gpu = cuda.mem_alloc(TOTAL * 4) + delta_gpu = cuda.mem_alloc(4) + action_gpu = cuda.mem_alloc(4) + obs_gpu = cuda.mem_alloc(4) + + # Dummy arrays + cuda.memset_d32(indx_gpu, 0, TOTAL) + cuda.memset_d32(delta_gpu, 0, 1) + cuda.memset_d32(action_gpu, 0, 1) + cuda.memset_d32(obs_gpu, 0, 1) + + block = (NT, 1, 1) + grid = (NX // NT, NY, NZ) + + # Initialize parabolic flow + init_fn(flag_gpu, ddf_gpu, block=block, grid=grid) + ctx.synchronize() + + # Overwrite flags with cylinder geometry + cuda.memcpy_htod(flag_gpu, flag_host) + + # Step loop + print(f" Running {n_steps} steps (NX={NX}, NY={NY}, NZ={NZ}, omega={omega:.4f})...") + t0 = time.time() + for i in range(n_steps): + step_fn(flag_gpu, ddf_gpu, temp_gpu, indx_gpu, delta_gpu, + action_gpu, obs_gpu, + block=block, grid=grid) + ddf_gpu, temp_gpu = temp_gpu, ddf_gpu + if (i + 1) % 5000 == 0: + ctx.synchronize() + elapsed = time.time() - t0 + mlups = TOTAL * (i + 1) / elapsed / 1e6 + print(f" step {i+1}/{n_steps} ({mlups:.1f} MLUPS)") + ctx.synchronize() + dt = time.time() - t0 + mlups = TOTAL * n_steps / dt / 1e6 + print(f" Done: {dt:.1f}s, {mlups:.1f} MLUPS") + + # Copy back + ddf = np.zeros(TOTAL * NQ, dtype=np.float32) + cuda.memcpy_dtoh(ddf, ddf_gpu) + return ddf + + finally: + ctx.pop() + + +# ━━━━━━━━━━━━━━━━━━━━━━━ Visualization ━━━━━━━━━━━━━━━━━━━━━━━ + +def plot_d3q19_cylinder(ddf, flag, Re, u0, out_path): + """4-panel figure at z=NZ/2: vel-mag, vorticity, streamlines, u(y) profile.""" + z_mid = NZ // 2 + rho, ux, uy, uz = extract_fields_3d(ddf, z_mid) + vel_mag = np.sqrt(ux**2 + uy**2 + uz**2) + + mask2d = ((flag.reshape(NZ, NY, NX)[z_mid] & (SOLID_FLAG | OBSTACLE_FLAG)) != 0) + vel_masked = np.ma.array(vel_mag, mask=mask2d) + + fig, axes = plt.subplots(2, 2, figsize=(16, 10)) + + # (a) Velocity magnitude + ax = axes[0, 0] + im = ax.imshow(vel_masked, origin='lower', aspect='auto', + cmap='jet', extent=[0, NX, 0, NY]) + circ = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7) + ax.add_patch(circ) + plt.colorbar(im, ax=ax, label='|u|') + ax.set_title(f'Velocity magnitude (z={z_mid})') + ax.set_xlabel('x'); ax.set_ylabel('y') + + # (b) Vorticity ω_z = ∂v/∂x − ∂u/∂y + ax = axes[0, 1] + dvdx = np.gradient(uy, axis=1) + dudy = np.gradient(ux, axis=0) + omega_z = dvdx - dudy + om_masked = np.ma.array(omega_z, mask=mask2d) + vmax = np.percentile(np.abs(omega_z[~mask2d]), 99) if np.any(~mask2d) else 1e-3 + im = ax.imshow(om_masked, origin='lower', aspect='auto', + cmap='RdBu_r', extent=[0, NX, 0, NY], + vmin=-vmax, vmax=vmax) + circ2 = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7) + ax.add_patch(circ2) + plt.colorbar(im, ax=ax, label='ω_z') + ax.set_title('Vorticity ω_z') + ax.set_xlabel('x'); ax.set_ylabel('y') + + # (c) Streamlines + ax = axes[1, 0] + X, Y = np.meshgrid(np.arange(NX), np.arange(NY)) + ux_c = ux.copy(); ux_c[mask2d] = 0 + uy_c = uy.copy(); uy_c[mask2d] = 0 + speed = np.sqrt(ux_c**2 + uy_c**2) + ax.streamplot(X, Y, ux_c, uy_c, color=speed, cmap='jet', + density=2.5, linewidth=0.7) + circ3 = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7) + ax.add_patch(circ3) + ax.set_xlim(0, NX); ax.set_ylim(0, NY) + ax.set_aspect('auto') + ax.set_title('Streamlines') + ax.set_xlabel('x'); ax.set_ylabel('y') + + # (d) u_x(y) profiles at different x stations + ax = axes[1, 1] + y_arr = np.arange(NY) + # Analytical parabolic inlet + yy = (y_arr - 0.5 * (NY - 1)) / (NY - 2.0) + u_ana = u0 * 1.5 * (1 - 4 * yy**2) + + x_stations = [NX // 8, NX // 4, NX // 2, 3 * NX // 4] + for xs in x_stations: + ax.plot(ux[:, xs], y_arr, label=f'x={xs}') + ax.plot(u_ana, y_arr, 'k--', lw=1.5, label='Analytical inlet') + ax.set_xlabel('u_x'); ax.set_ylabel('y') + ax.set_title('u_x(y) profiles') + ax.legend(fontsize=8) + ax.grid(True, alpha=0.3) + + fig.suptitle(f'D3Q19 SRT Cylinder — Re={Re:.0f}, D={2*CYL_R:.0f}, ' + f'Grid={NX}×{NY}×{NZ}', fontsize=13) + fig.tight_layout() + fig.savefig(out_path, dpi=150) + print(f" Saved: {out_path}") + plt.close(fig) + + +# ━━━━━━━━━━━━━━━━━━━━━━━ Main ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ + +def main(): + parser = argparse.ArgumentParser(description='D3Q19 SRT Cylinder Flow') + parser.add_argument('--device', type=int, default=0, + help='CUDA device ID (default: 0)') + parser.add_argument('--re', type=float, default=100.0, + help='Reynolds number based on diameter (default: 100)') + parser.add_argument('--u0', type=float, default=0.04, + help='Inlet characteristic velocity (default: 0.04)') + parser.add_argument('--steps', type=int, default=30000, + help='Number of LBM steps (default: 30000)') + args = parser.parse_args() + + out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output') + os.makedirs(out_dir, exist_ok=True) + + D = 2 * CYL_R + vis, omega = compute_vis_omega(args.re, D, args.u0) + print(f"\n===== D3Q19 SRT Cylinder Flow =====") + print(f" Re = {args.re:.0f}, D = {D:.0f}, U0 = {args.u0}") + print(f" ν = {vis:.6f}, ω = {omega:.4f}") + if omega > 1.95: + print(f" WARNING: omega={omega:.4f} is close to 2.0, stability may be poor.") + print(f" Consider reducing U0 or Re.") + + flag = build_cylinder_3d(CYL_CX, CYL_CY, CYL_R) + n_solid = np.sum(flag == SOLID_FLAG) + n_fluid = np.sum(flag == FLUID_FLAG) + print(f" Grid: {NX}×{NY}×{NZ} = {TOTAL} cells (fluid: {n_fluid}, solid: {n_solid})") + print(f" Memory: ~{2 * TOTAL * NQ * 4 / 1e6:.0f} MB for double-buffer DDF") + + ddf = run_d3q19(args.device, args.steps, vis, args.u0, flag) + + plot_d3q19_cylinder(ddf, flag, args.re, args.u0, + os.path.join(out_dir, f'cylinder_d3q19_re{int(args.re)}.png')) + + print("\nDone.") + + +if __name__ == '__main__': + main()