Frank14f/CelerisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

improve LES stability

Browse Source
This commit is contained in:
Frank14f 2026-03-29 22:16:20 +08:00
parent f3e2e557d4
commit 5076e5d789
13 changed files with 1163 additions and 141 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -69,3 +69,4 @@ venv.bak/
# reference:
ref/
output/

BIN
output/poiseuille_d2q9.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 191 KiB

After

Width:  |  Height:  |  Size: 122 KiB

8
src/CelerisLab/cuda/compiler.py
View File

@ -95,7 +95,9 @@ 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):
use_ddf_shifting: int = 0,
use_les: int = 0,
les_cs: float = 0.16):
"""Configure macros.h for the new modular kernel architecture.
Args:
@ -103,6 +105,8 @@ def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig,
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
use_les: 0=off (default), 1=Smagorinsky LES
les_cs: Smagorinsky constant C_s
"""
# First apply legacy config
config_kernal(config_cuda, config_field)
@ -113,6 +117,8 @@ def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig,
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)
lines = modify_macro(lines, "USE_LES", use_les)
lines = modify_macro(lines, "LES_CS", f"{les_cs:.6f}f")
write_lines(kernel_path("macros.h"), lines)
def config_object(n_obj: int):

85
src/CelerisLab/lbm/driver.py
View File

@ -23,6 +23,13 @@ class FlowField:
field_config: utils.FlowFieldConfig,
cuda_config: utils.CudaConfig,
device_id: Union[int, List[int]] = None,
use_kernel_v2: bool = True,
collision_model: int = 0,
streaming_model: int = 0,
store_precision: int = 0,
use_ddf_shifting: int = 0,
use_les: int = 0,
les_cs: float = 0.16,
):
self.field_config = field_config
self.cuda_config = cuda_config
@ -49,10 +56,26 @@ class FlowField:
utils.check_cuda_capability(field_config, cuda_config, device_id)
# Config kernel
compiler.config_kernal(cuda_config, field_config)
compiler.config_object(int(0))
# compiler.config_sensor(int(0))
self.use_kernel_v2 = bool(use_kernel_v2)
self.collision_model = int(collision_model)
self.streaming_model = int(streaming_model)
self.store_precision = int(store_precision)
self.use_ddf_shifting = int(use_ddf_shifting)
self.use_les = int(use_les)
self.les_cs = float(les_cs)
if self.collision_model not in (0, 1, 2):
raise ValueError("collision_model must be 0(SRT), 1(TRT), or 2(MRT).")
if self.streaming_model not in (0, 1):
raise ValueError("streaming_model must be 0(double-buffer) or 1(esopull).")
if self.store_precision not in (0, 1, 2):
raise ValueError("store_precision must be 0(FP32), 1(FP16S), or 2(FP16C).")
if self.use_ddf_shifting not in (0, 1):
raise ValueError("use_ddf_shifting must be 0 or 1.")
if self.use_les not in (0, 1):
raise ValueError("use_les must be 0 or 1.")
if not (0.0 < self.les_cs < 1.0):
raise ValueError("les_cs must be in (0, 1).")
# Set constants
if field_config.data_type == "FP32":
@ -84,11 +107,10 @@ class FlowField:
f"Unsupported lattice type {field_config.lattice} in {field_config.dimensionality} dimensions."
)
# Compile kernel
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
initflow = self.ptx.get_function("InitTubeFlow")
self.objects = {}
# Compile and load kernel
self._rebuild_kernel()
# Initialize memory
self.ddf = np.zeros(self.FIELD_SIZE * self.LATTICE, dtype=self.DATA_TYPE)
@ -105,11 +127,10 @@ class FlowField:
self.delta_gpu = cuda.mem_alloc(1)
self.vortex_gpu = cuda.mem_alloc(self.vortex_config.nbytes)
self.objects = {}
self.action = np.zeros(0, dtype=self.DATA_TYPE)
self.obs = np.zeros(0, dtype=self.DATA_TYPE)
initflow(
self.initflow(
self.flag_gpu,
self.ddf_gpu,
block=(self.cuda_config.threads_per_block, 1, 1),
@ -122,6 +143,38 @@ class FlowField:
cuda.memcpy_dtoh(self.flag, self.flag_gpu)
cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
def _configure_kernel(self):
if self.use_kernel_v2:
compiler.config_kernal_v2(
self.cuda_config,
self.field_config,
collision_model=self.collision_model,
streaming_model=self.streaming_model,
store_precision=self.store_precision,
use_ddf_shifting=self.use_ddf_shifting,
use_les=self.use_les,
les_cs=self.les_cs,
)
else:
compiler.config_kernal(self.cuda_config, self.field_config)
def _compile_and_load_kernel(self):
if self.use_kernel_v2:
compiler.compile_kernel_v2()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
self.step = self.ptx.get_function("OneStep")
self.initflow = self.ptx.get_function("InitTubeFlow_v2")
else:
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
self.initflow = self.ptx.get_function("InitTubeFlow")
def _rebuild_kernel(self):
self._configure_kernel()
compiler.config_object(len(self.objects))
self._compile_and_load_kernel()
def add_cylinder(self, center: Tuple[float, float, float], radius: float, id_obj: Optional[int] = None):
x_c, y_c, z_c = center
@ -193,10 +246,7 @@ class FlowField:
cuda.memcpy_htod(self.flag_gpu, self.flag)
cuda.memcpy_htod(self.indx_gpu, self.indx)
compiler.config_object(len(self.objects))
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
self._rebuild_kernel()
def add_sensor(self, center: Tuple[float, float, float], radius: float):
x_c, y_c, z_c = center
@ -235,10 +285,7 @@ class FlowField:
cuda.memcpy_htod(self.flag_gpu, self.flag)
cuda.memcpy_htod(self.indx_gpu, self.indx)
compiler.config_object(len(self.objects))
compiler.compile_kernel()
self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
self.step = self.ptx.get_function("OneStep")
self._rebuild_kernel()
def add_vortex(self, center: Tuple[float, float, float], radius: float, strength: float, direction: float, type: str):
x_c, y_c, z_c = center

151
src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
View File

@ -15,6 +15,31 @@
#ifndef CELERIS_BOUNDARY_INLET_OUTLET_CUH
#define CELERIS_BOUNDARY_INLET_OUTLET_CUH
#ifndef INLET_PROFILE
#define INLET_PROFILE 1
#endif
#ifndef OUTLET_MODE
#define OUTLET_MODE 0
#endif
#ifndef OUTLET_BACKFLOW_CLAMP
#define OUTLET_BACKFLOW_CLAMP 1
#endif
#ifndef OUTLET_BLEND_ALPHA
#define OUTLET_BLEND_ALPHA 0.70f
#endif
__device__ __forceinline__ float inlet_target_u(float y_coord) {
#if INLET_PROFILE == 0
return U0;
#else
float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
return U0 * 1.5f * (1.0f - 4.0f * yy * yy);
#endif
}
#if NQ == 9
// ---------------------------------------------------------------------------
@ -33,37 +58,21 @@ __device__ inline void apply_parabolic_inlet(float* __restrict__ f,
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;
float rho_neb, u_neb, v_neb;
compute_rho_u(f_neb, rho_neb, u_neb, v_neb);
// 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 u_target = inlet_target_u(y_coord);
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;
float feq_tar[9], feq_neb[9];
compute_feq(rho_neb, u_target, v_target, feq_tar);
compute_feq(rho_neb, u_neb, v_neb, feq_neb);
// 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;
f[1] = f_neb[1] - feq_neb[1] + feq_tar[1];
f[5] = f_neb[5] - feq_neb[5] + feq_tar[5];
f[7] = f_neb[7] - feq_neb[7] + feq_tar[7];
}
// ---------------------------------------------------------------------------
@ -76,35 +85,39 @@ __device__ inline void apply_pressure_outlet(float* __restrict__ f,
const float* __restrict__ f_neb,
float y_coord)
{
float p_out = 0.0f;
(void)y_coord;
// 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;
#if OUTLET_MODE == 1
// Simple zero-gradient copy for unknown incoming directions at outlet.
f[2] = f_neb[2];
f[8] = f_neb[8];
f[6] = f_neb[6];
#else
// 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;
// Convective non-reflecting style: extrapolate outlet density from neighbor
// and keep neighbor velocity.
float rho_neb, u_neb, v_neb;
compute_rho_u(f_neb, rho_neb, u_neb, v_neb);
#if OUTLET_BACKFLOW_CLAMP
u_neb = fmaxf(u_neb, 0.0f);
#endif
float rho_out = rho_neb;
// 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;
float feq_tar[9], feq_neb[9];
compute_feq(rho_out, u_neb, v_neb, feq_tar);
compute_feq(rho_neb, u_neb, v_neb, feq_neb);
// 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;
#if OUTLET_MODE == 2
const float a = OUTLET_BLEND_ALPHA;
f[2] = a * (f_neb[2] - feq_neb[2] + feq_tar[2]) + (1.0f - a) * f_neb[2];
f[8] = a * (f_neb[8] - feq_neb[8] + feq_tar[8]) + (1.0f - a) * f_neb[8];
f[6] = a * (f_neb[6] - feq_neb[6] + feq_tar[6]) + (1.0f - a) * f_neb[6];
#else
f[2] = f_neb[2] - feq_neb[2] + feq_tar[2];
f[8] = f_neb[8] - feq_neb[8] + feq_tar[8];
f[6] = f_neb[6] - feq_neb[6] + feq_tar[6];
#endif
#endif
}
#endif // NQ == 9
@ -128,8 +141,7 @@ __device__ inline void apply_parabolic_inlet_3d(float* __restrict__ f,
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);
float u_tar = inlet_target_u(y_coord);
// feq arrays
float feq_tar[19], feq_neb[19];
@ -148,24 +160,47 @@ __device__ inline void apply_pressure_outlet_3d(float* __restrict__ f,
const float* __restrict__ f_neb,
float y_coord)
{
(void)y_coord;
#if OUTLET_MODE == 1
// Simple zero-gradient copy for unknown incoming directions at outlet.
f[2] = f_neb[2];
f[8] = f_neb[8];
f[10] = f_neb[10];
f[14] = f_neb[14];
f[16] = f_neb[16];
#else
// Neighbor macros
float rho_neb, un, vn, wn;
compute_rho_u(f_neb, rho_neb, un, vn, wn);
#if OUTLET_BACKFLOW_CLAMP
un = fmaxf(un, 0.0f);
#endif
// 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);
// Convective non-reflecting style: extrapolate outlet density from neighbor
// and keep neighbor velocity.
float rho_out = rho_neb;
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);
compute_feq(rho_out, un, vn, wn, feq_tar);
compute_feq(rho_neb, un, vn, wn, feq_neb);
// Reconstruct cx<0 directions: i = 2, 8, 10, 14, 16
#if OUTLET_MODE == 2
const float a = OUTLET_BLEND_ALPHA;
f[2] = a * (f_neb[2] - feq_neb[2] + feq_tar[2]) + (1.0f - a) * f_neb[2];
f[8] = a * (f_neb[8] - feq_neb[8] + feq_tar[8]) + (1.0f - a) * f_neb[8];
f[10] = a * (f_neb[10] - feq_neb[10] + feq_tar[10]) + (1.0f - a) * f_neb[10];
f[14] = a * (f_neb[14] - feq_neb[14] + feq_tar[14]) + (1.0f - a) * f_neb[14];
f[16] = a * (f_neb[16] - feq_neb[16] + feq_tar[16]) + (1.0f - a) * f_neb[16];
#else
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
#endif
}
#endif // NQ == 19

33
src/CelerisLab/lbm/kernels/kernel_v2.cu
View File

@ -39,6 +39,7 @@
#include "operators/collision_trt.cuh"
#include "operators/collision_mrt.cuh"
#include "operators/forcing_guo.cuh"
#include "operators/turbulence_smag.cuh"
// ---------------------------------------------------------------------------
// Layer 3: Streaming
@ -208,18 +209,31 @@ __global__ void OneStep(
}
}
// Obstacle nodes: full-way bounce-back (same swap as wall nodes)
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
if (fl & LEGACY_FLUID) {
float feq[NQ], Fin[NQ];
compute_feq(rho_n, ux, uy, feq);
zero_forcing(Fin);
float omega_col = d_params.omega;
#if USE_LES
omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
#endif
omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
#if COLLISION_MODEL == 0
collide_srt(f, feq, Fin, d_params.omega);
collide_srt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 1
collide_trt(f, feq, Fin, d_params.omega);
collide_trt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 2
collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
#endif
}
@ -289,7 +303,18 @@ __global__ void OneStep(
float feq[NQ], Fin[NQ];
compute_feq(rho_n, ux, uy, uz, feq);
zero_forcing(Fin);
collide_srt(f, feq, Fin, d_params.omega);
float omega_col = d_params.omega;
#if USE_LES
omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
#endif
omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
#if COLLISION_MODEL == 0
collide_srt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 1
collide_trt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 2
collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
#endif
}
stream_pull_store(k, f, fi_out);

40
src/CelerisLab/lbm/kernels/macros.h
View File

@ -61,3 +61,43 @@
#ifndef USE_DDF_SHIFTING
#define USE_DDF_SHIFTING 0
#endif
// LES model: 0=off, 1=Smagorinsky
#ifndef USE_LES
#define USE_LES 0
#endif
// Smagorinsky constant C_s
#ifndef LES_CS
#define LES_CS 0.16f
#endif
// Inlet profile: 1=parabolic (channel), 0=uniform (external flow)
#ifndef INLET_PROFILE
#define INLET_PROFILE 1
#endif
// Outlet mode: 0=non-equilibrium extrapolation, 1=zero-gradient copy (more dissipative)
#ifndef OUTLET_MODE
#define OUTLET_MODE 0
#endif
// Outlet blend factor for damped outlet mode (OUTLET_MODE=2):
// f_out = a*(non-eq extrapolation) + (1-a)*(zero-gradient copy)
#ifndef OUTLET_BLEND_ALPHA
#define OUTLET_BLEND_ALPHA 0.70f
#endif
// Outlet backflow clamp: 0=off, 1=force non-negative streamwise velocity at outlet target
#ifndef OUTLET_BACKFLOW_CLAMP
#define OUTLET_BACKFLOW_CLAMP 1
#endif
// Global collision omega guardrails
#ifndef OMEGA_COLLISION_MIN
#define OMEGA_COLLISION_MIN 0.01f
#endif
#ifndef OMEGA_COLLISION_MAX
#define OMEGA_COLLISION_MAX 1.999f
#endif

123
src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh
View File

@ -31,19 +31,19 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
float omega)
{
// ----- Relaxation rates -----
// s_rho = s_jx = s_jy = 1.0 (conserved moments → overwrite to eq)
// s_rho = s_jx = s_jy = 0.0 (strictly conserved in collision)
// 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_rho = 0.0f; // conserved moments
const float s_e = 1.2f;
const float s_eps = 1.2f;
const float s_jx = 1.0f;
const float s_jx = 0.0f;
const float s_q = 1.2f;
const float s_jy = 1.0f;
const float s_jy = 0.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;
// ----- Pressure from local density -----
const float p = rho * (1.0f / 3.0f);
// ----- Forward transform: m = M * g (new paired ordering) -----
float m[9];
@ -61,14 +61,14 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
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
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];
@ -112,21 +112,106 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
#elif NQ == 19
// ---------------------------------------------------------------------------
// D3Q19 MRT (placeholder – use SRT or TRT until fully validated)
// D3Q19 MRT (tensor-projected multi-mode relaxation)
//
// This implementation performs a true multi-relaxation update by splitting
// non-equilibrium content into three subspaces:
// 1) deviatoric 2nd-order stress modes -> s_nu (shear, controlled by omega)
// 2) isotropic 2nd-order stress mode -> s_bulk (bulk, controlled by omega_bulk)
// 3) higher-order residual modes -> s_high (kinetic damping)
//
// It is formulated in Hermite/tensor space, avoiding explicit 19x19 matrices
// while still relaxing different mode families at different rates.
// ---------------------------------------------------------------------------
__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;
// Relaxation rates
const float s_nu = omega;
const float s_bulk = (d_params.omega_bulk > 0.0f) ? d_params.omega_bulk : 1.2f;
const float s_high = 1.4f;
const float one_minus_s_nu = 1.0f - s_nu;
const float one_minus_s_bulk = 1.0f - s_bulk;
const float one_minus_s_high = 1.0f - s_high;
// Non-equilibrium populations
float neq[19];
#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]));
neq[i] = g[i] - feq[i];
}
// Build non-equilibrium 2nd-order tensor Πneq = Σ (ci ci) (fi-feqi)
float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
float pixy = 0.0f, pixz = 0.0f, piyz = 0.0f;
#pragma unroll
for (int i = 0; i < 19; i++) {
const float ci_x = (float)d_cx[i];
const float ci_y = (float)d_cy[i];
const float ci_z = (float)d_cz[i];
const float fneq = neq[i];
pixx += fneq * ci_x * ci_x;
piyy += fneq * ci_y * ci_y;
pizz += fneq * ci_z * ci_z;
pixy += fneq * ci_x * ci_y;
pixz += fneq * ci_x * ci_z;
piyz += fneq * ci_y * ci_z;
}
// Isotropic / deviatoric split of Πneq
const float tr_pi = pixx + piyy + pizz;
const float tr_pi_thrd = tr_pi * (1.0f / 3.0f);
const float dev_xx = pixx - tr_pi_thrd;
const float dev_yy = piyy - tr_pi_thrd;
const float dev_zz = pizz - tr_pi_thrd;
const float dev_xy = pixy;
const float dev_xz = pixz;
const float dev_yz = piyz;
// cs^2 = 1/3 => 2*cs^4 = 2/9, inverse = 4.5
// Project neq onto second-order Hermite basis and relax mode families.
const float proj_pref = 4.5f;
const float c_tau = 1.0f - 0.5f * omega;
#pragma unroll
for (int i = 0; i < 19; i++) {
const float ci_x = (float)d_cx[i];
const float ci_y = (float)d_cy[i];
const float ci_z = (float)d_cz[i];
// Q = ci ci - cs^2 I
const float q_xx = ci_x * ci_x - CS2;
const float q_yy = ci_y * ci_y - CS2;
const float q_zz = ci_z * ci_z - CS2;
const float q_xy = ci_x * ci_y;
const float q_xz = ci_x * ci_z;
const float q_yz = ci_y * ci_z;
// Contractions Q:Πdev and Q:Πiso
const float q_dot_dev =
q_xx * dev_xx + q_yy * dev_yy + q_zz * dev_zz
+ 2.0f * (q_xy * dev_xy + q_xz * dev_xz + q_yz * dev_yz);
const float q_trace = q_xx + q_yy + q_zz; // = |ci|^2 - 1
const float q_dot_iso = q_trace * tr_pi_thrd;
const float neq_dev = d_w[i] * proj_pref * q_dot_dev;
const float neq_iso = d_w[i] * proj_pref * q_dot_iso;
const float neq_h = neq[i] - neq_dev - neq_iso;
g[i] = feq[i]
+ one_minus_s_nu * neq_dev
+ one_minus_s_bulk * neq_iso
+ one_minus_s_high * neq_h
+ c_tau * Fin[i];
}
}

111
src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh Normal file
View File

@ -0,0 +1,111 @@
// CelerisLab – operators/turbulence_smag.cuh
// Smagorinsky LES: compute effective omega from local non-equilibrium stress.
//
// Model:
// nu_t = (C_s * Delta)^2 * |S|, Delta = 1 (lattice unit)
// omega_eff = 1 / (3*(nu + nu_t) + 0.5)
//
// |S| is estimated from non-equilibrium second-order moments using tau0.
// ============================================================================
#ifndef CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
#define CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
#ifndef USE_LES
#define USE_LES 0
#endif
#ifndef LES_CS
#define LES_CS 0.16f
#endif
__device__ __forceinline__ float clamp_omega(float w) {
return fminf(1.999f, fmaxf(0.01f, w));
}
#if NQ == 9
__device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
const float* __restrict__ feq,
float rho,
float omega0)
{
const float rho_safe = fmaxf(rho, 1.0e-12f);
const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
// Πneq = Σ (ci ci)(fi-feqi)
float pixx = 0.0f, piyy = 0.0f, pixy = 0.0f;
#pragma unroll
for (int i = 0; i < 9; i++) {
const float fneq = f[i] - feq[i];
const float cx = (float)d_cx[i];
const float cy = (float)d_cy[i];
pixx += fneq * cx * cx;
piyy += fneq * cy * cy;
pixy += fneq * cx * cy;
}
const float denom = 2.0f * rho_safe * CS2 * tau0;
const float sxx = -pixx / denom;
const float syy = -piyy / denom;
const float sxy = -pixy / denom;
// |S| = sqrt(2 Sij Sij)
const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + 2.0f*sxy*sxy));
const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
const float nut = (LES_CS * LES_CS) * s_mag;
const float nue = nu0 + nut;
return clamp_omega(1.0f / (3.0f * nue + 0.5f));
}
#elif NQ == 19
__device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
const float* __restrict__ feq,
float rho,
float omega0)
{
const float rho_safe = fmaxf(rho, 1.0e-12f);
const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
// Πneq = Σ (ci ci)(fi-feqi)
float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
float pixy = 0.0f, pixz = 0.0f, piyz = 0.0f;
#pragma unroll
for (int i = 0; i < 19; i++) {
const float fneq = f[i] - feq[i];
const float cx = (float)d_cx[i];
const float cy = (float)d_cy[i];
const float cz = (float)d_cz[i];
pixx += fneq * cx * cx;
piyy += fneq * cy * cy;
pizz += fneq * cz * cz;
pixy += fneq * cx * cy;
pixz += fneq * cx * cz;
piyz += fneq * cy * cz;
}
const float denom = 2.0f * rho_safe * CS2 * tau0;
const float sxx = -pixx / denom;
const float syy = -piyy / denom;
const float szz = -pizz / denom;
const float sxy = -pixy / denom;
const float sxz = -pixz / denom;
const float syz = -piyz / denom;
// |S| = sqrt(2 Sij Sij)
const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + szz*szz
+ 2.0f*(sxy*sxy + sxz*sxz + syz*syz)));
const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
const float nut = (LES_CS * LES_CS) * s_mag;
const float nue = nu0 + nut;
return clamp_omega(1.0f / (3.0f * nue + 0.5f));
}
#endif // NQ
#endif // CELERIS_OPERATORS_TURBULENCE_SMAG_CUH

51
src/CelerisLab/lbm/kernels/step/one_step_double.cu
View File

@ -69,7 +69,8 @@ __global__ void StreamCollideDouble(
// ----- Boundary conditions (inlet/outlet/wall) -----
#if NQ == 9
if (fl & LEGACY_SOLID) {
if (x == 0) {
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++) {
@ -77,18 +78,31 @@ __global__ void StreamCollideDouble(
}
apply_parabolic_inlet(f, f_neb, (float)y);
}
else if (x == NX - 1) {
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 {
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
}
}
}
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;
}
}
#elif NQ == 19
if (fl & LEGACY_SOLID) {
if (x == 0) {
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, z);
for (int i = 0; i < NQ; i++) {
@ -96,13 +110,18 @@ __global__ void StreamCollideDouble(
}
apply_parabolic_inlet_3d(f, f_neb, (float)y);
}
else if (x == NX - 1) {
else if (x == (unsigned int)(NX - 1) && interior_y) {
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 {
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
}
}
}
#endif
@ -126,20 +145,36 @@ __global__ void StreamCollideDouble(
#if NQ == 9
compute_feq(rho_n, ux, uy, feq);
zero_forcing(Fin);
float omega_col = d_params.omega;
#if USE_LES
omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
#endif
omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
// ----- Collision -----
#if COLLISION_MODEL == 0
collide_srt(f, feq, Fin, d_params.omega);
collide_srt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 1
collide_trt(f, feq, Fin, d_params.omega);
collide_trt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 2
collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
#endif
#elif NQ == 19
compute_feq(rho_n, ux, uy, uz, feq);
zero_forcing(Fin);
collide_srt(f, feq, Fin, d_params.omega);
float omega_col = d_params.omega;
#if USE_LES
omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
#endif
omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
#if COLLISION_MODEL == 0
collide_srt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 1
collide_trt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 2
collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
#endif
#endif
}

77
src/CelerisLab/lbm/kernels/step/one_step_esopull.cu
View File

@ -70,7 +70,8 @@ __global__ void StreamCollideEsoPull(
// ----- Boundary conditions -----
#if NQ == 9
if (fl & LEGACY_SOLID) {
if (x == 0) {
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);
load_f_esopull(k_neb, f_neb, fi, j, t); // approximate: reuse j
@ -80,19 +81,67 @@ __global__ void StreamCollideEsoPull(
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)) {
else if (x == (unsigned int)(NX - 1) && interior_y) {
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);
} else {
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
}
}
}
if (fl & LEGACY_OBSTACLE) {
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
}
}
if (y == 1 || y == (unsigned int)(NY - 2)) {
apply_wall_bb_d2q9(y, f);
}
#elif NQ == 19
if (fl & LEGACY_SOLID) {
bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
if (x == 0 && interior_y) {
unsigned long k_neb = linear_index(x + 1u, y, z);
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_parabolic_inlet_3d(f, f_neb, (float)y);
}
else if (x == (unsigned int)(NX - 1) && interior_y) {
unsigned long k_neb = linear_index(x - 1u, y, z);
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_3d(f, f_neb, (float)y);
} else {
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
}
}
}
if (fl & LEGACY_OBSTACLE) {
#pragma unroll
for (int i = 1; i < NQ; i += 2) {
float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
}
}
if (y == 1 || y == (unsigned int)(NY - 2)) {
apply_wall_bb_d3q19_y(y, f);
}
#endif
// ----- Forcing -----
@ -123,13 +172,18 @@ __global__ void StreamCollideEsoPull(
float feq[NQ];
#if NQ == 9
compute_feq(rho_n, ux, uy, feq);
float omega_col = d_params.omega;
#if USE_LES
omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
#endif
omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
#if COLLISION_MODEL == 0
collide_srt(f, feq, Fin, d_params.omega);
collide_srt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 1
collide_trt(f, feq, Fin, d_params.omega);
collide_trt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 2
collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
#endif
#elif NQ == 19
@ -137,7 +191,18 @@ __global__ void StreamCollideEsoPull(
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);
float omega_col = d_params.omega;
#if USE_LES
omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
#endif
omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
#if COLLISION_MODEL == 0
collide_srt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 1
collide_trt(f, feq, Fin, omega_col);
#elif COLLISION_MODEL == 2
collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
#endif
#endif
}

71
tests/test_d2q9_visual.py
View File

@ -2,12 +2,13 @@
"""
D2Q9 Regression Test — Poiseuille Channel + Cylinder Flow
==========================================================
Uses the ORIGINAL kernel.cu with same grid / BCs.
Uses kernel_v2.cu by default (legacy kernel.cu remains optional fallback).
Produces matplotlib figures for visual validation.
Usage:
python tests/test_d2q9_visual.py --device 2
python tests/test_d2q9_visual.py --device 2 --cylinder
python tests/test_d2q9_visual.py --device 2 --legacy
"""
import sys, os, argparse, time
@ -47,7 +48,7 @@ INTERFACE_FLAG = 0b00001000
# ━━━━━━━━━━━━━━━━━━━━━━━ Helpers ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
def configure_macros(n_objs=0):
"""Write macros.h to match original kernel settings."""
"""Write macros.h for D2Q9 regression (v2 default, legacy optional)."""
lines = compiler.read_lines(compiler.kernel_path("macros.h"))
defs = {
'MULT_GPU': 'False', 'NT': NT,
@ -58,23 +59,31 @@ def configure_macros(n_objs=0):
'DIM': DIM, 'NQ': NQ,
'VIS': VIS, '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, # keep unshifted for v2 defaults
}
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.
"""
def extract_fields(ddf_host, use_ddf_shifting=False):
"""Compute rho, u, v from host DDF in original D2Q9 direction ordering."""
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
if use_ddf_shifting:
rho = np.sum(f, axis=0) + RHO
denom = np.full_like(rho, RHO)
else:
rho = np.sum(f, axis=0)
denom = rho
denom_safe = np.where(np.abs(denom) > 1e-12, denom, 1.0)
u = (f[1] + f[5] + f[8] - f[3] - f[6] - f[7]) / denom_safe
v = (f[2] + f[5] + f[6] - f[4] - f[7] - f[8]) / denom_safe
u = np.where(np.abs(denom) > 1e-12, u, 0.0)
v = np.where(np.abs(denom) > 1e-12, v, 0.0)
return rho, u, v
@ -132,7 +141,7 @@ def build_cylinder_data(cx, cy, radius):
# ━━━━━━━━━━━━━━━━━━━━━━━ Simulation ━━━━━━━━━━━━━━━━━━━━━━━━━━
def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host):
def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host, use_legacy=False):
"""Compile kernel, run LBM, return DDF on host."""
cuda.init()
dev = cuda.Device(device_id)
@ -141,11 +150,17 @@ def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host)
try:
configure_macros(n_objs)
if use_legacy:
compiler.compile_kernel()
ptx_path = compiler.kernel_path("kernel.ptx")
init_name = "InitTubeFlow"
else:
compiler.compile_kernel_v2()
ptx_path = compiler.kernel_path("kernel_v2.ptx")
init_name = "InitTubeFlow_v2"
mod = cuda.module_from_file(ptx_path)
step_fn = mod.get_function("OneStep")
init_fn = mod.get_function("InitTubeFlow")
init_fn = mod.get_function(init_name)
nbytes_ddf = TOTAL * NQ * 4
ddf_gpu = cuda.mem_alloc(nbytes_ddf)
@ -198,9 +213,9 @@ def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host)
# ━━━━━━━━━━━━━━━━━━━━━━━ Visualization ━━━━━━━━━━━━━━━━━━━━━━━
def plot_poiseuille(ddf, flag, out_path):
def plot_poiseuille(ddf, flag, out_path, use_ddf_shifting=False):
"""3-panel figure: velocity mag, u(y) profile, pressure along centerline."""
rho, u, v = extract_fields(ddf)
rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
vel_mag = np.sqrt(u**2 + v**2)
# Mask solid cells for display
@ -244,9 +259,9 @@ def plot_poiseuille(ddf, flag, out_path):
plt.close(fig)
def plot_cylinder(ddf, flag, cx, cy, radius, out_path):
def plot_cylinder(ddf, flag, cx, cy, radius, out_path, use_ddf_shifting=False):
"""3-panel figure: velocity magnitude (zoom), vorticity, streamlines."""
rho, u, v = extract_fields(ddf)
rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
vel_mag = np.sqrt(u**2 + v**2)
mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool)
@ -313,6 +328,8 @@ 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('--legacy', action='store_true',
help='Use legacy kernel.cu path (default uses kernel_v2.cu)')
parser.add_argument('--cylinder', action='store_true',
help='Also run cylinder flow test')
parser.add_argument('--steps-pois', type=int, default=5000,
@ -324,6 +341,10 @@ def main():
out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output')
os.makedirs(out_dir, exist_ok=True)
use_ddf_shifting = bool(args.legacy)
mode = 'legacy kernel.cu' if args.legacy else 'kernel_v2.cu'
print(f"\n[Mode] {mode}")
# ---- Test 1: Poiseuille ----
print("\n===== Test 1: Poiseuille Channel Flow =====")
flag_pois = np.ones(TOTAL, dtype=np.uint8)
@ -331,11 +352,13 @@ def main():
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'))
flag_pois, indx_pois, delta_pois,
use_legacy=args.legacy)
plot_poiseuille(ddf, flag, os.path.join(out_dir, 'poiseuille_d2q9.png'),
use_ddf_shifting=use_ddf_shifting)
# Error metric
rho, u, v = extract_fields(ddf)
rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
y_arr = np.arange(NY, dtype=float)
u_ana = analytical_poiseuille(y_arr)
x_mid = NX // 2
@ -351,9 +374,11 @@ def main():
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)
flag_cyl, indx_cyl, delta_cyl,
use_legacy=args.legacy)
plot_cylinder(ddf2, flag2, cyl_cx, cyl_cy, cyl_r,
os.path.join(out_dir, 'cylinder_d2q9.png'))
os.path.join(out_dir, 'cylinder_d2q9.png'),
use_ddf_shifting=use_ddf_shifting)
print("\nDone.")

547
tests/test_high_re_validation.py Normal file
View File

@ -0,0 +1,547 @@
#!/usr/bin/env python3
"""
High-Re Validation (kernel_v2)
==============================
Unified validation script for high-Re runs with optional LES.
Default targets:
- 2D D2Q9: Re=5000
- 3D D3Q19: Re=3000
The script configures macros.h temporarily, compiles kernel_v2, runs the case,
and restores macros.h automatically.
"""
import argparse
import json
import os
import struct
import sys
import time
sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src"))
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pycuda.driver as cuda
from CelerisLab.cuda import compiler
FLUID_FLAG = 0x01
SOLID_FLAG = 0x02
OBSTACLE_FLAG = 0x04
def collision_name(model):
return {0: "SRT", 1: "TRT", 2: "MRT"}.get(model, f"M{model}")
def make_case_tag(cfg):
les_tag = "LES" if cfg["use_les"] else "NoLES"
return (
f"{cfg['name']}_Re{int(cfg['target_re'])}_"
f"{collision_name(cfg['collision_model'])}_{les_tag}_"
f"OM{int(cfg['outlet_mode'])}_WMAX{cfg['omega_collision_max']:.3f}"
)
def validate_case(rho):
nan_count = int(np.isnan(rho).sum())
if nan_count > 0:
return False, "NaN detected"
rho_min = float(np.min(rho))
rho_max = float(np.max(rho))
if rho_min <= 0.0:
return False, "Non-positive density"
if rho_max >= 2.0:
return False, "Density blow-up"
return True, "OK"
def plot_case(cfg, host_ddf, out_dir):
nq = cfg["nq"]
nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
flag = cfg["flag"]
tag = make_case_tag(cfg)
out_path = os.path.join(out_dir, f"{tag}.png")
if nq == 9:
f = host_ddf.reshape(nq, ny, nx)
rho = np.sum(f, axis=0)
ux = np.zeros_like(rho)
uy = np.zeros_like(rho)
cx = [0, 1, -1, 0, 0, 1, -1, 1, -1]
cy = [0, 0, 0, 1, -1, 1, -1, -1, 1]
for i in range(nq):
ux += cx[i] * f[i]
uy += cy[i] * f[i]
rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
ux /= rho_safe
uy /= rho_safe
vel = np.sqrt(ux * ux + uy * uy)
mask = flag.reshape(ny, nx) != FLUID_FLAG
vel_m = np.ma.array(vel, mask=mask)
vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
vort_m = np.ma.array(vort, mask=mask)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
im0 = axes[0].imshow(vel_m, origin="lower", aspect="auto", cmap="turbo")
plt.colorbar(im0, ax=axes[0], label="|u|")
axes[0].set_title("Velocity Magnitude")
vmax = np.percentile(np.abs(vort[~mask]), 99) if np.any(~mask) else 1e-6
vmax = max(vmax, 1.0e-6)
im1 = axes[1].imshow(vort_m, origin="lower", aspect="auto", cmap="RdBu_r", vmin=-vmax, vmax=vmax)
plt.colorbar(im1, ax=axes[1], label="vorticity")
axes[1].set_title("Vorticity")
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
ux_s = np.ma.array(ux, mask=mask)
uy_s = np.ma.array(uy, mask=mask)
speed = np.ma.sqrt(ux_s * ux_s + uy_s * uy_s)
axes[2].streamplot(X, Y, ux_s, uy_s, color=speed, cmap="viridis", density=2.0, linewidth=0.7)
axes[2].set_xlim(0, nx)
axes[2].set_ylim(0, ny)
axes[2].set_title("Streamlines")
fig.suptitle(tag)
fig.tight_layout()
fig.savefig(out_path, dpi=150)
plt.close(fig)
return out_path
# D3Q19: visualize mid-z slice
f = host_ddf.reshape(nq, nz, ny, nx)
z0 = nz // 2
fs = f[:, z0, :, :]
rho = np.sum(fs, axis=0)
ux = np.zeros_like(rho)
uy = np.zeros_like(rho)
uz = np.zeros_like(rho)
cx = np.array([0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0])
cy = np.array([0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1])
cz = np.array([0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1])
for i in range(nq):
ux += cx[i] * fs[i]
uy += cy[i] * fs[i]
uz += cz[i] * fs[i]
rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
ux /= rho_safe
uy /= rho_safe
uz /= rho_safe
vel = np.sqrt(ux * ux + uy * uy + uz * uz)
mask3 = flag.reshape(nz, ny, nx)[z0] != FLUID_FLAG
vel_m = np.ma.array(vel, mask=mask3)
vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
vort_m = np.ma.array(vort, mask=mask3)
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
im0 = axes[0].imshow(vel_m, origin="lower", aspect="auto", cmap="turbo")
plt.colorbar(im0, ax=axes[0], label="|u|")
axes[0].set_title("Velocity Magnitude (z-mid)")
vmax = np.percentile(np.abs(vort[~mask3]), 99) if np.any(~mask3) else 1.0e-6
vmax = max(vmax, 1.0e-6)
im1 = axes[1].imshow(vort_m, origin="lower", aspect="auto", cmap="RdBu_r", vmin=-vmax, vmax=vmax)
plt.colorbar(im1, ax=axes[1], label="vorticity")
axes[1].set_title("Vorticity (z-mid)")
X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
ux_s = np.ma.array(ux, mask=mask3)
uy_s = np.ma.array(uy, mask=mask3)
speed = np.ma.sqrt(ux_s * ux_s + uy_s * uy_s)
axes[2].streamplot(X, Y, ux_s, uy_s, color=speed, cmap="viridis", density=2.0, linewidth=0.7)
axes[2].set_xlim(0, nx)
axes[2].set_ylim(0, ny)
axes[2].set_title("Streamlines (z-mid)")
fig.suptitle(tag)
fig.tight_layout()
fig.savefig(out_path, dpi=150)
plt.close(fig)
return out_path
def compute_vis_omega(reynolds, diameter, u0):
vis = u0 * diameter / reynolds
omega = 1.0 / (3.0 * vis + 0.5)
return vis, omega
def set_macros(nx, ny, nz, dim, nq, vis, u0, collision_model, use_les, les_cs,
outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
omega_collision_max):
lines = compiler.read_lines(compiler.kernel_path("macros.h"))
defs = {
"MULT_GPU": "False",
"NT": 128,
"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": "1.0",
"U0": u0,
"N_OBJS": 0,
"COLLISION_MODEL": collision_model,
"STREAMING_MODEL": 0,
"STORE_PRECISION": 0,
"USE_DDF_SHIFTING": 0,
"USE_LES": int(use_les),
"LES_CS": f"{les_cs:.6f}f",
"INLET_PROFILE": 0,
"OUTLET_MODE": int(outlet_mode),
"OUTLET_BACKFLOW_CLAMP": int(outlet_backflow_clamp),
"OUTLET_BLEND_ALPHA": f"{float(outlet_blend_alpha):.3f}f",
"OMEGA_COLLISION_MAX": f"{float(omega_collision_max):.3f}f",
}
for name, value in defs.items():
lines = compiler.modify_macro(lines, name, value)
compiler.write_lines(compiler.kernel_path("macros.h"), lines)
def build_flags_2d(nx, ny, cx, cy, radius):
n = nx * ny
flag = np.ones(n, dtype=np.uint8) * FLUID_FLAG
for y in range(ny):
for x in range(nx):
k = y * nx + x
if y == 0 or y == ny - 1 or x == 0 or x == nx - 1:
flag[k] = SOLID_FLAG
elif (x - cx) ** 2 + (y - cy) ** 2 < radius ** 2:
flag[k] = OBSTACLE_FLAG
return flag
def build_flags_3d(nx, ny, nz, cx, cy, radius):
n = nx * ny * nz
flag = np.ones(n, 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
if y == 0 or y == ny - 1 or x == 0 or x == nx - 1:
flag[k] = SOLID_FLAG
elif (x - cx) ** 2 + (y - cy) ** 2 < radius ** 2:
flag[k] = OBSTACLE_FLAG
return flag
def run_case(device_id, cfg):
nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
dim, nq = cfg["dim"], cfg["nq"]
n = nx * ny * nz
set_macros(
nx=nx,
ny=ny,
nz=nz,
dim=dim,
nq=nq,
vis=cfg["vis"],
u0=cfg["u0"],
collision_model=cfg["collision_model"],
use_les=cfg["use_les"],
les_cs=cfg["les_cs"],
outlet_mode=cfg["outlet_mode"],
outlet_backflow_clamp=cfg["outlet_backflow_clamp"],
outlet_blend_alpha=cfg["outlet_blend_alpha"],
omega_collision_max=cfg["omega_collision_max"],
)
compiler.compile_kernel_v2()
cuda.init()
dev = cuda.Device(device_id)
ctx = dev.make_context()
try:
mod = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
init_fn = mod.get_function("InitTubeFlow_v2")
step_fn = mod.get_function("OneStep")
params_ptr, params_size = mod.get_global("d_params")
params_data = struct.pack(
"IIIQfffffffI",
nx,
ny,
nz,
n,
cfg["omega"],
1.1,
0.0,
0.0,
0.0,
1.0,
cfg["u0"],
0,
)
if len(params_data) < params_size:
params_data += b"\x00" * (params_size - len(params_data))
cuda.memcpy_htod(params_ptr, params_data)
fsize = n * nq * 4
d_fi = cuda.mem_alloc(fsize)
d_fi2 = cuda.mem_alloc(fsize)
d_flag = cuda.mem_alloc(n)
d_indx = cuda.mem_alloc(n * 4)
d_delta = cuda.mem_alloc(4)
d_action = cuda.mem_alloc(4)
d_obs = cuda.mem_alloc(4)
cuda.memset_d32(d_indx, 0, n)
cuda.memset_d32(d_delta, 0, 1)
cuda.memset_d32(d_action, 0, 1)
cuda.memset_d32(d_obs, 0, 1)
block = (128, 1, 1)
grid = ((nx + 127) // 128, ny, nz)
init_fn(d_flag, d_fi, block=block, grid=grid)
cuda.memcpy_dtod(d_fi2, d_fi, fsize)
cuda.memcpy_htod(d_flag, cfg["flag"])
t0 = time.time()
for step in range(cfg["steps"]):
step_fn(d_flag, d_fi, d_fi2, d_indx, d_delta, d_action, d_obs, block=block, grid=grid)
d_fi, d_fi2 = d_fi2, d_fi
if (step + 1) % cfg["report_every"] == 0:
cuda.Context.synchronize()
host = np.empty(n * nq, dtype=np.float32)
cuda.memcpy_dtoh(host, d_fi)
if nq == 9:
rho = host.reshape(nq, ny, nx).sum(axis=0)
c = float(rho[ny // 2, nx // 2])
else:
rho = host.reshape(nq, nz, ny, nx).sum(axis=0)
c = float(rho[nz // 2, ny // 2, nx // 2])
nan_count = int(np.isnan(rho).sum())
print(f" step {step+1:7d}: rho_center={c:.6f}, nan={nan_count}")
if nan_count > 0:
break
cuda.Context.synchronize()
elapsed = time.time() - t0
host = np.empty(n * nq, dtype=np.float32)
cuda.memcpy_dtoh(host, d_fi)
if nq == 9:
rho = host.reshape(nq, ny, nx).sum(axis=0)
center = float(rho[ny // 2, nx // 2])
else:
rho = host.reshape(nq, nz, ny, nx).sum(axis=0)
center = float(rho[nz // 2, ny // 2, nx // 2])
ok, reason = validate_case(rho)
plot_path = None
if cfg.get("save_plot", True):
plot_path = plot_case(cfg, host, cfg["out_dir"])
return {
"case_tag": make_case_tag(cfg),
"name": cfg["name"],
"target_re": cfg["target_re"],
"steps": cfg["steps"],
"mlups": float(n * cfg["steps"] / elapsed / 1e6),
"nan_count": int(np.isnan(rho).sum()),
"rho_center": center,
"rho_min": float(np.nanmin(rho)),
"rho_max": float(np.nanmax(rho)),
"omega": cfg["omega"],
"vis": cfg["vis"],
"collision_model": cfg["collision_model"],
"use_les": bool(cfg["use_les"]),
"les_cs": float(cfg["les_cs"]),
"outlet_mode": int(cfg["outlet_mode"]),
"outlet_backflow_clamp": int(cfg["outlet_backflow_clamp"]),
"outlet_blend_alpha": float(cfg["outlet_blend_alpha"]),
"omega_collision_max": float(cfg["omega_collision_max"]),
"pass": bool(ok),
"reason": reason,
"plot_path": plot_path,
}
finally:
ctx.pop()
def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
omega_collision_max):
nx, ny, nz = 512, 256, 1
cx, cy, radius = 128.0, 128.0, 24.0
u0 = 0.03
vis, omega = compute_vis_omega(re2d, 2.0 * radius, u0)
return {
"name": "2D_D2Q9_highRe",
"dim": 2,
"nq": 9,
"nx": nx,
"ny": ny,
"nz": nz,
"flag": build_flags_2d(nx, ny, cx, cy, radius),
"u0": u0,
"vis": vis,
"omega": omega,
"steps": steps2d,
"report_every": max(steps2d // 10, 1),
"collision_model": collision_model,
"use_les": use_les,
"les_cs": les_cs,
"outlet_mode": int(outlet_mode),
"outlet_backflow_clamp": int(outlet_backflow_clamp),
"outlet_blend_alpha": float(outlet_blend_alpha),
"omega_collision_max": float(omega_collision_max),
"target_re": re2d,
"save_plot": True,
"out_dir": out_dir,
}
def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
omega_collision_max):
nx, ny, nz = 256, 128, 32
cx, cy, radius = 64.0, 64.0, 12.0
u0 = 0.04
vis, omega = compute_vis_omega(re3d, 2.0 * radius, u0)
return {
"name": "3D_D3Q19_highRe",
"dim": 3,
"nq": 19,
"nx": nx,
"ny": ny,
"nz": nz,
"flag": build_flags_3d(nx, ny, nz, cx, cy, radius),
"u0": u0,
"vis": vis,
"omega": omega,
"steps": steps3d,
"report_every": max(steps3d // 10, 1),
"collision_model": collision_model,
"use_les": use_les,
"les_cs": les_cs,
"outlet_mode": int(outlet_mode),
"outlet_backflow_clamp": int(outlet_backflow_clamp),
"outlet_blend_alpha": float(outlet_blend_alpha),
"omega_collision_max": float(omega_collision_max),
"target_re": re3d,
"save_plot": True,
"out_dir": out_dir,
}
def build_comprehensive_cases(args, out_dir):
cases = []
# Coverage matrix at moderate Re to verify all changed pathways.
for cm in (0, 1, 2):
for les in (0, 1):
cases.append(build_case_2d(re2d=200.0, steps2d=args.matrix_steps2d,
collision_model=cm, use_les=les,
les_cs=args.les_cs, out_dir=out_dir,
outlet_mode=args.outlet_mode,
outlet_backflow_clamp=1,
outlet_blend_alpha=args.outlet_blend_alpha,
omega_collision_max=args.omega_collision_max))
cases.append(build_case_3d(re3d=200.0, steps3d=args.matrix_steps3d,
collision_model=cm, use_les=les,
les_cs=args.les_cs, out_dir=out_dir,
outlet_mode=args.outlet_mode,
outlet_backflow_clamp=1,
outlet_blend_alpha=args.outlet_blend_alpha,
omega_collision_max=args.omega_collision_max))
return cases
def main():
parser = argparse.ArgumentParser(description="High-Re validation for kernel_v2")
parser.add_argument("--device", type=int, default=0)
parser.add_argument("--re2d", type=float, default=5000.0)
parser.add_argument("--re3d", type=float, default=3000.0)
parser.add_argument("--steps2d", type=int, default=10000)
parser.add_argument("--steps3d", type=int, default=20000)
parser.add_argument("--collision", type=int, default=1, choices=[0, 1, 2],
help="0=SRT, 1=TRT, 2=MRT")
parser.add_argument("--use-les", action="store_true", default=True,
help="Enable Smagorinsky LES")
parser.add_argument("--no-les", action="store_false", dest="use_les")
parser.add_argument("--les-cs", type=float, default=0.16)
parser.add_argument("--outlet-mode", type=int, default=0, choices=[0, 1, 2],
help="0=non-equilibrium extrapolation, 1=zero-gradient copy, 2=damped blend")
parser.add_argument("--outlet-blend-alpha", type=float, default=0.70,
help="Blend alpha for outlet-mode 2")
parser.add_argument("--omega-collision-max", type=float, default=1.999,
help="Upper clamp for collision omega")
parser.add_argument("--only", choices=["2d", "3d", "both"], default="both")
parser.add_argument("--comprehensive", action="store_true",
help="Run coverage matrix: SRT/TRT/MRT x LES on/off for 2D and 3D")
parser.add_argument("--matrix-steps2d", type=int, default=1000)
parser.add_argument("--matrix-steps3d", type=int, default=600)
args = parser.parse_args()
macro_path = compiler.kernel_path("macros.h")
macro_backup = compiler.read_lines(macro_path)
out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "output")
os.makedirs(out_dir, exist_ok=True)
out_json = os.path.join(out_dir, "high_re_validation_summary.json")
try:
results = []
if args.only in ("2d", "both"):
c2 = build_case_2d(args.re2d, args.steps2d, args.collision, args.use_les,
args.les_cs, out_dir, args.outlet_mode, 1,
args.outlet_blend_alpha, args.omega_collision_max)
print("\n=== Running 2D high-Re case ===")
print(f" target Re={args.re2d:.1f}, vis={c2['vis']:.6e}, omega={c2['omega']:.6f}")
results.append(run_case(args.device, c2))
if args.only in ("3d", "both"):
c3 = build_case_3d(args.re3d, args.steps3d, args.collision, args.use_les,
args.les_cs, out_dir, args.outlet_mode, 1,
args.outlet_blend_alpha, args.omega_collision_max)
print("\n=== Running 3D high-Re case ===")
print(f" target Re={args.re3d:.1f}, vis={c3['vis']:.6e}, omega={c3['omega']:.6f}")
results.append(run_case(args.device, c3))
if args.comprehensive:
print("\n=== Running comprehensive coverage matrix ===")
for cfg in build_comprehensive_cases(args, out_dir):
print(f" {cfg['name']} Re={cfg['target_re']:.1f} "
f"{collision_name(cfg['collision_model'])} LES={int(cfg['use_les'])}")
results.append(run_case(args.device, cfg))
with open(out_json, "w", encoding="utf-8") as f:
json.dump(results, f, indent=2)
print("\n=== Summary ===")
n_pass = 0
for r in results:
if r["pass"]:
n_pass += 1
print(f"{r['name']}: nan={r['nan_count']}, rho_center={r['rho_center']:.6f}, "
f"rho[min,max]=[{r['rho_min']:.6f}, {r['rho_max']:.6f}], "
f"MLUPS={r['mlups']:.1f}, pass={r['pass']} ({r['reason']})")
if r.get("plot_path"):
print(f" plot: {r['plot_path']}")
print(f"Pass rate: {n_pass}/{len(results)}")
print(f"Saved: {out_json}")
finally:
compiler.write_lines(macro_path, macro_backup)
if __name__ == "__main__":
main()