From 926ca27407c26ca11283282a6e53249925725376 Mon Sep 17 00:00:00 2001 From: TRON Date: Sun, 18 Oct 2026 09:16:25 -0700 Subject: [PATCH] experiment: isolate full-state regularized west velocity inlet --- src/CelerisLab/config.py | 6 ++++ .../lbm/kernels/boundary/inlet/common.cuh | 2 +- .../boundary/inlet/experimental_full9.cuh | 36 +++++++++++++++++++ .../lbm/kernels/boundary/inlet_outlet.cuh | 5 +++ 4 files changed, 48 insertions(+), 1 deletion(-) create mode 100644 src/CelerisLab/lbm/kernels/boundary/inlet/experimental_full9.cuh diff --git a/src/CelerisLab/config.py b/src/CelerisLab/config.py index bff1c45..7b79f04 100644 --- a/src/CelerisLab/config.py +++ b/src/CelerisLab/config.py @@ -32,6 +32,7 @@ INLET_SCHEME_MAP = { "channel_stabilized": 1, "equilibrium": 2, "regularized": 3, + "experimental_regularized_full9": 4, } OUTLET_MAP = {"neq_extrap": 0, "zero_gradient": 1, "blended": 2} Y_WALL_BC_MAP = {"bounce_back": 0, "free_slip": 1, "on_site_velocity": 2} @@ -162,6 +163,11 @@ class LBMConfig: raise ValueError("wall_slot_width requires on_site_velocity and velocity=0") if self.top_wall_bc == "reservoir" and (self.y_wall_bc != "on_site_velocity" or self.velocity != 0): raise ValueError("reservoir requires on_site_velocity and velocity=0") + if self.inlet_scheme == "experimental_regularized_full9": + if (not self.is_d2q9 or self.collision != "SRT" + or self.streaming != "double_buffer" or self.store_precision != "FP32" + or self.les_enabled or self.y_wall_bc != "bounce_back"): + raise ValueError("experimental_regularized_full9 requires D2Q9/SRT/FP32/double_buffer, bounce_back walls, no LES") if self.y_wall_bc == "on_site_velocity": if (not self.is_d2q9 or self.collision != "SRT" or self.streaming != "double_buffer" or self.store_precision != "FP32" diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh index fab9a69..38f6d33 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/common.cuh @@ -74,7 +74,7 @@ __device__ __forceinline__ float west_velocity_rho_closure_d3q19( __device__ __forceinline__ bool inlet_scheme_uses_post_collision_ghost() { -#if INLET_SCHEME == 0 +#if INLET_SCHEME == 0 || INLET_SCHEME == 4 return true; #else return false; diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet/experimental_full9.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet/experimental_full9.cuh new file mode 100644 index 0000000..54a5d82 --- /dev/null +++ b/src/CelerisLab/lbm/kernels/boundary/inlet/experimental_full9.cuh @@ -0,0 +1,36 @@ +// Experimental D2Q9 west velocity boundary: full-state stress regularization. +// Palabos boundaryTemplates::compute_PiNeq / regularized velocity completion. +// This returns PRE-collision populations. The existing step kernel collides once. +#ifndef CELERIS_EXPERIMENTAL_FULL9_INLET_CUH +#define CELERIS_EXPERIMENTAL_FULL9_INLET_CUH + +#if LATTICE_MODEL == LATTICE_D2Q9 +__device__ inline void apply_experimental_full9_inlet_d2q9( + float* __restrict__ f, float ux_target, float uy_target) +{ + const float rho = west_velocity_rho_closure_d2q9(f, ux_target); + float eq[9]; + compute_feq(rho, ux_target, uy_target, eq); + // Only known pulled directions are read; unknown wrapped q1/5/7 are ignored. + // f and eq use the same shifting convention, so weights cancel in K. + const float k2 = f[2] - eq[2]; + const float k3 = f[3] - eq[3]; + const float k4 = f[4] - eq[4]; + const float k6 = f[6] - eq[6]; + const float k8 = f[8] - eq[8]; + const float pxx = 2.0f * (k2 + k6 + k8); + const float pyy = k3 + k4 + 2.0f * (k6 + k8); + const float pxy = 2.0f * (k6 - k8); + const int cx[9] = {0, 1,-1, 0, 0, 1,-1, 1,-1}; + const int cy[9] = {0, 0, 0, 1,-1, 1,-1,-1, 1}; + for (int q = 0; q < 9; ++q) { + const float weight = q == 0 ? W0_VAL : (q < 5 ? WS_VAL : WE_VAL); + const float stress = weight * ( + 4.5f * (((float)(cx[q]*cx[q]) - 1.0f/3.0f) * pxx + + ((float)(cy[q]*cy[q]) - 1.0f/3.0f) * pyy) + + 9.0f * (float)(cx[q]*cy[q]) * pxy); + f[q] = eq[q] + stress; + } +} +#endif +#endif diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh index 0add5ba..906899d 100644 --- a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh +++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh @@ -24,6 +24,7 @@ #include "inlet/channel_stabilized.cuh" #include "inlet/equilibrium.cuh" #include "inlet/regularized.cuh" +#include "inlet/experimental_full9.cuh" #include "outlet/pressure_neq.cuh" #if DIM == 2 @@ -45,6 +46,8 @@ __device__ __forceinline__ void apply_inlet_pull_d2q9( #if INLET_SCHEME == 0 apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target); +#elif INLET_SCHEME == 4 + apply_experimental_full9_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); @@ -143,6 +146,8 @@ __device__ __forceinline__ void apply_inlet_esopull_d2q9( #if INLET_SCHEME == 0 apply_zou_he_left_velocity_inlet_d2q9(f, u_target, v_target); +#elif INLET_SCHEME == 4 + apply_experimental_full9_inlet_d2q9(f, u_target, v_target); #elif INLET_SCHEME == 1 || INLET_SCHEME == 3 #if INLET_SCHEME == 1 apply_channel_stabilized_inlet_d2q9(f, f_neb, (float)y);