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

Add TRT-LES diagnostics and plotting functionality

Browse Source 
- Introduced new functions for computing and plotting TRT-LES fields.
- Enhanced case diagnostics to include mass drift and inlet variance metrics.
- Updated configuration to support inlet profile selection and TRT magic parameter.
- Modified existing functions to accommodate new diagnostic calculations.
- Improved case tagging to include inlet profile and TRT parameters.
- Added checks for fluid dynamics and diagnostics in the run_case function.
This commit is contained in:
Frank14f 2026-04-06 15:57:44 +08:00
parent 5076e5d789
commit b339f302d3
16 changed files with 632 additions and 773 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -28,6 +28,7 @@ wheels/
# IDE
.vscode/
.github/
.idea/
*.swp
*.swo

BIN
output/cylinder_d2q9.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 513 KiB

BIN
output/cylinder_d3q19_re100.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 633 KiB

BIN
output/cylinder_d3q19_re20.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 581 KiB

BIN
output/poiseuille_d2q9.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 122 KiB

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

@ -194,11 +194,16 @@ class FlowField:
else:
raise ValueError(f"Unsupported data type {self.DATA_TYPE}.")
# Ensure host-side DDF mirrors current device state before local edits.
cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
for x in range(int(x_c - radius) - 1, int(x_c + radius) + 1):
for y in range(int(y_c - radius) - 1, int(y_c + radius) + 1):
if (x - x_c) ** 2 + (y - y_c) ** 2 < radius**2:
k = x + y * self.FIELD_SHAPE[0]
self.flag[k] = SOLID
for i in range(self.LATTICE):
self.ddf[k + i * self.FIELD_SIZE] = self.WW[i]
delta_temp = np.zeros(11, dtype=self.DATA_TYPE)
delta_temp[0] = id_object.view(self.DATA_TYPE)
for i in range(self.LATTICE):
@ -245,6 +250,8 @@ class FlowField:
cuda.memcpy_htod(self.delta_gpu, self.delta_curve)
cuda.memcpy_htod(self.flag_gpu, self.flag)
cuda.memcpy_htod(self.indx_gpu, self.indx)
cuda.memcpy_htod(self.ddf_gpu, self.ddf)
cuda.memcpy_htod(self.temp_gpu, self.ddf)
self._rebuild_kernel()

19
src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh
View File

@ -43,6 +43,25 @@ __device__ inline void apply_bounce_back(float* __restrict__ f,
// 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_y_pull(unsigned int y,
float* __restrict__ f,
const fpxx* __restrict__ fi_in,
unsigned long k)
{
if (y == 1) {
// Directions sourced from y=0 wall in pull step: +y, +x+y, -x+y.
f[3] = load_ddf(fi_in, index_f(k, 4u));
f[5] = load_ddf(fi_in, index_f(k, 6u));
f[8] = load_ddf(fi_in, index_f(k, 7u));
}
else if (y == (unsigned int)(NY - 2)) {
// Directions sourced from y=NY-1 wall in pull step: -y, -x-y, +x-y.
f[4] = load_ddf(fi_in, index_f(k, 3u));
f[6] = load_ddf(fi_in, index_f(k, 5u));
f[7] = load_ddf(fi_in, index_f(k, 8u));
}
}
__device__ inline void apply_wall_bb_d2q9(unsigned int y,
float* __restrict__ f)
{

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

@ -31,6 +31,14 @@
#define OUTLET_BLEND_ALPHA 0.70f
#endif
#ifndef OUTLET_SRT_NEQ_DAMP
#define OUTLET_SRT_NEQ_DAMP 0.50f
#endif
#ifndef INLET_TRT_NEQ_DAMP
#define INLET_TRT_NEQ_DAMP 0.50f
#endif
__device__ __forceinline__ float inlet_target_u(float y_coord) {
#if INLET_PROFILE == 0
return U0;
@ -69,10 +77,22 @@ __device__ inline void apply_parabolic_inlet(float* __restrict__ f,
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] - 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];
#if COLLISION_MODEL == 1
// TRT path: reconstruct full population set at inlet using damped donor
// non-equilibrium transport. This follows the high-Re stable family that
// replaces all boundary-node populations, reducing odd-mode contamination.
const float beta = INLET_TRT_NEQ_DAMP;
#pragma unroll
for (int i = 0; i < 9; i++) {
const float fneq = f_neb[i] - feq_neb[i];
f[i] = feq_tar[i] + beta * fneq;
}
#else
const float beta = 1.0f;
f[1] = feq_tar[1] + beta * (f_neb[1] - feq_neb[1]);
f[5] = feq_tar[5] + beta * (f_neb[5] - feq_neb[5]);
f[7] = feq_tar[7] + beta * (f_neb[7] - feq_neb[7]);
#endif
}
// ---------------------------------------------------------------------------
@ -94,20 +114,30 @@ __device__ inline void apply_pressure_outlet(float* __restrict__ f,
f[6] = f_neb[6];
#else
// Convective non-reflecting style: extrapolate outlet density from neighbor
// and keep neighbor velocity.
// Prescribed-pressure outlet: keep neighbor velocity and impose reference
// outlet density for NEQ reconstruction.
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;
float rho_out = RHO;
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);
#if OUTLET_MODE == 2
#if COLLISION_MODEL == 0
// SRT path: use full-population damped NEQ reconstruction at outlet to
// suppress checkerboard/grid noise from high-frequency donor content.
// Reference: high-Re outlet regularization rationale and NEQ decomposition.
const float beta = OUTLET_SRT_NEQ_DAMP;
#pragma unroll
for (int i = 0; i < 9; i++) {
const float fneq = f_neb[i] - feq_neb[i];
f[i] = feq_tar[i] + beta * fneq;
}
#elif 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];
@ -178,15 +208,22 @@ __device__ inline void apply_pressure_outlet_3d(float* __restrict__ f,
un = fmaxf(un, 0.0f);
#endif
// Convective non-reflecting style: extrapolate outlet density from neighbor
// and keep neighbor velocity.
float rho_out = rho_neb;
// Prescribed-pressure outlet: keep neighbor velocity and impose reference
// outlet density for NEQ reconstruction.
float rho_out = RHO;
float feq_tar[19], feq_neb[19];
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
#if COLLISION_MODEL == 0
const float beta = OUTLET_SRT_NEQ_DAMP;
#pragma unroll
for (int i = 0; i < 19; i++) {
const float fneq = f_neb[i] - feq_neb[i];
f[i] = feq_tar[i] + beta * fneq;
}
#elif 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];

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

@ -91,9 +91,8 @@ __global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi)
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)));
// Keep startup field consistent with runtime inlet configuration.
float u_init = inlet_target_u((float)y);
if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
flag[k] = LEGACY_SOLID;
@ -116,9 +115,8 @@ __global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi)
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)));
// Keep startup field consistent with runtime inlet configuration.
float u_init = inlet_target_u((float)y);
if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
flag[k] = LEGACY_SOLID;
@ -217,9 +215,14 @@ __global__ void OneStep(
}
}
if ((fl & LEGACY_FLUID) && (y == 1u || y == (unsigned int)(NY - 2))) {
apply_wall_bb_d2q9_y_pull(y, f, fi_in, k);
}
// Collision
if (fl & LEGACY_FLUID) {
float feq[NQ], Fin[NQ];
compute_rho_u(f, rho_n, ux, uy);
compute_feq(rho_n, ux, uy, feq);
zero_forcing(Fin);
float omega_col = d_params.omega;
@ -298,9 +301,14 @@ __global__ void OneStep(
}
}
if ((fl & LEGACY_FLUID) && (y == 1 || y == NY - 2)) {
apply_wall_bb_d3q19_y_pull(y, f, fi_in, k);
}
// ---- Collision (fluid only) ----
if (fl & LEGACY_FLUID) {
float feq[NQ], Fin[NQ];
compute_rho_u(f, rho_n, ux, uy, uz);
compute_feq(rho_n, ux, uy, uz, feq);
zero_forcing(Fin);
float omega_col = d_params.omega;

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

@ -101,3 +101,8 @@
#ifndef OMEGA_COLLISION_MAX
#define OMEGA_COLLISION_MAX 1.999f
#endif
// TRT magic parameter Lambda used to map omega+ -> omega-
#ifndef TRT_MAGIC_PARAM
#define TRT_MAGIC_PARAM 0.1875f
#endif

5
src/CelerisLab/lbm/kernels/operators/collision_trt.cuh
View File

@ -13,8 +13,11 @@
#ifndef CELERIS_OPERATORS_COLLISION_TRT_CUH
#define CELERIS_OPERATORS_COLLISION_TRT_CUH
// Magic parameter Λ = 3/16 (optimal for porous-media / bounce-back wall location)
// Magic parameter Λ for TRT wall/collision coupling.
// Default keeps previous behavior; can be overridden in macros.h / tests.
#ifndef TRT_MAGIC_PARAM
#define TRT_MAGIC_PARAM (0.1875f)
#endif
__device__ __forceinline__ float compute_omega_minus(float omega_plus) {
return 1.0f / (TRT_MAGIC_PARAM / (1.0f / omega_plus - 0.5f) + 0.5f);

103
src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
View File

@ -19,8 +19,21 @@
#define LES_CS 0.16f
#endif
#ifndef LES_FP_ITERS
#define LES_FP_ITERS 3
#endif
#ifndef LES_FP_RELAX
#define LES_FP_RELAX 0.70f
#endif
#ifndef LES_NUT_MAX_RATIO
#define LES_NUT_MAX_RATIO 20.0f
#endif
__device__ __forceinline__ float clamp_omega(float w) {
return fminf(1.999f, fmaxf(0.01f, w));
return fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, w));
}
#if NQ == 9
@ -32,6 +45,7 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
{
const float rho_safe = fmaxf(rho, 1.0e-12f);
const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
// Πneq = Σ (ci ci)(fi-feqi)
float pixx = 0.0f, piyy = 0.0f, pixy = 0.0f;
@ -46,18 +60,36 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
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;
// Self-consistent tau iteration for LES/TRT coupling.
// S is estimated from Pi_neq with current tau_eff, then nu_t updates tau_eff.
float tau_eff = tau0;
const float tau_max = 1.0f / fmaxf(OMEGA_COLLISION_MIN, 1.0e-6f);
// |S| = sqrt(2 Sij Sij)
const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + 2.0f*sxy*sxy));
#pragma unroll
for (int it = 0; it < LES_FP_ITERS; ++it) {
const float denom = 2.0f * rho_safe * CS2 * fmaxf(tau_eff, 0.500001f);
const float sxx = -pixx / denom;
const float syy = -piyy / denom;
const float sxy = -pixy / denom;
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));
// Use deviatoric strain to reduce SGS sensitivity to bulk/pressure-wave content.
const float tr = 0.5f * (sxx + syy);
const float sxx_dev = sxx - tr;
const float syy_dev = syy - tr;
const float sxy_dev = sxy;
// |S_dev| = sqrt(2 S_dev,ij S_dev,ij)
const float s_mag = sqrtf(2.0f * (sxx_dev*sxx_dev + syy_dev*syy_dev + 2.0f*sxy_dev*sxy_dev));
float nut = (LES_CS * LES_CS) * s_mag;
const float nut_cap = fmaxf(0.0f, LES_NUT_MAX_RATIO * nu0);
nut = fminf(fmaxf(nut, 0.0f), nut_cap);
const float tau_new = fminf(tau_max, fmaxf(0.500001f, 0.5f + 3.0f * (nu0 + nut)));
tau_eff = LES_FP_RELAX * tau_new + (1.0f - LES_FP_RELAX) * tau_eff;
}
return clamp_omega(1.0f / tau_eff);
}
#elif NQ == 19
@ -69,6 +101,7 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
{
const float rho_safe = fmaxf(rho, 1.0e-12f);
const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
// Πneq = Σ (ci ci)(fi-feqi)
float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
@ -88,22 +121,42 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
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;
// Self-consistent tau iteration for LES/TRT coupling.
float tau_eff = tau0;
const float tau_max = 1.0f / fmaxf(OMEGA_COLLISION_MIN, 1.0e-6f);
// |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)));
#pragma unroll
for (int it = 0; it < LES_FP_ITERS; ++it) {
const float denom = 2.0f * rho_safe * CS2 * fmaxf(tau_eff, 0.500001f);
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;
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));
// Use deviatoric strain to reduce SGS sensitivity to bulk/pressure-wave content.
const float tr = (sxx + syy + szz) / 3.0f;
const float sxx_dev = sxx - tr;
const float syy_dev = syy - tr;
const float szz_dev = szz - tr;
const float sxy_dev = sxy;
const float sxz_dev = sxz;
const float syz_dev = syz;
// |S_dev| = sqrt(2 S_dev,ij S_dev,ij)
const float s_mag = sqrtf(2.0f * (sxx_dev*sxx_dev + syy_dev*syy_dev + szz_dev*szz_dev
+ 2.0f*(sxy_dev*sxy_dev + sxz_dev*sxz_dev + syz_dev*syz_dev)));
float nut = (LES_CS * LES_CS) * s_mag;
const float nut_cap = fmaxf(0.0f, LES_NUT_MAX_RATIO * nu0);
nut = fminf(fmaxf(nut, 0.0f), nut_cap);
const float tau_new = fminf(tau_max, fmaxf(0.500001f, 0.5f + 3.0f * (nu0 + nut)));
tau_eff = LES_FP_RELAX * tau_new + (1.0f - LES_FP_RELAX) * tau_eff;
}
return clamp_omega(1.0f / tau_eff);
}
#endif // NQ

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

@ -128,12 +128,12 @@ __global__ void StreamCollideDouble(
// ----- Wall bounce-back (top/bottom walls) -----
#if NQ == 9
if (y == 1 || y == NY - 2) {
apply_wall_bb_d2q9(y, f);
if ((fl & LEGACY_FLUID) && (y == 1u || y == (unsigned int)(NY - 2))) {
apply_wall_bb_d2q9_y_pull(y, f, fi_in, k);
}
#elif NQ == 19
if (y == 1 || y == NY - 2) {
apply_wall_bb_d3q19_y(y, f);
if ((fl & LEGACY_FLUID) && (y == 1u || y == (unsigned int)(NY - 2))) {
apply_wall_bb_d3q19_y_pull(y, f, fi_in, k);
}
#endif
@ -143,6 +143,7 @@ __global__ void StreamCollideDouble(
float Fin[NQ];
#if NQ == 9
compute_rho_u(f, rho_n, ux, uy);
compute_feq(rho_n, ux, uy, feq);
zero_forcing(Fin);
float omega_col = d_params.omega;
@ -161,6 +162,7 @@ __global__ void StreamCollideDouble(
#endif
#elif NQ == 19
compute_rho_u(f, rho_n, ux, uy, uz);
compute_feq(rho_n, ux, uy, uz, feq);
zero_forcing(Fin);
float omega_col = d_params.omega;

387
tests/test_d2q9_visual.py
View File

@ -1,387 +0,0 @@
#!/usr/bin/env python3
"""
D2Q9 Regression Test — Poiseuille Channel + Cylinder Flow
==========================================================
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
# 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 for D2Q9 regression (v2 default, legacy optional)."""
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,
'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, use_ddf_shifting=False):
"""Compute rho, u, v from host DDF in original D2Q9 direction ordering."""
f = ddf_host.reshape(NQ, NY, NX)
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
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, use_legacy=False):
"""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)
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(init_name)
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, use_ddf_shifting=False):
"""3-panel figure: velocity mag, u(y) profile, pressure along centerline."""
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
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, use_ddf_shifting=False):
"""3-panel figure: velocity magnitude (zoom), vorticity, streamlines."""
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)
# 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('--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,
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)
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)
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,
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, use_ddf_shifting=use_ddf_shifting)
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,
use_legacy=args.legacy)
plot_cylinder(ddf2, flag2, cyl_cx, cyl_cy, cyl_r,
os.path.join(out_dir, 'cylinder_d2q9.png'),
use_ddf_shifting=use_ddf_shifting)
print("\nDone.")
if __name__ == '__main__':
main()

327
tests/test_d3q19_cylinder.py
View File

@ -1,327 +0,0 @@
#!/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()

460
tests/test_high_re_validation.py
View File

@ -33,6 +33,11 @@ FLUID_FLAG = 0x01
SOLID_FLAG = 0x02
OBSTACLE_FLAG = 0x04
OMEGA_COLLISION_MIN_DEFAULT = 0.01
LES_POST_FP_ITERS = 3
LES_POST_FP_RELAX = 0.70
LES_POST_NUT_MAX_RATIO = 20.0
def collision_name(model):
return {0: "SRT", 1: "TRT", 2: "MRT"}.get(model, f"M{model}")
@ -40,10 +45,13 @@ def collision_name(model):
def make_case_tag(cfg):
les_tag = "LES" if cfg["use_les"] else "NoLES"
inlet_tag = f"IP{int(cfg.get('inlet_profile', 0))}"
lam_tag = f"LAM{float(cfg.get('trt_magic_param', 0.1875)):.4f}"
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}"
f"OM{int(cfg['outlet_mode'])}_{inlet_tag}_{lam_tag}_"
f"WMAX{cfg['omega_collision_max']:.3f}"
)
@ -61,6 +69,366 @@ def validate_case(rho):
return True, "OK"
def lattice_weights(nq):
if nq == 9:
return np.array(
[4.0 / 9.0] + [1.0 / 9.0] * 4 + [1.0 / 36.0] * 4,
dtype=np.float32,
)
if nq == 19:
return np.array(
[1.0 / 3.0] + [1.0 / 18.0] * 6 + [1.0 / 36.0] * 12,
dtype=np.float32,
)
raise ValueError(f"Unsupported nq={nq}")
def impose_rest_state_on_nonfluid(cfg, host_ddf):
nq = cfg["nq"]
nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
w = lattice_weights(nq)
if nq == 9:
f = host_ddf.reshape(nq, ny, nx)
nonfluid = cfg["flag"].reshape(ny, nx) != FLUID_FLAG
for i in range(nq):
f[i, nonfluid] = w[i]
return host_ddf
f = host_ddf.reshape(nq, nz, ny, nx)
nonfluid = cfg["flag"].reshape(nz, ny, nx) != FLUID_FLAG
for i in range(nq):
f[i, nonfluid] = w[i]
return host_ddf
def compute_case_diagnostics(cfg, host_ddf):
nq = cfg["nq"]
nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
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)
fluid = cfg["flag"].reshape(ny, nx) == FLUID_FLAG
mass = float(np.nansum(rho[fluid]))
x0 = 1
x1 = min(nx - 1, 33)
inlet_window = np.zeros_like(fluid)
inlet_window[1:ny - 1, x0:x1] = True
win_mask = fluid & inlet_window
inlet_var = float(np.nanvar(vel[win_mask])) if np.any(win_mask) else float("nan")
# Quantify how well the first interior column follows the designed inlet profile.
x_probe = 1
line_mask = fluid[:, x_probe]
line_u = ux[:, x_probe]
y = np.arange(ny, dtype=np.float32)
if int(cfg.get("inlet_profile", 0)) == 0:
target = np.full(ny, float(cfg["u0"]), dtype=np.float32)
else:
yy = (y - 0.5 * (ny - 1)) / (ny - 2.0)
target = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy * yy)
if np.any(line_mask):
diff = line_u[line_mask] - target[line_mask]
t_ref = float(np.max(np.abs(target[line_mask])))
denom = max(t_ref, 1.0e-8)
inlet_rel_l2 = float(np.sqrt(np.mean(diff * diff)) / denom)
inlet_rel_linf = float(np.max(np.abs(diff)) / denom)
else:
inlet_rel_l2 = float("nan")
inlet_rel_linf = float("nan")
# Inlet-plane-wave indicator: streamwise oscillation of column-averaged
# macros in the pre-obstacle region.
cx = float(cfg.get("cx", 0.25 * nx))
radius = float(cfg.get("radius", ny / 12.0))
x_pre0 = 1
x_pre1 = min(nx - 2, max(x_pre0 + 4, int(cx - 1.5 * radius)))
col_u = []
col_r = []
for xp in range(x_pre0, x_pre1):
col_mask = fluid[1:ny - 1, xp]
if np.any(col_mask):
col_u.append(float(np.mean(ux[1:ny - 1, xp][col_mask])))
col_r.append(float(np.mean(rho[1:ny - 1, xp][col_mask])))
if int(cfg.get("inlet_profile", 0)) == 0:
u_target_mean = float(cfg["u0"])
else:
y_int = np.arange(1, ny - 1, dtype=np.float32)
yy_int = (y_int - 0.5 * (ny - 1)) / (ny - 2.0)
target_int = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy_int * yy_int)
u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])
if len(col_u) >= 4:
col_u_arr = np.array(col_u, dtype=np.float64)
col_r_arr = np.array(col_r, dtype=np.float64)
inlet_wave_ux_rel = float(np.std(col_u_arr - u_target_mean) / max(abs(u_target_mean), 1.0e-8))
rho_ref = max(abs(float(np.mean(col_r_arr))), 1.0e-8)
inlet_wave_rho_rel = float(np.std(col_r_arr) / rho_ref)
else:
inlet_wave_ux_rel = float("nan")
inlet_wave_rho_rel = float("nan")
# TRT checker/grid-noise indicator: odd-even imbalance in wake ux field.
xw0 = min(nx - 3, max(2, int(cx + 2.0 * radius)))
xw1 = min(nx - 2, max(xw0 + 4, int(cx + 12.0 * radius)))
yw0, yw1 = 2, ny - 2
if xw1 > xw0 + 2 and yw1 > yw0 + 2:
reg = ux[yw0:yw1, xw0:xw1].astype(np.float64)
reg_mask = fluid[yw0:yw1, xw0:xw1]
if np.any(reg_mask):
valid = reg[reg_mask]
m = float(np.mean(valid))
centered = np.where(reg_mask, reg - m, np.nan)
rms = float(np.sqrt(np.mean((valid - m) * (valid - m))))
yy_i, xx_i = np.indices(centered.shape)
even_vals = centered[((xx_i + yy_i) & 1) == 0]
odd_vals = centered[((xx_i + yy_i) & 1) == 1]
even_vals = even_vals[np.isfinite(even_vals)]
odd_vals = odd_vals[np.isfinite(odd_vals)]
if rms > 1.0e-12 and even_vals.size > 8 and odd_vals.size > 8:
wake_checker_rel = float(abs(np.mean(even_vals) - np.mean(odd_vals)) / rms)
else:
wake_checker_rel = float("nan")
pair_mask = reg_mask[:, :-1] & reg_mask[:, 1:]
a = centered[:, :-1][pair_mask]
b = centered[:, 1:][pair_mask]
if a.size > 16 and np.std(a) > 1.0e-12 and np.std(b) > 1.0e-12:
corr = float(np.corrcoef(a, b)[0, 1])
wake_checker_anti_corr_x = float(max(0.0, -corr))
else:
wake_checker_anti_corr_x = float("nan")
else:
wake_checker_rel = float("nan")
wake_checker_anti_corr_x = float("nan")
else:
wake_checker_rel = float("nan")
wake_checker_anti_corr_x = float("nan")
return {
"mass": mass,
"inlet_var": inlet_var,
"inlet_line_rel_l2": inlet_rel_l2,
"inlet_line_rel_linf": inlet_rel_linf,
"inlet_wave_ux_rel": inlet_wave_ux_rel,
"inlet_wave_rho_rel": inlet_wave_rho_rel,
"wake_checker_rel": wake_checker_rel,
"wake_checker_anti_corr_x": wake_checker_anti_corr_x,
}
f = host_ddf.reshape(nq, nz, ny, nx)
rho = np.sum(f, 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] * f[i]
uy += cy[i] * f[i]
uz += cz[i] * f[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)
fluid = cfg["flag"].reshape(nz, ny, nx) == FLUID_FLAG
mass = float(np.nansum(rho[fluid]))
x0 = 1
x1 = min(nx - 1, 17)
inlet_window = np.zeros_like(fluid)
inlet_window[:, 1:ny - 1, x0:x1] = True
win_mask = fluid & inlet_window
inlet_var = float(np.nanvar(vel[win_mask])) if np.any(win_mask) else float("nan")
cx = float(cfg.get("cx", 0.25 * nx))
radius = float(cfg.get("radius", ny / 12.0))
x_pre0 = 1
x_pre1 = min(nx - 2, max(x_pre0 + 4, int(cx - 1.5 * radius)))
col_u = []
col_r = []
for xp in range(x_pre0, x_pre1):
col_mask = fluid[:, 1:ny - 1, xp]
if np.any(col_mask):
col_u.append(float(np.mean(ux[:, 1:ny - 1, xp][col_mask])))
col_r.append(float(np.mean(rho[:, 1:ny - 1, xp][col_mask])))
if int(cfg.get("inlet_profile", 0)) == 0:
u_target_mean = float(cfg["u0"])
else:
y_int = np.arange(1, ny - 1, dtype=np.float32)
yy_int = (y_int - 0.5 * (ny - 1)) / (ny - 2.0)
target_int = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy_int * yy_int)
u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])
if len(col_u) >= 4:
col_u_arr = np.array(col_u, dtype=np.float64)
col_r_arr = np.array(col_r, dtype=np.float64)
inlet_wave_ux_rel = float(np.std(col_u_arr - u_target_mean) / max(abs(u_target_mean), 1.0e-8))
rho_ref = max(abs(float(np.mean(col_r_arr))), 1.0e-8)
inlet_wave_rho_rel = float(np.std(col_r_arr) / rho_ref)
else:
inlet_wave_ux_rel = float("nan")
inlet_wave_rho_rel = float("nan")
return {
"mass": mass,
"inlet_var": inlet_var,
"inlet_wave_ux_rel": inlet_wave_ux_rel,
"inlet_wave_rho_rel": inlet_wave_rho_rel,
"wake_checker_rel": float("nan"),
"wake_checker_anti_corr_x": float("nan"),
"inlet_line_rel_l2": float("nan"),
"inlet_line_rel_linf": float("nan"),
}
def compute_trt_les_fields_2d(cfg, host_ddf):
if cfg["nq"] != 9:
return None
nx, ny = cfg["nx"], cfg["ny"]
f = host_ddf.reshape(9, ny, nx).astype(np.float64)
cx = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float64).reshape(9, 1, 1)
cy = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float64).reshape(9, 1, 1)
w = lattice_weights(9).astype(np.float64).reshape(9, 1, 1)
rho = np.sum(f, axis=0)
rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
ux = np.sum(f * cx, axis=0) / rho_safe
uy = np.sum(f * cy, axis=0) / rho_safe
u2 = ux * ux + uy * uy
cu = 3.0 * (cx * ux[None, :, :] + cy * uy[None, :, :])
feq = w * rho[None, :, :] * (1.0 + cu + 0.5 * cu * cu - 1.5 * u2[None, :, :])
fneq = f - feq
pixx = np.sum(fneq * cx * cx, axis=0)
piyy = np.sum(fneq * cy * cy, axis=0)
pixy = np.sum(fneq * cx * cy, axis=0)
omega0 = float(cfg["omega"])
omega_min = float(cfg.get("omega_collision_min", OMEGA_COLLISION_MIN_DEFAULT))
omega_max = float(cfg["omega_collision_max"])
tau0 = max(1.0 / max(omega0, 1.0e-6), 0.500001)
nu0 = (tau0 - 0.5) * (1.0 / 3.0)
tau_max = 1.0 / max(omega_min, 1.0e-6)
tau_eff = np.full((ny, nx), tau0, dtype=np.float64)
nut = np.zeros((ny, nx), dtype=np.float64)
rho_ref = np.maximum(rho_safe, 1.0e-12)
cs2 = 1.0 / 3.0
nut_cap = max(0.0, LES_POST_NUT_MAX_RATIO * nu0)
cs = float(cfg.get("les_cs", 0.16))
for _ in range(LES_POST_FP_ITERS):
denom = 2.0 * rho_ref * cs2 * np.maximum(tau_eff, 0.500001)
sxx = -pixx / denom
syy = -piyy / denom
sxy = -pixy / denom
tr = 0.5 * (sxx + syy)
sxx_dev = sxx - tr
syy_dev = syy - tr
sxy_dev = sxy
s_mag = np.sqrt(np.maximum(0.0, 2.0 * (sxx_dev * sxx_dev + syy_dev * syy_dev + 2.0 * sxy_dev * sxy_dev)))
nut = np.clip((cs * cs) * s_mag, 0.0, nut_cap)
tau_new = np.clip(0.5 + 3.0 * (nu0 + nut), 0.500001, tau_max)
tau_eff = LES_POST_FP_RELAX * tau_new + (1.0 - LES_POST_FP_RELAX) * tau_eff
omega_plus = np.clip(1.0 / tau_eff, omega_min, omega_max)
denom_odd = np.maximum(1.0 / omega_plus - 0.5, 1.0e-9)
lam = float(cfg.get("trt_magic_param", 0.1875))
omega_minus = 1.0 / (lam / denom_odd + 0.5)
fluid = cfg["flag"].reshape(ny, nx) == FLUID_FLAG
rho_mean = float(np.mean(rho[fluid])) if np.any(fluid) else float(np.mean(rho))
rho_prime = rho - rho_mean
return {
"fluid": fluid,
"omega_plus": omega_plus,
"omega_minus": omega_minus,
"nut": nut,
"rho_prime": rho_prime,
}
def plot_trt_les_maps(cfg, host_ddf, out_dir):
if cfg["nq"] != 9 or cfg["collision_model"] != 1 or not cfg["use_les"]:
return None
fields = compute_trt_les_fields_2d(cfg, host_ddf)
if fields is None:
return None
fluid = fields["fluid"]
tag = make_case_tag(cfg)
out_path = os.path.join(out_dir, f"{tag}_trt_les_fields.png")
fig, axes = plt.subplots(2, 2, figsize=(12, 9))
panels = [
("omega_plus", r"$\omega_+$", "viridis", False),
("omega_minus", r"$\omega_-$", "viridis", False),
("nut", r"$\nu_t$", "magma", False),
("rho_prime", r"$\rho-\overline{\rho}$", "RdBu_r", True),
]
for ax, (key, title, cmap, signed) in zip(axes.ravel(), panels):
arr = fields[key]
arr_m = np.ma.array(arr, mask=~fluid)
finite_vals = arr[fluid]
if finite_vals.size == 0:
vmin, vmax = 0.0, 1.0
elif signed:
span = np.percentile(np.abs(finite_vals), 99)
span = max(float(span), 1.0e-12)
vmin, vmax = -span, span
else:
q1 = float(np.percentile(finite_vals, 1))
q99 = float(np.percentile(finite_vals, 99))
if not np.isfinite(q1) or not np.isfinite(q99) or q99 <= q1:
q1 = float(np.min(finite_vals))
q99 = float(np.max(finite_vals))
if q99 <= q1:
q99 = q1 + 1.0e-12
vmin, vmax = q1, q99
im = ax.imshow(arr_m, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax)
plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
ax.set_title(title)
ax.set_xlabel("x")
ax.set_ylabel("y")
fig.suptitle(f"{tag} TRT-LES diagnostics")
fig.tight_layout()
fig.savefig(out_path, dpi=160)
plt.close(fig)
return out_path
def plot_case(cfg, host_ddf, out_dir):
nq = cfg["nq"]
nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
@ -176,7 +544,7 @@ def compute_vis_omega(reynolds, diameter, u0):
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):
omega_collision_max, inlet_profile, trt_magic_param):
lines = compiler.read_lines(compiler.kernel_path("macros.h"))
defs = {
"MULT_GPU": "False",
@ -203,11 +571,12 @@ def set_macros(nx, ny, nz, dim, nq, vis, u0, collision_model, use_les, les_cs,
"USE_DDF_SHIFTING": 0,
"USE_LES": int(use_les),
"LES_CS": f"{les_cs:.6f}f",
"INLET_PROFILE": 0,
"INLET_PROFILE": int(inlet_profile),
"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",
"TRT_MAGIC_PARAM": f"{float(trt_magic_param):.6f}f",
}
for name, value in defs.items():
lines = compiler.modify_macro(lines, name, value)
@ -261,6 +630,8 @@ def run_case(device_id, cfg):
outlet_backflow_clamp=cfg["outlet_backflow_clamp"],
outlet_blend_alpha=cfg["outlet_blend_alpha"],
omega_collision_max=cfg["omega_collision_max"],
inlet_profile=cfg["inlet_profile"],
trt_magic_param=cfg["trt_magic_param"],
)
compiler.compile_kernel_v2()
@ -310,10 +681,18 @@ def run_case(device_id, cfg):
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"])
host0 = np.empty(n * nq, dtype=np.float32)
cuda.memcpy_dtoh(host0, d_fi)
host0 = impose_rest_state_on_nonfluid(cfg, host0)
cuda.memcpy_htod(d_fi, host0)
cuda.memcpy_htod(d_fi2, host0)
diag0 = compute_case_diagnostics(cfg, host0)
mass0 = diag0["mass"]
inlet_var0 = diag0["inlet_var"]
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)
@ -347,10 +726,26 @@ def run_case(device_id, cfg):
center = float(rho[nz // 2, ny // 2, nx // 2])
ok, reason = validate_case(rho)
diag_end = compute_case_diagnostics(cfg, host)
if np.isfinite(mass0) and mass0 != 0.0 and np.isfinite(diag_end["mass"]):
mass_drift = abs(diag_end["mass"] - mass0) / abs(mass0)
else:
mass_drift = float("nan")
if np.isfinite(inlet_var0) and inlet_var0 > 0.0 and np.isfinite(diag_end["inlet_var"]):
inlet_var_ratio_to_init = diag_end["inlet_var"] / inlet_var0
else:
inlet_var_ratio_to_init = float("nan")
plot_path = None
if cfg.get("save_plot", True):
plot_path = plot_case(cfg, host, cfg["out_dir"])
trt_les_map_path = None
if cfg.get("save_trt_les_maps", True):
trt_les_map_path = plot_trt_les_maps(cfg, host, cfg["out_dir"])
return {
"case_tag": make_case_tag(cfg),
"name": cfg["name"],
@ -361,6 +756,18 @@ def run_case(device_id, cfg):
"rho_center": center,
"rho_min": float(np.nanmin(rho)),
"rho_max": float(np.nanmax(rho)),
"mass0": float(mass0),
"mass_end": float(diag_end["mass"]),
"mass_drift": float(mass_drift),
"inlet_var0": float(inlet_var0),
"inlet_var_end": float(diag_end["inlet_var"]),
"inlet_var_ratio_to_init": float(inlet_var_ratio_to_init),
"inlet_line_rel_l2": float(diag_end.get("inlet_line_rel_l2", float("nan"))),
"inlet_line_rel_linf": float(diag_end.get("inlet_line_rel_linf", float("nan"))),
"inlet_wave_ux_rel": float(diag_end.get("inlet_wave_ux_rel", float("nan"))),
"inlet_wave_rho_rel": float(diag_end.get("inlet_wave_rho_rel", float("nan"))),
"wake_checker_rel": float(diag_end.get("wake_checker_rel", float("nan"))),
"wake_checker_anti_corr_x": float(diag_end.get("wake_checker_anti_corr_x", float("nan"))),
"omega": cfg["omega"],
"vis": cfg["vis"],
"collision_model": cfg["collision_model"],
@ -369,10 +776,13 @@ def run_case(device_id, cfg):
"outlet_mode": int(cfg["outlet_mode"]),
"outlet_backflow_clamp": int(cfg["outlet_backflow_clamp"]),
"outlet_blend_alpha": float(cfg["outlet_blend_alpha"]),
"inlet_profile": int(cfg["inlet_profile"]),
"omega_collision_max": float(cfg["omega_collision_max"]),
"trt_magic_param": float(cfg["trt_magic_param"]),
"pass": bool(ok),
"reason": reason,
"plot_path": plot_path,
"trt_les_map_path": trt_les_map_path,
}
finally:
ctx.pop()
@ -380,7 +790,7 @@ def run_case(device_id, cfg):
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):
omega_collision_max, inlet_profile, trt_magic_param):
nx, ny, nz = 512, 256, 1
cx, cy, radius = 128.0, 128.0, 24.0
u0 = 0.03
@ -392,6 +802,9 @@ def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
"nx": nx,
"ny": ny,
"nz": nz,
"cx": cx,
"cy": cy,
"radius": radius,
"flag": build_flags_2d(nx, ny, cx, cy, radius),
"u0": u0,
"vis": vis,
@ -404,15 +817,18 @@ def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
"outlet_mode": int(outlet_mode),
"outlet_backflow_clamp": int(outlet_backflow_clamp),
"outlet_blend_alpha": float(outlet_blend_alpha),
"inlet_profile": int(inlet_profile),
"omega_collision_max": float(omega_collision_max),
"trt_magic_param": float(trt_magic_param),
"target_re": re2d,
"save_plot": True,
"save_trt_les_maps": 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):
omega_collision_max, inlet_profile, trt_magic_param):
nx, ny, nz = 256, 128, 32
cx, cy, radius = 64.0, 64.0, 12.0
u0 = 0.04
@ -424,6 +840,9 @@ def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
"nx": nx,
"ny": ny,
"nz": nz,
"cx": cx,
"cy": cy,
"radius": radius,
"flag": build_flags_3d(nx, ny, nz, cx, cy, radius),
"u0": u0,
"vis": vis,
@ -436,9 +855,12 @@ def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
"outlet_mode": int(outlet_mode),
"outlet_backflow_clamp": int(outlet_backflow_clamp),
"outlet_blend_alpha": float(outlet_blend_alpha),
"inlet_profile": int(inlet_profile),
"omega_collision_max": float(omega_collision_max),
"trt_magic_param": float(trt_magic_param),
"target_re": re3d,
"save_plot": True,
"save_trt_les_maps": True,
"out_dir": out_dir,
}
@ -454,14 +876,18 @@ def build_comprehensive_cases(args, 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))
omega_collision_max=args.omega_collision_max,
inlet_profile=args.inlet_profile,
trt_magic_param=args.trt_magic_param))
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))
omega_collision_max=args.omega_collision_max,
inlet_profile=args.inlet_profile,
trt_magic_param=args.trt_magic_param))
return cases
@ -480,10 +906,14 @@ def main():
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("--inlet-profile", type=int, default=1, choices=[0, 1],
help="0=uniform inlet, 1=parabolic inlet")
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("--trt-magic-param", type=float, default=0.002,
help="TRT magic parameter Lambda used in omega- mapping")
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")
@ -504,7 +934,8 @@ def main():
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)
args.outlet_blend_alpha, args.omega_collision_max,
args.inlet_profile, args.trt_magic_param)
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))
@ -512,7 +943,8 @@ def main():
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)
args.outlet_blend_alpha, args.omega_collision_max,
args.inlet_profile, args.trt_magic_param)
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))
@ -534,9 +966,15 @@ def main():
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"mass_drift={r['mass_drift']:.3e}, inlet_var_end={r['inlet_var_end']:.3e}, "
f"inlet_relL2={r['inlet_line_rel_l2']:.3e}, inlet_relLinf={r['inlet_line_rel_linf']:.3e}, "
f"waveUx={r['inlet_wave_ux_rel']:.3e}, waveRho={r['inlet_wave_rho_rel']:.3e}, "
f"chkRel={r['wake_checker_rel']:.3e}, chkAntiX={r['wake_checker_anti_corr_x']:.3e}, "
f"MLUPS={r['mlups']:.1f}, pass={r['pass']} ({r['reason']})")
if r.get("plot_path"):
print(f" plot: {r['plot_path']}")
if r.get("trt_les_map_path"):
print(f" trt_les_map: {r['trt_les_map_path']}")
print(f"Pass rate: {n_pass}/{len(results)}")
print(f"Saved: {out_json}")
finally: