diff --git a/.gitignore b/.gitignore
index ee1fdd0..2982506 100644
--- a/.gitignore
+++ b/.gitignore
@@ -68,4 +68,5 @@ venv.bak/
 *.log
 
 # reference:
-ref/
\ No newline at end of file
+ref/
+output/
\ No newline at end of file
diff --git a/output/poiseuille_d2q9.png b/output/poiseuille_d2q9.png
index 752f3a1..36c5a39 100644
Binary files a/output/poiseuille_d2q9.png and b/output/poiseuille_d2q9.png differ
diff --git a/src/CelerisLab/cuda/compiler.py b/src/CelerisLab/cuda/compiler.py
index 2164b71..18eb528 100644
--- a/src/CelerisLab/cuda/compiler.py
+++ b/src/CelerisLab/cuda/compiler.py
@@ -95,7 +95,9 @@ def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig,
                      collision_model: int = 2,
                      streaming_model: int = 0,
                      store_precision: int = 0,
-                     use_ddf_shifting: int = 0):
+                     use_ddf_shifting: int = 0,
+                     use_les: int = 0,
+                     les_cs: float = 0.16):
     """Configure macros.h for the new modular kernel architecture.
 
     Args:
@@ -103,6 +105,8 @@ def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig,
         streaming_model:  0=double-buffer (default), 1=Esoteric-Pull
         store_precision:  0=FP32 (default), 1=FP16S, 2=FP16C
         use_ddf_shifting: 0=off (default), 1=on
+        use_les:          0=off (default), 1=Smagorinsky LES
+        les_cs:           Smagorinsky constant C_s
     """
     # First apply legacy config
     config_kernal(config_cuda, config_field)
@@ -113,6 +117,8 @@ def config_kernal_v2(config_cuda: CudaConfig, config_field: FlowFieldConfig,
     lines = modify_macro(lines, "STREAMING_MODEL", streaming_model)
     lines = modify_macro(lines, "STORE_PRECISION", store_precision)
     lines = modify_macro(lines, "USE_DDF_SHIFTING", use_ddf_shifting)
+    lines = modify_macro(lines, "USE_LES", use_les)
+    lines = modify_macro(lines, "LES_CS", f"{les_cs:.6f}f")
     write_lines(kernel_path("macros.h"), lines)
 
 def config_object(n_obj: int):
diff --git a/src/CelerisLab/lbm/driver.py b/src/CelerisLab/lbm/driver.py
index 94649f2..d33bd17 100644
--- a/src/CelerisLab/lbm/driver.py
+++ b/src/CelerisLab/lbm/driver.py
@@ -23,6 +23,13 @@ class FlowField:
         field_config: utils.FlowFieldConfig,
         cuda_config: utils.CudaConfig,
         device_id: Union[int, List[int]] = None,
+        use_kernel_v2: bool = True,
+        collision_model: int = 0,
+        streaming_model: int = 0,
+        store_precision: int = 0,
+        use_ddf_shifting: int = 0,
+        use_les: int = 0,
+        les_cs: float = 0.16,
     ):
         self.field_config = field_config
         self.cuda_config = cuda_config
@@ -49,10 +56,26 @@ class FlowField:
 
         utils.check_cuda_capability(field_config, cuda_config, device_id)
 
-        # Config kernel
-        compiler.config_kernal(cuda_config, field_config)
-        compiler.config_object(int(0))
-        # compiler.config_sensor(int(0))
+        self.use_kernel_v2 = bool(use_kernel_v2)
+        self.collision_model = int(collision_model)
+        self.streaming_model = int(streaming_model)
+        self.store_precision = int(store_precision)
+        self.use_ddf_shifting = int(use_ddf_shifting)
+        self.use_les = int(use_les)
+        self.les_cs = float(les_cs)
+
+        if self.collision_model not in (0, 1, 2):
+            raise ValueError("collision_model must be 0(SRT), 1(TRT), or 2(MRT).")
+        if self.streaming_model not in (0, 1):
+            raise ValueError("streaming_model must be 0(double-buffer) or 1(esopull).")
+        if self.store_precision not in (0, 1, 2):
+            raise ValueError("store_precision must be 0(FP32), 1(FP16S), or 2(FP16C).")
+        if self.use_ddf_shifting not in (0, 1):
+            raise ValueError("use_ddf_shifting must be 0 or 1.")
+        if self.use_les not in (0, 1):
+            raise ValueError("use_les must be 0 or 1.")
+        if not (0.0 < self.les_cs < 1.0):
+            raise ValueError("les_cs must be in (0, 1).")
 
         # Set constants
         if field_config.data_type == "FP32":
@@ -84,11 +107,10 @@ class FlowField:
                 f"Unsupported lattice type {field_config.lattice} in {field_config.dimensionality} dimensions."
             )
 
-        # Compile kernel
-        compiler.compile_kernel()
-        self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
-        self.step = self.ptx.get_function("OneStep")
-        initflow = self.ptx.get_function("InitTubeFlow")
+        self.objects = {}
+
+        # Compile and load kernel
+        self._rebuild_kernel()
 
         # Initialize memory
         self.ddf = np.zeros(self.FIELD_SIZE * self.LATTICE, dtype=self.DATA_TYPE)
@@ -105,11 +127,10 @@ class FlowField:
         self.delta_gpu = cuda.mem_alloc(1)
         self.vortex_gpu = cuda.mem_alloc(self.vortex_config.nbytes)
 
-        self.objects = {}
         self.action = np.zeros(0, dtype=self.DATA_TYPE)
         self.obs = np.zeros(0, dtype=self.DATA_TYPE)
 
-        initflow(
+        self.initflow(
             self.flag_gpu,
             self.ddf_gpu,
             block=(self.cuda_config.threads_per_block, 1, 1),
@@ -122,6 +143,38 @@ class FlowField:
         cuda.memcpy_dtoh(self.flag, self.flag_gpu)
         cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
 
+    def _configure_kernel(self):
+        if self.use_kernel_v2:
+            compiler.config_kernal_v2(
+                self.cuda_config,
+                self.field_config,
+                collision_model=self.collision_model,
+                streaming_model=self.streaming_model,
+                store_precision=self.store_precision,
+                use_ddf_shifting=self.use_ddf_shifting,
+                use_les=self.use_les,
+                les_cs=self.les_cs,
+            )
+        else:
+            compiler.config_kernal(self.cuda_config, self.field_config)
+
+    def _compile_and_load_kernel(self):
+        if self.use_kernel_v2:
+            compiler.compile_kernel_v2()
+            self.ptx = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
+            self.step = self.ptx.get_function("OneStep")
+            self.initflow = self.ptx.get_function("InitTubeFlow_v2")
+        else:
+            compiler.compile_kernel()
+            self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
+            self.step = self.ptx.get_function("OneStep")
+            self.initflow = self.ptx.get_function("InitTubeFlow")
+
+    def _rebuild_kernel(self):
+        self._configure_kernel()
+        compiler.config_object(len(self.objects))
+        self._compile_and_load_kernel()
+
     def add_cylinder(self, center: Tuple[float, float, float], radius: float, id_obj: Optional[int] = None):
         x_c, y_c, z_c = center
 
@@ -193,10 +246,7 @@ class FlowField:
         cuda.memcpy_htod(self.flag_gpu, self.flag)
         cuda.memcpy_htod(self.indx_gpu, self.indx)
 
-        compiler.config_object(len(self.objects))
-        compiler.compile_kernel()
-        self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
-        self.step = self.ptx.get_function("OneStep")
+        self._rebuild_kernel()
 
     def add_sensor(self, center: Tuple[float, float, float], radius: float):
         x_c, y_c, z_c = center
@@ -235,10 +285,7 @@ class FlowField:
         cuda.memcpy_htod(self.flag_gpu, self.flag)
         cuda.memcpy_htod(self.indx_gpu, self.indx)
 
-        compiler.config_object(len(self.objects))
-        compiler.compile_kernel()
-        self.ptx = cuda.module_from_file(compiler.kernel_path("kernel.ptx"))
-        self.step = self.ptx.get_function("OneStep")
+        self._rebuild_kernel()
 
     def add_vortex(self, center: Tuple[float, float, float], radius: float, strength: float, direction: float, type: str):
         x_c, y_c, z_c = center
diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
index 3a87a5a..1fde9ca 100644
--- a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
+++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
@@ -15,6 +15,31 @@
 #ifndef CELERIS_BOUNDARY_INLET_OUTLET_CUH
 #define CELERIS_BOUNDARY_INLET_OUTLET_CUH
 
+#ifndef INLET_PROFILE
+#define INLET_PROFILE 1
+#endif
+
+#ifndef OUTLET_MODE
+#define OUTLET_MODE 0
+#endif
+
+#ifndef OUTLET_BACKFLOW_CLAMP
+#define OUTLET_BACKFLOW_CLAMP 1
+#endif
+
+#ifndef OUTLET_BLEND_ALPHA
+#define OUTLET_BLEND_ALPHA 0.70f
+#endif
+
+__device__ __forceinline__ float inlet_target_u(float y_coord) {
+#if INLET_PROFILE == 0
+    return U0;
+#else
+    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
+    return U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+#endif
+}
+
 #if NQ == 9
 
 // ---------------------------------------------------------------------------
@@ -33,37 +58,21 @@ __device__ inline void apply_parabolic_inlet(float* __restrict__ f,
                                              float y_coord)
 {
     // Neighbor macros
-    float p_neb = (f_neb[0]+f_neb[1]+f_neb[2]+f_neb[3]+f_neb[4]
-                  +f_neb[5]+f_neb[6]+f_neb[7]+f_neb[8]) / 3.0f;
+    float rho_neb, u_neb, v_neb;
+    compute_rho_u(f_neb, rho_neb, u_neb, v_neb);
 
     // Target velocity (parabolic profile)
-    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
-    float u_target = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+    float u_target = inlet_target_u(y_coord);
     float v_target = 0.0f;
 
-    // Neighbor velocity
-    float u_neb = (f_neb[1]-f_neb[2]+f_neb[5]-f_neb[6]+f_neb[7]-f_neb[8]) / RHO;
-    float v_neb = (f_neb[3]-f_neb[4]+f_neb[5]-f_neb[6]-f_neb[7]+f_neb[8]) / RHO;
-
-    // feq for direction i=1 (cx=1, cy=0), w=1/9:
-    //   feq = (2p + RHO*(2u² + 2u - v²)) / 6
-    float feq1_target = (2.0f*p_neb + RHO*(2.0f*u_target*u_target + 2.0f*u_target - v_target*v_target)) / 6.0f;
-    float feq1_neb    = (2.0f*p_neb + RHO*(2.0f*u_neb*u_neb + 2.0f*u_neb - v_neb*v_neb)) / 6.0f;
-
-    // feq for direction i=5 (cx=1, cy=1), w=1/36:
-    //   feq = (p + RHO*(u² + 3uv + u + v² + v)) / 12
-    float feq5_target = (p_neb + RHO*(u_target*u_target + 3.0f*u_target*v_target + u_target + v_target*v_target + v_target)) / 12.0f;
-    float feq5_neb    = (p_neb + RHO*(u_neb*u_neb + 3.0f*u_neb*v_neb + u_neb + v_neb*v_neb + v_neb)) / 12.0f;
-
-    // feq for direction i=7 (cx=1, cy=-1), w=1/36:
-    //   feq = (p + RHO*(u² - 3uv + u + v² - v)) / 12
-    float feq7_target = (p_neb + RHO*(u_target*u_target - 3.0f*u_target*v_target + u_target + v_target*v_target - v_target)) / 12.0f;
-    float feq7_neb    = (p_neb + RHO*(u_neb*u_neb - 3.0f*u_neb*v_neb + u_neb + v_neb*v_neb - v_neb)) / 12.0f;
+    float feq_tar[9], feq_neb[9];
+    compute_feq(rho_neb, u_target, v_target, feq_tar);
+    compute_feq(rho_neb, u_neb, v_neb, feq_neb);
 
     // Non-equilibrium extrapolation
-    f[1] = f_neb[1] - feq1_neb + feq1_target;
-    f[5] = f_neb[5] - feq5_neb + feq5_target;
-    f[7] = f_neb[7] - feq7_neb + feq7_target;
+    f[1] = f_neb[1] - feq_neb[1] + feq_tar[1];
+    f[5] = f_neb[5] - feq_neb[5] + feq_tar[5];
+    f[7] = f_neb[7] - feq_neb[7] + feq_tar[7];
 }
 
 // ---------------------------------------------------------------------------
@@ -76,35 +85,39 @@ __device__ inline void apply_pressure_outlet(float* __restrict__ f,
                                              const float* __restrict__ f_neb,
                                              float y_coord)
 {
-    float p_out = 0.0f;
+    (void)y_coord;
 
-    // Target velocity (parabolic, same as inlet for consistency)
-    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
-    float u_target = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
-    float v_target = 0.0f;
+#if OUTLET_MODE == 1
+    // Simple zero-gradient copy for unknown incoming directions at outlet.
+    f[2] = f_neb[2];
+    f[8] = f_neb[8];
+    f[6] = f_neb[6];
+#else
 
-    // Neighbor velocity
-    float u_neb = (f_neb[1]-f_neb[2]+f_neb[5]-f_neb[6]+f_neb[7]-f_neb[8]) / RHO;
-    float v_neb = (f_neb[3]-f_neb[4]+f_neb[5]-f_neb[6]-f_neb[7]+f_neb[8]) / RHO;
+    // Convective non-reflecting style: extrapolate outlet density from neighbor
+    // and keep neighbor velocity.
+    float rho_neb, u_neb, v_neb;
+    compute_rho_u(f_neb, rho_neb, u_neb, v_neb);
+#if OUTLET_BACKFLOW_CLAMP
+    u_neb = fmaxf(u_neb, 0.0f);
+#endif
+    float rho_out = rho_neb;
 
-    // feq for direction i=2 (cx=-1, cy=0), w=1/9:
-    //   feq = (2p - RHO*(-2u² + 2u + v²)) / 6
-    float feq2_target = (2.0f*p_out - RHO*(-2.0f*u_target*u_target + 2.0f*u_target + v_target*v_target)) / 6.0f;
-    float feq2_neb    = (2.0f*p_out - RHO*(-2.0f*u_neb*u_neb + 2.0f*u_neb + v_neb*v_neb)) / 6.0f;
+    float feq_tar[9], feq_neb[9];
+    compute_feq(rho_out, u_neb, v_neb, feq_tar);
+    compute_feq(rho_neb, u_neb, v_neb, feq_neb);
 
-    // feq for direction i=8 (cx=-1, cy=1), w=1/36:
-    //   feq = (p + RHO*(u² - 3uv - u + v² + v)) / 12
-    float feq8_target = (p_out + RHO*(u_target*u_target - 3.0f*u_target*v_target - u_target + v_target*v_target + v_target)) / 12.0f;
-    float feq8_neb    = (p_out + RHO*(u_neb*u_neb - 3.0f*u_neb*v_neb - u_neb + v_neb*v_neb + v_neb)) / 12.0f;
-
-    // feq for direction i=6 (cx=-1, cy=-1), w=1/36:
-    //   feq = (p + RHO*(u² + 3uv - u + v² - v)) / 12
-    float feq6_target = (p_out + RHO*(u_target*u_target + 3.0f*u_target*v_target - u_target + v_target*v_target - v_target)) / 12.0f;
-    float feq6_neb    = (p_out + RHO*(u_neb*u_neb + 3.0f*u_neb*v_neb - u_neb + v_neb*v_neb - v_neb)) / 12.0f;
-
-    f[2] = f_neb[2] - feq2_neb + feq2_target;
-    f[8] = f_neb[8] - feq8_neb + feq8_target;
-    f[6] = f_neb[6] - feq6_neb + feq6_target;
+#if OUTLET_MODE == 2
+    const float a = OUTLET_BLEND_ALPHA;
+    f[2] = a * (f_neb[2] - feq_neb[2] + feq_tar[2]) + (1.0f - a) * f_neb[2];
+    f[8] = a * (f_neb[8] - feq_neb[8] + feq_tar[8]) + (1.0f - a) * f_neb[8];
+    f[6] = a * (f_neb[6] - feq_neb[6] + feq_tar[6]) + (1.0f - a) * f_neb[6];
+#else
+    f[2] = f_neb[2] - feq_neb[2] + feq_tar[2];
+    f[8] = f_neb[8] - feq_neb[8] + feq_tar[8];
+    f[6] = f_neb[6] - feq_neb[6] + feq_tar[6];
+#endif
+#endif
 }
 
 #endif // NQ == 9
@@ -128,8 +141,7 @@ __device__ inline void apply_parabolic_inlet_3d(float* __restrict__ f,
     compute_rho_u(f_neb, rho_neb, un, vn, wn);
 
     // Target velocity (parabolic in y, uniform in z)
-    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
-    float u_tar = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
+    float u_tar = inlet_target_u(y_coord);
 
     // feq arrays
     float feq_tar[19], feq_neb[19];
@@ -148,24 +160,47 @@ __device__ inline void apply_pressure_outlet_3d(float* __restrict__ f,
                                                 const float* __restrict__ f_neb,
                                                 float y_coord)
 {
+    (void)y_coord;
+
+#if OUTLET_MODE == 1
+    // Simple zero-gradient copy for unknown incoming directions at outlet.
+    f[2]  = f_neb[2];
+    f[8]  = f_neb[8];
+    f[10] = f_neb[10];
+    f[14] = f_neb[14];
+    f[16] = f_neb[16];
+#else
+
     // Neighbor macros
     float rho_neb, un, vn, wn;
     compute_rho_u(f_neb, rho_neb, un, vn, wn);
+#if OUTLET_BACKFLOW_CLAMP
+    un = fmaxf(un, 0.0f);
+#endif
 
-    // Target: p_out = 0 gauge → rho = 1.0  (using neighbor velocity)
-    float yy = (y_coord - 0.5f * (NY - 1)) / (NY - 2.0f);
-    float u_tar = U0 * 1.5f * (1.0f - 4.0f * yy * yy);
-
+    // Convective non-reflecting style: extrapolate outlet density from neighbor
+    // and keep neighbor velocity.
+    float rho_out = rho_neb;
     float feq_tar[19], feq_neb[19];
-    compute_feq(RHO, u_tar, 0.0f, 0.0f, feq_tar);
-    compute_feq(RHO, un, vn, wn, feq_neb);
+    compute_feq(rho_out, un, vn, wn, feq_tar);
+    compute_feq(rho_neb, un, vn, wn, feq_neb);
 
     // Reconstruct cx<0 directions: i = 2, 8, 10, 14, 16
+#if OUTLET_MODE == 2
+    const float a = OUTLET_BLEND_ALPHA;
+    f[2]  = a * (f_neb[2]  - feq_neb[2]  + feq_tar[2])  + (1.0f - a) * f_neb[2];
+    f[8]  = a * (f_neb[8]  - feq_neb[8]  + feq_tar[8])  + (1.0f - a) * f_neb[8];
+    f[10] = a * (f_neb[10] - feq_neb[10] + feq_tar[10]) + (1.0f - a) * f_neb[10];
+    f[14] = a * (f_neb[14] - feq_neb[14] + feq_tar[14]) + (1.0f - a) * f_neb[14];
+    f[16] = a * (f_neb[16] - feq_neb[16] + feq_tar[16]) + (1.0f - a) * f_neb[16];
+#else
     f[2]  = f_neb[2]  - feq_neb[2]  + feq_tar[2];
     f[8]  = f_neb[8]  - feq_neb[8]  + feq_tar[8];
     f[10] = f_neb[10] - feq_neb[10] + feq_tar[10];
     f[14] = f_neb[14] - feq_neb[14] + feq_tar[14];
     f[16] = f_neb[16] - feq_neb[16] + feq_tar[16];
+#endif
+#endif
 }
 
 #endif // NQ == 19
diff --git a/src/CelerisLab/lbm/kernels/kernel_v2.cu b/src/CelerisLab/lbm/kernels/kernel_v2.cu
index 199eefb..4c041ab 100644
--- a/src/CelerisLab/lbm/kernels/kernel_v2.cu
+++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu
@@ -39,6 +39,7 @@
 #include "operators/collision_trt.cuh"
 #include "operators/collision_mrt.cuh"
 #include "operators/forcing_guo.cuh"
+#include "operators/turbulence_smag.cuh"
 
 // ---------------------------------------------------------------------------
 // Layer 3: Streaming
@@ -208,18 +209,31 @@ __global__ void OneStep(
         }
     }
 
+    // Obstacle nodes: full-way bounce-back (same swap as wall nodes)
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+        }
+    }
+
     // Collision
     if (fl & LEGACY_FLUID) {
         float feq[NQ], Fin[NQ];
         compute_feq(rho_n, ux, uy, feq);
         zero_forcing(Fin);
+        float omega_col = d_params.omega;
+    #if USE_LES
+        omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+    #endif
+        omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
 
 #if COLLISION_MODEL == 0
-        collide_srt(f, feq, Fin, d_params.omega);
+        collide_srt(f, feq, Fin, omega_col);
 #elif COLLISION_MODEL == 1
-        collide_trt(f, feq, Fin, d_params.omega);
+        collide_trt(f, feq, Fin, omega_col);
 #elif COLLISION_MODEL == 2
-        collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
+        collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
 #endif
     }
 
@@ -289,7 +303,18 @@ __global__ void OneStep(
         float feq[NQ], Fin[NQ];
         compute_feq(rho_n, ux, uy, uz, feq);
         zero_forcing(Fin);
-        collide_srt(f, feq, Fin, d_params.omega);
+        float omega_col = d_params.omega;
+#if USE_LES
+        omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+#endif
+        omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
+    #if COLLISION_MODEL == 0
+        collide_srt(f, feq, Fin, omega_col);
+    #elif COLLISION_MODEL == 1
+        collide_trt(f, feq, Fin, omega_col);
+    #elif COLLISION_MODEL == 2
+        collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
+    #endif
     }
 
     stream_pull_store(k, f, fi_out);
diff --git a/src/CelerisLab/lbm/kernels/macros.h b/src/CelerisLab/lbm/kernels/macros.h
index dda58fc..d784f34 100644
--- a/src/CelerisLab/lbm/kernels/macros.h
+++ b/src/CelerisLab/lbm/kernels/macros.h
@@ -60,4 +60,44 @@
 // DDF-shifting: 0=off, 1=on
 #ifndef USE_DDF_SHIFTING
 #define USE_DDF_SHIFTING 0
+#endif
+
+// LES model: 0=off, 1=Smagorinsky
+#ifndef USE_LES
+#define USE_LES 0
+#endif
+
+// Smagorinsky constant C_s
+#ifndef LES_CS
+#define LES_CS 0.16f
+#endif
+
+// Inlet profile: 1=parabolic (channel), 0=uniform (external flow)
+#ifndef INLET_PROFILE
+#define INLET_PROFILE 1
+#endif
+
+// Outlet mode: 0=non-equilibrium extrapolation, 1=zero-gradient copy (more dissipative)
+#ifndef OUTLET_MODE
+#define OUTLET_MODE 0
+#endif
+
+// Outlet blend factor for damped outlet mode (OUTLET_MODE=2):
+// f_out = a*(non-eq extrapolation) + (1-a)*(zero-gradient copy)
+#ifndef OUTLET_BLEND_ALPHA
+#define OUTLET_BLEND_ALPHA 0.70f
+#endif
+
+// Outlet backflow clamp: 0=off, 1=force non-negative streamwise velocity at outlet target
+#ifndef OUTLET_BACKFLOW_CLAMP
+#define OUTLET_BACKFLOW_CLAMP 1
+#endif
+
+// Global collision omega guardrails
+#ifndef OMEGA_COLLISION_MIN
+#define OMEGA_COLLISION_MIN 0.01f
+#endif
+
+#ifndef OMEGA_COLLISION_MAX
+#define OMEGA_COLLISION_MAX 1.999f
 #endif
\ No newline at end of file
diff --git a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh
index 5aab699..e8e3d80 100644
--- a/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh
+++ b/src/CelerisLab/lbm/kernels/operators/collision_mrt.cuh
@@ -31,19 +31,19 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
                                             float omega)
 {
     // ----- Relaxation rates -----
-    //   s_rho = s_jx = s_jy = 1.0  (conserved moments → overwrite to eq)
+    //   s_rho = s_jx = s_jy = 0.0  (strictly conserved in collision)
     //   s_e   = s_eps = s_q  = 1.2  (kinetic transport)
     //   s_nu  = omega               (viscosity-related)
-    const float s_rho = 1.0f;   // conserved moment relaxation
+    const float s_rho = 0.0f;   // conserved moments
     const float s_e   = 1.2f;
     const float s_eps = 1.2f;
-    const float s_jx  = 1.0f;
+    const float s_jx  = 0.0f;
     const float s_q   = 1.2f;
-    const float s_jy  = 1.0f;
+    const float s_jy  = 0.0f;
     const float s_nu  = omega;
 
-    // ----- Pressure (from density) -----
-    const float p = (g[0]+g[1]+g[2]+g[3]+g[4]+g[5]+g[6]+g[7]+g[8]) / 3.0f;
+    // ----- Pressure from local density -----
+    const float p = rho * (1.0f / 3.0f);
 
     // ----- Forward transform: m = M * g  (new paired ordering) -----
     float m[9];
@@ -61,14 +61,14 @@ __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
     const float u2 = ux * ux + uy * uy;
     float meq[9];
     meq[0] = 3.0f * p;                        // ρ
-    meq[1] = -6.0f * p + 3.0f * RHO * u2;     // e
-    meq[2] =  3.0f * p - 3.0f * RHO * u2;     // ε
-    meq[3] = RHO * ux;                         // jx
-    meq[4] = -RHO * ux;                        // qx
-    meq[5] = RHO * uy;                         // jy
-    meq[6] = -RHO * uy;                        // qy
-    meq[7] = RHO * (ux * ux - uy * uy);        // pxx
-    meq[8] = RHO * ux * uy;                     // pxy
+    meq[1] = -6.0f * p + 3.0f * rho * u2;     // e
+    meq[2] =  3.0f * p - 3.0f * rho * u2;     // ε
+    meq[3] = rho * ux;                         // jx
+    meq[4] = -rho * ux;                        // qx
+    meq[5] = rho * uy;                         // jy
+    meq[6] = -rho * uy;                        // qy
+    meq[7] = rho * (ux * ux - uy * uy);        // pxx
+    meq[8] = rho * ux * uy;                     // pxy
 
     // ----- Relaxation:  delta_m[i] = s_i * (meq[i] - m[i]) -----
     float dm[9];
@@ -112,21 +112,106 @@ __device__ __forceinline__ void collide_mrt_no_force(float* __restrict__ g,
 #elif NQ == 19
 
 // ---------------------------------------------------------------------------
-// D3Q19 MRT  (placeholder – use SRT or TRT until fully validated)
+// D3Q19 MRT (tensor-projected multi-mode relaxation)
+//
+// This implementation performs a true multi-relaxation update by splitting
+// non-equilibrium content into three subspaces:
+//   1) deviatoric 2nd-order stress modes  -> s_nu   (shear, controlled by omega)
+//   2) isotropic 2nd-order stress mode    -> s_bulk (bulk, controlled by omega_bulk)
+//   3) higher-order residual modes        -> s_high (kinetic damping)
+//
+// It is formulated in Hermite/tensor space, avoiding explicit 19x19 matrices
+// while still relaxing different mode families at different rates.
 // ---------------------------------------------------------------------------
 __device__ __forceinline__ void collide_mrt(float* __restrict__ g,
                                             float rho, float ux, float uy, float uz,
                                             const float* __restrict__ Fin,
                                             float omega)
 {
-    // TODO: implement D3Q19 MRT with 19-moment Gram–Schmidt basis
-    // Fall back to SRT for now
     float feq[19];
     compute_feq(rho, ux, uy, uz, feq);
-    const float c_tau = 1.0f - 0.5f * omega;
+
+    // Relaxation rates
+    const float s_nu   = omega;
+    const float s_bulk = (d_params.omega_bulk > 0.0f) ? d_params.omega_bulk : 1.2f;
+    const float s_high = 1.4f;
+
+    const float one_minus_s_nu   = 1.0f - s_nu;
+    const float one_minus_s_bulk = 1.0f - s_bulk;
+    const float one_minus_s_high = 1.0f - s_high;
+
+    // Non-equilibrium populations
+    float neq[19];
     #pragma unroll
     for (int i = 0; i < 19; i++) {
-        g[i] = fmaf(1.0f - omega, g[i], fmaf(omega, feq[i], c_tau * Fin[i]));
+        neq[i] = g[i] - feq[i];
+    }
+
+    // Build non-equilibrium 2nd-order tensor Πneq = Σ (ci ci) (fi-feqi)
+    float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
+    float pixy = 0.0f, pixz = 0.0f, piyz = 0.0f;
+
+    #pragma unroll
+    for (int i = 0; i < 19; i++) {
+        const float ci_x = (float)d_cx[i];
+        const float ci_y = (float)d_cy[i];
+        const float ci_z = (float)d_cz[i];
+        const float fneq = neq[i];
+        pixx += fneq * ci_x * ci_x;
+        piyy += fneq * ci_y * ci_y;
+        pizz += fneq * ci_z * ci_z;
+        pixy += fneq * ci_x * ci_y;
+        pixz += fneq * ci_x * ci_z;
+        piyz += fneq * ci_y * ci_z;
+    }
+
+    // Isotropic / deviatoric split of Πneq
+    const float tr_pi      = pixx + piyy + pizz;
+    const float tr_pi_thrd = tr_pi * (1.0f / 3.0f);
+
+    const float dev_xx = pixx - tr_pi_thrd;
+    const float dev_yy = piyy - tr_pi_thrd;
+    const float dev_zz = pizz - tr_pi_thrd;
+    const float dev_xy = pixy;
+    const float dev_xz = pixz;
+    const float dev_yz = piyz;
+
+    // cs^2 = 1/3 => 2*cs^4 = 2/9, inverse = 4.5
+    // Project neq onto second-order Hermite basis and relax mode families.
+    const float proj_pref = 4.5f;
+    const float c_tau = 1.0f - 0.5f * omega;
+
+    #pragma unroll
+    for (int i = 0; i < 19; i++) {
+        const float ci_x = (float)d_cx[i];
+        const float ci_y = (float)d_cy[i];
+        const float ci_z = (float)d_cz[i];
+
+        // Q = ci ci - cs^2 I
+        const float q_xx = ci_x * ci_x - CS2;
+        const float q_yy = ci_y * ci_y - CS2;
+        const float q_zz = ci_z * ci_z - CS2;
+        const float q_xy = ci_x * ci_y;
+        const float q_xz = ci_x * ci_z;
+        const float q_yz = ci_y * ci_z;
+
+        // Contractions Q:Πdev and Q:Πiso
+        const float q_dot_dev =
+            q_xx * dev_xx + q_yy * dev_yy + q_zz * dev_zz
+            + 2.0f * (q_xy * dev_xy + q_xz * dev_xz + q_yz * dev_yz);
+
+        const float q_trace = q_xx + q_yy + q_zz;  // = |ci|^2 - 1
+        const float q_dot_iso = q_trace * tr_pi_thrd;
+
+        const float neq_dev = d_w[i] * proj_pref * q_dot_dev;
+        const float neq_iso = d_w[i] * proj_pref * q_dot_iso;
+        const float neq_h   = neq[i] - neq_dev - neq_iso;
+
+        g[i] = feq[i]
+             + one_minus_s_nu   * neq_dev
+             + one_minus_s_bulk * neq_iso
+             + one_minus_s_high * neq_h
+             + c_tau * Fin[i];
     }
 }
 
diff --git a/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh b/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
new file mode 100644
index 0000000..bc8b00b
--- /dev/null
+++ b/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
@@ -0,0 +1,111 @@
+// CelerisLab – operators/turbulence_smag.cuh
+// Smagorinsky LES: compute effective omega from local non-equilibrium stress.
+//
+// Model:
+//   nu_t = (C_s * Delta)^2 * |S|,  Delta = 1 (lattice unit)
+//   omega_eff = 1 / (3*(nu + nu_t) + 0.5)
+//
+// |S| is estimated from non-equilibrium second-order moments using tau0.
+// ============================================================================
+
+#ifndef CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
+#define CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
+
+#ifndef USE_LES
+#define USE_LES 0
+#endif
+
+#ifndef LES_CS
+#define LES_CS 0.16f
+#endif
+
+__device__ __forceinline__ float clamp_omega(float w) {
+    return fminf(1.999f, fmaxf(0.01f, w));
+}
+
+#if NQ == 9
+
+__device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
+                                                    const float* __restrict__ feq,
+                                                    float rho,
+                                                    float omega0)
+{
+    const float rho_safe = fmaxf(rho, 1.0e-12f);
+    const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
+
+    // Πneq = Σ (ci ci)(fi-feqi)
+    float pixx = 0.0f, piyy = 0.0f, pixy = 0.0f;
+
+    #pragma unroll
+    for (int i = 0; i < 9; i++) {
+        const float fneq = f[i] - feq[i];
+        const float cx = (float)d_cx[i];
+        const float cy = (float)d_cy[i];
+        pixx += fneq * cx * cx;
+        piyy += fneq * cy * cy;
+        pixy += fneq * cx * cy;
+    }
+
+    const float denom = 2.0f * rho_safe * CS2 * tau0;
+    const float sxx = -pixx / denom;
+    const float syy = -piyy / denom;
+    const float sxy = -pixy / denom;
+
+    // |S| = sqrt(2 Sij Sij)
+    const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + 2.0f*sxy*sxy));
+
+    const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
+    const float nut = (LES_CS * LES_CS) * s_mag;
+    const float nue = nu0 + nut;
+    return clamp_omega(1.0f / (3.0f * nue + 0.5f));
+}
+
+#elif NQ == 19
+
+__device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
+                                                    const float* __restrict__ feq,
+                                                    float rho,
+                                                    float omega0)
+{
+    const float rho_safe = fmaxf(rho, 1.0e-12f);
+    const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
+
+    // Πneq = Σ (ci ci)(fi-feqi)
+    float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
+    float pixy = 0.0f, pixz = 0.0f, piyz = 0.0f;
+
+    #pragma unroll
+    for (int i = 0; i < 19; i++) {
+        const float fneq = f[i] - feq[i];
+        const float cx = (float)d_cx[i];
+        const float cy = (float)d_cy[i];
+        const float cz = (float)d_cz[i];
+        pixx += fneq * cx * cx;
+        piyy += fneq * cy * cy;
+        pizz += fneq * cz * cz;
+        pixy += fneq * cx * cy;
+        pixz += fneq * cx * cz;
+        piyz += fneq * cy * cz;
+    }
+
+    const float denom = 2.0f * rho_safe * CS2 * tau0;
+    const float sxx = -pixx / denom;
+    const float syy = -piyy / denom;
+    const float szz = -pizz / denom;
+    const float sxy = -pixy / denom;
+    const float sxz = -pixz / denom;
+    const float syz = -piyz / denom;
+
+    // |S| = sqrt(2 Sij Sij)
+    const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + szz*szz
+                                   + 2.0f*(sxy*sxy + sxz*sxz + syz*syz)));
+
+    const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
+    const float nut = (LES_CS * LES_CS) * s_mag;
+    const float nue = nu0 + nut;
+    return clamp_omega(1.0f / (3.0f * nue + 0.5f));
+}
+
+#endif // NQ
+
+#endif // CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
diff --git a/src/CelerisLab/lbm/kernels/step/one_step_double.cu b/src/CelerisLab/lbm/kernels/step/one_step_double.cu
index 960c0f8..b6cf1a6 100644
--- a/src/CelerisLab/lbm/kernels/step/one_step_double.cu
+++ b/src/CelerisLab/lbm/kernels/step/one_step_double.cu
@@ -69,7 +69,8 @@ __global__ void StreamCollideDouble(
     // ----- Boundary conditions (inlet/outlet/wall) -----
 #if NQ == 9
     if (fl & LEGACY_SOLID) {
-        if (x == 0) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
             float f_neb[NQ];
             unsigned long k_neb = linear_index(x + 1u, y);
             for (int i = 0; i < NQ; i++) {
@@ -77,18 +78,31 @@ __global__ void StreamCollideDouble(
             }
             apply_parabolic_inlet(f, f_neb, (float)y);
         }
-        else if (x == NX - 1) {
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
             float f_neb[NQ];
             unsigned long k_neb = linear_index(x - 1u, y);
             for (int i = 0; i < NQ; i++) {
                 f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
             }
             apply_pressure_outlet(f, f_neb, (float)y);
+        } else {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+            }
+        }
+    }
+
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
         }
     }
 #elif NQ == 19
     if (fl & LEGACY_SOLID) {
-        if (x == 0) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
             float f_neb[NQ];
             unsigned long k_neb = linear_index(x + 1u, y, z);
             for (int i = 0; i < NQ; i++) {
@@ -96,13 +110,18 @@ __global__ void StreamCollideDouble(
             }
             apply_parabolic_inlet_3d(f, f_neb, (float)y);
         }
-        else if (x == NX - 1) {
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
             float f_neb[NQ];
             unsigned long k_neb = linear_index(x - 1u, y, z);
             for (int i = 0; i < NQ; i++) {
                 f_neb[i] = load_ddf(fi_in, index_f(k_neb, (unsigned int)i));
             }
             apply_pressure_outlet_3d(f, f_neb, (float)y);
+        } else {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float t = f[i]; f[i] = f[i+1]; f[i+1] = t;
+            }
         }
     }
 #endif
@@ -126,20 +145,36 @@ __global__ void StreamCollideDouble(
 #if NQ == 9
         compute_feq(rho_n, ux, uy, feq);
         zero_forcing(Fin);
+    float omega_col = d_params.omega;
+#if USE_LES
+    omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+#endif
+    omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
 
         // ----- Collision -----
 #if COLLISION_MODEL == 0
-        collide_srt(f, feq, Fin, d_params.omega);
+    collide_srt(f, feq, Fin, omega_col);
 #elif COLLISION_MODEL == 1
-        collide_trt(f, feq, Fin, d_params.omega);
+    collide_trt(f, feq, Fin, omega_col);
 #elif COLLISION_MODEL == 2
-        collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
+    collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
 #endif
 
 #elif NQ == 19
         compute_feq(rho_n, ux, uy, uz, feq);
         zero_forcing(Fin);
-        collide_srt(f, feq, Fin, d_params.omega);
+    float omega_col = d_params.omega;
+#if USE_LES
+    omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+#endif
+    omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
+#if COLLISION_MODEL == 0
+    collide_srt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 1
+    collide_trt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 2
+    collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
+#endif
 #endif
     }
 
diff --git a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu
index 907fd63..d93b094 100644
--- a/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu
+++ b/src/CelerisLab/lbm/kernels/step/one_step_esopull.cu
@@ -70,7 +70,8 @@ __global__ void StreamCollideEsoPull(
     // ----- Boundary conditions -----
 #if NQ == 9
     if (fl & LEGACY_SOLID) {
-        if (x == 0) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
             float f_neb[NQ];
             unsigned long k_neb = linear_index(x + 1u, y);
             load_f_esopull(k_neb, f_neb, fi, j, t);  // approximate: reuse j
@@ -80,19 +81,67 @@ __global__ void StreamCollideEsoPull(
             load_f_esopull(k_neb, f_neb, fi, j_neb, t);
             apply_parabolic_inlet(f, f_neb, (float)y);
         }
-        else if (x == (unsigned int)(NX - 1)) {
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
             unsigned long k_neb = linear_index(x - 1u, y);
             unsigned long j_neb[NQ];
             compute_neighbors(k_neb, j_neb);
             float f_neb[NQ];
             load_f_esopull(k_neb, f_neb, fi, j_neb, t);
             apply_pressure_outlet(f, f_neb, (float)y);
+        } else {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
+            }
+        }
+    }
+
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
         }
     }
 
     if (y == 1 || y == (unsigned int)(NY - 2)) {
         apply_wall_bb_d2q9(y, f);
     }
+#elif NQ == 19
+    if (fl & LEGACY_SOLID) {
+        bool interior_y = (y > 0u) && (y < (unsigned int)(NY - 1));
+        if (x == 0 && interior_y) {
+            unsigned long k_neb = linear_index(x + 1u, y, z);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_parabolic_inlet_3d(f, f_neb, (float)y);
+        }
+        else if (x == (unsigned int)(NX - 1) && interior_y) {
+            unsigned long k_neb = linear_index(x - 1u, y, z);
+            unsigned long j_neb[NQ];
+            compute_neighbors(k_neb, j_neb);
+            float f_neb[NQ];
+            load_f_esopull(k_neb, f_neb, fi, j_neb, t);
+            apply_pressure_outlet_3d(f, f_neb, (float)y);
+        } else {
+            #pragma unroll
+            for (int i = 1; i < NQ; i += 2) {
+                float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
+            }
+        }
+    }
+
+    if (fl & LEGACY_OBSTACLE) {
+        #pragma unroll
+        for (int i = 1; i < NQ; i += 2) {
+            float ttmp = f[i]; f[i] = f[i+1]; f[i+1] = ttmp;
+        }
+    }
+
+    if (y == 1 || y == (unsigned int)(NY - 2)) {
+        apply_wall_bb_d3q19_y(y, f);
+    }
 #endif
 
     // ----- Forcing -----
@@ -123,13 +172,18 @@ __global__ void StreamCollideEsoPull(
         float feq[NQ];
 #if NQ == 9
         compute_feq(rho_n, ux, uy, feq);
+    float omega_col = d_params.omega;
+#if USE_LES
+    omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+#endif
+    omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
 
 #if COLLISION_MODEL == 0
-        collide_srt(f, feq, Fin, d_params.omega);
+    collide_srt(f, feq, Fin, omega_col);
 #elif COLLISION_MODEL == 1
-        collide_trt(f, feq, Fin, d_params.omega);
+    collide_trt(f, feq, Fin, omega_col);
 #elif COLLISION_MODEL == 2
-        collide_mrt(f, rho_n, ux, uy, Fin, d_params.omega);
+    collide_mrt(f, rho_n, ux, uy, Fin, omega_col);
 #endif
 
 #elif NQ == 19
@@ -137,7 +191,18 @@ __global__ void StreamCollideEsoPull(
         if (force_field != nullptr) fzn += force_field[2*TOTAL_CELLS + k];
         apply_guo_velocity_correction(ux, uy, uz, fxn, fyn, fzn, rho_n);
         compute_feq(rho_n, ux, uy, uz, feq);
-        collide_srt(f, feq, Fin, d_params.omega);
+    float omega_col = d_params.omega;
+#if USE_LES
+    omega_col = compute_omega_smag(f, feq, rho_n, omega_col);
+#endif
+    omega_col = fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, omega_col));
+#if COLLISION_MODEL == 0
+    collide_srt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 1
+    collide_trt(f, feq, Fin, omega_col);
+#elif COLLISION_MODEL == 2
+    collide_mrt(f, rho_n, ux, uy, uz, Fin, omega_col);
+#endif
 #endif
     }
 
diff --git a/tests/test_d2q9_visual.py b/tests/test_d2q9_visual.py
index 8c96c84..9f60b66 100644
--- a/tests/test_d2q9_visual.py
+++ b/tests/test_d2q9_visual.py
@@ -2,12 +2,13 @@
 """
 D2Q9 Regression Test — Poiseuille Channel + Cylinder Flow
 ==========================================================
-Uses the ORIGINAL kernel.cu with same grid / BCs.
+Uses kernel_v2.cu by default (legacy kernel.cu remains optional fallback).
 Produces matplotlib figures for visual validation.
 
 Usage:
     python tests/test_d2q9_visual.py --device 2
     python tests/test_d2q9_visual.py --device 2 --cylinder
+    python tests/test_d2q9_visual.py --device 2 --legacy
 """
 
 import sys, os, argparse, time
@@ -47,7 +48,7 @@ INTERFACE_FLAG = 0b00001000
 # ━━━━━━━━━━━━━━━━━━━━━━━  Helpers  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
 
 def configure_macros(n_objs=0):
-    """Write macros.h to match original kernel settings."""
+    """Write macros.h for D2Q9 regression (v2 default, legacy optional)."""
     lines = compiler.read_lines(compiler.kernel_path("macros.h"))
     defs = {
         'MULT_GPU': 'False', 'NT': NT,
@@ -58,23 +59,31 @@ def configure_macros(n_objs=0):
         'DIM': DIM, 'NQ': NQ,
         'VIS': VIS, 'RHO': f'{RHO}', 'U0': U0,
         'N_OBJS': n_objs,
+        'COLLISION_MODEL': 0,     # SRT
+        'STREAMING_MODEL': 0,     # double-buffer
+        'STORE_PRECISION': 0,     # FP32
+        'USE_DDF_SHIFTING': 0,    # keep unshifted for v2 defaults
     }
     for name, val in defs.items():
         lines = compiler.modify_macro(lines, name, val)
     compiler.write_lines(compiler.kernel_path("macros.h"), lines)
 
 
-def extract_fields(ddf_host):
-    """Compute rho, u, v from host DDF (original const.h ordering).
-
-    The original kernel uses DDF-shifting: stores f_shifted = f - w_i*RHO.
-    So sum(f_shifted) = rho - RHO  (~0 for incompressible flow),
-    and momentum = sum(e_x * f_shifted) works because sum(w_i * e_x) = 0.
-    """
+def extract_fields(ddf_host, use_ddf_shifting=False):
+    """Compute rho, u, v from host DDF in original D2Q9 direction ordering."""
     f = ddf_host.reshape(NQ, NY, NX)
-    rho = np.sum(f, axis=0) + RHO          # un-shift: rho = sum(f_shifted) + RHO
-    u = (f[1] + f[5] + f[8] - f[3] - f[6] - f[7]) / RHO
-    v = (f[2] + f[5] + f[6] - f[4] - f[7] - f[8]) / RHO
+    if use_ddf_shifting:
+        rho = np.sum(f, axis=0) + RHO
+        denom = np.full_like(rho, RHO)
+    else:
+        rho = np.sum(f, axis=0)
+        denom = rho
+
+    denom_safe = np.where(np.abs(denom) > 1e-12, denom, 1.0)
+    u = (f[1] + f[5] + f[8] - f[3] - f[6] - f[7]) / denom_safe
+    v = (f[2] + f[5] + f[6] - f[4] - f[7] - f[8]) / denom_safe
+    u = np.where(np.abs(denom) > 1e-12, u, 0.0)
+    v = np.where(np.abs(denom) > 1e-12, v, 0.0)
     return rho, u, v
 
 
@@ -132,7 +141,7 @@ def build_cylinder_data(cx, cy, radius):
 
 # ━━━━━━━━━━━━━━━━━━━━━━━  Simulation  ━━━━━━━━━━━━━━━━━━━━━━━━━━
 
-def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host):
+def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host, use_legacy=False):
     """Compile kernel, run LBM, return DDF on host."""
     cuda.init()
     dev = cuda.Device(device_id)
@@ -141,11 +150,17 @@ def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host)
 
     try:
         configure_macros(n_objs)
-        compiler.compile_kernel()
-        ptx_path = compiler.kernel_path("kernel.ptx")
+        if use_legacy:
+            compiler.compile_kernel()
+            ptx_path = compiler.kernel_path("kernel.ptx")
+            init_name = "InitTubeFlow"
+        else:
+            compiler.compile_kernel_v2()
+            ptx_path = compiler.kernel_path("kernel_v2.ptx")
+            init_name = "InitTubeFlow_v2"
         mod = cuda.module_from_file(ptx_path)
         step_fn = mod.get_function("OneStep")
-        init_fn = mod.get_function("InitTubeFlow")
+        init_fn = mod.get_function(init_name)
 
         nbytes_ddf = TOTAL * NQ * 4
         ddf_gpu   = cuda.mem_alloc(nbytes_ddf)
@@ -198,9 +213,9 @@ def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host)
 
 # ━━━━━━━━━━━━━━━━━━━━━━━  Visualization  ━━━━━━━━━━━━━━━━━━━━━━━
 
-def plot_poiseuille(ddf, flag, out_path):
+def plot_poiseuille(ddf, flag, out_path, use_ddf_shifting=False):
     """3-panel figure: velocity mag, u(y) profile, pressure along centerline."""
-    rho, u, v = extract_fields(ddf)
+    rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
     vel_mag = np.sqrt(u**2 + v**2)
 
     # Mask solid cells for display
@@ -244,9 +259,9 @@ def plot_poiseuille(ddf, flag, out_path):
     plt.close(fig)
 
 
-def plot_cylinder(ddf, flag, cx, cy, radius, out_path):
+def plot_cylinder(ddf, flag, cx, cy, radius, out_path, use_ddf_shifting=False):
     """3-panel figure: velocity magnitude (zoom), vorticity, streamlines."""
-    rho, u, v = extract_fields(ddf)
+    rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
     vel_mag = np.sqrt(u**2 + v**2)
     mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool)
 
@@ -313,6 +328,8 @@ def main():
     parser = argparse.ArgumentParser(description='D2Q9 Regression Test')
     parser.add_argument('--device', type=int, default=2,
                         help='CUDA device ID (default: 2)')
+    parser.add_argument('--legacy', action='store_true',
+                        help='Use legacy kernel.cu path (default uses kernel_v2.cu)')
     parser.add_argument('--cylinder', action='store_true',
                         help='Also run cylinder flow test')
     parser.add_argument('--steps-pois', type=int, default=5000,
@@ -324,6 +341,10 @@ def main():
     out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output')
     os.makedirs(out_dir, exist_ok=True)
 
+    use_ddf_shifting = bool(args.legacy)
+    mode = 'legacy kernel.cu' if args.legacy else 'kernel_v2.cu'
+    print(f"\n[Mode] {mode}")
+
     # ---- Test 1: Poiseuille ----
     print("\n===== Test 1: Poiseuille Channel Flow =====")
     flag_pois = np.ones(TOTAL, dtype=np.uint8)
@@ -331,11 +352,13 @@ def main():
     delta_pois = np.zeros(1, dtype=np.float32)
 
     ddf, flag = run_simulation(args.device, args.steps_pois, 0,
-                               flag_pois, indx_pois, delta_pois)
-    plot_poiseuille(ddf, flag, os.path.join(out_dir, 'poiseuille_d2q9.png'))
+                               flag_pois, indx_pois, delta_pois,
+                               use_legacy=args.legacy)
+    plot_poiseuille(ddf, flag, os.path.join(out_dir, 'poiseuille_d2q9.png'),
+                    use_ddf_shifting=use_ddf_shifting)
 
     # Error metric
-    rho, u, v = extract_fields(ddf)
+    rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
     y_arr = np.arange(NY, dtype=float)
     u_ana = analytical_poiseuille(y_arr)
     x_mid = NX // 2
@@ -351,9 +374,11 @@ def main():
         flag_cyl, indx_cyl, delta_cyl = build_cylinder_data(cyl_cx, cyl_cy, cyl_r)
 
         ddf2, flag2 = run_simulation(args.device, args.steps_cyl, 1,
-                                     flag_cyl, indx_cyl, delta_cyl)
+                         flag_cyl, indx_cyl, delta_cyl,
+                         use_legacy=args.legacy)
         plot_cylinder(ddf2, flag2, cyl_cx, cyl_cy, cyl_r,
-                      os.path.join(out_dir, 'cylinder_d2q9.png'))
+                  os.path.join(out_dir, 'cylinder_d2q9.png'),
+                  use_ddf_shifting=use_ddf_shifting)
 
     print("\nDone.")
 
diff --git a/tests/test_high_re_validation.py b/tests/test_high_re_validation.py
new file mode 100644
index 0000000..9b6fcd3
--- /dev/null
+++ b/tests/test_high_re_validation.py
@@ -0,0 +1,547 @@
+#!/usr/bin/env python3
+"""
+High-Re Validation (kernel_v2)
+==============================
+Unified validation script for high-Re runs with optional LES.
+
+Default targets:
+- 2D D2Q9: Re=5000
+- 3D D3Q19: Re=3000
+
+The script configures macros.h temporarily, compiles kernel_v2, runs the case,
+and restores macros.h automatically.
+"""
+
+import argparse
+import json
+import os
+import struct
+import sys
+import time
+
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "src"))
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import pycuda.driver as cuda
+
+from CelerisLab.cuda import compiler
+
+FLUID_FLAG = 0x01
+SOLID_FLAG = 0x02
+OBSTACLE_FLAG = 0x04
+
+
+def collision_name(model):
+    return {0: "SRT", 1: "TRT", 2: "MRT"}.get(model, f"M{model}")
+
+
+def make_case_tag(cfg):
+    les_tag = "LES" if cfg["use_les"] else "NoLES"
+    return (
+        f"{cfg['name']}_Re{int(cfg['target_re'])}_"
+        f"{collision_name(cfg['collision_model'])}_{les_tag}_"
+        f"OM{int(cfg['outlet_mode'])}_WMAX{cfg['omega_collision_max']:.3f}"
+    )
+
+
+def validate_case(rho):
+    nan_count = int(np.isnan(rho).sum())
+    if nan_count > 0:
+        return False, "NaN detected"
+
+    rho_min = float(np.min(rho))
+    rho_max = float(np.max(rho))
+    if rho_min <= 0.0:
+        return False, "Non-positive density"
+    if rho_max >= 2.0:
+        return False, "Density blow-up"
+    return True, "OK"
+
+
+def plot_case(cfg, host_ddf, out_dir):
+    nq = cfg["nq"]
+    nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
+    flag = cfg["flag"]
+    tag = make_case_tag(cfg)
+    out_path = os.path.join(out_dir, f"{tag}.png")
+
+    if nq == 9:
+        f = host_ddf.reshape(nq, ny, nx)
+        rho = np.sum(f, axis=0)
+        ux = np.zeros_like(rho)
+        uy = np.zeros_like(rho)
+        cx = [0, 1, -1, 0, 0, 1, -1, 1, -1]
+        cy = [0, 0, 0, 1, -1, 1, -1, -1, 1]
+        for i in range(nq):
+            ux += cx[i] * f[i]
+            uy += cy[i] * f[i]
+        rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
+        ux /= rho_safe
+        uy /= rho_safe
+        vel = np.sqrt(ux * ux + uy * uy)
+
+        mask = flag.reshape(ny, nx) != FLUID_FLAG
+        vel_m = np.ma.array(vel, mask=mask)
+        vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
+        vort_m = np.ma.array(vort, mask=mask)
+
+        fig, axes = plt.subplots(1, 3, figsize=(16, 5))
+
+        im0 = axes[0].imshow(vel_m, origin="lower", aspect="auto", cmap="turbo")
+        plt.colorbar(im0, ax=axes[0], label="|u|")
+        axes[0].set_title("Velocity Magnitude")
+
+        vmax = np.percentile(np.abs(vort[~mask]), 99) if np.any(~mask) else 1e-6
+        vmax = max(vmax, 1.0e-6)
+        im1 = axes[1].imshow(vort_m, origin="lower", aspect="auto", cmap="RdBu_r", vmin=-vmax, vmax=vmax)
+        plt.colorbar(im1, ax=axes[1], label="vorticity")
+        axes[1].set_title("Vorticity")
+
+        X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
+        ux_s = np.ma.array(ux, mask=mask)
+        uy_s = np.ma.array(uy, mask=mask)
+        speed = np.ma.sqrt(ux_s * ux_s + uy_s * uy_s)
+        axes[2].streamplot(X, Y, ux_s, uy_s, color=speed, cmap="viridis", density=2.0, linewidth=0.7)
+        axes[2].set_xlim(0, nx)
+        axes[2].set_ylim(0, ny)
+        axes[2].set_title("Streamlines")
+
+        fig.suptitle(tag)
+        fig.tight_layout()
+        fig.savefig(out_path, dpi=150)
+        plt.close(fig)
+        return out_path
+
+    # D3Q19: visualize mid-z slice
+    f = host_ddf.reshape(nq, nz, ny, nx)
+    z0 = nz // 2
+    fs = f[:, z0, :, :]
+    rho = np.sum(fs, axis=0)
+    ux = np.zeros_like(rho)
+    uy = np.zeros_like(rho)
+    uz = np.zeros_like(rho)
+    cx = np.array([0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0])
+    cy = np.array([0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1])
+    cz = np.array([0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1])
+    for i in range(nq):
+        ux += cx[i] * fs[i]
+        uy += cy[i] * fs[i]
+        uz += cz[i] * fs[i]
+    rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
+    ux /= rho_safe
+    uy /= rho_safe
+    uz /= rho_safe
+    vel = np.sqrt(ux * ux + uy * uy + uz * uz)
+
+    mask3 = flag.reshape(nz, ny, nx)[z0] != FLUID_FLAG
+    vel_m = np.ma.array(vel, mask=mask3)
+    vort = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
+    vort_m = np.ma.array(vort, mask=mask3)
+
+    fig, axes = plt.subplots(1, 3, figsize=(16, 5))
+    im0 = axes[0].imshow(vel_m, origin="lower", aspect="auto", cmap="turbo")
+    plt.colorbar(im0, ax=axes[0], label="|u|")
+    axes[0].set_title("Velocity Magnitude (z-mid)")
+
+    vmax = np.percentile(np.abs(vort[~mask3]), 99) if np.any(~mask3) else 1.0e-6
+    vmax = max(vmax, 1.0e-6)
+    im1 = axes[1].imshow(vort_m, origin="lower", aspect="auto", cmap="RdBu_r", vmin=-vmax, vmax=vmax)
+    plt.colorbar(im1, ax=axes[1], label="vorticity")
+    axes[1].set_title("Vorticity (z-mid)")
+
+    X, Y = np.meshgrid(np.arange(nx), np.arange(ny))
+    ux_s = np.ma.array(ux, mask=mask3)
+    uy_s = np.ma.array(uy, mask=mask3)
+    speed = np.ma.sqrt(ux_s * ux_s + uy_s * uy_s)
+    axes[2].streamplot(X, Y, ux_s, uy_s, color=speed, cmap="viridis", density=2.0, linewidth=0.7)
+    axes[2].set_xlim(0, nx)
+    axes[2].set_ylim(0, ny)
+    axes[2].set_title("Streamlines (z-mid)")
+
+    fig.suptitle(tag)
+    fig.tight_layout()
+    fig.savefig(out_path, dpi=150)
+    plt.close(fig)
+    return out_path
+
+
+def compute_vis_omega(reynolds, diameter, u0):
+    vis = u0 * diameter / reynolds
+    omega = 1.0 / (3.0 * vis + 0.5)
+    return vis, omega
+
+
+def set_macros(nx, ny, nz, dim, nq, vis, u0, collision_model, use_les, les_cs,
+               outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
+               omega_collision_max):
+    lines = compiler.read_lines(compiler.kernel_path("macros.h"))
+    defs = {
+        "MULT_GPU": "False",
+        "NT": 128,
+        "X_1U": nx,
+        "Y_1U": ny,
+        "Z_1U": nz,
+        "LBtype": "float",
+        "UX": 1,
+        "UY": 1,
+        "UZ": 1,
+        "NX": nx,
+        "NY": ny,
+        "NZ": nz,
+        "DIM": dim,
+        "NQ": nq,
+        "VIS": f"{vis:.10f}",
+        "RHO": "1.0",
+        "U0": u0,
+        "N_OBJS": 0,
+        "COLLISION_MODEL": collision_model,
+        "STREAMING_MODEL": 0,
+        "STORE_PRECISION": 0,
+        "USE_DDF_SHIFTING": 0,
+        "USE_LES": int(use_les),
+        "LES_CS": f"{les_cs:.6f}f",
+        "INLET_PROFILE": 0,
+        "OUTLET_MODE": int(outlet_mode),
+        "OUTLET_BACKFLOW_CLAMP": int(outlet_backflow_clamp),
+        "OUTLET_BLEND_ALPHA": f"{float(outlet_blend_alpha):.3f}f",
+        "OMEGA_COLLISION_MAX": f"{float(omega_collision_max):.3f}f",
+    }
+    for name, value in defs.items():
+        lines = compiler.modify_macro(lines, name, value)
+    compiler.write_lines(compiler.kernel_path("macros.h"), lines)
+
+
+def build_flags_2d(nx, ny, cx, cy, radius):
+    n = nx * ny
+    flag = np.ones(n, dtype=np.uint8) * FLUID_FLAG
+    for y in range(ny):
+        for x in range(nx):
+            k = y * nx + x
+            if y == 0 or y == ny - 1 or x == 0 or x == nx - 1:
+                flag[k] = SOLID_FLAG
+            elif (x - cx) ** 2 + (y - cy) ** 2 < radius ** 2:
+                flag[k] = OBSTACLE_FLAG
+    return flag
+
+
+def build_flags_3d(nx, ny, nz, cx, cy, radius):
+    n = nx * ny * nz
+    flag = np.ones(n, dtype=np.uint8) * FLUID_FLAG
+    for z in range(nz):
+        for y in range(ny):
+            for x in range(nx):
+                k = z * ny * nx + y * nx + x
+                if y == 0 or y == ny - 1 or x == 0 or x == nx - 1:
+                    flag[k] = SOLID_FLAG
+                elif (x - cx) ** 2 + (y - cy) ** 2 < radius ** 2:
+                    flag[k] = OBSTACLE_FLAG
+    return flag
+
+
+def run_case(device_id, cfg):
+    nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
+    dim, nq = cfg["dim"], cfg["nq"]
+    n = nx * ny * nz
+
+    set_macros(
+        nx=nx,
+        ny=ny,
+        nz=nz,
+        dim=dim,
+        nq=nq,
+        vis=cfg["vis"],
+        u0=cfg["u0"],
+        collision_model=cfg["collision_model"],
+        use_les=cfg["use_les"],
+        les_cs=cfg["les_cs"],
+        outlet_mode=cfg["outlet_mode"],
+        outlet_backflow_clamp=cfg["outlet_backflow_clamp"],
+        outlet_blend_alpha=cfg["outlet_blend_alpha"],
+        omega_collision_max=cfg["omega_collision_max"],
+    )
+    compiler.compile_kernel_v2()
+
+    cuda.init()
+    dev = cuda.Device(device_id)
+    ctx = dev.make_context()
+    try:
+        mod = cuda.module_from_file(compiler.kernel_path("kernel_v2.ptx"))
+        init_fn = mod.get_function("InitTubeFlow_v2")
+        step_fn = mod.get_function("OneStep")
+
+        params_ptr, params_size = mod.get_global("d_params")
+        params_data = struct.pack(
+            "IIIQfffffffI",
+            nx,
+            ny,
+            nz,
+            n,
+            cfg["omega"],
+            1.1,
+            0.0,
+            0.0,
+            0.0,
+            1.0,
+            cfg["u0"],
+            0,
+        )
+        if len(params_data) < params_size:
+            params_data += b"\x00" * (params_size - len(params_data))
+        cuda.memcpy_htod(params_ptr, params_data)
+
+        fsize = n * nq * 4
+        d_fi = cuda.mem_alloc(fsize)
+        d_fi2 = cuda.mem_alloc(fsize)
+        d_flag = cuda.mem_alloc(n)
+        d_indx = cuda.mem_alloc(n * 4)
+        d_delta = cuda.mem_alloc(4)
+        d_action = cuda.mem_alloc(4)
+        d_obs = cuda.mem_alloc(4)
+
+        cuda.memset_d32(d_indx, 0, n)
+        cuda.memset_d32(d_delta, 0, 1)
+        cuda.memset_d32(d_action, 0, 1)
+        cuda.memset_d32(d_obs, 0, 1)
+
+        block = (128, 1, 1)
+        grid = ((nx + 127) // 128, ny, nz)
+
+        init_fn(d_flag, d_fi, block=block, grid=grid)
+        cuda.memcpy_dtod(d_fi2, d_fi, fsize)
+
+        cuda.memcpy_htod(d_flag, cfg["flag"])
+
+        t0 = time.time()
+        for step in range(cfg["steps"]):
+            step_fn(d_flag, d_fi, d_fi2, d_indx, d_delta, d_action, d_obs, block=block, grid=grid)
+            d_fi, d_fi2 = d_fi2, d_fi
+
+            if (step + 1) % cfg["report_every"] == 0:
+                cuda.Context.synchronize()
+                host = np.empty(n * nq, dtype=np.float32)
+                cuda.memcpy_dtoh(host, d_fi)
+                if nq == 9:
+                    rho = host.reshape(nq, ny, nx).sum(axis=0)
+                    c = float(rho[ny // 2, nx // 2])
+                else:
+                    rho = host.reshape(nq, nz, ny, nx).sum(axis=0)
+                    c = float(rho[nz // 2, ny // 2, nx // 2])
+                nan_count = int(np.isnan(rho).sum())
+                print(f"  step {step+1:7d}: rho_center={c:.6f}, nan={nan_count}")
+                if nan_count > 0:
+                    break
+
+        cuda.Context.synchronize()
+        elapsed = time.time() - t0
+
+        host = np.empty(n * nq, dtype=np.float32)
+        cuda.memcpy_dtoh(host, d_fi)
+        if nq == 9:
+            rho = host.reshape(nq, ny, nx).sum(axis=0)
+            center = float(rho[ny // 2, nx // 2])
+        else:
+            rho = host.reshape(nq, nz, ny, nx).sum(axis=0)
+            center = float(rho[nz // 2, ny // 2, nx // 2])
+
+        ok, reason = validate_case(rho)
+        plot_path = None
+        if cfg.get("save_plot", True):
+            plot_path = plot_case(cfg, host, cfg["out_dir"])
+
+        return {
+            "case_tag": make_case_tag(cfg),
+            "name": cfg["name"],
+            "target_re": cfg["target_re"],
+            "steps": cfg["steps"],
+            "mlups": float(n * cfg["steps"] / elapsed / 1e6),
+            "nan_count": int(np.isnan(rho).sum()),
+            "rho_center": center,
+            "rho_min": float(np.nanmin(rho)),
+            "rho_max": float(np.nanmax(rho)),
+            "omega": cfg["omega"],
+            "vis": cfg["vis"],
+            "collision_model": cfg["collision_model"],
+            "use_les": bool(cfg["use_les"]),
+            "les_cs": float(cfg["les_cs"]),
+            "outlet_mode": int(cfg["outlet_mode"]),
+            "outlet_backflow_clamp": int(cfg["outlet_backflow_clamp"]),
+            "outlet_blend_alpha": float(cfg["outlet_blend_alpha"]),
+            "omega_collision_max": float(cfg["omega_collision_max"]),
+            "pass": bool(ok),
+            "reason": reason,
+            "plot_path": plot_path,
+        }
+    finally:
+        ctx.pop()
+
+
+def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
+                  outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
+                  omega_collision_max):
+    nx, ny, nz = 512, 256, 1
+    cx, cy, radius = 128.0, 128.0, 24.0
+    u0 = 0.03
+    vis, omega = compute_vis_omega(re2d, 2.0 * radius, u0)
+    return {
+        "name": "2D_D2Q9_highRe",
+        "dim": 2,
+        "nq": 9,
+        "nx": nx,
+        "ny": ny,
+        "nz": nz,
+        "flag": build_flags_2d(nx, ny, cx, cy, radius),
+        "u0": u0,
+        "vis": vis,
+        "omega": omega,
+        "steps": steps2d,
+        "report_every": max(steps2d // 10, 1),
+        "collision_model": collision_model,
+        "use_les": use_les,
+        "les_cs": les_cs,
+        "outlet_mode": int(outlet_mode),
+        "outlet_backflow_clamp": int(outlet_backflow_clamp),
+        "outlet_blend_alpha": float(outlet_blend_alpha),
+        "omega_collision_max": float(omega_collision_max),
+        "target_re": re2d,
+        "save_plot": True,
+        "out_dir": out_dir,
+    }
+
+def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
+                  outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
+                  omega_collision_max):
+    nx, ny, nz = 256, 128, 32
+    cx, cy, radius = 64.0, 64.0, 12.0
+    u0 = 0.04
+    vis, omega = compute_vis_omega(re3d, 2.0 * radius, u0)
+    return {
+        "name": "3D_D3Q19_highRe",
+        "dim": 3,
+        "nq": 19,
+        "nx": nx,
+        "ny": ny,
+        "nz": nz,
+        "flag": build_flags_3d(nx, ny, nz, cx, cy, radius),
+        "u0": u0,
+        "vis": vis,
+        "omega": omega,
+        "steps": steps3d,
+        "report_every": max(steps3d // 10, 1),
+        "collision_model": collision_model,
+        "use_les": use_les,
+        "les_cs": les_cs,
+        "outlet_mode": int(outlet_mode),
+        "outlet_backflow_clamp": int(outlet_backflow_clamp),
+        "outlet_blend_alpha": float(outlet_blend_alpha),
+        "omega_collision_max": float(omega_collision_max),
+        "target_re": re3d,
+        "save_plot": True,
+        "out_dir": out_dir,
+    }
+
+
+def build_comprehensive_cases(args, out_dir):
+    cases = []
+    # Coverage matrix at moderate Re to verify all changed pathways.
+    for cm in (0, 1, 2):
+        for les in (0, 1):
+            cases.append(build_case_2d(re2d=200.0, steps2d=args.matrix_steps2d,
+                                       collision_model=cm, use_les=les,
+                                       les_cs=args.les_cs, out_dir=out_dir,
+                                       outlet_mode=args.outlet_mode,
+                                       outlet_backflow_clamp=1,
+                                       outlet_blend_alpha=args.outlet_blend_alpha,
+                                       omega_collision_max=args.omega_collision_max))
+            cases.append(build_case_3d(re3d=200.0, steps3d=args.matrix_steps3d,
+                                       collision_model=cm, use_les=les,
+                                       les_cs=args.les_cs, out_dir=out_dir,
+                                       outlet_mode=args.outlet_mode,
+                                       outlet_backflow_clamp=1,
+                                       outlet_blend_alpha=args.outlet_blend_alpha,
+                                       omega_collision_max=args.omega_collision_max))
+    return cases
+
+
+def main():
+    parser = argparse.ArgumentParser(description="High-Re validation for kernel_v2")
+    parser.add_argument("--device", type=int, default=0)
+    parser.add_argument("--re2d", type=float, default=5000.0)
+    parser.add_argument("--re3d", type=float, default=3000.0)
+    parser.add_argument("--steps2d", type=int, default=10000)
+    parser.add_argument("--steps3d", type=int, default=20000)
+    parser.add_argument("--collision", type=int, default=1, choices=[0, 1, 2],
+                        help="0=SRT, 1=TRT, 2=MRT")
+    parser.add_argument("--use-les", action="store_true", default=True,
+                        help="Enable Smagorinsky LES")
+    parser.add_argument("--no-les", action="store_false", dest="use_les")
+    parser.add_argument("--les-cs", type=float, default=0.16)
+    parser.add_argument("--outlet-mode", type=int, default=0, choices=[0, 1, 2],
+                        help="0=non-equilibrium extrapolation, 1=zero-gradient copy, 2=damped blend")
+    parser.add_argument("--outlet-blend-alpha", type=float, default=0.70,
+                        help="Blend alpha for outlet-mode 2")
+    parser.add_argument("--omega-collision-max", type=float, default=1.999,
+                        help="Upper clamp for collision omega")
+    parser.add_argument("--only", choices=["2d", "3d", "both"], default="both")
+    parser.add_argument("--comprehensive", action="store_true",
+                        help="Run coverage matrix: SRT/TRT/MRT x LES on/off for 2D and 3D")
+    parser.add_argument("--matrix-steps2d", type=int, default=1000)
+    parser.add_argument("--matrix-steps3d", type=int, default=600)
+    args = parser.parse_args()
+
+    macro_path = compiler.kernel_path("macros.h")
+    macro_backup = compiler.read_lines(macro_path)
+
+    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "output")
+    os.makedirs(out_dir, exist_ok=True)
+    out_json = os.path.join(out_dir, "high_re_validation_summary.json")
+
+    try:
+        results = []
+
+        if args.only in ("2d", "both"):
+            c2 = build_case_2d(args.re2d, args.steps2d, args.collision, args.use_les,
+                               args.les_cs, out_dir, args.outlet_mode, 1,
+                               args.outlet_blend_alpha, args.omega_collision_max)
+            print("\n=== Running 2D high-Re case ===")
+            print(f"  target Re={args.re2d:.1f}, vis={c2['vis']:.6e}, omega={c2['omega']:.6f}")
+            results.append(run_case(args.device, c2))
+
+        if args.only in ("3d", "both"):
+            c3 = build_case_3d(args.re3d, args.steps3d, args.collision, args.use_les,
+                               args.les_cs, out_dir, args.outlet_mode, 1,
+                               args.outlet_blend_alpha, args.omega_collision_max)
+            print("\n=== Running 3D high-Re case ===")
+            print(f"  target Re={args.re3d:.1f}, vis={c3['vis']:.6e}, omega={c3['omega']:.6f}")
+            results.append(run_case(args.device, c3))
+
+        if args.comprehensive:
+            print("\n=== Running comprehensive coverage matrix ===")
+            for cfg in build_comprehensive_cases(args, out_dir):
+                print(f"  {cfg['name']} Re={cfg['target_re']:.1f} "
+                      f"{collision_name(cfg['collision_model'])} LES={int(cfg['use_les'])}")
+                results.append(run_case(args.device, cfg))
+
+        with open(out_json, "w", encoding="utf-8") as f:
+            json.dump(results, f, indent=2)
+
+        print("\n=== Summary ===")
+        n_pass = 0
+        for r in results:
+            if r["pass"]:
+                n_pass += 1
+            print(f"{r['name']}: nan={r['nan_count']}, rho_center={r['rho_center']:.6f}, "
+                  f"rho[min,max]=[{r['rho_min']:.6f}, {r['rho_max']:.6f}], "
+                  f"MLUPS={r['mlups']:.1f}, pass={r['pass']} ({r['reason']})")
+            if r.get("plot_path"):
+                print(f"  plot: {r['plot_path']}")
+        print(f"Pass rate: {n_pass}/{len(results)}")
+        print(f"Saved: {out_json}")
+    finally:
+        compiler.write_lines(macro_path, macro_backup)
+
+
+if __name__ == "__main__":
+    main()