diff --git a/README.md b/README.md index 44a6b72..cc8bed3 100644 --- a/README.md +++ b/README.md @@ -97,7 +97,7 @@ The on-disk schema matches `src/CelerisLab/configs/config_lbm.json` (nested sect "ddf_shifting": false, "les": { "enabled": false, "cs": 0.16, "closed_form": true }, "trt": { "magic_param": 0.1875 }, - "inlet": { "profile": "parabolic", "trt_neq_damp": 0.5 }, + "inlet": { "profile": "parabolic", "scheme": "zou_he_local", "trt_neq_damp": 0.5 }, "outlet": { "mode": "neq_extrap", "backflow_clamp": true, diff --git a/src/CelerisLab/common/checkpoint.py b/src/CelerisLab/common/checkpoint.py index 21ab15d..9b4d584 100644 --- a/src/CelerisLab/common/checkpoint.py +++ b/src/CelerisLab/common/checkpoint.py @@ -85,6 +85,7 @@ def save_checkpoint(field, stepper, lbm_cfg, bodies, path=None): "les_closed_form": cfg.les_closed_form, "trt_magic_param": cfg.trt_magic_param, "inlet_profile": cfg.inlet_profile, + "inlet_scheme": cfg.inlet_scheme, "outlet_mode": cfg.outlet_mode, "omega_min": cfg.omega_min, "omega_max": cfg.omega_max, } diff --git a/src/CelerisLab/config.py b/src/CelerisLab/config.py index 42ac837..4edb505 100644 --- a/src/CelerisLab/config.py +++ b/src/CelerisLab/config.py @@ -26,6 +26,12 @@ COLLISION_MAP = {"SRT": 0, "TRT": 1, "MRT": 2} STREAMING_MAP = {"double_buffer": 0, "esopull": 1} PRECISION_MAP = {"FP32": 0, "FP16S": 1, "FP16C": 2} INLET_MAP = {"uniform": 0, "parabolic": 1} +INLET_SCHEME_MAP = { + "zou_he_local": 0, + "channel_stabilized": 1, + "equilibrium": 2, + "regularized": 3, +} OUTLET_MAP = {"neq_extrap": 0, "zero_gradient": 1, "blended": 2} Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1} DTYPE_MAP = {"FP32": "float", "FP64": "double"} @@ -66,10 +72,12 @@ class LBMConfig: les_closed_form: bool = True trt_magic_param: float = 0.1875 inlet_profile: str = "parabolic" + inlet_scheme: str = "zou_he_local" outlet_mode: str = "neq_extrap" outlet_blend_alpha: float = 0.7 outlet_backflow_clamp: bool = True inlet_trt_neq_damp: float = 0.5 + inlet_regularized_neq_damp: float = 0.5 outlet_srt_neq_damp: float = 0.5 y_wall_bc: str = "bounce_back" omega_min: float = 0.01 @@ -106,6 +114,8 @@ class LBMConfig: assert self.data_type in DTYPE_MAP, f"Unknown data_type: {self.data_type}" assert self.store_precision in PRECISION_MAP, ( f"Unknown store_precision: {self.store_precision}") + assert self.inlet_scheme in INLET_SCHEME_MAP, ( + f"Unknown inlet_scheme: {self.inlet_scheme}") if self.store_precision == "FP16C": raise ValueError( "store_precision='FP16C' is not supported in the current runtime path. " @@ -119,6 +129,12 @@ class LBMConfig: assert self.nx > 0 and self.ny > 0 and self.nz > 0 if self.dim == 2: assert self.nz == 1, "nz must be 1 for 2D" + if not (0.0 <= self.inlet_trt_neq_damp <= 1.0): + raise ValueError("method.inlet.trt_neq_damp must lie in [0, 1].") + if not (0.0 <= self.inlet_regularized_neq_damp <= 1.0): + raise ValueError("method.inlet.regularized_neq_damp must lie in [0, 1].") + if not (0.0 <= self.outlet_srt_neq_damp <= 1.0): + raise ValueError("method.outlet.srt_neq_damp must lie in [0, 1].") if self.dim == 3 and self.y_wall_bc == "free_slip": raise ValueError("y_wall_bc='free_slip' is currently supported for D2Q9 only.") if self.omega_max >= 2.0: @@ -156,6 +172,7 @@ class LBMConfig: "LES_CS": f"{self.les_cs:.6f}f", "LES_CLOSED_FORM": int(self.les_closed_form), "INLET_PROFILE": INLET_MAP[self.inlet_profile], + "INLET_SCHEME": INLET_SCHEME_MAP[self.inlet_scheme], "OUTLET_MODE": OUTLET_MAP[self.outlet_mode], "OUTLET_BLEND_ALPHA": f"{self.outlet_blend_alpha:.3f}f", "OUTLET_BACKFLOW_CLAMP": int(self.outlet_backflow_clamp), @@ -164,6 +181,7 @@ class LBMConfig: "OMEGA_COLLISION_MAX": f"{self.omega_max:.3f}f", "TRT_MAGIC_PARAM": f"{self.trt_magic_param:.6f}f", "INLET_TRT_NEQ_DAMP": f"{self.inlet_trt_neq_damp:.4f}f", + "INLET_REGULARIZED_NEQ_DAMP": f"{self.inlet_regularized_neq_damp:.4f}f", "OUTLET_SRT_NEQ_DAMP": f"{self.outlet_srt_neq_damp:.4f}f", } @@ -216,10 +234,12 @@ def load_lbm_config(path: Optional[str] = None) -> LBMConfig: les_closed_form=m["les"].get("closed_form", True), trt_magic_param=m["trt"]["magic_param"], inlet_profile=m["inlet"]["profile"], + inlet_scheme=m["inlet"].get("scheme", "zou_he_local"), outlet_mode=m["outlet"]["mode"], outlet_blend_alpha=m["outlet"]["blend_alpha"], outlet_backflow_clamp=m["outlet"]["backflow_clamp"], inlet_trt_neq_damp=m["inlet"].get("trt_neq_damp", 0.5), + inlet_regularized_neq_damp=m["inlet"].get("regularized_neq_damp", 0.5), outlet_srt_neq_damp=m["outlet"].get("srt_neq_damp", 0.5), y_wall_bc=m.get("y_wall_bc", "bounce_back"), omega_min=m["omega_guard"]["min"], @@ -235,4 +255,4 @@ def load_body_config(path: Optional[str] = None) -> BodyConfig: fpath = _find_config("config_body.json", path) with open(fpath) as f: d = json.load(f) - return BodyConfig(objects=d.get("objects", [])) + return BodyConfig(objects=d.get("objects", [])) \ No newline at end of file diff --git a/src/CelerisLab/configs/CONFIG.md b/src/CelerisLab/configs/CONFIG.md index 5cba1db..6c87f32 100644 --- a/src/CelerisLab/configs/CONFIG.md +++ b/src/CelerisLab/configs/CONFIG.md @@ -35,8 +35,10 @@ Python `config.py` 只负责读取和校验，不是配置位置。 | `les.cs` | float | 0.16 | | Smagorinsky 常数 | | `les.closed_form` | bool | true | | 闭合形式 τ_eff（vs 迭代） | | `trt.magic_param` | float | 0.1875 | | TRT Λ 参数，高 Re 建议 0.001 | -| `inlet.profile` | string | `"parabolic"` | `uniform`, `parabolic` | 入口速度剖面 | -| `inlet.trt_neq_damp` | float | 0.5 | [0, 1] | TRT 入口 NEQ donor 阻尼；更小更平滑、精度略降 | +| `inlet.profile` | string | `"parabolic"` | `uniform`, `parabolic` | 入口速度剖面（物理目标速度，与 scheme 独立） | +| `inlet.scheme` | string | `"zou_he_local"` | `zou_he_local`, `channel_stabilized`, `equilibrium`, `regularized` | 入口数值闭合。`zou_he_local` 为本地 Zou-He，适合研究或 MRT 路径；`channel_stabilized` 为 donor NEQ 稳定化入口，适合高阻塞或更保守的量产路径；`equilibrium` 直接写入 `feq` 源态，适合 ghost-source 架构下的稳健 SRT 基线；`regularized` 使用本地宏量加 incoming donor NEQ 阻尼，是介于 `equilibrium` 与 `channel_stabilized` 之间的实验入口 | +| `inlet.trt_neq_damp` | float | 0.5 | [0, 1] | 仅 `channel_stabilized`：TRT 入口 donor NEQ 阻尼；更小更平滑、精度略降 | +| `inlet.regularized_neq_damp` | float | 0.5 | [0, 1] | 仅 `regularized`：incoming 方向 donor NEQ 阻尼；0 退化到 unknown 方向仅平衡态，1 为 unknown 方向全 donor NEQ | | `outlet.mode` | string | `"neq_extrap"` | `neq_extrap`, `zero_gradient`, `blended` | 出口条件 | | `outlet.backflow_clamp` | bool | true | | 出口回流钳位 | | `outlet.blend_alpha` | float | 0.7 | | `blended` 下对未知入域方向的混合系数（**所有碰撞模型**共用同一路径） | diff --git a/src/CelerisLab/configs/config_lbm.json b/src/CelerisLab/configs/config_lbm.json index 32b6b8a..0c12fa8 100644 --- a/src/CelerisLab/configs/config_lbm.json +++ b/src/CelerisLab/configs/config_lbm.json @@ -33,9 +33,13 @@ }, "inlet": { "profile": "parabolic", + "scheme": "zou_he_local", "_profile": "uniform | parabolic", + "_scheme": "zou_he_local | channel_stabilized | equilibrium | regularized", "trt_neq_damp": 0.5, - "_trt_neq_damp": "TRT inlet NEQ donor damping [0,1]. Lower = smoother inlet, less accurate." + "_trt_neq_damp": "channel_stabilized only. TRT donor NEQ damping [0,1]. Lower = smoother inlet, less accurate.", + "regularized_neq_damp": 0.5, + "_regularized_neq_damp": "regularized only. Incoming-direction donor NEQ damping [0,1]. 0 = equilibrium-only on unknowns, 1 = full donor NEQ on unknowns." }, "outlet": { "mode": "neq_extrap", @@ -58,4 +62,4 @@ "compute_capability": "auto", "_compute_capability": "auto | e.g. 7.0" } -} +} \ No newline at end of file diff --git a/src/CelerisLab/cuda/compiler_v2.py b/src/CelerisLab/cuda/compiler_v2.py index d595ff0..2624ddb 100644 --- a/src/CelerisLab/cuda/compiler_v2.py +++ b/src/CelerisLab/cuda/compiler_v2.py @@ -115,6 +115,7 @@ def generate_config(cfg: LBMConfig, n_objects: int = 0): #define LES_CLOSED_FORM {m['LES_CLOSED_FORM']} #define INLET_PROFILE {m['INLET_PROFILE']} +#define INLET_SCHEME {m['INLET_SCHEME']} #define OUTLET_MODE {m['OUTLET_MODE']} #define OUTLET_BLEND_ALPHA {m['OUTLET_BLEND_ALPHA']} #define OUTLET_BACKFLOW_CLAMP {m['OUTLET_BACKFLOW_CLAMP']} @@ -125,10 +126,12 @@ def generate_config(cfg: LBMConfig, n_objects: int = 0): #define TRT_MAGIC_PARAM {m['TRT_MAGIC_PARAM']} // NEQ damping coefficients for inlet/outlet BC reconstruction. -// TRT inlet: damps donor non-equilibrium to reduce inlet noise at high Re. -// SRT outlet: damps donor non-equilibrium to suppress checkerboard noise. -#define INLET_TRT_NEQ_DAMP {m['INLET_TRT_NEQ_DAMP']} -#define OUTLET_SRT_NEQ_DAMP {m['OUTLET_SRT_NEQ_DAMP']} +// TRT inlet: donor damping used by the channel_stabilized inlet. +// Regularized inlet: damping on incoming-direction donor NEQ. +// SRT outlet: damped outlet reconstruction to suppress checkerboard noise. +#define INLET_TRT_NEQ_DAMP {m['INLET_TRT_NEQ_DAMP']} +#define INLET_REGULARIZED_NEQ_DAMP {m['INLET_REGULARIZED_NEQ_DAMP']} +#define OUTLET_SRT_NEQ_DAMP {m['OUTLET_SRT_NEQ_DAMP']} #endif """) @@ -219,4 +222,3 @@ def load_module(ptx_path: Optional[str] = None) -> cuda.Module: def compile_kernel_v2(arch: str = "sm_70") -> str: """Alias for compile_kernel() kept for test-script compatibility.""" return compile_kernel(arch=arch) - diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/channel_stabilized.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/channel_stabilized.cuh new file mode 100644 index 0000000..13a964e --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/channel_stabilized.cuh @@ -0,0 +1,110 @@ +// CelerisLab – boundary/inlet/channel_stabilized.cuh +// Donor-based west velocity inlet closures designed for robust ghost-source use. +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_INLET_CHANNEL_STABILIZED_CUH +#define CELERIS_BOUNDARY_INLET_CHANNEL_STABILIZED_CUH + +#if LATTICE_MODEL == LATTICE_D2Q9 +__device__ inline void apply_channel_stabilized_inlet_d2q9( + float* __restrict__ f, + const float* __restrict__ f_neb, + float y_coord) +{ + float rho_neb, u_neb, v_neb; + compute_rho_u(f_neb, rho_neb, u_neb, v_neb); + + const float u_target = inlet_target_u(y_coord); + const float v_target = 0.0f; + const float rho_in = west_velocity_rho_closure_d2q9(f, u_target); + + float feq_tar[9], feq_neb[9]; + compute_feq(rho_in, u_target, v_target, feq_tar); + compute_feq(rho_neb, u_neb, v_neb, feq_neb); + +#if COLLISION_MODEL == 1 + const float beta_n = INLET_TRT_NEQ_DAMP; +#else + const float beta_n = 1.0f; +#endif + + const float f1_try = feq_tar[1] + beta_n * (f_neb[1] - feq_neb[1]); + const float known_sum = f[0] + f[2] + f[3] + f[4] + f[6] + f[8]; + float pair_diff = rho_in * v_target + (f[4] - f[3]) + (f[6] - f[8]); + const float f1_hi = fmaxf(0.0f, rho_in - known_sum - fabsf(pair_diff)); + const float f1 = fminf(fmaxf(f1_try, 0.0f), f1_hi); + float pair_sum = rho_in - known_sum - f1; + + if (fabsf(pair_diff) > pair_sum) { + pair_diff = copysignf(pair_sum, pair_diff); + } + + f[1] = f1; + f[5] = 0.5f * (pair_sum + pair_diff); + f[7] = 0.5f * (pair_sum - pair_diff); +} +#endif + +#if LATTICE_MODEL == LATTICE_D3Q19 +__device__ inline void apply_channel_stabilized_inlet_d3q19( + float* __restrict__ f, + const float* __restrict__ f_neb, + float y_coord) +{ + float rho_neb, un, vn, wn; + compute_rho_u(f_neb, rho_neb, un, vn, wn); + + const float u_tar = inlet_target_u(y_coord); + const float v_tar = 0.0f; + const float w_tar = 0.0f; + const float rho_in = west_velocity_rho_closure_d3q19(f, u_tar); + + float feq_tar[19], feq_neb[19]; + compute_feq(rho_in, u_tar, v_tar, w_tar, feq_tar); + compute_feq(rho_neb, un, vn, wn, feq_neb); + +#if COLLISION_MODEL == 1 + const float beta_n = INLET_TRT_NEQ_DAMP; +#else + const float beta_n = 1.0f; +#endif + + const float f1_try = feq_tar[1] + beta_n * (f_neb[1] - feq_neb[1]); + const float zsum_try = (feq_tar[9] + feq_tar[15]) + + beta_n * ((f_neb[9] - feq_neb[9]) + + (f_neb[15] - feq_neb[15])); + + const float known_sum = f[0] + f[2] + f[3] + f[4] + f[5] + f[6] + + f[8] + f[10] + f[11] + f[12] + f[14] + + f[16] + f[17] + f[18]; + const float rem_total = rho_in - known_sum; + + float y_diff = rho_in * v_tar + - (f[3] - f[4] - f[8] + f[11] - f[12] + f[14] + f[17] - f[18]); + float z_diff = rho_in * w_tar + - (f[5] - f[6] - f[10] + f[11] - f[12] + f[16] + f[18] - f[17]); + + const float f1_hi = fmaxf(0.0f, rem_total - fabsf(y_diff)); + const float f1 = fminf(fmaxf(f1_try, 0.0f), f1_hi); + + const float zsum_hi = fmaxf(0.0f, rem_total - f1 - fabsf(y_diff)); + const float z_sum = fminf(fmaxf(zsum_try, 0.0f), zsum_hi); + if (fabsf(z_diff) > z_sum) { + z_diff = copysignf(z_sum, z_diff); + } + + float y_sum = rem_total - f1 - z_sum; + y_sum = fmaxf(y_sum, 0.0f); + if (fabsf(y_diff) > y_sum) { + y_diff = copysignf(y_sum, y_diff); + } + + f[1] = f1; + f[9] = 0.5f * (z_sum + z_diff); + f[15] = 0.5f * (z_sum - z_diff); + f[7] = 0.5f * (y_sum + y_diff); + f[13] = 0.5f * (y_sum - y_diff); +} +#endif + +#endif // CELERIS_BOUNDARY_INLET_CHANNEL_STABILIZED_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh new file mode 100644 index 0000000..5f53658 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh @@ -0,0 +1,72 @@ +// CelerisLab – boundary/inlet/common.cuh +// Shared helpers for west velocity inlet source-state generation. +// +// Important semantic contract in this solver: +// - x=0 inlet cells are SOLID|BC_INLET ghost source nodes +// - they generate the state later pulled by the first interior fluid column +// - inlet methods therefore construct a source state, not a textbook fluid-node +// boundary update in isolation +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_INLET_COMMON_CUH +#define CELERIS_BOUNDARY_INLET_COMMON_CUH + +#ifndef INLET_PROFILE +#define INLET_PROFILE 1 +#endif + +#ifndef INLET_SCHEME +#define INLET_SCHEME 0 +#endif + +#ifndef INLET_TRT_NEQ_DAMP +#define INLET_TRT_NEQ_DAMP 0.50f +#endif + +#ifndef INLET_REGULARIZED_NEQ_DAMP +#define INLET_REGULARIZED_NEQ_DAMP 0.50f +#endif + +__device__ __forceinline__ float inlet_target_u(float y_coord) { +#if INLET_PROFILE == 0 + return U0; +#else + const float y_clamped = fminf((float)(NY - 2), fmaxf(1.0f, y_coord)); + const float H = fmaxf((float)(NY - 2), 1.0f); + const float eta = (y_clamped - 0.5f) / H; + const float shape = fmaxf(0.0f, 4.0f * eta * (1.0f - eta)); + return U0 * 1.5f * shape; +#endif +} + +#if LATTICE_MODEL == LATTICE_D2Q9 +__device__ __forceinline__ float west_velocity_rho_closure_d2q9( + const float* __restrict__ f, + float ux_target) +{ + return (f[0] + f[3] + f[4] + 2.0f * (f[2] + f[6] + f[8])) + / (1.0f - ux_target); +} +#endif + +#if LATTICE_MODEL == LATTICE_D3Q19 +__device__ __forceinline__ float west_velocity_rho_closure_d3q19( + const float* __restrict__ f, + float ux_target) +{ + return (f[0] + f[3] + f[4] + f[5] + f[6] + f[11] + f[12] + f[17] + f[18] + + 2.0f * (f[2] + f[8] + f[10] + f[14] + f[16])) + / (1.0f - ux_target); +} +#endif + +__device__ __forceinline__ bool inlet_scheme_uses_post_collision_ghost() +{ +#if INLET_SCHEME == 0 + return true; +#else + return false; +#endif +} + +#endif // CELERIS_BOUNDARY_INLET_COMMON_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/equilibrium.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/equilibrium.cuh new file mode 100644 index 0000000..c95ae78 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/equilibrium.cuh @@ -0,0 +1,48 @@ +// CelerisLab – boundary/inlet/equilibrium.cuh +// West velocity inlet source state built from full equilibrium. +// +// This method recovers rho from local west-boundary mass closure and then +// overwrites the full ghost-node state with feq(rho, u_target). +// It is robust for ghost-source architectures because it injects no boundary +// non-equilibrium content. +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_INLET_EQUILIBRIUM_CUH +#define CELERIS_BOUNDARY_INLET_EQUILIBRIUM_CUH + +#if LATTICE_MODEL == LATTICE_D2Q9 +__device__ inline void apply_equilibrium_left_velocity_inlet_d2q9( + float* __restrict__ f, + float ux_target, + float uy_target) +{ + const float rho = west_velocity_rho_closure_d2q9(f, ux_target); + float feq[9]; + compute_feq(rho, ux_target, uy_target, feq); + + #pragma unroll + for (int i = 0; i < 9; i++) { + f[i] = feq[i]; + } +} +#endif + +#if LATTICE_MODEL == LATTICE_D3Q19 +__device__ inline void apply_equilibrium_left_velocity_inlet_d3q19( + float* __restrict__ f, + float ux_target, + float uy_target, + float uz_target) +{ + const float rho = west_velocity_rho_closure_d3q19(f, ux_target); + float feq[19]; + compute_feq(rho, ux_target, uy_target, uz_target, feq); + + #pragma unroll + for (int i = 0; i < 19; i++) { + f[i] = feq[i]; + } +} +#endif + +#endif // CELERIS_BOUNDARY_INLET_EQUILIBRIUM_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/regularized.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/regularized.cuh new file mode 100644 index 0000000..7775c96 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/regularized.cuh @@ -0,0 +1,59 @@ +// CelerisLab – boundary/inlet/regularized.cuh +// West velocity inlet with local macro state and donor-damped incoming NEQ. +// +// This method keeps the target macro state from local west-boundary closure but +// avoids a full local algebraic source state by injecting only damped donor NEQ +// on incoming directions. +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_INLET_REGULARIZED_CUH +#define CELERIS_BOUNDARY_INLET_REGULARIZED_CUH + +#if LATTICE_MODEL == LATTICE_D2Q9 +__device__ inline void apply_regularized_left_velocity_inlet_d2q9( + float* __restrict__ f, + const float* __restrict__ f_neb, + float ux_target, + float uy_target) +{ + const float rho = west_velocity_rho_closure_d2q9(f, ux_target); + float rho_neb, u_neb, v_neb; + compute_rho_u(f_neb, rho_neb, u_neb, v_neb); + + float feq_tar[9], feq_neb[9]; + compute_feq(rho, ux_target, uy_target, feq_tar); + compute_feq(rho_neb, u_neb, v_neb, feq_neb); + + const float beta = INLET_REGULARIZED_NEQ_DAMP; + 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 + +#if LATTICE_MODEL == LATTICE_D3Q19 +__device__ inline void apply_regularized_left_velocity_inlet_d3q19( + float* __restrict__ f, + const float* __restrict__ f_neb, + float ux_target, + float uy_target, + float uz_target) +{ + const float rho = west_velocity_rho_closure_d3q19(f, ux_target); + float rho_neb, u_neb, v_neb, w_neb; + compute_rho_u(f_neb, rho_neb, u_neb, v_neb, w_neb); + + float feq_tar[19], feq_neb[19]; + compute_feq(rho, ux_target, uy_target, uz_target, feq_tar); + compute_feq(rho_neb, u_neb, v_neb, w_neb, feq_neb); + + const float beta = INLET_REGULARIZED_NEQ_DAMP; + f[1] = feq_tar[1] + beta * (f_neb[1] - feq_neb[1]); + f[7] = feq_tar[7] + beta * (f_neb[7] - feq_neb[7]); + f[9] = feq_tar[9] + beta * (f_neb[9] - feq_neb[9]); + f[13] = feq_tar[13] + beta * (f_neb[13] - feq_neb[13]); + f[15] = feq_tar[15] + beta * (f_neb[15] - feq_neb[15]); +} +#endif + +#endif // CELERIS_BOUNDARY_INLET_REGULARIZED_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh new file mode 100644 index 0000000..36c38b4 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/zou_he_local.cuh @@ -0,0 +1,90 @@ +// CelerisLab – boundary/inlet/zou_he_local.cuh +// Local on-site west velocity inlet closures. +// +// D2Q9 and D3Q19 variants recover rho from local mass closure and reconstruct +// unknown incoming populations from local target velocity constraints. +// +// In the current ghost-source architecture these reconstructed states are later +// used as pull sources. Therefore this method requires post-BC ghost collision +// in the step kernel. +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_INLET_ZOU_HE_LOCAL_CUH +#define CELERIS_BOUNDARY_INLET_ZOU_HE_LOCAL_CUH + +#if LATTICE_MODEL == LATTICE_D2Q9 + +// Free-slip y-walls: at inlet rows y=1 and y=NY-2, pull can source wall nodes for +// some known directions. Copy those from stored DDF at (x=1, same y) only. +__device__ inline void repair_zou_he_west_knowns_d2q9( + float* __restrict__ f, + const fpxx* __restrict__ fi_in, + unsigned int x, + unsigned int y) +{ + if (x != 0u) return; + + const unsigned long k_int = linear_index(x + 1u, y); + + if (y == 1u) { + f[3] = load_ddf(fi_in, index_f(k_int, 3u)); + f[8] = load_ddf(fi_in, index_f(k_int, 8u)); + } else if (y == (unsigned int)(NY - 2)) { + f[4] = load_ddf(fi_in, index_f(k_int, 4u)); + f[6] = load_ddf(fi_in, index_f(k_int, 6u)); + } +} + +__device__ inline void apply_zou_he_left_velocity_inlet_d2q9( + float* __restrict__ f, + float ux_target, + float uy_target) +{ + const float rho = west_velocity_rho_closure_d2q9(f, ux_target); + + f[1] = f[2] + (2.0f / 3.0f) * rho * ux_target; + + f[5] = f[6] + + 0.5f * (f[4] - f[3]) + + (1.0f / 6.0f) * rho * ux_target + + 0.5f * rho * uy_target; + + f[7] = f[8] + + 0.5f * (f[3] - f[4]) + + (1.0f / 6.0f) * rho * ux_target + - 0.5f * rho * uy_target; +} + +#endif // LATTICE_MODEL == LATTICE_D2Q9 + +#if LATTICE_MODEL == LATTICE_D3Q19 + +// Hecht-Harting style D3Q19 on-site west velocity closure adapted to the +// codebase paired ordering: +// 7 = (+x,+y), 8 = (-x,-y) +// 13 = (+x,-y), 14 = (-x,+y) +// 9 = (+x,+z), 10 = (-x,-z) +// 15 = (+x,-z), 16 = (-x,+z) +__device__ inline void apply_zou_he_left_velocity_inlet_d3q19( + float* __restrict__ f, + float ux_target, + float uy_target, + float uz_target) +{ + const float rho = west_velocity_rho_closure_d3q19(f, ux_target); + + const float Nyx = 0.5f * (f[3] + f[11] + f[17] - (f[4] + f[12] + f[18])) + - (1.0f / 3.0f) * rho * uy_target; + const float Nzx = 0.5f * (f[5] + f[11] + f[18] - (f[6] + f[17] + f[12])) + - (1.0f / 3.0f) * rho * uz_target; + + f[1] = f[2] + (1.0f / 3.0f) * rho * ux_target; + f[7] = f[8] + (1.0f / 6.0f) * rho * (ux_target + uy_target) - Nyx; + f[13] = f[14] + (1.0f / 6.0f) * rho * (ux_target - uy_target) + Nyx; + f[9] = f[10] + (1.0f / 6.0f) * rho * (ux_target + uz_target) - Nzx; + f[15] = f[16] + (1.0f / 6.0f) * rho * (ux_target - uz_target) + Nzx; +} + +#endif // LATTICE_MODEL == LATTICE_D3Q19 + +#endif // CELERIS_BOUNDARY_INLET_ZOU_HE_LOCAL_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh index 52ceea9..13d00d8 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh @@ -1,359 +1,225 @@ // CelerisLab – boundary/inlet_outlet.cuh -// Inlet and outlet boundary conditions (D2Q9). +// Inlet and outlet dispatch layer. // -// Parabolic inlet (non-equilibrium extrapolation, Zou-He style): -// Left wall (x=0): reconstruct cx>0 populations (i=1,5,7) +// This file contains only profile helpers, method includes, and compile-time +// switching. Each concrete inlet or outlet implementation lives in its own file +// under boundary/inlet/ or boundary/outlet/. // -// Pressure outlet (non-equilibrium extrapolation): -// Right wall (x=NX-1): reconstruct cx<0 populations (i=2,6,8) -// -// New paired D2Q9 ordering: -// cx = {0, 1,-1, 0, 0, 1,-1, 1,-1} -// cy = {0, 0, 0, 1,-1, 1,-1,-1, 1} +// Inlet scheme IDs: +// 0 = zou_he_local +// 1 = channel_stabilized +// 2 = equilibrium +// 3 = regularized // ============================================================================ #ifndef CELERIS_BOUNDARY_INLET_OUTLET_CUH #define CELERIS_BOUNDARY_INLET_OUTLET_CUH -#ifndef INLET_PROFILE -#define INLET_PROFILE 1 +#ifndef INLET_SCHEME +#define INLET_SCHEME 0 #endif -#ifndef OUTLET_MODE -#define OUTLET_MODE 0 -#endif +#include "inlet/common.cuh" +#include "inlet/zou_he_local.cuh" +#include "inlet/channel_stabilized.cuh" +#include "inlet/equilibrium.cuh" +#include "inlet/regularized.cuh" +#include "outlet/pressure_neq.cuh" -#ifndef OUTLET_BACKFLOW_CLAMP -#define OUTLET_BACKFLOW_CLAMP 1 -#endif - -#ifndef OUTLET_BLEND_ALPHA -#define OUTLET_BLEND_ALPHA 0.70f -#endif - -// OUTLET_SRT_NEQ_DAMP and INLET_TRT_NEQ_DAMP are injected by config_method.h. -// These fallback defaults are only active if building outside the normal -// Python compile pipeline (e.g. standalone nvcc tests). -#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 - // Uniform profile: U0 is the imposed streamwise velocity everywhere on the - // inlet fluid band. - return U0; -#else - // Parabolic profile on the fluid-node band y in [1, NY-2]. - // - // U0 is treated here as the mean streamwise inlet velocity. The returned - // peak centerline velocity is 1.5 * U0, matching the discrete Poiseuille - // profile used throughout initialization and boundary reconstruction. - // Keep this convention aligned with case setup and validation scripts. - const float y_clamped = fminf((float)(NY - 2), fmaxf(1.0f, y_coord)); - const float H = fmaxf((float)(NY - 2), 1.0f); - const float eta = (y_clamped - 0.5f) / H; // first and last fluid rows map near 0 and 1 - const float shape = fmaxf(0.0f, 4.0f * eta * (1.0f - eta)); - return U0 * 1.5f * shape; -#endif -} - -#if LATTICE_MODEL == LATTICE_D2Q9 - -// --------------------------------------------------------------------------- -// Parabolic inlet (x = 0, non-equilibrium extrapolation) -// -// f, f_neb are local DDF arrays: -// f = populations at the boundary node (x=0) -// f_neb = populations at the interior neighbor (x=1) -// y = y-coordinate of the boundary node -// -// Reconstructs f[1], f[5], f[7] (cx > 0 directions in new ordering). -// -// Velocity convention: -// uniform -> U0 is the imposed inlet velocity -// parabolic -> U0 is the mean inlet velocity, so inlet_target_u() returns a -// centerline peak of 1.5 * U0 -// -// Reconstruction keeps the Zou-He mass closure but does not copy all three -// unknown-direction NEQ parts from the donor. -// -// Why this split form: -// - In narrow high-blockage channels, the donor diagonals f_neb[5], f_neb[7] -// are strongly contaminated by near-wall shear and tend to inject a spurious -// negative shift into the first interior column. Empirically this is not just -// a corner-node artifact: the whole x=1 profile can be biased low when the -// channel becomes very narrow. -// - Removing NEQ entirely hurts stability, so the streamwise unknown f[1] keeps -// donor NEQ information. -// - The diagonal unknowns are instead reconstructed from local density and -// transverse-velocity constraints. A positivity limiter is applied only if -// the local constraints are mutually inconsistent, preferring exact rho and -// streamwise flux over exact v_target at that node. -// --------------------------------------------------------------------------- -__device__ inline void apply_parabolic_inlet(float* __restrict__ f, - const float* __restrict__ f_neb, - float y_coord) +#if DIM == 2 +__device__ __forceinline__ void apply_inlet_pull_d2q9( + float* __restrict__ f, + unsigned int x, + unsigned int y, + const fpxx* __restrict__ fi_in) { - // Donor macros from the first interior fluid column. - float rho_neb, u_neb, v_neb; - compute_rho_u(f_neb, rho_neb, u_neb, v_neb); - - // Target inlet velocity. - const float u_target = inlet_target_u(y_coord); + const float u_target = inlet_target_u((float)y); const float v_target = 0.0f; - // Zou-He mass closure at the west boundary. - // Known (after pull and any wall pre-repair): f[0],f[2],f[3],f[4],f[6],f[8] - // Unknown to reconstruct: f[1],f[5],f[7] - const float rho_in = (f[0] + f[3] + f[4] + 2.0f * (f[2] + f[6] + f[8])) - / (1.0f - u_target); - - float feq_tar[9], feq_neb[9]; - compute_feq(rho_in, u_target, v_target, feq_tar); - compute_feq(rho_neb, u_neb, v_neb, feq_neb); - -#if COLLISION_MODEL == 1 - const float beta_n = INLET_TRT_NEQ_DAMP; -#else - const float beta_n = 1.0f; +#if INLET_SCHEME == 0 +#if Y_WALL_BC == 1 + repair_zou_he_west_knowns_d2q9(f, fi_in, x, y); #endif - - // Keep donor NEQ only for the streamwise incoming population. This retains - // the stabilizing normal-flux information without feeding both diagonal - // donor modes back into the inlet every step. - const float f1_try = feq_tar[1] + beta_n * (f_neb[1] - feq_neb[1]); - - // Known-part density contribution. - const float known_sum = f[0] + f[2] + f[3] + f[4] + f[6] + f[8]; - - // From uy = (f3 - f4 + f5 - f6 - f7 + f8) / rho, so with v_target = 0: - // f5 - f7 = rho*v_target + (f4 - f3) + (f6 - f8) - float pair_diff = rho_in * v_target + (f[4] - f[3]) + (f[6] - f[8]); - - // Density fixes f5 + f7 once f1 is chosen: - // f5 + f7 = rho - known_sum - f1 - // To keep both diagonals non-negative we need pair_sum >= |pair_diff|, - // hence f1 <= rho - known_sum - |pair_diff|. - const float f1_hi = fmaxf(0.0f, rho_in - known_sum - fabsf(pair_diff)); - const float f1 = fminf(fmaxf(f1_try, 0.0f), f1_hi); - - float pair_sum = rho_in - known_sum - f1; - - // If the local constraints are still inconsistent because of roundoff or an - // extremely distorted incoming state, clip the transverse difference rather - // than emitting negative diagonal populations. - if (fabsf(pair_diff) > pair_sum) { - pair_diff = copysignf(pair_sum, pair_diff); + apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target); +#elif INLET_SCHEME == 1 || INLET_SCHEME == 3 + float f_neb[NQ]; + const 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)); } - - f[1] = f1; - f[5] = 0.5f * (pair_sum + pair_diff); - f[7] = 0.5f * (pair_sum - pair_diff); + #if INLET_SCHEME == 1 + apply_channel_stabilized_inlet_d2q9(f, f_neb, (float)y); + #else + apply_regularized_left_velocity_inlet_d2q9(f, f_neb, u_target, v_target); + #endif +#elif INLET_SCHEME == 2 + apply_equilibrium_left_velocity_inlet_d2q9(f, u_target, v_target); +#else + #error "Unsupported INLET_SCHEME for D2Q9" +#endif } -// --------------------------------------------------------------------------- -// Pressure outlet (x = NX-1, non-equilibrium extrapolation) -// -// Reconstructs f[2], f[6], f[8] (cx < 0 directions in new ordering) -// p_out = 0 (gauge pressure), uses velocity from neighbor. -// --------------------------------------------------------------------------- -__device__ inline void apply_pressure_outlet(float* __restrict__ f, - const float* __restrict__ f_neb, - float y_coord) +__device__ __forceinline__ void apply_outlet_pull_d2q9( + float* __restrict__ f, + unsigned int x, + unsigned int y, + const fpxx* __restrict__ fi_in) { - (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[6] = f_neb[6]; -#else - - // 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; - - 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 COLLISION_MODEL == 0 || COLLISION_MODEL == 1 - // SRT and TRT path: use full-population damped NEQ reconstruction at - // outlet to suppress checkerboard and boundary-source noise. - 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; + float f_neb[NQ]; + const 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)); } -#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]; - 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 + apply_pressure_outlet_d2q9(f, f_neb); } +#endif -#endif // LATTICE_D2Q9 - -// ============================================================================ -// D3Q19 inlet / outlet (non-equilibrium extrapolation) -// -// Parabolic inlet (x=0): reconstruct cx>0 populations i=1,7,9,13,15 -// Pressure outlet (x=NX-1): reconstruct cx<0 populations i=2,8,10,14,16 -// -// Uses generic feq computation from macro.cuh to avoid hand-expanded formulas. -// ============================================================================ -#if LATTICE_MODEL == LATTICE_D3Q19 - -__device__ inline void apply_parabolic_inlet_3d(float* __restrict__ f, - const float* __restrict__ f_neb, - float y_coord) +#if DIM == 3 +__device__ __forceinline__ void apply_inlet_pull_d3q19( + float* __restrict__ f, + unsigned int x, + unsigned int y, + unsigned int z, + const fpxx* __restrict__ fi_in) { - // Donor macros from the first interior fluid column. - float rho_neb, un, vn, wn; - compute_rho_u(f_neb, rho_neb, un, vn, wn); + const float u_target = inlet_target_u((float)y); + const float v_target = 0.0f; + const float w_target = 0.0f; - // Target velocity: parabolic in y, uniform in z. - const float u_tar = inlet_target_u(y_coord); - const float v_tar = 0.0f; - const float w_tar = 0.0f; - - // Zou-He mass balance at the west boundary. - // Unknown cx>0 populations are i = 1, 7, 9, 13, 15. - const float rho_in = (f[0] + f[3] + f[4] + f[5] + f[6] + f[11] + f[12] + f[17] + f[18] - + 2.0f * (f[2] + f[8] + f[10] + f[14] + f[16])) - / (1.0f - u_tar); - - float feq_tar[19], feq_neb[19]; - compute_feq(rho_in, u_tar, v_tar, w_tar, feq_tar); - compute_feq(rho_neb, un, vn, wn, feq_neb); - -#if COLLISION_MODEL == 1 - const float beta_n = INLET_TRT_NEQ_DAMP; +#if INLET_SCHEME == 0 + apply_zou_he_left_velocity_inlet_d3q19(f, u_target, v_target, w_target); +#elif INLET_SCHEME == 1 || INLET_SCHEME == 3 + float f_neb[NQ]; + const 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)); + } + #if INLET_SCHEME == 1 + apply_channel_stabilized_inlet_d3q19(f, f_neb, (float)y); + #else + apply_regularized_left_velocity_inlet_d3q19(f, f_neb, u_target, v_target, w_target); + #endif +#elif INLET_SCHEME == 2 + apply_equilibrium_left_velocity_inlet_d3q19(f, u_target, v_target, w_target); #else - const float beta_n = 1.0f; + #error "Unsupported INLET_SCHEME for D3Q19" #endif - - // D3Q19 counterpart of the D2Q9 split strategy: - // - keep donor NEQ on the pure streamwise incoming population f[1] - // - retain the x-z pair only through its total incoming mass, because the - // present channel setup has y-walls but no z-wall BC path - // - reconstruct the y-coupled diagonals from local rho/uy constraints to - // avoid feeding wall-shear contamination back into the inlet every step - const float f1_try = feq_tar[1] + beta_n * (f_neb[1] - feq_neb[1]); - const float zsum_try = (feq_tar[9] + feq_tar[15]) - + beta_n * ((f_neb[9] - feq_neb[9]) - + (f_neb[15] - feq_neb[15])); - - const float known_sum = f[0] + f[2] + f[3] + f[4] + f[5] + f[6] - + f[8] + f[10] + f[11] + f[12] + f[14] - + f[16] + f[17] + f[18]; - const float rem_total = rho_in - known_sum; - - // uy constraint: - // uy = (f3 - f4 + f7 - f8 + f11 - f12 + f14 - f13 + f17 - f18) / rho - float y_diff = rho_in * v_tar - - (f[3] - f[4] - f[8] + f[11] - f[12] + f[14] + f[17] - f[18]); - - // uz constraint: - // uz = (f5 - f6 + f9 - f10 + f11 - f12 + f16 - f15 + f18 - f17) / rho - float z_diff = rho_in * w_tar - - (f[5] - f[6] - f[10] + f[11] - f[12] + f[16] + f[18] - f[17]); - - // Reserve enough total mass for the y-diagonal pair to satisfy positivity. - const float f1_hi = fmaxf(0.0f, rem_total - fabsf(y_diff)); - const float f1 = fminf(fmaxf(f1_try, 0.0f), f1_hi); - - const float zsum_hi = fmaxf(0.0f, rem_total - f1 - fabsf(y_diff)); - const float z_sum = fminf(fmaxf(zsum_try, 0.0f), zsum_hi); - if (fabsf(z_diff) > z_sum) { - z_diff = copysignf(z_sum, z_diff); - } - - float y_sum = rem_total - f1 - z_sum; - y_sum = fmaxf(y_sum, 0.0f); - if (fabsf(y_diff) > y_sum) { - y_diff = copysignf(y_sum, y_diff); - } - - f[1] = f1; - f[9] = 0.5f * (z_sum + z_diff); - f[15] = 0.5f * (z_sum - z_diff); - f[7] = 0.5f * (y_sum + y_diff); - f[13] = 0.5f * (y_sum - y_diff); } -__device__ inline void apply_pressure_outlet_3d(float* __restrict__ f, - const float* __restrict__ f_neb, - float y_coord) +__device__ __forceinline__ void apply_outlet_pull_d3q19( + float* __restrict__ f, + unsigned int x, + unsigned int y, + unsigned int z, + const fpxx* __restrict__ fi_in) { - (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 - - // 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 COLLISION_MODEL == 0 || COLLISION_MODEL == 1 - 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; + float f_neb[NQ]; + const 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)); } -#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]; - 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]; + apply_pressure_outlet_d3q19(f, f_neb); +} #endif + +#if DIM == 2 +__device__ __forceinline__ void apply_inlet_esopull_d2q9( + float* __restrict__ f, + unsigned int x, + unsigned int y, + const fpxx* __restrict__ fi, + unsigned long t) +{ + const float u_target = inlet_target_u((float)y); + const float v_target = 0.0f; + +#if INLET_SCHEME == 0 +#if Y_WALL_BC == 1 + repair_zou_he_west_knowns_d2q9(f, fi, x, y); +#endif + apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target); +#elif INLET_SCHEME == 1 || INLET_SCHEME == 3 + const unsigned long k_neb = linear_index(x + 1u, y); + unsigned long j_neb[NQ]; + compute_neighbors(x + 1u, y, j_neb); + float f_neb[NQ]; + load_f_esopull(k_neb, f_neb, fi, j_neb, t); + #if INLET_SCHEME == 1 + apply_channel_stabilized_inlet_d2q9(f, f_neb, (float)y); + #else + apply_regularized_left_velocity_inlet_d2q9(f, f_neb, u_target, v_target); + #endif +#elif INLET_SCHEME == 2 + apply_equilibrium_left_velocity_inlet_d2q9(f, u_target, v_target); +#else + #error "Unsupported INLET_SCHEME for D2Q9" #endif } -#endif // LATTICE_D3Q19 +__device__ __forceinline__ void apply_outlet_esopull_d2q9( + float* __restrict__ f, + unsigned int x, + unsigned int y, + const fpxx* __restrict__ fi, + unsigned long t) +{ + const unsigned long k_neb = linear_index(x - 1u, y); + unsigned long j_neb[NQ]; + compute_neighbors(x - 1u, y, j_neb); + float f_neb[NQ]; + load_f_esopull(k_neb, f_neb, fi, j_neb, t); + apply_pressure_outlet_d2q9(f, f_neb); +} +#endif -#endif // CELERIS_BOUNDARY_INLET_OUTLET_CUH \ No newline at end of file +#if DIM == 3 +__device__ __forceinline__ void apply_inlet_esopull_d3q19( + float* __restrict__ f, + unsigned int x, + unsigned int y, + unsigned int z, + const fpxx* __restrict__ fi, + unsigned long t) +{ + const float u_target = inlet_target_u((float)y); + const float v_target = 0.0f; + const float w_target = 0.0f; + +#if INLET_SCHEME == 0 + apply_zou_he_left_velocity_inlet_d3q19(f, u_target, v_target, w_target); +#elif INLET_SCHEME == 1 || INLET_SCHEME == 3 + const 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); + #if INLET_SCHEME == 1 + apply_channel_stabilized_inlet_d3q19(f, f_neb, (float)y); + #else + apply_regularized_left_velocity_inlet_d3q19(f, f_neb, u_target, v_target, w_target); + #endif +#elif INLET_SCHEME == 2 + apply_equilibrium_left_velocity_inlet_d3q19(f, u_target, v_target, w_target); +#else + #error "Unsupported INLET_SCHEME for D3Q19" +#endif +} + +__device__ __forceinline__ void apply_outlet_esopull_d3q19( + float* __restrict__ f, + unsigned int x, + unsigned int y, + unsigned int z, + const fpxx* __restrict__ fi, + unsigned long t) +{ + const 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_d3q19(f, f_neb); +} +#endif + +#endif // CELERIS_BOUNDARY_INLET_OUTLET_CUH diff --git a/src/CelerisLab/lbm/kernels/boundary/outlet/pressure_neq.cuh b/src/CelerisLab/lbm/kernels/boundary/outlet/pressure_neq.cuh new file mode 100644 index 0000000..043960d --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/outlet/pressure_neq.cuh @@ -0,0 +1,114 @@ +// CelerisLab – boundary/outlet/pressure_neq.cuh +// Pressure outlet and zero-gradient outlet closures. +// ============================================================================ + +#ifndef CELERIS_BOUNDARY_OUTLET_PRESSURE_NEQ_CUH +#define CELERIS_BOUNDARY_OUTLET_PRESSURE_NEQ_CUH + +#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 + +#ifndef OUTLET_SRT_NEQ_DAMP +#define OUTLET_SRT_NEQ_DAMP 0.50f +#endif + +#if LATTICE_MODEL == LATTICE_D2Q9 +__device__ inline void apply_pressure_outlet_d2q9( + float* __restrict__ f, + const float* __restrict__ f_neb) +{ +#if OUTLET_MODE == 1 + f[2] = f_neb[2]; + f[8] = f_neb[8]; + f[6] = f_neb[6]; +#else + 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 + const 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 COLLISION_MODEL == 0 || COLLISION_MODEL == 1 + 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]; + 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 + +#if LATTICE_MODEL == LATTICE_D3Q19 +__device__ inline void apply_pressure_outlet_d3q19( + float* __restrict__ f, + const float* __restrict__ f_neb) +{ +#if OUTLET_MODE == 1 + 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 + 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 + const 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); + +#if COLLISION_MODEL == 0 || COLLISION_MODEL == 1 + 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]; + 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 + +#endif // CELERIS_BOUNDARY_OUTLET_PRESSURE_NEQ_CUH diff --git a/src/CelerisLab/lbm/kernels/config/config_grid.h b/src/CelerisLab/lbm/kernels/config/config_grid.h index b6ea941..bc7e935 100644 --- a/src/CelerisLab/lbm/kernels/config/config_grid.h +++ b/src/CelerisLab/lbm/kernels/config/config_grid.h @@ -6,8 +6,8 @@ #define NT 256 #define MULT_GPU 0 -#define NX 1351 -#define NY 601 +#define NX 401 +#define NY 201 #define NZ 1 // ---- Lattice model (single source of truth) ---- diff --git a/src/CelerisLab/lbm/kernels/config/config_method.h b/src/CelerisLab/lbm/kernels/config/config_method.h index 0ae2977..e07cec1 100644 --- a/src/CelerisLab/lbm/kernels/config/config_method.h +++ b/src/CelerisLab/lbm/kernels/config/config_method.h @@ -13,19 +13,22 @@ #define LES_CLOSED_FORM 1 #define INLET_PROFILE 0 +#define INLET_SCHEME 3 #define OUTLET_MODE 0 #define OUTLET_BLEND_ALPHA 0.700f #define OUTLET_BACKFLOW_CLAMP 1 -#define Y_WALL_BC 1 +#define Y_WALL_BC 0 #define OMEGA_COLLISION_MIN 0.01f #define OMEGA_COLLISION_MAX 1.960f #define TRT_MAGIC_PARAM 0.187500f // NEQ damping coefficients for inlet/outlet BC reconstruction. -// TRT inlet: damps donor non-equilibrium to reduce inlet noise at high Re. -// SRT outlet: damps donor non-equilibrium to suppress checkerboard noise. -#define INLET_TRT_NEQ_DAMP 0.5000f -#define OUTLET_SRT_NEQ_DAMP 0.5000f +// TRT inlet: donor damping used by the channel_stabilized inlet. +// Regularized inlet: damping on incoming-direction donor NEQ. +// SRT outlet: damped outlet reconstruction to suppress checkerboard noise. +#define INLET_TRT_NEQ_DAMP 0.5000f +#define INLET_REGULARIZED_NEQ_DAMP 0.5000f +#define OUTLET_SRT_NEQ_DAMP 0.5000f #endif diff --git a/src/CelerisLab/lbm/kernels/config/config_objects.h b/src/CelerisLab/lbm/kernels/config/config_objects.h index 8357ada..e7cad65 100644 --- a/src/CelerisLab/lbm/kernels/config/config_objects.h +++ b/src/CelerisLab/lbm/kernels/config/config_objects.h @@ -3,6 +3,6 @@ #ifndef CELERIS_CONFIG_OBJECTS_H #define CELERIS_CONFIG_OBJECTS_H -#define N_OBJS 1 +#define N_OBJS 0 #endif diff --git a/src/CelerisLab/lbm/kernels/config/config_physics.h b/src/CelerisLab/lbm/kernels/config/config_physics.h index 788bfc8..a08cc9b 100644 --- a/src/CelerisLab/lbm/kernels/config/config_physics.h +++ b/src/CelerisLab/lbm/kernels/config/config_physics.h @@ -4,7 +4,7 @@ #define CELERIS_CONFIG_PHYSICS_H #define LBtype float -#define VIS 0.0056250000 +#define VIS 0.0090000000 #define RHO 1.0 #define U0 0.03 diff --git a/src/CelerisLab/lbm/kernels/step/one_step_double.cu b/src/CelerisLab/lbm/kernels/step/one_step_double.cu index bc20ddc..285a77a 100644 --- a/src/CelerisLab/lbm/kernels/step/one_step_double.cu +++ b/src/CelerisLab/lbm/kernels/step/one_step_double.cu @@ -22,18 +22,10 @@ __device__ __forceinline__ void apply_boundary_pull( bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); if (is_inlet(fl) && interior_y) { - float f_neb[NQ]; - unsigned long k_neb = linear_index(x + 1u, y); - for (int i = 0; i < NQ; i++) - f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i)); - apply_parabolic_inlet(f, f_neb, (float)y); + apply_inlet_pull_d2q9(f, x, y, fi_in); } else if (is_outlet(fl) && 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); + apply_outlet_pull_d2q9(f, x, y, fi_in); } else { bounce_back_swap(f); @@ -52,18 +44,10 @@ __device__ __forceinline__ void apply_boundary_pull_3d( bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); if (is_inlet(fl) && 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_parabolic_inlet_3d(f, f_neb, (float)y); + apply_inlet_pull_d3q19(f, x, y, z, fi_in); } else if (is_outlet(fl) && 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); + apply_outlet_pull_d3q19(f, x, y, z, fi_in); } else { bounce_back_swap(f); @@ -85,6 +69,7 @@ void OneStep( uint16_t fl = flag[k]; unsigned long j[NQ]; + const bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); compute_neighbors(x, y, j); float f[NQ]; @@ -108,8 +93,15 @@ void OneStep( #endif } - // Collision (fluid only) - if (is_fluid(fl)) { + // Collision + // Normal path: fluid nodes collide. + // Path A repair: when using the local Zou-He inlet, the west inlet nodes are + // still tagged SOLID|BC_INLET ghost nodes, but their post-BC state is later + // used as a pull source for x=1. Leaving that ghost state uncollided was the + // main donor/ghost semantic bug behind the observed inlet blow-ups. + const bool collide_inlet_ghost = is_inlet(fl) && interior_y + && inlet_scheme_uses_post_collision_ghost(); + if (is_fluid(fl) || collide_inlet_ghost) { float rho_n, ux, uy; compute_rho_u(f, rho_n, ux, uy); collide_dispatch(f, rho_n, ux, uy); @@ -124,6 +116,7 @@ void OneStep( uint16_t fl = flag[k]; unsigned long j[NQ]; + const bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); compute_neighbors(k, j); float f[NQ]; @@ -138,7 +131,9 @@ void OneStep( if (is_fluid(fl) && (y == 1u || y == (unsigned int)(NY - 2))) apply_wall_bb_d3q19_y_pull(y, f, fi_in, k); - if (is_fluid(fl)) { + const bool collide_inlet_ghost = is_inlet(fl) && interior_y + && inlet_scheme_uses_post_collision_ghost(); + if (is_fluid(fl) || collide_inlet_ghost) { float rho_n, ux, uy, uz; compute_rho_u(f, rho_n, ux, uy, uz); collide_dispatch(f, rho_n, ux, uy, uz); diff --git a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu index 2b215f0..c5eefbf 100644 --- a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu +++ b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu @@ -29,20 +29,10 @@ __device__ __forceinline__ void apply_boundary_esopull( bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); if (is_inlet(fl) && interior_y) { - unsigned long k_neb = linear_index(x + 1u, y); - unsigned long j_neb[NQ]; - compute_neighbors(x + 1u, y, j_neb); - float f_neb[NQ]; - load_f_esopull(k_neb, f_neb, fi, j_neb, t); - apply_parabolic_inlet(f, f_neb, (float)y); + apply_inlet_esopull_d2q9(f, x, y, fi, t); } else if (is_outlet(fl) && interior_y) { - unsigned long k_neb = linear_index(x - 1u, y); - unsigned long j_neb[NQ]; - compute_neighbors(x - 1u, y, 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); + apply_outlet_esopull_d2q9(f, x, y, fi, t); } else { bounce_back_swap(f); @@ -62,20 +52,10 @@ __device__ __forceinline__ void apply_boundary_esopull_3d( bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); if (is_inlet(fl) && 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); + apply_inlet_esopull_d3q19(f, x, y, z, fi, t); } else if (is_outlet(fl) && 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); + apply_outlet_esopull_d3q19(f, x, y, z, fi, t); } else { bounce_back_swap(f); @@ -99,6 +79,7 @@ void EsoPullStep( uint16_t fl = flag[k]; unsigned long j[NQ]; + const bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); compute_neighbors(x, y, j); float f[NQ]; @@ -118,8 +99,13 @@ void EsoPullStep( if (is_obstacle(fl) && !is_curved(fl)) bounce_back_swap(f); - // Collision (fluid only) - if (is_fluid(fl)) { + // Collision + // Same donor/ghost warning as the double-buffer path: for the local Zou-He + // inlet, the ghost-node state must be regularized before the next pull uses + // it as the source state for the first interior column. + const bool collide_inlet_ghost = is_inlet(fl) && interior_y + && inlet_scheme_uses_post_collision_ghost(); + if (is_fluid(fl) || collide_inlet_ghost) { float rho_n, ux, uy; compute_rho_u(f, rho_n, ux, uy); collide_dispatch(f, rho_n, ux, uy); @@ -134,6 +120,7 @@ void EsoPullStep( uint16_t fl = flag[k]; unsigned long j[NQ]; + const bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1)); compute_neighbors(k, j); float f[NQ]; @@ -151,7 +138,9 @@ void EsoPullStep( if (is_obstacle(fl) && !is_curved(fl)) bounce_back_swap(f); - if (is_fluid(fl)) { + const bool collide_inlet_ghost = is_inlet(fl) && interior_y + && inlet_scheme_uses_post_collision_ghost(); + if (is_fluid(fl) || collide_inlet_ghost) { float rho_n, ux, uy, uz; compute_rho_u(f, rho_n, ux, uy, uz); collide_dispatch(f, rho_n, ux, uy, uz); diff --git a/tests/Rotating_cylinder_validation_plan.md b/tests/Kan99b_validation.md similarity index 100% rename from tests/Rotating_cylinder_validation_plan.md rename to tests/Kan99b_validation.md diff --git a/tests/Sah04_St_validation_matrix.md b/tests/Sah04_validation.md similarity index 100% rename from tests/Sah04_St_validation_matrix.md rename to tests/Sah04_validation.md diff --git a/tests/inlet_plan.md b/tests/inlet_plan.md new file mode 100644 index 0000000..cbfd430 --- /dev/null +++ b/tests/inlet_plan.md @@ -0,0 +1,75 @@ +## Inlet module refactor plan + +The inlet module should be treated as a source-state generator for west ghost nodes, not as a grab-bag of formulas inside one boundary file. In the current solver, inlet cells are `SOLID | BC_INLET` ghost nodes whose state is later pulled by the first interior column. That semantic should stay explicit in the code structure. + +## Current target structure + +- `boundary/inlet/common` + - shared profile logic such as `inlet_target_u(y)` + - shared rho-closure helpers for west velocity inlet + - helper stating whether a scheme requires post-BC ghost collision +- `boundary/inlet/zou_he_local` + - local on-site Zou-He source-state reconstruction + - D2Q9 and D3Q19 west inlet versions +- `boundary/inlet/channel_stabilized` + - donor-based stabilized inlet for high blockage or conservative production runs + - D2Q9 and D3Q19 west inlet versions +- `boundary/inlet/equilibrium` + - full `feq` source-state construction from local rho closure and target velocity + - D2Q9 and D3Q19 versions +- `boundary/inlet/regularized` + - local macro state plus damped donor NEQ on incoming directions + - D2Q9 and D3Q19 versions +- `boundary/outlet/pressure_neq` + - pressure outlet and zero-gradient outlet implementations +- `boundary/inlet_outlet` + - compile-time dispatch only + - no long method bodies + - streaming-specific donor assembly and method selection + +## Design rule for each inlet scheme + +Each scheme should answer the same question: + +- given a west ghost node after pull loading, what source state should be stored there for the next interior pull + +That makes the interface stable even when methods differ in how much donor information they use. + +## Scheme meanings + +| Scheme | Main idea | Best fit | Main caution | +|---|---|---|---| +| `zou_he_local` | textbook local algebraic closure | MRT, research comparisons, clean local baseline | in ghost-source semantics it requires post-BC ghost collision and can be noisy for high-omega SRT | +| `channel_stabilized` | donor-based stabilized inlet | high blockage, production robustness, conservative benchmark work | less pure as a local boundary method | +| `equilibrium` | write full `feq` source state | robust SRT baseline, simple ghost-source compatibility | may suppress inlet NEQ too aggressively for some validation targets | +| `regularized` | local macro state plus damped incoming donor NEQ | middle ground between `equilibrium` and donor-heavy methods | still an experimental family and may need tuning | + +## Collision policy + +Post-BC ghost collision must be owned by the scheme, not hard-coded as a general inlet rule. + +Current policy: + +- `zou_he_local` requires post-BC ghost collision +- `channel_stabilized` does not +- `equilibrium` does not +- `regularized` does not + +This should remain encoded through a helper such as `inlet_scheme_uses_post_collision_ghost()` rather than scattered `INLET_SCHEME == ...` checks. + +## Why this split matters + +The earlier instability work showed that the main difficulty was not a single formula error. The real issue was mixing methods that assume different node semantics: + +- local fluid-boundary formulas such as Zou-He +- ghost-source node architecture in the solver +- different collision-model tolerances, especially SRT versus MRT + +Keeping each inlet method in its own file makes those assumptions visible and lowers the chance of mixing donor and ghost semantics by accident. + +## Next cleanups worth doing later + +1. Split outlet schemes into separate files as more outlet variants are added. +2. If inlet junction handling grows, move row-specific or corner-specific logic into dedicated helpers instead of embedding it inside the main schemes. +3. When validation settles, add a small test matrix document mapping recommended schemes to benchmark families and collision models. +4. If the solver later moves away from ghost-source inlet nodes, keep this folder layout but replace the per-scheme internals rather than rebuilding the whole dispatch layer. diff --git a/tests/run_inlet_channel_diagnostic.py b/tests/run_inlet_channel_diagnostic.py new file mode 100644 index 0000000..0f88a3f --- /dev/null +++ b/tests/run_inlet_channel_diagnostic.py @@ -0,0 +1,574 @@ +# CelerisLab/tests/run_inlet_channel_diagnostic.py +"""Empty-channel inlet diagnostic: field snapshots and line profiles. + +Runs no-cylinder channel flows to isolate inlet / wall / outlet effects before +adding a body. See user matrix in module docstring sections A–C. + +Usage:: + + conda run -n pycuda_3_10 python tests/run_inlet_channel_diagnostic.py --part all + conda run -n pycuda_3_10 python tests/run_inlet_channel_diagnostic.py --part a + conda run -n pycuda_3_10 python tests/run_inlet_channel_diagnostic.py --part b --nx 401 --ny 201 + +Output (default ``tests/output/inlet_channel_diag/``):: + + A_baseline/{SRT,MRT}/snapshots/step_XXXXXX.{npz,png} + B_matrix/caseNN_.../snapshots/... + B_matrix/caseNN_.../lines/ux_lines_stepXXXXXX.png + summary.csv +""" + +from __future__ import annotations + +import argparse +import csv +import json +import os +import sys +import tempfile +from dataclasses import dataclass +from typing import Any, Dict, Iterable, List, Optional, Sequence, Tuple + +import numpy as np + +_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json") + +# Default snapshot steps for parts A and B. +DEFAULT_SNAPSHOT_STEPS: Tuple[int, ...] = (100, 500, 1000, 1500, 2000) + + +@dataclass(frozen=True) +class CaseSpec: + """One empty-channel configuration.""" + + case_id: str + label: str + inlet_scheme: str + y_wall_bc: str + outlet_mode: str + collision: str = "MRT" + + +def _load_json(path: str) -> dict: + with open(path, "r", encoding="utf-8") as f: + return json.load(f) + + +def _write_json(path: str, payload: dict) -> None: + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + with open(path, "w", encoding="utf-8") as f: + json.dump(payload, f, indent=2) + + +def vorticity_z(ux: np.ndarray, uy: np.ndarray) -> np.ndarray: + """ωz = ∂uy/∂x − ∂ux/∂y on (ny, nx) arrays.""" + ux = np.asarray(ux, dtype=np.float64) + uy = np.asarray(uy, dtype=np.float64) + return np.gradient(uy, axis=1) - np.gradient(ux, axis=0) + + +def _line_y_indices(ny: int) -> List[Tuple[int, str]]: + return [ + (1, "y1"), + (ny // 2, f"y{ny // 2}"), + (ny - 2, f"y{ny - 2}"), + ] + + +def _build_cfg( + base: dict, + *, + nx: int, + ny: int, + collision: str, + inlet_scheme: str, + inlet_profile: str, + y_wall_bc: str, + outlet_mode: str, + velocity: float, + viscosity: float, +) -> dict: + cfg = json.loads(json.dumps(base)) + cfg["grid"]["nx"] = int(nx) + cfg["grid"]["ny"] = int(ny) + cfg["grid"]["nz"] = 1 + cfg["physics"]["velocity"] = float(velocity) + cfg["physics"]["viscosity"] = float(viscosity) + cfg["physics"]["rho"] = 1.0 + cfg["method"]["collision"] = str(collision).upper() + cfg["method"]["streaming"] = "double_buffer" + cfg["method"]["store_precision"] = "FP32" + cfg["method"]["les"]["enabled"] = False + cfg["method"]["inlet"]["profile"] = str(inlet_profile) + cfg["method"]["inlet"]["scheme"] = str(inlet_scheme) + cfg["method"]["y_wall_bc"] = str(y_wall_bc) + cfg["method"]["outlet"]["mode"] = str(outlet_mode) + return cfg + + +def _field_stats(rho: np.ndarray, ux: np.ndarray, vort: np.ndarray) -> Dict[str, float]: + def _f(a: np.ndarray) -> float: + b = a[np.isfinite(a)] + return float("nan") if b.size == 0 else float(np.max(np.abs(b))) + + return { + "rho_min": float(np.nanmin(rho)) if np.isfinite(rho).any() else float("nan"), + "rho_max": float(np.nanmax(rho)) if np.isfinite(rho).any() else float("nan"), + "ux_max": _f(ux), + "vort_max": _f(vort), + "finite": bool(np.isfinite(rho).all() and np.isfinite(ux).all()), + } + + +def _save_snapshot_npz( + path: str, + *, + step: int, + rho: np.ndarray, + ux: np.ndarray, + uy: np.ndarray, + vort: np.ndarray, + meta: dict, +) -> None: + os.makedirs(os.path.dirname(path), exist_ok=True) + np.savez_compressed( + path, + rho=np.asarray(rho, dtype=np.float32), + ux=np.asarray(ux, dtype=np.float32), + uy=np.asarray(uy, dtype=np.float32), + vort=np.asarray(vort, dtype=np.float32), + step=np.int32(step), + meta_json=np.asarray(json.dumps(meta)), + ) + + +def _save_field_pngs( + out_dir: str, + prefix: str, + *, + rho: np.ndarray, + ux: np.ndarray, + vort: np.ndarray, + title: str, +) -> List[str]: + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return [] + + os.makedirs(out_dir, exist_ok=True) + ny, nx = rho.shape + extent = (0, nx - 1, 0, ny - 1) + paths: List[str] = [] + + def _one(arr: np.ndarray, name: str, cmap: str, sym: bool) -> None: + a = np.asarray(arr, dtype=np.float64) + fin = a[np.isfinite(a)] + if fin.size == 0: + vmin, vmax = -1.0, 1.0 + elif sym: + v = float(np.percentile(np.abs(fin), 99.5)) or 1.0 + vmin, vmax = -v, v + else: + vmin = float(np.percentile(fin, 0.5)) + vmax = float(np.percentile(fin, 99.5)) + if vmax <= vmin: + vmax = vmin + 1.0 + fig, ax = plt.subplots(figsize=(min(16.0, max(8.0, nx / 80.0)), min(8.0, max(3.0, ny / 50.0)))) + im = ax.imshow(a, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax, extent=extent) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_title(f"{title} — {name}") + fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) + fig.tight_layout() + p = os.path.join(out_dir, f"{prefix}_{name}.png") + fig.savefig(p, dpi=140, bbox_inches="tight") + plt.close(fig) + paths.append(p) + + _one(rho, "rho", "viridis", sym=False) + _one(ux, "ux", "RdBu_r", sym=True) + _one(vort, "vort", "RdBu_r", sym=True) + return paths + + +def _save_line_plots( + path: str, + *, + rho: np.ndarray, + ux: np.ndarray, + step: int, + case_label: str, + y_rows: Sequence[Tuple[int, str]], +) -> None: + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return + + ny, nx = rho.shape + x = np.arange(nx, dtype=np.float64) + + fig, axes = plt.subplots(2, 1, figsize=(min(14.0, max(8.0, nx / 60.0)), 7.0), sharex=True) + for y_idx, y_lab in y_rows: + y_idx = int(np.clip(y_idx, 0, ny - 1)) + axes[0].plot(x, ux[y_idx, :], label=y_lab, linewidth=1.2) + axes[1].plot(x, rho[y_idx, :], label=y_lab, linewidth=1.2) + + axes[0].set_ylabel("u_x") + axes[0].legend(loc="best", fontsize=8) + axes[0].grid(True, alpha=0.3) + axes[1].set_ylabel("rho") + axes[1].set_xlabel("x (lattice)") + axes[1].legend(loc="best", fontsize=8) + axes[1].grid(True, alpha=0.3) + fig.suptitle(f"{case_label} — line profiles at step {step}") + fig.tight_layout() + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + + +def _run_channel( + case: CaseSpec, + *, + base_cfg: dict, + nx: int, + ny: int, + velocity: float, + viscosity: float, + out_root: str, + snapshot_steps: Sequence[int], + max_step: int, + save_png: bool, + save_lines: bool, + stop_on_nan: bool, +) -> List[Dict[str, Any]]: + """Run one case; write snapshots and optional line plots. Return summary rows.""" + sys.path.insert(0, os.path.join(_REPO, "src")) + from CelerisLab import Simulation # noqa: WPS433 + + cfg = _build_cfg( + base_cfg, + nx=nx, + ny=ny, + collision=case.collision, + inlet_scheme=case.inlet_scheme, + inlet_profile="uniform", + y_wall_bc=case.y_wall_bc, + outlet_mode=case.outlet_mode, + velocity=velocity, + viscosity=viscosity, + ) + bdoc = {"objects": []} + + run_dir = os.path.join(out_root, case.case_id) + snap_dir = os.path.join(run_dir, "snapshots") + line_dir = os.path.join(run_dir, "lines") + os.makedirs(snap_dir, exist_ok=True) + if save_lines: + os.makedirs(line_dir, exist_ok=True) + + tmpd = tempfile.mkdtemp(prefix="inlet_diag_") + lbm_tmp = os.path.join(tmpd, "config_lbm.json") + body_tmp = os.path.join(tmpd, "config_body.json") + _write_json(lbm_tmp, cfg) + _write_json(body_tmp, bdoc) + + meta_base = { + "case_id": case.case_id, + "label": case.label, + "nx": nx, + "ny": ny, + "inlet_scheme": case.inlet_scheme, + "inlet_profile": "uniform", + "y_wall_bc": case.y_wall_bc, + "outlet_mode": case.outlet_mode, + "collision": case.collision, + "velocity": velocity, + "viscosity": viscosity, + } + _write_json(os.path.join(run_dir, "case_meta.json"), meta_base) + + sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp) + sim.initialize() + + y_rows = _line_y_indices(ny) + want_steps = sorted({int(s) for s in snapshot_steps if 0 < int(s) <= max_step}) + next_snap = 0 + rows: List[Dict[str, Any]] = [] + + for step in range(1, max_step + 1): + sim.step(1) + if next_snap < len(want_steps) and step == want_steps[next_snap]: + macro = sim.get_macroscopic() + rho = np.asarray(macro["rho"], dtype=np.float64) + ux = np.asarray(macro["ux"], dtype=np.float64) + uy = np.asarray(macro["uy"], dtype=np.float64) + vort = vorticity_z(ux, uy) + stats = _field_stats(rho, ux, vort) + stem = f"step_{step:06d}" + meta = {**meta_base, "step": step, **stats} + npz_path = os.path.join(snap_dir, f"{stem}.npz") + _save_snapshot_npz( + npz_path, + step=step, + rho=rho, + ux=ux, + uy=uy, + vort=vort, + meta=meta, + ) + if save_png: + _save_field_pngs( + snap_dir, + stem, + rho=rho, + ux=ux, + vort=vort, + title=f"{case.label} step {step}", + ) + if save_lines: + line_png = os.path.join(line_dir, f"lines_{stem}.png") + _save_line_plots( + line_png, + rho=rho, + ux=ux, + step=step, + case_label=case.label, + y_rows=y_rows, + ) + # Also save raw 1D data for replotting. + line_npz = os.path.join(line_dir, f"lines_{stem}.npz") + payload = {"x": np.arange(nx, dtype=np.float32)} + for y_idx, y_lab in y_rows: + payload[f"ux_{y_lab}"] = ux[y_idx, :].astype(np.float32) + payload[f"rho_{y_lab}"] = rho[y_idx, :].astype(np.float32) + payload["step"] = np.int32(step) + np.savez_compressed(line_npz, **payload) + + rows.append( + { + "case_id": case.case_id, + "label": case.label, + "collision": case.collision, + "inlet_scheme": case.inlet_scheme, + "y_wall_bc": case.y_wall_bc, + "outlet_mode": case.outlet_mode, + "step": step, + **stats, + "npz": npz_path, + } + ) + if stop_on_nan and not stats["finite"]: + sim.close() + rows[-1]["stopped_early"] = True + return rows + next_snap += 1 + + sim.close() + return rows + + +def _part_a_cases() -> List[CaseSpec]: + # Kan99b-style baseline: zou_he + free_slip + neq_extrap; SRT and MRT. + base = CaseSpec( + case_id="", + label="", + inlet_scheme="zou_he_local", + y_wall_bc="free_slip", + outlet_mode="neq_extrap", + ) + out: List[CaseSpec] = [] + for coll in ("SRT", "MRT"): + cid = f"A_{coll.lower()}_zouhe_fs_neq" + out.append( + CaseSpec( + case_id=cid, + label=f"A baseline {coll} zou_he / free_slip / neq_extrap", + inlet_scheme=base.inlet_scheme, + y_wall_bc=base.y_wall_bc, + outlet_mode=base.outlet_mode, + collision=coll, + ) + ) + return out + + +def _part_b_cases() -> List[CaseSpec]: + return [ + CaseSpec( + case_id="B_case01_zouhe_fs_neq", + label="1 zou_he / free_slip / neq_extrap", + inlet_scheme="zou_he_local", + y_wall_bc="free_slip", + outlet_mode="neq_extrap", + collision="MRT", + ), + CaseSpec( + case_id="B_case02_zouhe_bb_neq", + label="2 zou_he / bounce_back / neq_extrap", + inlet_scheme="zou_he_local", + y_wall_bc="bounce_back", + outlet_mode="neq_extrap", + collision="MRT", + ), + CaseSpec( + case_id="B_case03_zouhe_fs_zgrad", + label="3 zou_he / free_slip / zero_gradient", + inlet_scheme="zou_he_local", + y_wall_bc="free_slip", + outlet_mode="zero_gradient", + collision="MRT", + ), + CaseSpec( + case_id="B_case04_stab_fs_neq", + label="4 channel_stabilized / free_slip / neq_extrap", + inlet_scheme="channel_stabilized", + y_wall_bc="free_slip", + outlet_mode="neq_extrap", + collision="MRT", + ), + ] + + +def _write_summary_csv(path: str, rows: Sequence[Dict[str, Any]]) -> None: + if not rows: + return + keys = [ + "case_id", + "label", + "collision", + "inlet_scheme", + "y_wall_bc", + "outlet_mode", + "step", + "finite", + "rho_min", + "rho_max", + "ux_max", + "vort_max", + "stopped_early", + ] + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + with open(path, "w", encoding="utf-8", newline="") as f: + w = csv.DictWriter(f, fieldnames=keys, extrasaction="ignore") + w.writeheader() + for r in rows: + w.writerow(r) + + +def main() -> int: + ap = argparse.ArgumentParser(description="Empty-channel inlet diagnostic (fields + line profiles)") + ap.add_argument( + "--part", + choices=("a", "b", "all"), + default="all", + help="A=baseline SRT/MRT; B=4-case matrix; all=both", + ) + ap.add_argument("--nx", type=int, default=401, help="Channel length (lattice)") + ap.add_argument("--ny", type=int, default=201, help="Channel height (lattice)") + ap.add_argument("--velocity", type=float, default=0.03, help="Uniform inlet U0") + ap.add_argument("--viscosity", type=float, default=0.009, help="Kinematic viscosity") + ap.add_argument( + "--out-dir", + type=str, + default=os.path.join(_REPO, "tests", "output", "inlet_channel_diag"), + ) + ap.add_argument( + "--steps", + type=str, + default="", + help="Comma-separated snapshot steps (default: 100,500,1000,1500,2000)", + ) + ap.add_argument("--no-png", action="store_true", help="Skip rho/ux/vort PNG maps") + ap.add_argument("--no-lines", action="store_true", help="Skip ux/rho line-profile plots") + ap.add_argument( + "--continue-on-nan", + action="store_true", + help="Keep stepping after non-finite fields (default: stop case early)", + ) + args = ap.parse_args() + + if not os.path.isfile(_DEFAULT_LBM): + print(f"Missing config: {_DEFAULT_LBM}", file=sys.stderr) + return 2 + + if args.steps.strip(): + snap_steps = tuple(int(s.strip()) for s in args.steps.split(",") if s.strip()) + else: + snap_steps = DEFAULT_SNAPSHOT_STEPS + max_step = max(snap_steps) + + base_cfg = _load_json(_DEFAULT_LBM) + out_dir = os.path.abspath(args.out_dir) + os.makedirs(out_dir, exist_ok=True) + + cases: List[CaseSpec] = [] + if args.part in ("a", "all"): + cases.extend(_part_a_cases()) + if args.part in ("b", "all"): + cases.extend(_part_b_cases()) + + save_png = not args.no_png + # Part A: field maps only; Part B: fields + line plots. + all_rows: List[Dict[str, Any]] = [] + + for case in cases: + part = "A_baseline" if case.case_id.startswith("A_") else "B_matrix" + root = os.path.join(out_dir, part) + save_lines = not args.no_lines and part == "B_matrix" + print(f"--- {case.case_id}: {case.label} ({case.collision}) ---", flush=True) + rows = _run_channel( + case, + base_cfg=base_cfg, + nx=args.nx, + ny=args.ny, + velocity=args.velocity, + viscosity=args.viscosity, + out_root=root, + snapshot_steps=snap_steps, + max_step=max_step, + save_png=save_png, + save_lines=save_lines, + stop_on_nan=not args.continue_on_nan, + ) + all_rows.extend(rows) + for r in rows: + fin = "OK" if r.get("finite") else "NONFINITE" + print( + f" step {r['step']:5d} {fin} rho=[{r['rho_min']:.4f},{r['rho_max']:.4f}] " + f"ux_max={r['ux_max']:.4f} vort_max={r['vort_max']:.4f}", + flush=True, + ) + if rows and rows[-1].get("stopped_early"): + print(" (stopped early: non-finite fields)", flush=True) + + summary_path = os.path.join(out_dir, "summary.csv") + _write_summary_csv(summary_path, all_rows) + manifest = { + "snapshot_steps": list(snap_steps), + "nx": args.nx, + "ny": args.ny, + "velocity": args.velocity, + "viscosity": args.viscosity, + "line_y_indices": [{"y": y, "label": lab} for y, lab in _line_y_indices(args.ny)], + "cases": [c.case_id for c in cases], + } + _write_json(os.path.join(out_dir, "manifest.json"), manifest) + + print(f"Wrote: {summary_path}", flush=True) + print(f"Wrote: {os.path.join(out_dir, 'manifest.json')}", flush=True) + print(f"Output root: {out_dir}", flush=True) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/run_inlet_ghost_timing_experiment.py b/tests/run_inlet_ghost_timing_experiment.py new file mode 100644 index 0000000..d8ecc0e --- /dev/null +++ b/tests/run_inlet_ghost_timing_experiment.py @@ -0,0 +1,460 @@ +# CelerisLab/tests/run_inlet_ghost_timing_experiment.py +"""Minimal ghost-inlet timing experiments (post field-diagnostic). + +Experiment 1 — DDF time series at inlet center vs first interior column: + (x=0, y=NY/2) and (x=1, y=NY/2), steps 1..N. + Populations f1,f2,f5,f6,f7,f8 plus rho, ux (macro and sum f). + +Experiment 2 — Same channel with ``inlet.collide=false`` vs ``true``: + When collide is on, inlet ghost nodes undergo collision after Zou-He BC. + Compare rho_max / ux_max vs step to test ghost-source timing hypothesis. + +Usage:: + + conda run -n pycuda_3_10 python tests/run_inlet_ghost_timing_experiment.py + conda run -n pycuda_3_10 python tests/run_inlet_ghost_timing_experiment.py --exp 1 --steps 50 + conda run -n pycuda_3_10 python tests/run_inlet_ghost_timing_experiment.py --exp 2 --steps 500 +""" + +from __future__ import annotations + +import argparse +import csv +import json +import os +import sys +import tempfile +from typing import Any, Dict, List, Optional, Sequence, Tuple + +import numpy as np + +_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json") + +# D2Q9 indices logged (see zou_he_velocity.cuh). +POP_IDX: Tuple[int, ...] = (1, 2, 5, 6, 7, 8) +POP_NAMES: Tuple[str, ...] = tuple(f"f{i}" for i in POP_IDX) + +# cx, cy for macroscopic ux, uy from local f (matches descriptors.cuh D2Q9). +_CX = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float64) +_CY = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float64) + + +def _load_json(path: str) -> dict: + with open(path, "r", encoding="utf-8") as f: + return json.load(f) + + +def _write_json(path: str, payload: dict) -> None: + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + with open(path, "w", encoding="utf-8") as f: + json.dump(payload, f, indent=2) + + +def _build_cfg( + base: dict, + *, + nx: int, + ny: int, + collision: str, + inlet_collide: bool, + velocity: float, + viscosity: float, +) -> dict: + cfg = json.loads(json.dumps(base)) + cfg["grid"]["nx"] = int(nx) + cfg["grid"]["ny"] = int(ny) + cfg["grid"]["nz"] = 1 + cfg["physics"]["velocity"] = float(velocity) + cfg["physics"]["viscosity"] = float(viscosity) + cfg["physics"]["rho"] = 1.0 + cfg["method"]["collision"] = str(collision).upper() + cfg["method"]["streaming"] = "double_buffer" + cfg["method"]["store_precision"] = "FP32" + cfg["method"]["les"]["enabled"] = False + cfg["method"]["inlet"]["profile"] = "uniform" + cfg["method"]["inlet"]["scheme"] = "zou_he_local" + cfg["method"]["inlet"]["collide"] = bool(inlet_collide) + cfg["method"]["y_wall_bc"] = "free_slip" + cfg["method"]["outlet"]["mode"] = "neq_extrap" + return cfg + + +def _macro_from_f(f: np.ndarray) -> Tuple[float, float, float]: + f = np.asarray(f, dtype=np.float64) + rho = float(np.sum(f)) + if abs(rho) < 1e-14: + return rho, 0.0, 0.0 + ux = float(np.dot(f, _CX) / rho) + uy = float(np.dot(f, _CY) / rho) + return rho, ux, uy + + +def _sample_node(ddf_qnyx: np.ndarray, x: int, y: int) -> Dict[str, float]: + f = ddf_qnyx[:, y, x].astype(np.float64) + rho_m, ux_m, uy_m = _macro_from_f(f) + out: Dict[str, float] = { + "rho_sum": rho_m, + "ux_macro": ux_m, + "uy_macro": uy_m, + } + for i, name in zip(POP_IDX, POP_NAMES): + out[name] = float(f[i]) + return out + + +def _run_steps( + cfg: dict, + *, + steps: int, + y_mid: int, + probe_x: Sequence[int] = (0, 1), +) -> Tuple[List[Dict[str, Any]], Dict[str, Any]]: + sys.path.insert(0, os.path.join(_REPO, "src")) + from CelerisLab import Simulation # noqa: WPS433 + + bdoc = {"objects": []} + tmpd = tempfile.mkdtemp(prefix="ghost_timing_") + lbm_tmp = os.path.join(tmpd, "config_lbm.json") + body_tmp = os.path.join(tmpd, "config_body.json") + with open(lbm_tmp, "w", encoding="utf-8") as f: + json.dump(cfg, f) + with open(body_tmp, "w", encoding="utf-8") as f: + json.dump(bdoc, f) + + sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp) + sim.initialize() + nx, ny = sim.lbm_cfg.nx, sim.lbm_cfg.ny + y_mid = int(np.clip(y_mid, 1, ny - 2)) + + rows: List[Dict[str, Any]] = [] + for step in range(1, int(steps) + 1): + sim.step(1) + sim.field.download_ddf() + farr = sim.field.ddf.reshape(sim.lbm_cfg.nq, ny, nx) + + rec: Dict[str, Any] = {"step": step} + for x in probe_x: + tag = "inlet" if x == 0 else "interior" + s = _sample_node(farr, int(x), y_mid) + for k, v in s.items(): + rec[f"{tag}_{k}"] = v + + # Pull semantics at interior: f[2] is read from stored f[2] at x=0 (same link index). + rec["cross_f2_match"] = abs(rec["interior_f2"] - rec["inlet_f2"]) < 1e-5 + rec["cross_f1_inlet_to_int_pull"] = float( + farr[1, y_mid, 1] + ) # after step, what x=1 holds in f1 + rec["delta_inlet_f1"] = ( + float("nan") if step == 1 else rec["inlet_f1"] - rows[-1]["inlet_f1"] + ) + rec["delta_inlet_ux"] = ( + float("nan") + if step == 1 + else rec["inlet_ux_macro"] - rows[-1]["inlet_ux_macro"] + ) + + macro = sim.get_macroscopic() + rho_f = np.asarray(macro["rho"], dtype=np.float64) + ux_f = np.asarray(macro["ux"], dtype=np.float64) + rec["domain_rho_max"] = float(np.nanmax(rho_f)) + rec["domain_rho_min"] = float(np.nanmin(rho_f)) + rec["domain_ux_max"] = float(np.nanmax(np.abs(ux_f))) + rec["finite"] = bool(np.isfinite(rho_f).all() and np.isfinite(ux_f).all()) + + rows.append(rec) + + meta = { + "nx": nx, + "ny": ny, + "y_mid": y_mid, + "probe_x": list(probe_x), + "inlet_collide": bool(cfg["method"]["inlet"].get("collide", False)), + "collision": cfg["method"]["collision"], + "inlet_scheme": cfg["method"]["inlet"]["scheme"], + "y_wall_bc": cfg["method"]["y_wall_bc"], + "outlet_mode": cfg["method"]["outlet"]["mode"], + "velocity": cfg["physics"]["velocity"], + "viscosity": cfg["physics"]["viscosity"], + } + sim.close() + return rows, meta + + +def _write_csv(path: str, rows: Sequence[Dict[str, Any]]) -> None: + if not rows: + return + keys: List[str] = [] + for r in rows: + for k in r: + if k not in keys: + keys.append(k) + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + with open(path, "w", encoding="utf-8", newline="") as f: + w = csv.DictWriter(f, fieldnames=keys) + w.writeheader() + w.writerows(rows) + + +def _plot_exp1(out_dir: str, rows: Sequence[Dict[str, Any]], y_mid: int) -> List[str]: + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return [] + + steps = [int(r["step"]) for r in rows] + paths: List[str] = [] + + def _ts(key: str, label: str, ax, **kw): + ax.plot(steps, [r[key] for r in rows], label=label, **kw) + + # Populations + fig, axes = plt.subplots(2, 1, figsize=(11, 7), sharex=True) + for name in POP_NAMES: + _ts(f"inlet_{name}", f"inlet {name}", axes[0], linewidth=1.0) + axes[0].set_ylabel("f at x=0") + axes[0].legend(ncol=3, fontsize=7, loc="best") + axes[0].grid(True, alpha=0.3) + for name in POP_NAMES: + _ts(f"interior_{name}", f"int {name}", axes[1], linewidth=1.0) + axes[1].set_ylabel("f at x=1") + axes[1].set_xlabel("step") + axes[1].legend(ncol=3, fontsize=7, loc="best") + axes[1].grid(True, alpha=0.3) + fig.suptitle(f"Exp1 populations (y={y_mid})") + fig.tight_layout() + p1 = os.path.join(out_dir, "exp1_populations.png") + fig.savefig(p1, dpi=150, bbox_inches="tight") + plt.close(fig) + paths.append(p1) + + # rho / ux + step-to-step deltas + fig, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True) + _ts("inlet_ux_macro", "inlet ux", axes[0]) + _ts("interior_ux_macro", "interior ux", axes[0]) + axes[0].axhline(0.03, color="k", ls="--", lw=0.8, label="U0") + axes[0].set_ylabel("u_x") + axes[0].legend(fontsize=8) + axes[0].grid(True, alpha=0.3) + _ts("inlet_rho_sum", "inlet rho", axes[1]) + _ts("interior_rho_sum", "interior rho", axes[1]) + axes[1].set_ylabel("rho") + axes[1].grid(True, alpha=0.3) + _ts("delta_inlet_f1", "|Δf1| inlet", axes[2]) + axes[2].set_ylabel("Δf1") + axes[2].set_xlabel("step") + axes[2].grid(True, alpha=0.3) + fig.suptitle(f"Exp1 macro / inlet f1 increment (y={y_mid})") + fig.tight_layout() + p2 = os.path.join(out_dir, "exp1_macro_delta.png") + fig.savefig(p2, dpi=150, bbox_inches="tight") + plt.close(fig) + paths.append(p2) + + return paths + + +def _plot_exp2(out_dir: str, rows_a: Sequence[Dict[str, Any]], rows_b: Sequence[Dict[str, Any]]) -> List[str]: + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return [] + + fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True) + for rows, lab, c in ( + (rows_a, "ghost (no collide)", "C0"), + (rows_b, "inlet collide", "C1"), + ): + steps = [int(r["step"]) for r in rows] + axes[0].plot(steps, [r["domain_rho_max"] for r in rows], label=lab, color=c) + axes[1].plot(steps, [r["domain_ux_max"] for r in rows], label=lab, color=c) + axes[0].set_ylabel("rho_max") + axes[0].legend() + axes[0].grid(True, alpha=0.3) + axes[1].set_ylabel("|ux|_max") + axes[1].set_xlabel("step") + axes[1].legend() + axes[1].grid(True, alpha=0.3) + fig.suptitle("Exp2: ghost vs inlet collide") + fig.tight_layout() + p = os.path.join(out_dir, "exp2_stability_compare.png") + fig.savefig(p, dpi=150, bbox_inches="tight") + plt.close(fig) + return [p] + + +def _oscillation_summary(rows: Sequence[Dict[str, Any]], *, last_n: int = 20) -> Dict[str, float]: + """High-frequency proxy: std of inlet f1 and ux over the last *last_n* steps.""" + if len(rows) < 2: + return {} + tail = rows[-min(last_n, len(rows)) :] + f1 = np.array([r["inlet_f1"] for r in tail], dtype=np.float64) + ux = np.array([r["inlet_ux_macro"] for r in tail], dtype=np.float64) + d1 = np.array([r["delta_inlet_f1"] for r in tail if np.isfinite(r["delta_inlet_f1"])], dtype=np.float64) + return { + "std_inlet_f1_last": float(np.std(f1)), + "std_inlet_ux_last": float(np.std(ux)), + "mean_abs_delta_f1_last": float(np.mean(np.abs(d1))) if d1.size else float("nan"), + } + + +def run_exp1( + base: dict, + *, + out_dir: str, + nx: int, + ny: int, + steps: int, + collision: str, + velocity: float, + viscosity: float, +) -> None: + y_mid = ny // 2 + cfg = _build_cfg( + base, + nx=nx, + ny=ny, + collision=collision, + inlet_collide=False, + velocity=velocity, + viscosity=viscosity, + ) + print(f"Exp1: zou_he ghost inlet, y_mid={y_mid}, steps={steps}", flush=True) + rows, meta = _run_steps(cfg, steps=steps, y_mid=y_mid) + meta["oscillation"] = _oscillation_summary(rows) + + exp_dir = os.path.join(out_dir, "exp1_ddf_timeseries") + os.makedirs(exp_dir, exist_ok=True) + _write_csv(os.path.join(exp_dir, "timeseries.csv"), rows) + _write_json(os.path.join(exp_dir, "meta.json"), meta) + plots = _plot_exp1(exp_dir, rows, y_mid) + for p in plots: + print(f" plot: {p}", flush=True) + + # Console summary for quick read + print(" last 5 steps (inlet center):", flush=True) + for r in rows[-5:]: + print( + f" step {r['step']:3d} f1={r['inlet_f1']:.6f} ux={r['inlet_ux_macro']:.6f} " + f"Δf1={r['delta_inlet_f1']:.2e} rho_max={r['domain_rho_max']:.4f} finite={r['finite']}", + flush=True, + ) + print(f" oscillation: {meta['oscillation']}", flush=True) + print(f"Wrote: {exp_dir}/timeseries.csv", flush=True) + + +def run_exp2( + base: dict, + *, + out_dir: str, + nx: int, + ny: int, + steps: int, + collision: str, + velocity: float, + viscosity: float, +) -> None: + y_mid = ny // 2 + exp_dir = os.path.join(out_dir, "exp2_inlet_collide") + os.makedirs(exp_dir, exist_ok=True) + + summaries: Dict[str, Any] = {} + all_rows: Dict[str, List[Dict[str, Any]]] = {} + + for collide, tag in ((False, "ghost_no_collide"), (True, "ghost_with_collide")): + cfg = _build_cfg( + base, + nx=nx, + ny=ny, + collision=collision, + inlet_collide=collide, + velocity=velocity, + viscosity=viscosity, + ) + print(f"Exp2 [{tag}]: inlet.collide={collide}, steps={steps}", flush=True) + rows, meta = _run_steps(cfg, steps=steps, y_mid=y_mid) + all_rows[tag] = rows + _write_csv(os.path.join(exp_dir, f"{tag}.csv"), rows) + + last_finite = next( + (int(r["step"]) for r in rows if not r.get("finite", True)), + None, + ) + summaries[tag] = { + **meta, + "first_nonfinite_step": last_finite, + "final_rho_max": rows[-1]["domain_rho_max"] if rows else None, + "final_finite": rows[-1].get("finite") if rows else None, + "oscillation": _oscillation_summary(rows), + } + print( + f" final rho_max={summaries[tag]['final_rho_max']:.4f} " + f"finite={summaries[tag]['final_finite']} " + f"first_nonfinite={summaries[tag]['first_nonfinite_step']}", + flush=True, + ) + + _write_json(os.path.join(exp_dir, "summary.json"), summaries) + plots = _plot_exp2(exp_dir, all_rows["ghost_no_collide"], all_rows["ghost_with_collide"]) + for p in plots: + print(f" plot: {p}", flush=True) + print(f"Wrote: {exp_dir}/summary.json", flush=True) + + +def main() -> int: + ap = argparse.ArgumentParser(description="Ghost inlet timing experiments") + ap.add_argument("--exp", choices=("1", "2", "all"), default="all") + ap.add_argument("--steps", type=int, default=50, help="Steps for exp1 (default 50)") + ap.add_argument("--steps2", type=int, default=500, help="Steps for exp2 (default 500)") + ap.add_argument("--nx", type=int, default=401) + ap.add_argument("--ny", type=int, default=201) + ap.add_argument("--collision", default="MRT", choices=("SRT", "TRT", "MRT")) + ap.add_argument("--velocity", type=float, default=0.03) + ap.add_argument("--viscosity", type=float, default=0.009) + ap.add_argument( + "--out-dir", + default=os.path.join(_REPO, "tests", "output", "inlet_ghost_timing"), + ) + args = ap.parse_args() + + if not os.path.isfile(_DEFAULT_LBM): + print(f"Missing {_DEFAULT_LBM}", file=sys.stderr) + return 2 + + base = _load_json(_DEFAULT_LBM) + out_dir = os.path.abspath(args.out_dir) + os.makedirs(out_dir, exist_ok=True) + + if args.exp in ("1", "all"): + run_exp1( + base, + out_dir=out_dir, + nx=args.nx, + ny=args.ny, + steps=args.steps, + collision=args.collision, + velocity=args.velocity, + viscosity=args.viscosity, + ) + if args.exp in ("2", "all"): + run_exp2( + base, + out_dir=out_dir, + nx=args.nx, + ny=args.ny, + steps=args.steps2, + collision=args.collision, + velocity=args.velocity, + viscosity=args.viscosity, + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/run_inlet_scenario_fields.py b/tests/run_inlet_scenario_fields.py new file mode 100644 index 0000000..4a092bc --- /dev/null +++ b/tests/run_inlet_scenario_fields.py @@ -0,0 +1,587 @@ +# CelerisLab/tests/run_inlet_scenario_fields.py +"""Three-scenario inlet field export: final or last-stable step at 5000 LBM steps. + +Scenarios: + - empty_channel: zou_he_local × {SRT,MRT} × {free_slip,bounce_back} + - empty_channel_inlet_matrix: all inlet schemes × {SRT,MRT}, bounce_back only + - kan99b: zou_he_local × {SRT,MRT} × {free_slip,bounce_back}, Re=100, domain M + - sah04_case9: channel_stabilized × {SRT,MRT}, high-blockage case 9 geometry + +Outputs per run (under ``--out-dir``): + fields/final_{rho,ux,vort}.png, fields/final.npz + lines/lines_ux_rho.png, lines/lines.npz + run_meta.json + +Usage:: + + conda run -n pycuda_3_10 python tests/run_inlet_scenario_fields.py + conda run -n pycuda_3_10 python tests/run_inlet_scenario_fields.py --scenario empty_channel + conda run -n pycuda_3_10 python tests/run_inlet_scenario_fields.py --scenario empty_channel_inlet_matrix +""" + +from __future__ import annotations + +import argparse +import csv +import json +import os +import sys +import tempfile +from dataclasses import dataclass, replace +from typing import Any, Dict, List, Optional, Sequence, Tuple + +import numpy as np + +_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) +_DEFAULT_LBM = os.path.join(_REPO, "src", "CelerisLab", "configs", "config_lbm.json") + +# Kan99b lattice contract +_KAN_U_INF = 0.03 +_KAN_D = 30.0 +_KAN_R = 15.0 +_KAN_RE = 100.0 +_KAN_ALPHA = 1.0 + +# Sah04 case 9 (high tier) +_SAH_D = 30 +_SAH_NX = 80 * _SAH_D + 2 +_SAH_NY = 35 +_SAH_CX = 40.0 * _SAH_D + 0.5 +_SAH_CY = 17.0 +_SAH_RE = 200.0 +_SAH_U_MAX = 0.1 + + +@dataclass(frozen=True) +class RunSpec: + """One simulation run specification.""" + + scenario: str + run_id: str + label: str + nx: int + ny: int + collision: str + inlet_scheme: str + inlet_profile: str + y_wall_bc: str + outlet_mode: str + velocity: float + viscosity: float + steps: int + has_cylinder: bool + cylinder_center: Tuple[float, float] = (0.0, 0.0) + cylinder_radius: float = 0.0 + cylinder_omega: float = 0.0 + + +def _load_json(path: str) -> dict: + with open(path, "r", encoding="utf-8") as f: + return json.load(f) + + +def _write_json(path: str, payload: dict) -> None: + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + with open(path, "w", encoding="utf-8") as f: + json.dump(payload, f, indent=2) + + +def vorticity_z(ux: np.ndarray, uy: np.ndarray) -> np.ndarray: + ux = np.asarray(ux, dtype=np.float64) + uy = np.asarray(uy, dtype=np.float64) + return np.gradient(uy, axis=1) - np.gradient(ux, axis=0) + + +def _line_y_indices(ny: int) -> List[Tuple[int, str]]: + return [(1, "y1"), (ny // 2, f"y{ny // 2}"), (ny - 2, f"y{ny - 2}")] + + +_INLET_SCHEMES = ( + "zou_he_local", + "channel_stabilized", + "equilibrium", + "regularized", +) + + +def _empty_channel_inlet_matrix_specs() -> List[RunSpec]: + """All inlet schemes on empty channel, bounce_back, SRT/MRT (5000-step field export).""" + specs: List[RunSpec] = [] + for coll in ("SRT", "MRT"): + for scheme in _INLET_SCHEMES: + run_id = f"{coll.lower()}_{scheme}" + specs.append( + RunSpec( + scenario="empty_channel_inlet_matrix", + run_id=run_id, + label=f"empty {scheme} {coll} bounce_back", + nx=401, + ny=201, + collision=coll, + inlet_scheme=scheme, + inlet_profile="uniform", + y_wall_bc="bounce_back", + outlet_mode="neq_extrap", + velocity=0.03, + viscosity=0.009, + steps=5000, + has_cylinder=False, + ) + ) + return specs + + +def _all_specs() -> List[RunSpec]: + specs: List[RunSpec] = [] + for coll in ("SRT", "MRT"): + for wall in ("free_slip", "bounce_back"): + wid = f"{coll.lower()}_{wall}" + specs.append( + RunSpec( + scenario="empty_channel", + run_id=wid, + label=f"empty zou_he {coll} {wall}", + nx=401, + ny=201, + collision=coll, + inlet_scheme="zou_he_local", + inlet_profile="uniform", + y_wall_bc=wall, + outlet_mode="neq_extrap", + velocity=0.03, + viscosity=0.009, + steps=5000, + has_cylinder=False, + ) + ) + dom_m = (1351, 601, (450.0, 300.0)) + nu_k = _KAN_U_INF * _KAN_D / _KAN_RE + omega = 2.0 * _KAN_ALPHA * _KAN_U_INF / _KAN_D + for coll in ("SRT", "MRT"): + for wall in ("free_slip", "bounce_back"): + wid = f"{coll.lower()}_{wall}" + specs.append( + RunSpec( + scenario="kan99b", + run_id=wid, + label=f"kan99b zou_he {coll} {wall}", + nx=dom_m[0], + ny=dom_m[1], + collision=coll, + inlet_scheme="zou_he_local", + inlet_profile="uniform", + y_wall_bc=wall, + outlet_mode="neq_extrap", + velocity=_KAN_U_INF, + viscosity=nu_k, + steps=5000, + has_cylinder=True, + cylinder_center=dom_m[2], + cylinder_radius=_KAN_R, + cylinder_omega=omega, + ) + ) + u0_mean = _SAH_U_MAX / 1.5 + nu_s = _SAH_U_MAX * _SAH_D / _SAH_RE + for coll in ("SRT", "MRT"): + wid = f"{coll.lower()}_channel_stab" + specs.append( + RunSpec( + scenario="sah04_case9", + run_id=wid, + label=f"sah04 case9 channel_stab {coll}", + nx=_SAH_NX, + ny=_SAH_NY, + collision=coll, + inlet_scheme="channel_stabilized", + inlet_profile="parabolic", + y_wall_bc="bounce_back", + outlet_mode="neq_extrap", + velocity=u0_mean, + viscosity=nu_s, + steps=5000, + has_cylinder=True, + cylinder_center=(_SAH_CX, _SAH_CY), + cylinder_radius=0.5 * _SAH_D, + cylinder_omega=0.0, + ) + ) + return specs + + +def _build_cfg(base: dict, spec: RunSpec) -> dict: + cfg = json.loads(json.dumps(base)) + cfg["grid"]["nx"] = spec.nx + cfg["grid"]["ny"] = spec.ny + cfg["grid"]["nz"] = 1 + cfg["physics"]["velocity"] = float(spec.velocity) + cfg["physics"]["viscosity"] = float(spec.viscosity) + cfg["physics"]["rho"] = 1.0 + cfg["method"]["collision"] = spec.collision.upper() + cfg["method"]["streaming"] = "double_buffer" + cfg["method"]["store_precision"] = "FP32" + cfg["method"]["les"]["enabled"] = False + cfg["method"]["inlet"]["profile"] = spec.inlet_profile + cfg["method"]["inlet"]["scheme"] = spec.inlet_scheme + cfg["method"]["y_wall_bc"] = spec.y_wall_bc + cfg["method"]["outlet"]["mode"] = spec.outlet_mode + return cfg + + +def _body_doc(spec: RunSpec) -> dict: + if not spec.has_cylinder: + return {"objects": []} + return { + "objects": [ + { + "type": "cylinder", + "center": [float(spec.cylinder_center[0]), float(spec.cylinder_center[1])], + "radius": float(spec.cylinder_radius), + "omega": float(spec.cylinder_omega), + } + ] + } + + +def _save_field_pngs( + out_dir: str, + prefix: str, + *, + rho: np.ndarray, + ux: np.ndarray, + vort: np.ndarray, + title: str, +) -> List[str]: + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return [] + + os.makedirs(out_dir, exist_ok=True) + ny, nx = rho.shape + extent = (0, nx - 1, 0, ny - 1) + paths: List[str] = [] + + def _one(arr: np.ndarray, name: str, cmap: str, sym: bool) -> None: + a = np.asarray(arr, dtype=np.float64) + fin = a[np.isfinite(a)] + if fin.size == 0: + vmin, vmax = -1.0, 1.0 + elif sym: + v = float(np.percentile(np.abs(fin), 99.5)) or 1.0 + vmin, vmax = -v, v + else: + vmin = float(np.percentile(fin, 0.5)) + vmax = float(np.percentile(fin, 99.5)) + if vmax <= vmin: + vmax = vmin + 1.0 + fw = min(18.0, max(8.0, nx / 70.0)) + fh = min(10.0, max(3.0, ny / 45.0)) + fig, ax = plt.subplots(figsize=(fw, fh)) + im = ax.imshow(a, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax, extent=extent) + ax.set_xlabel("x") + ax.set_ylabel("y") + ax.set_title(f"{title} — {name}") + fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04) + fig.tight_layout() + p = os.path.join(out_dir, f"{prefix}_{name}.png") + fig.savefig(p, dpi=150, bbox_inches="tight") + plt.close(fig) + paths.append(p) + + _one(rho, "rho", "viridis", sym=False) + _one(ux, "ux", "RdBu_r", sym=True) + _one(vort, "vort", "RdBu_r", sym=True) + return paths + + +def _save_line_plots( + path: str, + *, + rho: np.ndarray, + ux: np.ndarray, + step: int, + label: str, + y_rows: Sequence[Tuple[int, str]], +) -> None: + try: + import matplotlib + + matplotlib.use("Agg") + import matplotlib.pyplot as plt + except ImportError: + return + + ny, nx = rho.shape + x = np.arange(nx, dtype=np.float64) + fig, axes = plt.subplots(2, 1, figsize=(min(14.0, max(8.0, nx / 55.0)), 7.0), sharex=True) + for y_idx, y_lab in y_rows: + yi = int(np.clip(y_idx, 0, ny - 1)) + axes[0].plot(x, ux[yi, :], label=y_lab, linewidth=1.0) + axes[1].plot(x, rho[yi, :], label=y_lab, linewidth=1.0) + axes[0].set_ylabel("u_x") + axes[0].legend(loc="best", fontsize=8) + axes[0].grid(True, alpha=0.3) + axes[1].set_ylabel("rho") + axes[1].set_xlabel("x (lattice)") + axes[1].legend(loc="best", fontsize=8) + axes[1].grid(True, alpha=0.3) + fig.suptitle(f"{label} — ux/rho lines at step {step}") + fig.tight_layout() + os.makedirs(os.path.dirname(path) or ".", exist_ok=True) + fig.savefig(path, dpi=150, bbox_inches="tight") + plt.close(fig) + + +def _snapshot_from_sim(sim) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + macro = sim.get_macroscopic() + rho = np.asarray(macro["rho"], dtype=np.float64) + ux = np.asarray(macro["ux"], dtype=np.float64) + uy = np.asarray(macro["uy"], dtype=np.float64) + vort = vorticity_z(ux, uy) + return rho, ux, uy, vort + + +def _is_stable_fields( + rho: np.ndarray, + ux: np.ndarray, + *, + rho_lo: float = 0.85, + rho_hi: float = 1.25, + ux_cap: float = 0.15, +) -> bool: + """Finite fields within a physically plausible band (reject pre-blow-up states).""" + if not (np.isfinite(rho).all() and np.isfinite(ux).all()): + return False + r0 = float(np.min(rho)) + r1 = float(np.max(rho)) + umax = float(np.max(np.abs(ux))) + return (rho_lo <= r0) and (r1 <= rho_hi) and (umax <= ux_cap) + + +def run_one(spec: RunSpec, base_cfg: dict, out_root: str) -> Dict[str, Any]: + sys.path.insert(0, os.path.join(_REPO, "src")) + import pycuda.driver as cuda + from CelerisLab import Simulation # noqa: WPS433 + + run_dir = os.path.join(out_root, spec.scenario, spec.run_id) + field_dir = os.path.join(run_dir, "fields") + line_dir = os.path.join(run_dir, "lines") + os.makedirs(field_dir, exist_ok=True) + os.makedirs(line_dir, exist_ok=True) + + cfg = _build_cfg(base_cfg, spec) + tmpd = tempfile.mkdtemp(prefix="inlet_scenario_") + lbm_tmp = os.path.join(tmpd, "config_lbm.json") + body_tmp = os.path.join(tmpd, "config_body.json") + _write_json(lbm_tmp, cfg) + _write_json(body_tmp, _body_doc(spec)) + + sim = Simulation(lbm_config_path=lbm_tmp, body_config_path=body_tmp) + if spec.has_cylinder and spec.cylinder_omega != 0.0: + sim.bodies.get(0).state.omega = np.float32(spec.cylinder_omega) + sim.initialize() + + stream = cuda.Stream() + y_rows = _line_y_indices(spec.ny) + + last_good: Optional[Dict[str, Any]] = None + first_bad_step: Optional[int] = None + force_bad_step: Optional[int] = None + + print(f" [{spec.scenario}/{spec.run_id}] {spec.label} steps={spec.steps}", flush=True) + + for step in range(1, spec.steps + 1): + if spec.has_cylinder: + sim.bodies.zero_force_segment_async(stream) + sim.stepper.step( + 1, + action_gpu=sim.bodies.action_gpu, + obs_gpu=sim.bodies.obs_gpu, + stream=stream, + ) + if step % 100 == 0 or step == spec.steps: + stream.synchronize() + sim.bodies.download_obs_full_async(stream) + stream.synchronize() + fvec = sim.bodies.read_force(0) + if not (np.isfinite(fvec[0]) and np.isfinite(fvec[1])): + if force_bad_step is None: + force_bad_step = step + else: + sim.step(1) + + rho, ux, uy, vort = _snapshot_from_sim(sim) + if _is_stable_fields(rho, ux): + last_good = { + "step": step, + "rho": rho.copy(), + "ux": ux.copy(), + "uy": uy.copy(), + "vort": vort.copy(), + } + elif first_bad_step is None: + first_bad_step = step + + sim.close() + + if last_good is None: + raise RuntimeError(f"No finite snapshot for {spec.run_id}") + + out_step = int(last_good["step"]) + rho = last_good["rho"] + ux = last_good["ux"] + uy = last_good["uy"] + vort = last_good["vort"] + requested_final = spec.steps + used_last_stable = out_step < requested_final + + meta = { + "scenario": spec.scenario, + "run_id": spec.run_id, + "label": spec.label, + "nx": spec.nx, + "ny": spec.ny, + "collision": spec.collision, + "inlet_scheme": spec.inlet_scheme, + "inlet_profile": spec.inlet_profile, + "y_wall_bc": spec.y_wall_bc, + "outlet_mode": spec.outlet_mode, + "velocity": spec.velocity, + "viscosity": spec.viscosity, + "requested_steps": requested_final, + "output_step": out_step, + "used_last_stable": used_last_stable, + "first_nonfinite_step": first_bad_step, + "first_force_nonfinite_step": force_bad_step, + "rho_min": float(np.min(rho)), + "rho_max": float(np.max(rho)), + "ux_max": float(np.max(np.abs(ux))), + "vort_max": float(np.max(np.abs(vort[np.isfinite(vort)]))) if np.isfinite(vort).any() else float("nan"), + } + _write_json(os.path.join(run_dir, "run_meta.json"), meta) + + stem = f"step_{out_step:06d}" + np.savez_compressed( + os.path.join(field_dir, "final.npz"), + rho=rho.astype(np.float32), + ux=ux.astype(np.float32), + uy=uy.astype(np.float32), + vort=vort.astype(np.float32), + step=np.int32(out_step), + ) + + title = f"{spec.label} (step {out_step}" + (", last stable" if used_last_stable else ", final") + ")" + pngs = _save_field_pngs(field_dir, "final", rho=rho, ux=ux, vort=vort, title=title) + _save_line_plots( + os.path.join(line_dir, "lines_ux_rho.png"), + rho=rho, + ux=ux, + step=out_step, + label=spec.label, + y_rows=y_rows, + ) + line_payload: Dict[str, Any] = {"x": np.arange(spec.nx, dtype=np.float32), "step": np.int32(out_step)} + for y_idx, y_lab in y_rows: + yi = int(np.clip(y_idx, 0, spec.ny - 1)) + line_payload[f"ux_{y_lab}"] = ux[yi, :].astype(np.float32) + line_payload[f"rho_{y_lab}"] = rho[yi, :].astype(np.float32) + np.savez_compressed(os.path.join(line_dir, "lines.npz"), **line_payload) + + status = "last_stable" if used_last_stable else "final" + print( + f" -> {status} step {out_step} rho=[{meta['rho_min']:.4f},{meta['rho_max']:.4f}] " + f"ux_max={meta['ux_max']:.4f} force_bad={force_bad_step}", + flush=True, + ) + return {**meta, "field_pngs": pngs, "run_dir": run_dir} + + +def main() -> int: + ap = argparse.ArgumentParser(description="Three-scenario inlet field export (5000 steps)") + ap.add_argument( + "--scenario", + choices=( + "empty_channel", + "empty_channel_inlet_matrix", + "kan99b", + "sah04_case9", + "all", + ), + default="all", + ) + ap.add_argument( + "--collision", + default="", + help="Optional filter: SRT or MRT only", + ) + ap.add_argument("--steps", type=int, default=5000) + ap.add_argument( + "--out-dir", + default=os.path.join(_REPO, "tests", "output", "inlet_scenario_fields"), + ) + args = ap.parse_args() + + if not os.path.isfile(_DEFAULT_LBM): + print(f"Missing {_DEFAULT_LBM}", file=sys.stderr) + return 2 + + base = _load_json(_DEFAULT_LBM) + if args.scenario == "empty_channel_inlet_matrix": + specs = _empty_channel_inlet_matrix_specs() + elif args.scenario == "all": + specs = _all_specs() + _empty_channel_inlet_matrix_specs() + else: + specs = _all_specs() + if args.scenario != "all": + specs = [s for s in specs if s.scenario == args.scenario] + if args.collision.strip(): + coll = args.collision.strip().upper() + specs = [s for s in specs if s.collision.upper() == coll] + specs = [replace(s, steps=int(args.steps)) for s in specs] + + out_dir = os.path.abspath(args.out_dir) + os.makedirs(out_dir, exist_ok=True) + + rows: List[Dict[str, Any]] = [] + for spec in specs: + try: + row = run_one(spec, base, out_dir) + rows.append(row) + except Exception as e: # noqa: BLE001 + print(f" FAILED {spec.scenario}/{spec.run_id}: {e}", flush=True) + rows.append( + { + "scenario": spec.scenario, + "run_id": spec.run_id, + "label": spec.label, + "error": str(e), + } + ) + + summary_path = os.path.join(out_dir, "summary.csv") + if rows: + keys: List[str] = [] + for r in rows: + for k in r: + if k not in keys and k != "field_pngs": + keys.append(k) + with open(summary_path, "w", encoding="utf-8", newline="") as f: + w = csv.DictWriter(f, fieldnames=keys, extrasaction="ignore") + w.writeheader() + w.writerows(rows) + + _write_json( + os.path.join(out_dir, "manifest.json"), + {"steps": args.steps, "scenario_filter": args.scenario, "runs": [s.run_id for s in specs]}, + ) + print(f"Wrote: {summary_path}", flush=True) + print(f"Output: {out_dir}", flush=True) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/tests/run_sah04_st_matrix.py b/tests/run_sah04_st_matrix.py index 0266128..83675b9 100644 --- a/tests/run_sah04_st_matrix.py +++ b/tests/run_sah04_st_matrix.py @@ -46,7 +46,7 @@ from typing import Any, Dict, List, Optional, Sequence, Tuple import numpy as np import pycuda.driver as cuda -_PKG_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) +_PKG_ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..")) _DEFAULT_LBM = os.path.join(_PKG_ROOT, "src", "CelerisLab", "configs", "config_lbm.json") # D=30 fixed; Lx_fluid = 80D per Sah04 confined setup