diff --git a/.gitignore b/.gitignore
index 2982506..4928f4a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -28,6 +28,7 @@ wheels/
 
 # IDE
 .vscode/
+.github/
 .idea/
 *.swp
 *.swo
diff --git a/output/cylinder_d2q9.png b/output/cylinder_d2q9.png
deleted file mode 100644
index c75dfb4..0000000
Binary files a/output/cylinder_d2q9.png and /dev/null differ
diff --git a/output/cylinder_d3q19_re100.png b/output/cylinder_d3q19_re100.png
deleted file mode 100644
index 4617ce0..0000000
Binary files a/output/cylinder_d3q19_re100.png and /dev/null differ
diff --git a/output/cylinder_d3q19_re20.png b/output/cylinder_d3q19_re20.png
deleted file mode 100644
index 87355fd..0000000
Binary files a/output/cylinder_d3q19_re20.png and /dev/null differ
diff --git a/output/poiseuille_d2q9.png b/output/poiseuille_d2q9.png
deleted file mode 100644
index 36c5a39..0000000
Binary files a/output/poiseuille_d2q9.png and /dev/null differ
diff --git a/src/CelerisLab/lbm/driver.py b/src/CelerisLab/lbm/driver.py
index d33bd17..8f844b7 100644
--- a/src/CelerisLab/lbm/driver.py
+++ b/src/CelerisLab/lbm/driver.py
@@ -193,12 +193,17 @@ class FlowField:
             # max_id = max(self.objects.keys())
         else:
             raise ValueError(f"Unsupported data type {self.DATA_TYPE}.")
+
+        # Ensure host-side DDF mirrors current device state before local edits.
+        cuda.memcpy_dtoh(self.ddf, self.ddf_gpu)
         
         for x in range(int(x_c - radius) - 1, int(x_c + radius) + 1):
             for y in range(int(y_c - radius) - 1, int(y_c + radius) + 1):
                 if (x - x_c) ** 2 + (y - y_c) ** 2 < radius**2:
                     k = x + y * self.FIELD_SHAPE[0]
                     self.flag[k] = SOLID
+                    for i in range(self.LATTICE):
+                        self.ddf[k + i * self.FIELD_SIZE] = self.WW[i]
                     delta_temp = np.zeros(11, dtype=self.DATA_TYPE)
                     delta_temp[0] = id_object.view(self.DATA_TYPE)
                     for i in range(self.LATTICE):
@@ -245,6 +250,8 @@ class FlowField:
         cuda.memcpy_htod(self.delta_gpu, self.delta_curve)
         cuda.memcpy_htod(self.flag_gpu, self.flag)
         cuda.memcpy_htod(self.indx_gpu, self.indx)
+        cuda.memcpy_htod(self.ddf_gpu, self.ddf)
+        cuda.memcpy_htod(self.temp_gpu, self.ddf)
 
         self._rebuild_kernel()
 
diff --git a/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh b/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh
index eb9e21e..4ab49f8 100644
--- a/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh
+++ b/src/CelerisLab/lbm/kernels/boundary/bounce_back.cuh
@@ -43,6 +43,25 @@ __device__ inline void apply_bounce_back(float* __restrict__ f,
 // For nodes at y=1 (adjacent to y=0 wall) or y=NY-2 (adjacent to y=NY-1 wall)
 // ---------------------------------------------------------------------------
 #if NQ == 9
+__device__ inline void apply_wall_bb_d2q9_y_pull(unsigned int y,
+                                                  float* __restrict__ f,
+                                                  const fpxx* __restrict__ fi_in,
+                                                  unsigned long k)
+{
+    if (y == 1) {
+        // Directions sourced from y=0 wall in pull step: +y, +x+y, -x+y.
+        f[3] = load_ddf(fi_in, index_f(k, 4u));
+        f[5] = load_ddf(fi_in, index_f(k, 6u));
+        f[8] = load_ddf(fi_in, index_f(k, 7u));
+    }
+    else if (y == (unsigned int)(NY - 2)) {
+        // Directions sourced from y=NY-1 wall in pull step: -y, -x-y, +x-y.
+        f[4] = load_ddf(fi_in, index_f(k, 3u));
+        f[6] = load_ddf(fi_in, index_f(k, 5u));
+        f[7] = load_ddf(fi_in, index_f(k, 8u));
+    }
+}
+
 __device__ inline void apply_wall_bb_d2q9(unsigned int y,
                                            float* __restrict__ f)
 {
diff --git a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
index 1fde9ca..3b1e284 100644
--- a/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
+++ b/src/CelerisLab/lbm/kernels/boundary/inlet_outlet.cuh
@@ -31,6 +31,14 @@
 #define OUTLET_BLEND_ALPHA 0.70f
 #endif
 
+#ifndef OUTLET_SRT_NEQ_DAMP
+#define OUTLET_SRT_NEQ_DAMP 0.50f
+#endif
+
+#ifndef INLET_TRT_NEQ_DAMP
+#define INLET_TRT_NEQ_DAMP 0.50f
+#endif
+
 __device__ __forceinline__ float inlet_target_u(float y_coord) {
 #if INLET_PROFILE == 0
     return U0;
@@ -69,10 +77,22 @@ __device__ inline void apply_parabolic_inlet(float* __restrict__ f,
     compute_feq(rho_neb, u_target, v_target, feq_tar);
     compute_feq(rho_neb, u_neb, v_neb, feq_neb);
 
-    // Non-equilibrium extrapolation
-    f[1] = f_neb[1] - feq_neb[1] + feq_tar[1];
-    f[5] = f_neb[5] - feq_neb[5] + feq_tar[5];
-    f[7] = f_neb[7] - feq_neb[7] + feq_tar[7];
+#if COLLISION_MODEL == 1
+    // TRT path: reconstruct full population set at inlet using damped donor
+    // non-equilibrium transport. This follows the high-Re stable family that
+    // replaces all boundary-node populations, reducing odd-mode contamination.
+    const float beta = INLET_TRT_NEQ_DAMP;
+    #pragma unroll
+    for (int i = 0; i < 9; i++) {
+        const float fneq = f_neb[i] - feq_neb[i];
+        f[i] = feq_tar[i] + beta * fneq;
+    }
+#else
+    const float beta = 1.0f;
+    f[1] = feq_tar[1] + beta * (f_neb[1] - feq_neb[1]);
+    f[5] = feq_tar[5] + beta * (f_neb[5] - feq_neb[5]);
+    f[7] = feq_tar[7] + beta * (f_neb[7] - feq_neb[7]);
+#endif
 }
 
 // ---------------------------------------------------------------------------
@@ -94,20 +114,30 @@ __device__ inline void apply_pressure_outlet(float* __restrict__ f,
     f[6] = f_neb[6];
 #else
 
-    // Convective non-reflecting style: extrapolate outlet density from neighbor
-    // and keep neighbor velocity.
+    // Prescribed-pressure outlet: keep neighbor velocity and impose reference
+    // outlet density for NEQ reconstruction.
     float rho_neb, u_neb, v_neb;
     compute_rho_u(f_neb, rho_neb, u_neb, v_neb);
 #if OUTLET_BACKFLOW_CLAMP
     u_neb = fmaxf(u_neb, 0.0f);
 #endif
-    float rho_out = rho_neb;
+    float rho_out = RHO;
 
     float feq_tar[9], feq_neb[9];
     compute_feq(rho_out, u_neb, v_neb, feq_tar);
     compute_feq(rho_neb, u_neb, v_neb, feq_neb);
 
-#if OUTLET_MODE == 2
+#if COLLISION_MODEL == 0
+    // SRT path: use full-population damped NEQ reconstruction at outlet to
+    // suppress checkerboard/grid noise from high-frequency donor content.
+    // Reference: high-Re outlet regularization rationale and NEQ decomposition.
+    const float beta = OUTLET_SRT_NEQ_DAMP;
+    #pragma unroll
+    for (int i = 0; i < 9; i++) {
+        const float fneq = f_neb[i] - feq_neb[i];
+        f[i] = feq_tar[i] + beta * fneq;
+    }
+#elif OUTLET_MODE == 2
     const float a = OUTLET_BLEND_ALPHA;
     f[2] = a * (f_neb[2] - feq_neb[2] + feq_tar[2]) + (1.0f - a) * f_neb[2];
     f[8] = a * (f_neb[8] - feq_neb[8] + feq_tar[8]) + (1.0f - a) * f_neb[8];
@@ -178,15 +208,22 @@ __device__ inline void apply_pressure_outlet_3d(float* __restrict__ f,
     un = fmaxf(un, 0.0f);
 #endif
 
-    // Convective non-reflecting style: extrapolate outlet density from neighbor
-    // and keep neighbor velocity.
-    float rho_out = rho_neb;
+    // Prescribed-pressure outlet: keep neighbor velocity and impose reference
+    // outlet density for NEQ reconstruction.
+    float rho_out = RHO;
     float feq_tar[19], feq_neb[19];
     compute_feq(rho_out, un, vn, wn, feq_tar);
     compute_feq(rho_neb, un, vn, wn, feq_neb);
 
     // Reconstruct cx<0 directions: i = 2, 8, 10, 14, 16
-#if OUTLET_MODE == 2
+#if COLLISION_MODEL == 0
+    const float beta = OUTLET_SRT_NEQ_DAMP;
+    #pragma unroll
+    for (int i = 0; i < 19; i++) {
+        const float fneq = f_neb[i] - feq_neb[i];
+        f[i] = feq_tar[i] + beta * fneq;
+    }
+#elif OUTLET_MODE == 2
     const float a = OUTLET_BLEND_ALPHA;
     f[2]  = a * (f_neb[2]  - feq_neb[2]  + feq_tar[2])  + (1.0f - a) * f_neb[2];
     f[8]  = a * (f_neb[8]  - feq_neb[8]  + feq_tar[8])  + (1.0f - a) * f_neb[8];
diff --git a/src/CelerisLab/lbm/kernels/kernel_v2.cu b/src/CelerisLab/lbm/kernels/kernel_v2.cu
index 4c041ab..717e4b2 100644
--- a/src/CelerisLab/lbm/kernels/kernel_v2.cu
+++ b/src/CelerisLab/lbm/kernels/kernel_v2.cu
@@ -91,9 +91,8 @@ __global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi)
     index_from_thread(x, y, k);
     if (x >= (unsigned int)NX || y >= (unsigned int)NY) return;
 
-    float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1))
-                                             * ((float)y - 0.5f * (NY - 1))
-                                / ((float)(NY - 2) * (float)(NY - 2)));
+    // Keep startup field consistent with runtime inlet configuration.
+    float u_init = inlet_target_u((float)y);
 
     if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
         flag[k] = LEGACY_SOLID;
@@ -116,9 +115,8 @@ __global__ void InitTubeFlow_v2(uint8_t* flag, fpxx* fi)
     index_from_thread(x, y, z, k);
     if (x >= (unsigned int)NX || y >= (unsigned int)NY || z >= (unsigned int)NZ) return;
 
-    float u_init = U0 * 1.5f * (1.0f - 4.0f * ((float)y - 0.5f * (NY - 1))
-                                             * ((float)y - 0.5f * (NY - 1))
-                                / ((float)(NY - 2) * (float)(NY - 2)));
+    // Keep startup field consistent with runtime inlet configuration.
+    float u_init = inlet_target_u((float)y);
 
     if (y == 0 || y == NY - 1 || x == 0 || x == NX - 1) {
         flag[k] = LEGACY_SOLID;
@@ -217,9 +215,14 @@ __global__ void OneStep(
         }
     }
 
+    if ((fl & LEGACY_FLUID) && (y == 1u || y == (unsigned int)(NY - 2))) {
+        apply_wall_bb_d2q9_y_pull(y, f, fi_in, k);
+    }
+
     // Collision
     if (fl & LEGACY_FLUID) {
         float feq[NQ], Fin[NQ];
+        compute_rho_u(f, rho_n, ux, uy);
         compute_feq(rho_n, ux, uy, feq);
         zero_forcing(Fin);
         float omega_col = d_params.omega;
@@ -298,9 +301,14 @@ __global__ void OneStep(
         }
     }
 
+    if ((fl & LEGACY_FLUID) && (y == 1 || y == NY - 2)) {
+        apply_wall_bb_d3q19_y_pull(y, f, fi_in, k);
+    }
+
     // ---- Collision (fluid only) ----
     if (fl & LEGACY_FLUID) {
         float feq[NQ], Fin[NQ];
+        compute_rho_u(f, rho_n, ux, uy, uz);
         compute_feq(rho_n, ux, uy, uz, feq);
         zero_forcing(Fin);
         float omega_col = d_params.omega;
diff --git a/src/CelerisLab/lbm/kernels/macros.h b/src/CelerisLab/lbm/kernels/macros.h
index d784f34..0ed1b25 100644
--- a/src/CelerisLab/lbm/kernels/macros.h
+++ b/src/CelerisLab/lbm/kernels/macros.h
@@ -100,4 +100,9 @@
 
 #ifndef OMEGA_COLLISION_MAX
 #define OMEGA_COLLISION_MAX 1.999f
+#endif
+
+// TRT magic parameter Lambda used to map omega+ -> omega-
+#ifndef TRT_MAGIC_PARAM
+#define TRT_MAGIC_PARAM 0.1875f
 #endif
\ No newline at end of file
diff --git a/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh b/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh
index 1d0ad1e..7afdd6b 100644
--- a/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh
+++ b/src/CelerisLab/lbm/kernels/operators/collision_trt.cuh
@@ -13,8 +13,11 @@
 #ifndef CELERIS_OPERATORS_COLLISION_TRT_CUH
 #define CELERIS_OPERATORS_COLLISION_TRT_CUH
 
-// Magic parameter Λ = 3/16 (optimal for porous-media / bounce-back wall location)
+// Magic parameter Λ for TRT wall/collision coupling.
+// Default keeps previous behavior; can be overridden in macros.h / tests.
+#ifndef TRT_MAGIC_PARAM
 #define TRT_MAGIC_PARAM (0.1875f)
+#endif
 
 __device__ __forceinline__ float compute_omega_minus(float omega_plus) {
     return 1.0f / (TRT_MAGIC_PARAM / (1.0f / omega_plus - 0.5f) + 0.5f);
diff --git a/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh b/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
index bc8b00b..afb6080 100644
--- a/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
+++ b/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
@@ -19,8 +19,21 @@
 #define LES_CS 0.16f
 #endif
 
+#ifndef LES_FP_ITERS
+#define LES_FP_ITERS 3
+#endif
+
+#ifndef LES_FP_RELAX
+#define LES_FP_RELAX 0.70f
+#endif
+
+#ifndef LES_NUT_MAX_RATIO
+#define LES_NUT_MAX_RATIO 20.0f
+#endif
+
+
 __device__ __forceinline__ float clamp_omega(float w) {
-    return fminf(1.999f, fmaxf(0.01f, w));
+    return fminf(OMEGA_COLLISION_MAX, fmaxf(OMEGA_COLLISION_MIN, w));
 }
 
 #if NQ == 9
@@ -32,6 +45,7 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
 {
     const float rho_safe = fmaxf(rho, 1.0e-12f);
     const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
+    const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
 
     // Πneq = Σ (ci ci)(fi-feqi)
     float pixx = 0.0f, piyy = 0.0f, pixy = 0.0f;
@@ -46,18 +60,36 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
         pixy += fneq * cx * cy;
     }
 
-    const float denom = 2.0f * rho_safe * CS2 * tau0;
-    const float sxx = -pixx / denom;
-    const float syy = -piyy / denom;
-    const float sxy = -pixy / denom;
+    // Self-consistent tau iteration for LES/TRT coupling.
+    // S is estimated from Pi_neq with current tau_eff, then nu_t updates tau_eff.
+    float tau_eff = tau0;
+    const float tau_max = 1.0f / fmaxf(OMEGA_COLLISION_MIN, 1.0e-6f);
 
-    // |S| = sqrt(2 Sij Sij)
-    const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + 2.0f*sxy*sxy));
+    #pragma unroll
+    for (int it = 0; it < LES_FP_ITERS; ++it) {
+        const float denom = 2.0f * rho_safe * CS2 * fmaxf(tau_eff, 0.500001f);
+        const float sxx = -pixx / denom;
+        const float syy = -piyy / denom;
+        const float sxy = -pixy / denom;
 
-    const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
-    const float nut = (LES_CS * LES_CS) * s_mag;
-    const float nue = nu0 + nut;
-    return clamp_omega(1.0f / (3.0f * nue + 0.5f));
+        // Use deviatoric strain to reduce SGS sensitivity to bulk/pressure-wave content.
+        const float tr = 0.5f * (sxx + syy);
+        const float sxx_dev = sxx - tr;
+        const float syy_dev = syy - tr;
+        const float sxy_dev = sxy;
+
+        // |S_dev| = sqrt(2 S_dev,ij S_dev,ij)
+        const float s_mag = sqrtf(2.0f * (sxx_dev*sxx_dev + syy_dev*syy_dev + 2.0f*sxy_dev*sxy_dev));
+
+        float nut = (LES_CS * LES_CS) * s_mag;
+        const float nut_cap = fmaxf(0.0f, LES_NUT_MAX_RATIO * nu0);
+        nut = fminf(fmaxf(nut, 0.0f), nut_cap);
+
+        const float tau_new = fminf(tau_max, fmaxf(0.500001f, 0.5f + 3.0f * (nu0 + nut)));
+        tau_eff = LES_FP_RELAX * tau_new + (1.0f - LES_FP_RELAX) * tau_eff;
+    }
+
+    return clamp_omega(1.0f / tau_eff);
 }
 
 #elif NQ == 19
@@ -69,6 +101,7 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
 {
     const float rho_safe = fmaxf(rho, 1.0e-12f);
     const float tau0 = fmaxf(1.0f / fmaxf(omega0, 1.0e-6f), 0.500001f);
+    const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
 
     // Πneq = Σ (ci ci)(fi-feqi)
     float pixx = 0.0f, piyy = 0.0f, pizz = 0.0f;
@@ -88,22 +121,42 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
         piyz += fneq * cy * cz;
     }
 
-    const float denom = 2.0f * rho_safe * CS2 * tau0;
-    const float sxx = -pixx / denom;
-    const float syy = -piyy / denom;
-    const float szz = -pizz / denom;
-    const float sxy = -pixy / denom;
-    const float sxz = -pixz / denom;
-    const float syz = -piyz / denom;
+    // Self-consistent tau iteration for LES/TRT coupling.
+    float tau_eff = tau0;
+    const float tau_max = 1.0f / fmaxf(OMEGA_COLLISION_MIN, 1.0e-6f);
 
-    // |S| = sqrt(2 Sij Sij)
-    const float s_mag = sqrtf(2.0f * (sxx*sxx + syy*syy + szz*szz
-                                   + 2.0f*(sxy*sxy + sxz*sxz + syz*syz)));
+    #pragma unroll
+    for (int it = 0; it < LES_FP_ITERS; ++it) {
+        const float denom = 2.0f * rho_safe * CS2 * fmaxf(tau_eff, 0.500001f);
+        const float sxx = -pixx / denom;
+        const float syy = -piyy / denom;
+        const float szz = -pizz / denom;
+        const float sxy = -pixy / denom;
+        const float sxz = -pixz / denom;
+        const float syz = -piyz / denom;
 
-    const float nu0 = (tau0 - 0.5f) * (1.0f / 3.0f);
-    const float nut = (LES_CS * LES_CS) * s_mag;
-    const float nue = nu0 + nut;
-    return clamp_omega(1.0f / (3.0f * nue + 0.5f));
+        // Use deviatoric strain to reduce SGS sensitivity to bulk/pressure-wave content.
+        const float tr = (sxx + syy + szz) / 3.0f;
+        const float sxx_dev = sxx - tr;
+        const float syy_dev = syy - tr;
+        const float szz_dev = szz - tr;
+        const float sxy_dev = sxy;
+        const float sxz_dev = sxz;
+        const float syz_dev = syz;
+
+        // |S_dev| = sqrt(2 S_dev,ij S_dev,ij)
+        const float s_mag = sqrtf(2.0f * (sxx_dev*sxx_dev + syy_dev*syy_dev + szz_dev*szz_dev
+                           + 2.0f*(sxy_dev*sxy_dev + sxz_dev*sxz_dev + syz_dev*syz_dev)));
+
+        float nut = (LES_CS * LES_CS) * s_mag;
+        const float nut_cap = fmaxf(0.0f, LES_NUT_MAX_RATIO * nu0);
+        nut = fminf(fmaxf(nut, 0.0f), nut_cap);
+
+        const float tau_new = fminf(tau_max, fmaxf(0.500001f, 0.5f + 3.0f * (nu0 + nut)));
+        tau_eff = LES_FP_RELAX * tau_new + (1.0f - LES_FP_RELAX) * tau_eff;
+    }
+
+    return clamp_omega(1.0f / tau_eff);
 }
 
 #endif // NQ
diff --git a/src/CelerisLab/lbm/kernels/step/one_step_double.cu b/src/CelerisLab/lbm/kernels/step/one_step_double.cu
index b6cf1a6..8712dca 100644
--- a/src/CelerisLab/lbm/kernels/step/one_step_double.cu
+++ b/src/CelerisLab/lbm/kernels/step/one_step_double.cu
@@ -128,12 +128,12 @@ __global__ void StreamCollideDouble(
 
     // ----- Wall bounce-back (top/bottom walls) -----
 #if NQ == 9
-    if (y == 1 || y == NY - 2) {
-        apply_wall_bb_d2q9(y, f);
+    if ((fl & LEGACY_FLUID) && (y == 1u || y == (unsigned int)(NY - 2))) {
+        apply_wall_bb_d2q9_y_pull(y, f, fi_in, k);
     }
 #elif NQ == 19
-    if (y == 1 || y == NY - 2) {
-        apply_wall_bb_d3q19_y(y, f);
+    if ((fl & LEGACY_FLUID) && (y == 1u || y == (unsigned int)(NY - 2))) {
+        apply_wall_bb_d3q19_y_pull(y, f, fi_in, k);
     }
 #endif
 
@@ -143,6 +143,7 @@ __global__ void StreamCollideDouble(
         float Fin[NQ];
 
 #if NQ == 9
+    compute_rho_u(f, rho_n, ux, uy);
         compute_feq(rho_n, ux, uy, feq);
         zero_forcing(Fin);
     float omega_col = d_params.omega;
@@ -161,6 +162,7 @@ __global__ void StreamCollideDouble(
 #endif
 
 #elif NQ == 19
+    compute_rho_u(f, rho_n, ux, uy, uz);
         compute_feq(rho_n, ux, uy, uz, feq);
         zero_forcing(Fin);
     float omega_col = d_params.omega;
diff --git a/tests/test_d2q9_visual.py b/tests/test_d2q9_visual.py
deleted file mode 100644
index 9f60b66..0000000
--- a/tests/test_d2q9_visual.py
+++ /dev/null
@@ -1,387 +0,0 @@
-#!/usr/bin/env python3
-"""
-D2Q9 Regression Test — Poiseuille Channel + Cylinder Flow
-==========================================================
-Uses kernel_v2.cu by default (legacy kernel.cu remains optional fallback).
-Produces matplotlib figures for visual validation.
-
-Usage:
-    python tests/test_d2q9_visual.py --device 2
-    python tests/test_d2q9_visual.py --device 2 --cylinder
-    python tests/test_d2q9_visual.py --device 2 --legacy
-"""
-
-import sys, os, argparse, time
-
-# Ensure CelerisLab package is importable
-sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src'))
-
-import numpy as np
-import pycuda.driver as cuda
-
-import matplotlib
-matplotlib.use('Agg')
-import matplotlib.pyplot as plt
-from matplotlib.colors import Normalize
-
-from CelerisLab.cuda import compiler
-from CelerisLab.common import preprocess as preproc
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Configuration  ━━━━━━━━━━━━━━━━━━━━━━━
-NX, NY = 1280, 512
-NQ = 9
-NT = 128
-DIM = 2
-VIS = 0.002
-U0 = 0.01
-RHO = 1.0
-TOTAL = NX * NY
-
-# Original direction vectors (const.h ordering)
-E = np.array([[0,0],[1,0],[0,1],[-1,0],[0,-1],[1,1],[-1,1],[-1,-1],[1,-1]], dtype=np.int32)
-OPP = np.array([0, 3, 4, 1, 2, 7, 8, 5, 6], dtype=np.int32)
-
-FLUID_FLAG = 0b00000001
-SOLID_FLAG = 0b00000010
-INTERFACE_FLAG = 0b00001000
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Helpers  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
-
-def configure_macros(n_objs=0):
-    """Write macros.h for D2Q9 regression (v2 default, legacy optional)."""
-    lines = compiler.read_lines(compiler.kernel_path("macros.h"))
-    defs = {
-        'MULT_GPU': 'False', 'NT': NT,
-        'X_1U': 128, 'Y_1U': 32, 'Z_1U': 1,
-        'LBtype': 'float',
-        'UX': 10, 'UY': 16, 'UZ': 1,
-        'NX': NX, 'NY': NY, 'NZ': 1,
-        'DIM': DIM, 'NQ': NQ,
-        'VIS': VIS, 'RHO': f'{RHO}', 'U0': U0,
-        'N_OBJS': n_objs,
-        'COLLISION_MODEL': 0,     # SRT
-        'STREAMING_MODEL': 0,     # double-buffer
-        'STORE_PRECISION': 0,     # FP32
-        'USE_DDF_SHIFTING': 0,    # keep unshifted for v2 defaults
-    }
-    for name, val in defs.items():
-        lines = compiler.modify_macro(lines, name, val)
-    compiler.write_lines(compiler.kernel_path("macros.h"), lines)
-
-
-def extract_fields(ddf_host, use_ddf_shifting=False):
-    """Compute rho, u, v from host DDF in original D2Q9 direction ordering."""
-    f = ddf_host.reshape(NQ, NY, NX)
-    if use_ddf_shifting:
-        rho = np.sum(f, axis=0) + RHO
-        denom = np.full_like(rho, RHO)
-    else:
-        rho = np.sum(f, axis=0)
-        denom = rho
-
-    denom_safe = np.where(np.abs(denom) > 1e-12, denom, 1.0)
-    u = (f[1] + f[5] + f[8] - f[3] - f[6] - f[7]) / denom_safe
-    v = (f[2] + f[5] + f[6] - f[4] - f[7] - f[8]) / denom_safe
-    u = np.where(np.abs(denom) > 1e-12, u, 0.0)
-    v = np.where(np.abs(denom) > 1e-12, v, 0.0)
-    return rho, u, v
-
-
-def analytical_poiseuille(y_arr):
-    """Analytical parabolic profile matching InitTubeFlow."""
-    yy = (y_arr - 0.5 * (NY - 1)) / (NY - 2.0)
-    return U0 * 1.5 * (1 - 4 * yy**2)
-
-
-def build_cylinder_data(cx, cy, radius):
-    """Replicate driver.py add_cylinder logic for flag / delta / indx."""
-    flag = np.ones(TOTAL, dtype=np.uint8)   # init all FLUID
-    indx = np.zeros(TOTAL, dtype=np.int32)
-    delta_list = []
-    index_offset = 0
-
-    # Build Poiseuille flag first (walls + solid borders)
-    for y in range(NY):
-        for x in range(NX):
-            k = x + y * NX
-            if y == 0 or y == NY - 1 or x == 0 or x == NX - 1:
-                flag[k] = SOLID_FLAG
-
-    # Add cylinder
-    for x in range(int(cx - radius) - 1, int(cx + radius) + 1):
-        for y in range(int(cy - radius) - 1, int(cy + radius) + 1):
-            if (x - cx)**2 + (y - cy)**2 < radius**2:
-                k = x + y * NX
-                flag[k] = SOLID_FLAG
-                dt = np.zeros(11, dtype=np.float32)
-                dt[0] = np.int32(0).view(np.float32)  # id_object = 0
-                has_interface = False
-                for i in range(NQ):
-                    xn = x + E[i][0]
-                    yn = y + E[i][1]
-                    if (xn - cx)**2 + (yn - cy)**2 >= radius**2:
-                        has_interface = True
-                        xi, yi = preproc.find_circle_intersection(
-                            x, y, xn, yn, cx, cy, radius)
-                        d_neb = np.sqrt((xi - xn)**2 + (yi - yn)**2)
-                        e_len = np.sqrt(E[i][0]**2 + E[i][1]**2)
-                        if e_len > 0:
-                            dt[i] = d_neb / e_len
-                if has_interface:
-                    flag[k] |= INTERFACE_FLAG
-                    dt[9] = (cy - y) / radius
-                    dt[10] = (x - cx) / radius
-                    indx[k] = index_offset
-                    delta_list.append(dt)
-                    index_offset += 11
-
-    delta = np.concatenate(delta_list) if delta_list else np.zeros(1, dtype=np.float32)
-    return flag, indx, delta
-
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Simulation  ━━━━━━━━━━━━━━━━━━━━━━━━━━
-
-def run_simulation(device_id, n_steps, n_objs, flag_host, indx_host, delta_host, use_legacy=False):
-    """Compile kernel, run LBM, return DDF on host."""
-    cuda.init()
-    dev = cuda.Device(device_id)
-    ctx = dev.make_context()
-    print(f"[GPU {device_id}] {dev.name()}")
-
-    try:
-        configure_macros(n_objs)
-        if use_legacy:
-            compiler.compile_kernel()
-            ptx_path = compiler.kernel_path("kernel.ptx")
-            init_name = "InitTubeFlow"
-        else:
-            compiler.compile_kernel_v2()
-            ptx_path = compiler.kernel_path("kernel_v2.ptx")
-            init_name = "InitTubeFlow_v2"
-        mod = cuda.module_from_file(ptx_path)
-        step_fn = mod.get_function("OneStep")
-        init_fn = mod.get_function(init_name)
-
-        nbytes_ddf = TOTAL * NQ * 4
-        ddf_gpu   = cuda.mem_alloc(nbytes_ddf)
-        temp_gpu  = cuda.mem_alloc(nbytes_ddf)
-        flag_gpu  = cuda.mem_alloc(flag_host.nbytes)
-        indx_gpu  = cuda.mem_alloc(indx_host.nbytes)
-        delta_gpu = cuda.mem_alloc(max(delta_host.nbytes, 4))
-
-        action_host = np.zeros(max(n_objs, 1), dtype=np.float32)
-        obs_host    = np.zeros(max(n_objs * DIM, 1), dtype=np.float32)
-        action_gpu = cuda.mem_alloc(action_host.nbytes)
-        obs_gpu    = cuda.mem_alloc(obs_host.nbytes)
-
-        cuda.memcpy_htod(action_gpu, action_host)
-        cuda.memcpy_htod(obs_gpu, obs_host)
-
-        block = (NT, 1, 1)
-        grid  = (NX // NT, NY, 1)
-
-        # Init Poiseuille
-        init_fn(flag_gpu, ddf_gpu, block=block, grid=grid)
-        ctx.synchronize()
-
-        # Overwrite flag / indx / delta for cylinder case
-        cuda.memcpy_htod(flag_gpu, flag_host)
-        cuda.memcpy_htod(indx_gpu, indx_host)
-        cuda.memcpy_htod(delta_gpu, delta_host)
-
-        # Step loop
-        t0 = time.time()
-        for i in range(n_steps):
-            step_fn(flag_gpu, ddf_gpu, temp_gpu, indx_gpu, delta_gpu,
-                    action_gpu, obs_gpu, block=block, grid=grid)
-            ddf_gpu, temp_gpu = temp_gpu, ddf_gpu
-        ctx.synchronize()
-        dt = time.time() - t0
-        mlups = TOTAL * n_steps / dt / 1e6
-        print(f"  {n_steps} steps in {dt:.2f}s  ({mlups:.1f} MLUPS)")
-
-        # Copy back
-        ddf = np.zeros(TOTAL * NQ, dtype=np.float32)
-        cuda.memcpy_dtoh(ddf, ddf_gpu)
-        flag_out = np.zeros(TOTAL, dtype=np.uint8)
-        cuda.memcpy_dtoh(flag_out, flag_gpu)
-
-        return ddf, flag_out
-    finally:
-        ctx.pop()
-
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Visualization  ━━━━━━━━━━━━━━━━━━━━━━━
-
-def plot_poiseuille(ddf, flag, out_path, use_ddf_shifting=False):
-    """3-panel figure: velocity mag, u(y) profile, pressure along centerline."""
-    rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
-    vel_mag = np.sqrt(u**2 + v**2)
-
-    # Mask solid cells for display
-    mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool)
-    vel_mag_masked = np.ma.array(vel_mag, mask=mask)
-
-    fig, axes = plt.subplots(1, 3, figsize=(18, 5))
-
-    # (a) Velocity magnitude heatmap
-    ax = axes[0]
-    im = ax.imshow(vel_mag_masked, origin='lower', aspect='auto',
-                   cmap='jet', extent=[0, NX, 0, NY])
-    plt.colorbar(im, ax=ax, label='|u|')
-    ax.set_title('Velocity magnitude')
-    ax.set_xlabel('x'); ax.set_ylabel('y')
-
-    # (b) u(y) profile at x = NX/2  vs. analytical
-    ax = axes[1]
-    x_mid = NX // 2
-    y_arr = np.arange(NY, dtype=float)
-    ax.plot(u[:, x_mid], y_arr, 'b-', lw=2, label='LBM')
-    ax.plot(analytical_poiseuille(y_arr), y_arr, 'r--', lw=1.5, label='Analytical')
-    ax.set_xlabel('u_x'); ax.set_ylabel('y')
-    ax.set_title(f'u(y) at x={x_mid}')
-    ax.legend()
-    ax.grid(True, alpha=0.3)
-
-    # (c) Pressure along centerline y = NY/2
-    ax = axes[2]
-    y_mid = NY // 2
-    p = rho / 3.0
-    ax.plot(np.arange(NX), p[y_mid, :], 'g-', lw=1.5)
-    ax.set_xlabel('x'); ax.set_ylabel('p = ρ/3')
-    ax.set_title(f'Pressure along centerline (y={y_mid})')
-    ax.grid(True, alpha=0.3)
-
-    fig.suptitle(f'D2Q9 Poiseuille – NX={NX}, NY={NY}, VIS={VIS}, U0={U0}', fontsize=13)
-    fig.tight_layout()
-    fig.savefig(out_path, dpi=150)
-    print(f"  Saved: {out_path}")
-    plt.close(fig)
-
-
-def plot_cylinder(ddf, flag, cx, cy, radius, out_path, use_ddf_shifting=False):
-    """3-panel figure: velocity magnitude (zoom), vorticity, streamlines."""
-    rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
-    vel_mag = np.sqrt(u**2 + v**2)
-    mask = (flag.reshape(NY, NX) & SOLID_FLAG).astype(bool)
-
-    # Zoom window around cylinder
-    pad = int(radius * 8)
-    x0 = max(int(cx - pad), 0)
-    x1 = min(int(cx + pad * 2), NX)
-    y0 = max(int(cy - pad), 0)
-    y1 = min(int(cy + pad), NY)
-
-    fig, axes = plt.subplots(1, 3, figsize=(20, 6))
-
-    # (a) Velocity magnitude (zoomed)
-    ax = axes[0]
-    vm_z = np.ma.array(vel_mag[y0:y1, x0:x1], mask=mask[y0:y1, x0:x1])
-    im = ax.imshow(vm_z, origin='lower', aspect='equal',
-                   cmap='jet', extent=[x0, x1, y0, y1])
-    circ = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7)
-    ax.add_patch(circ)
-    plt.colorbar(im, ax=ax, label='|u|')
-    ax.set_title('Velocity magnitude')
-
-    # (b) Vorticity
-    ax = axes[1]
-    dvdx = np.gradient(v, axis=1)
-    dudy = np.gradient(u, axis=0)
-    omega = dvdx - dudy
-    om_z = np.ma.array(omega[y0:y1, x0:x1], mask=mask[y0:y1, x0:x1])
-    vmax = np.percentile(np.abs(omega[~mask]), 99)
-    im = ax.imshow(om_z, origin='lower', aspect='equal',
-                   cmap='RdBu_r', extent=[x0, x1, y0, y1],
-                   vmin=-vmax, vmax=vmax)
-    circ2 = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7)
-    ax.add_patch(circ2)
-    plt.colorbar(im, ax=ax, label='ω')
-    ax.set_title('Vorticity')
-
-    # (c) Streamlines
-    ax = axes[2]
-    X, Y = np.meshgrid(np.arange(x0, x1), np.arange(y0, y1))
-    u_z = u[y0:y1, x0:x1].copy()
-    v_z = v[y0:y1, x0:x1].copy()
-    u_z[mask[y0:y1, x0:x1]] = 0
-    v_z[mask[y0:y1, x0:x1]] = 0
-    speed = np.sqrt(u_z**2 + v_z**2)
-    ax.streamplot(X, Y, u_z, v_z, color=speed, cmap='jet',
-                  density=2.0, linewidth=0.8)
-    circ3 = plt.Circle((cx, cy), radius, fill=True, color='gray', alpha=0.7)
-    ax.add_patch(circ3)
-    ax.set_xlim(x0, x1); ax.set_ylim(y0, y1)
-    ax.set_aspect('equal')
-    ax.set_title('Streamlines')
-
-    fig.suptitle(f'D2Q9 Cylinder Flow – Re_D={U0*1.5*2*radius/VIS:.0f}, D={2*radius}', fontsize=13)
-    fig.tight_layout()
-    fig.savefig(out_path, dpi=150)
-    print(f"  Saved: {out_path}")
-    plt.close(fig)
-
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Main  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
-
-def main():
-    parser = argparse.ArgumentParser(description='D2Q9 Regression Test')
-    parser.add_argument('--device', type=int, default=2,
-                        help='CUDA device ID (default: 2)')
-    parser.add_argument('--legacy', action='store_true',
-                        help='Use legacy kernel.cu path (default uses kernel_v2.cu)')
-    parser.add_argument('--cylinder', action='store_true',
-                        help='Also run cylinder flow test')
-    parser.add_argument('--steps-pois', type=int, default=5000,
-                        help='Steps for Poiseuille (default: 5000)')
-    parser.add_argument('--steps-cyl', type=int, default=30000,
-                        help='Steps for cylinder (default: 30000)')
-    args = parser.parse_args()
-
-    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output')
-    os.makedirs(out_dir, exist_ok=True)
-
-    use_ddf_shifting = bool(args.legacy)
-    mode = 'legacy kernel.cu' if args.legacy else 'kernel_v2.cu'
-    print(f"\n[Mode] {mode}")
-
-    # ---- Test 1: Poiseuille ----
-    print("\n===== Test 1: Poiseuille Channel Flow =====")
-    flag_pois = np.ones(TOTAL, dtype=np.uint8)
-    indx_pois = np.zeros(TOTAL, dtype=np.int32)
-    delta_pois = np.zeros(1, dtype=np.float32)
-
-    ddf, flag = run_simulation(args.device, args.steps_pois, 0,
-                               flag_pois, indx_pois, delta_pois,
-                               use_legacy=args.legacy)
-    plot_poiseuille(ddf, flag, os.path.join(out_dir, 'poiseuille_d2q9.png'),
-                    use_ddf_shifting=use_ddf_shifting)
-
-    # Error metric
-    rho, u, v = extract_fields(ddf, use_ddf_shifting=use_ddf_shifting)
-    y_arr = np.arange(NY, dtype=float)
-    u_ana = analytical_poiseuille(y_arr)
-    x_mid = NX // 2
-    u_num = u[:, x_mid]
-    # Interior cells only (skip walls)
-    err = np.max(np.abs(u_num[2:-2] - u_ana[2:-2])) / np.max(np.abs(u_ana[2:-2]))
-    print(f"  L∞ relative error at x={x_mid}: {err:.2e}")
-
-    # ---- Test 2: Cylinder ----
-    if args.cylinder:
-        print("\n===== Test 2: Flow Around Cylinder =====")
-        cyl_cx, cyl_cy, cyl_r = 256.0, 256.0, 32.0
-        flag_cyl, indx_cyl, delta_cyl = build_cylinder_data(cyl_cx, cyl_cy, cyl_r)
-
-        ddf2, flag2 = run_simulation(args.device, args.steps_cyl, 1,
-                         flag_cyl, indx_cyl, delta_cyl,
-                         use_legacy=args.legacy)
-        plot_cylinder(ddf2, flag2, cyl_cx, cyl_cy, cyl_r,
-                  os.path.join(out_dir, 'cylinder_d2q9.png'),
-                  use_ddf_shifting=use_ddf_shifting)
-
-    print("\nDone.")
-
-
-if __name__ == '__main__':
-    main()
diff --git a/tests/test_d3q19_cylinder.py b/tests/test_d3q19_cylinder.py
deleted file mode 100644
index bc69c9a..0000000
--- a/tests/test_d3q19_cylinder.py
+++ /dev/null
@@ -1,327 +0,0 @@
-#!/usr/bin/env python3
-"""
-D3Q19 SRT — Cylinder Wake Flow (Periodic Z)
-=============================================
-Tests the v2 modular kernel (kernel_v2.cu) in 3D.
-Cylinder axis along z, parabolic inlet, pressure outlet, no-slip y-walls.
-Produces cross-section visualizations at z=NZ/2.
-
-Usage:
-    python tests/test_d3q19_cylinder.py --device 3
-    python tests/test_d3q19_cylinder.py --device 3 --re 200 --steps 50000
-"""
-
-import sys, os, argparse, time, struct
-
-sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'src'))
-
-import numpy as np
-import pycuda.driver as cuda
-
-import matplotlib
-matplotlib.use('Agg')
-import matplotlib.pyplot as plt
-
-from CelerisLab.cuda import compiler
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Configuration  ━━━━━━━━━━━━━━━━━━━━━━━
-# Reasonable 3D grid — fits in < 500 MB GPU memory
-NX, NY, NZ = 256, 128, 32
-NQ  = 19
-NT  = 128
-DIM = 3
-RHO = 1.0
-CYL_CX, CYL_CY = 64.0, 64.0   # Cylinder center (x,y)
-CYL_R = 12.0                   # Cylinder radius
-TOTAL = NX * NY * NZ
-
-# D3Q19 paired direction ordering (from descriptors.cuh)
-#   0:rest  (1,2)±x  (3,4)±y  (5,6)±z
-#   (7,8)±(x+y) (9,10)±(x+z) (11,12)±(y+z)
-#   (13,14)±(x-y) (15,16)±(x-z) (17,18)±(y-z)
-CX = np.array([0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0], dtype=np.int32)
-CY = np.array([0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1], dtype=np.int32)
-CZ = np.array([0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1], dtype=np.int32)
-W  = np.array([1/3] + [1/18]*6 + [1/36]*12, dtype=np.float32)
-
-FLUID_FLAG = 0x01
-SOLID_FLAG = 0x02
-OBSTACLE_FLAG = 0x04   # triggers half-way BB at adjacent fluid nodes
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Helpers  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
-
-def compute_vis_omega(Re, D, U0):
-    """Compute viscosity and omega from target Reynolds number."""
-    vis = U0 * D / Re
-    omega = 1.0 / (3.0 * vis + 0.5)
-    return vis, omega
-
-
-def configure_macros_3d(vis, u0, n_objs=0):
-    """Write macros.h for D3Q19 SRT."""
-    lines = compiler.read_lines(compiler.kernel_path("macros.h"))
-    defs = {
-        'MULT_GPU': 'False', 'NT': NT,
-        'X_1U': NX, 'Y_1U': NY, 'Z_1U': NZ,
-        'LBtype': 'float',
-        'UX': 1, 'UY': 1, 'UZ': 1,
-        'NX': NX, 'NY': NY, 'NZ': NZ,
-        'DIM': DIM, 'NQ': NQ,
-        'VIS': f'{vis:.10f}', 'RHO': f'{RHO}', 'U0': u0,
-        'N_OBJS': n_objs,
-        'COLLISION_MODEL': 0,     # SRT
-        'STREAMING_MODEL': 0,     # double-buffer
-        'STORE_PRECISION': 0,     # FP32
-        'USE_DDF_SHIFTING': 0,
-    }
-    for name, val in defs.items():
-        lines = compiler.modify_macro(lines, name, val)
-    compiler.write_lines(compiler.kernel_path("macros.h"), lines)
-
-
-def build_cylinder_3d(cx, cy, radius):
-    """Build flag array for 3D cylinder (axis along z, periodic z)."""
-    flag = np.ones(TOTAL, dtype=np.uint8) * FLUID_FLAG
-    for z in range(NZ):
-        for y in range(NY):
-            for x in range(NX):
-                k = z * NY * NX + y * NX + x
-                # Channel walls
-                if y == 0 or y == NY - 1 or x == 0 or x == NX - 1:
-                    flag[k] = SOLID_FLAG
-                # Cylinder body (obstacle — triggers BB at fluid neighbors)
-                elif (x - cx)**2 + (y - cy)**2 < radius**2:
-                    flag[k] = OBSTACLE_FLAG
-    return flag
-
-
-def extract_fields_3d(ddf_host, z_slice):
-    """Extract rho, u, v, w at a given z-slice from D3Q19 DDF (v2 ordering)."""
-    # DDF layout: f[i * TOTAL + k]  where k = z*NY*NX + y*NX + x
-    f = ddf_host.reshape(NQ, NZ, NY, NX)
-    fz = f[:, z_slice, :, :]  # shape (NQ, NY, NX)
-    rho = np.sum(fz, axis=0)
-    ux = np.zeros_like(rho)
-    uy = np.zeros_like(rho)
-    uz_field = np.zeros_like(rho)
-    for i in range(NQ):
-        ux += CX[i] * fz[i]
-        uy += CY[i] * fz[i]
-        uz_field += CZ[i] * fz[i]
-    ux /= rho
-    uy /= rho
-    uz_field /= rho
-    return rho, ux, uy, uz_field
-
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Simulation  ━━━━━━━━━━━━━━━━━━━━━━━━━━
-
-def run_d3q19(device_id, n_steps, vis, u0, flag_host):
-    """Compile v2 kernel, run D3Q19 SRT, return DDF."""
-    omega = 1.0 / (3.0 * vis + 0.5)
-    cuda.init()
-    dev = cuda.Device(device_id)
-    ctx = dev.make_context()
-    print(f"[GPU {device_id}] {dev.name()}")
-
-    try:
-        configure_macros_3d(vis, u0)
-        compiler.compile_kernel_v2()
-        ptx_path = compiler.kernel_path("kernel_v2.ptx")
-        mod = cuda.module_from_file(ptx_path)
-
-        # Get kernels (extern "C" entries from kernel_v2.cu)
-        init_fn = mod.get_function("InitTubeFlow_v2")
-        step_fn = mod.get_function("OneStep")
-
-        # Set d_params.omega via __constant__ memory
-        params_ptr, params_size = mod.get_global("d_params")
-        # LBMParams struct layout (see params.cuh):
-        # Nx(4) Ny(4) Nz(4) N(8) omega(4) omega_bulk(4) fx(4) fy(4) fz(4)
-        # rho_ref(4) u_inlet(4) n_objects(4)
-        # Pack: unsigned int Nx, Ny, Nz; unsigned long N; float omega, omega_bulk, fx, fy, fz, rho_ref, u_inlet; unsigned int n_objects
-        params_data = struct.pack('IIIQfffffffI',
-                                  NX, NY, NZ,
-                                  TOTAL,
-                                  omega, 0.0,      # omega, omega_bulk
-                                  0.0, 0.0, 0.0,   # fx, fy, fz
-                                  RHO, u0,          # rho_ref, u_inlet
-                                  0)                 # n_objects
-        # Pad to match struct size
-        if len(params_data) < params_size:
-            params_data += b'\x00' * (params_size - len(params_data))
-        cuda.memcpy_htod(params_ptr, params_data)
-
-        # Allocate
-        nbytes_ddf = TOTAL * NQ * 4
-        ddf_gpu  = cuda.mem_alloc(nbytes_ddf)
-        temp_gpu = cuda.mem_alloc(nbytes_ddf)
-        flag_gpu = cuda.mem_alloc(flag_host.nbytes)
-        indx_gpu = cuda.mem_alloc(TOTAL * 4)
-        delta_gpu = cuda.mem_alloc(4)
-        action_gpu = cuda.mem_alloc(4)
-        obs_gpu = cuda.mem_alloc(4)
-
-        # Dummy arrays
-        cuda.memset_d32(indx_gpu, 0, TOTAL)
-        cuda.memset_d32(delta_gpu, 0, 1)
-        cuda.memset_d32(action_gpu, 0, 1)
-        cuda.memset_d32(obs_gpu, 0, 1)
-
-        block = (NT, 1, 1)
-        grid  = (NX // NT, NY, NZ)
-
-        # Initialize parabolic flow
-        init_fn(flag_gpu, ddf_gpu, block=block, grid=grid)
-        ctx.synchronize()
-
-        # Overwrite flags with cylinder geometry
-        cuda.memcpy_htod(flag_gpu, flag_host)
-
-        # Step loop
-        print(f"  Running {n_steps} steps (NX={NX}, NY={NY}, NZ={NZ}, omega={omega:.4f})...")
-        t0 = time.time()
-        for i in range(n_steps):
-            step_fn(flag_gpu, ddf_gpu, temp_gpu, indx_gpu, delta_gpu,
-                    action_gpu, obs_gpu,
-                    block=block, grid=grid)
-            ddf_gpu, temp_gpu = temp_gpu, ddf_gpu
-            if (i + 1) % 5000 == 0:
-                ctx.synchronize()
-                elapsed = time.time() - t0
-                mlups = TOTAL * (i + 1) / elapsed / 1e6
-                print(f"    step {i+1}/{n_steps}  ({mlups:.1f} MLUPS)")
-        ctx.synchronize()
-        dt = time.time() - t0
-        mlups = TOTAL * n_steps / dt / 1e6
-        print(f"  Done: {dt:.1f}s, {mlups:.1f} MLUPS")
-
-        # Copy back
-        ddf = np.zeros(TOTAL * NQ, dtype=np.float32)
-        cuda.memcpy_dtoh(ddf, ddf_gpu)
-        return ddf
-
-    finally:
-        ctx.pop()
-
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Visualization  ━━━━━━━━━━━━━━━━━━━━━━━
-
-def plot_d3q19_cylinder(ddf, flag, Re, u0, out_path):
-    """4-panel figure at z=NZ/2: vel-mag, vorticity, streamlines, u(y) profile."""
-    z_mid = NZ // 2
-    rho, ux, uy, uz = extract_fields_3d(ddf, z_mid)
-    vel_mag = np.sqrt(ux**2 + uy**2 + uz**2)
-
-    mask2d = ((flag.reshape(NZ, NY, NX)[z_mid] & (SOLID_FLAG | OBSTACLE_FLAG)) != 0)
-    vel_masked = np.ma.array(vel_mag, mask=mask2d)
-
-    fig, axes = plt.subplots(2, 2, figsize=(16, 10))
-
-    # (a) Velocity magnitude
-    ax = axes[0, 0]
-    im = ax.imshow(vel_masked, origin='lower', aspect='auto',
-                   cmap='jet', extent=[0, NX, 0, NY])
-    circ = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7)
-    ax.add_patch(circ)
-    plt.colorbar(im, ax=ax, label='|u|')
-    ax.set_title(f'Velocity magnitude (z={z_mid})')
-    ax.set_xlabel('x'); ax.set_ylabel('y')
-
-    # (b) Vorticity  ω_z = ∂v/∂x − ∂u/∂y
-    ax = axes[0, 1]
-    dvdx = np.gradient(uy, axis=1)
-    dudy = np.gradient(ux, axis=0)
-    omega_z = dvdx - dudy
-    om_masked = np.ma.array(omega_z, mask=mask2d)
-    vmax = np.percentile(np.abs(omega_z[~mask2d]), 99) if np.any(~mask2d) else 1e-3
-    im = ax.imshow(om_masked, origin='lower', aspect='auto',
-                   cmap='RdBu_r', extent=[0, NX, 0, NY],
-                   vmin=-vmax, vmax=vmax)
-    circ2 = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7)
-    ax.add_patch(circ2)
-    plt.colorbar(im, ax=ax, label='ω_z')
-    ax.set_title('Vorticity ω_z')
-    ax.set_xlabel('x'); ax.set_ylabel('y')
-
-    # (c) Streamlines
-    ax = axes[1, 0]
-    X, Y = np.meshgrid(np.arange(NX), np.arange(NY))
-    ux_c = ux.copy(); ux_c[mask2d] = 0
-    uy_c = uy.copy(); uy_c[mask2d] = 0
-    speed = np.sqrt(ux_c**2 + uy_c**2)
-    ax.streamplot(X, Y, ux_c, uy_c, color=speed, cmap='jet',
-                  density=2.5, linewidth=0.7)
-    circ3 = plt.Circle((CYL_CX, CYL_CY), CYL_R, fill=True, color='gray', alpha=0.7)
-    ax.add_patch(circ3)
-    ax.set_xlim(0, NX); ax.set_ylim(0, NY)
-    ax.set_aspect('auto')
-    ax.set_title('Streamlines')
-    ax.set_xlabel('x'); ax.set_ylabel('y')
-
-    # (d) u_x(y) profiles at different x stations
-    ax = axes[1, 1]
-    y_arr = np.arange(NY)
-    # Analytical parabolic inlet
-    yy = (y_arr - 0.5 * (NY - 1)) / (NY - 2.0)
-    u_ana = u0 * 1.5 * (1 - 4 * yy**2)
-
-    x_stations = [NX // 8, NX // 4, NX // 2, 3 * NX // 4]
-    for xs in x_stations:
-        ax.plot(ux[:, xs], y_arr, label=f'x={xs}')
-    ax.plot(u_ana, y_arr, 'k--', lw=1.5, label='Analytical inlet')
-    ax.set_xlabel('u_x'); ax.set_ylabel('y')
-    ax.set_title('u_x(y) profiles')
-    ax.legend(fontsize=8)
-    ax.grid(True, alpha=0.3)
-
-    fig.suptitle(f'D3Q19 SRT Cylinder — Re={Re:.0f}, D={2*CYL_R:.0f}, '
-                 f'Grid={NX}×{NY}×{NZ}', fontsize=13)
-    fig.tight_layout()
-    fig.savefig(out_path, dpi=150)
-    print(f"  Saved: {out_path}")
-    plt.close(fig)
-
-
-# ━━━━━━━━━━━━━━━━━━━━━━━  Main  ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
-
-def main():
-    parser = argparse.ArgumentParser(description='D3Q19 SRT Cylinder Flow')
-    parser.add_argument('--device', type=int, default=0,
-                        help='CUDA device ID (default: 0)')
-    parser.add_argument('--re', type=float, default=100.0,
-                        help='Reynolds number based on diameter (default: 100)')
-    parser.add_argument('--u0', type=float, default=0.04,
-                        help='Inlet characteristic velocity (default: 0.04)')
-    parser.add_argument('--steps', type=int, default=30000,
-                        help='Number of LBM steps (default: 30000)')
-    args = parser.parse_args()
-
-    out_dir = os.path.join(os.path.dirname(os.path.abspath(__file__)), '..', 'output')
-    os.makedirs(out_dir, exist_ok=True)
-
-    D = 2 * CYL_R
-    vis, omega = compute_vis_omega(args.re, D, args.u0)
-    print(f"\n===== D3Q19 SRT Cylinder Flow =====")
-    print(f"  Re = {args.re:.0f},  D = {D:.0f},  U0 = {args.u0}")
-    print(f"  ν  = {vis:.6f},  ω = {omega:.4f}")
-    if omega > 1.95:
-        print(f"  WARNING: omega={omega:.4f} is close to 2.0, stability may be poor.")
-        print(f"  Consider reducing U0 or Re.")
-
-    flag = build_cylinder_3d(CYL_CX, CYL_CY, CYL_R)
-    n_solid = np.sum(flag == SOLID_FLAG)
-    n_fluid = np.sum(flag == FLUID_FLAG)
-    print(f"  Grid: {NX}×{NY}×{NZ} = {TOTAL} cells  (fluid: {n_fluid}, solid: {n_solid})")
-    print(f"  Memory: ~{2 * TOTAL * NQ * 4 / 1e6:.0f} MB for double-buffer DDF")
-
-    ddf = run_d3q19(args.device, args.steps, vis, args.u0, flag)
-
-    plot_d3q19_cylinder(ddf, flag, args.re, args.u0,
-                        os.path.join(out_dir, f'cylinder_d3q19_re{int(args.re)}.png'))
-
-    print("\nDone.")
-
-
-if __name__ == '__main__':
-    main()
diff --git a/tests/test_high_re_validation.py b/tests/test_high_re_validation.py
index 9b6fcd3..8920894 100644
--- a/tests/test_high_re_validation.py
+++ b/tests/test_high_re_validation.py
@@ -33,6 +33,11 @@ FLUID_FLAG = 0x01
 SOLID_FLAG = 0x02
 OBSTACLE_FLAG = 0x04
 
+OMEGA_COLLISION_MIN_DEFAULT = 0.01
+LES_POST_FP_ITERS = 3
+LES_POST_FP_RELAX = 0.70
+LES_POST_NUT_MAX_RATIO = 20.0
+
 
 def collision_name(model):
     return {0: "SRT", 1: "TRT", 2: "MRT"}.get(model, f"M{model}")
@@ -40,10 +45,13 @@ def collision_name(model):
 
 def make_case_tag(cfg):
     les_tag = "LES" if cfg["use_les"] else "NoLES"
+    inlet_tag = f"IP{int(cfg.get('inlet_profile', 0))}"
+    lam_tag = f"LAM{float(cfg.get('trt_magic_param', 0.1875)):.4f}"
     return (
         f"{cfg['name']}_Re{int(cfg['target_re'])}_"
         f"{collision_name(cfg['collision_model'])}_{les_tag}_"
-        f"OM{int(cfg['outlet_mode'])}_WMAX{cfg['omega_collision_max']:.3f}"
+        f"OM{int(cfg['outlet_mode'])}_{inlet_tag}_{lam_tag}_"
+        f"WMAX{cfg['omega_collision_max']:.3f}"
     )
 
 
@@ -61,6 +69,366 @@ def validate_case(rho):
     return True, "OK"
 
 
+def lattice_weights(nq):
+    if nq == 9:
+        return np.array(
+            [4.0 / 9.0] + [1.0 / 9.0] * 4 + [1.0 / 36.0] * 4,
+            dtype=np.float32,
+        )
+    if nq == 19:
+        return np.array(
+            [1.0 / 3.0] + [1.0 / 18.0] * 6 + [1.0 / 36.0] * 12,
+            dtype=np.float32,
+        )
+    raise ValueError(f"Unsupported nq={nq}")
+
+
+def impose_rest_state_on_nonfluid(cfg, host_ddf):
+    nq = cfg["nq"]
+    nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
+    w = lattice_weights(nq)
+
+    if nq == 9:
+        f = host_ddf.reshape(nq, ny, nx)
+        nonfluid = cfg["flag"].reshape(ny, nx) != FLUID_FLAG
+        for i in range(nq):
+            f[i, nonfluid] = w[i]
+        return host_ddf
+
+    f = host_ddf.reshape(nq, nz, ny, nx)
+    nonfluid = cfg["flag"].reshape(nz, ny, nx) != FLUID_FLAG
+    for i in range(nq):
+        f[i, nonfluid] = w[i]
+    return host_ddf
+
+
+def compute_case_diagnostics(cfg, host_ddf):
+    nq = cfg["nq"]
+    nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
+
+    if nq == 9:
+        f = host_ddf.reshape(nq, ny, nx)
+        rho = np.sum(f, axis=0)
+        ux = np.zeros_like(rho)
+        uy = np.zeros_like(rho)
+        cx = [0, 1, -1, 0, 0, 1, -1, 1, -1]
+        cy = [0, 0, 0, 1, -1, 1, -1, -1, 1]
+        for i in range(nq):
+            ux += cx[i] * f[i]
+            uy += cy[i] * f[i]
+        rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
+        ux /= rho_safe
+        uy /= rho_safe
+        vel = np.sqrt(ux * ux + uy * uy)
+
+        fluid = cfg["flag"].reshape(ny, nx) == FLUID_FLAG
+        mass = float(np.nansum(rho[fluid]))
+
+        x0 = 1
+        x1 = min(nx - 1, 33)
+        inlet_window = np.zeros_like(fluid)
+        inlet_window[1:ny - 1, x0:x1] = True
+        win_mask = fluid & inlet_window
+        inlet_var = float(np.nanvar(vel[win_mask])) if np.any(win_mask) else float("nan")
+
+        # Quantify how well the first interior column follows the designed inlet profile.
+        x_probe = 1
+        line_mask = fluid[:, x_probe]
+        line_u = ux[:, x_probe]
+        y = np.arange(ny, dtype=np.float32)
+        if int(cfg.get("inlet_profile", 0)) == 0:
+            target = np.full(ny, float(cfg["u0"]), dtype=np.float32)
+        else:
+            yy = (y - 0.5 * (ny - 1)) / (ny - 2.0)
+            target = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy * yy)
+
+        if np.any(line_mask):
+            diff = line_u[line_mask] - target[line_mask]
+            t_ref = float(np.max(np.abs(target[line_mask])))
+            denom = max(t_ref, 1.0e-8)
+            inlet_rel_l2 = float(np.sqrt(np.mean(diff * diff)) / denom)
+            inlet_rel_linf = float(np.max(np.abs(diff)) / denom)
+        else:
+            inlet_rel_l2 = float("nan")
+            inlet_rel_linf = float("nan")
+
+        # Inlet-plane-wave indicator: streamwise oscillation of column-averaged
+        # macros in the pre-obstacle region.
+        cx = float(cfg.get("cx", 0.25 * nx))
+        radius = float(cfg.get("radius", ny / 12.0))
+        x_pre0 = 1
+        x_pre1 = min(nx - 2, max(x_pre0 + 4, int(cx - 1.5 * radius)))
+        col_u = []
+        col_r = []
+        for xp in range(x_pre0, x_pre1):
+            col_mask = fluid[1:ny - 1, xp]
+            if np.any(col_mask):
+                col_u.append(float(np.mean(ux[1:ny - 1, xp][col_mask])))
+                col_r.append(float(np.mean(rho[1:ny - 1, xp][col_mask])))
+
+        if int(cfg.get("inlet_profile", 0)) == 0:
+            u_target_mean = float(cfg["u0"])
+        else:
+            y_int = np.arange(1, ny - 1, dtype=np.float32)
+            yy_int = (y_int - 0.5 * (ny - 1)) / (ny - 2.0)
+            target_int = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy_int * yy_int)
+            u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])
+
+        if len(col_u) >= 4:
+            col_u_arr = np.array(col_u, dtype=np.float64)
+            col_r_arr = np.array(col_r, dtype=np.float64)
+            inlet_wave_ux_rel = float(np.std(col_u_arr - u_target_mean) / max(abs(u_target_mean), 1.0e-8))
+            rho_ref = max(abs(float(np.mean(col_r_arr))), 1.0e-8)
+            inlet_wave_rho_rel = float(np.std(col_r_arr) / rho_ref)
+        else:
+            inlet_wave_ux_rel = float("nan")
+            inlet_wave_rho_rel = float("nan")
+
+        # TRT checker/grid-noise indicator: odd-even imbalance in wake ux field.
+        xw0 = min(nx - 3, max(2, int(cx + 2.0 * radius)))
+        xw1 = min(nx - 2, max(xw0 + 4, int(cx + 12.0 * radius)))
+        yw0, yw1 = 2, ny - 2
+        if xw1 > xw0 + 2 and yw1 > yw0 + 2:
+            reg = ux[yw0:yw1, xw0:xw1].astype(np.float64)
+            reg_mask = fluid[yw0:yw1, xw0:xw1]
+            if np.any(reg_mask):
+                valid = reg[reg_mask]
+                m = float(np.mean(valid))
+                centered = np.where(reg_mask, reg - m, np.nan)
+                rms = float(np.sqrt(np.mean((valid - m) * (valid - m))))
+                yy_i, xx_i = np.indices(centered.shape)
+                even_vals = centered[((xx_i + yy_i) & 1) == 0]
+                odd_vals = centered[((xx_i + yy_i) & 1) == 1]
+                even_vals = even_vals[np.isfinite(even_vals)]
+                odd_vals = odd_vals[np.isfinite(odd_vals)]
+
+                if rms > 1.0e-12 and even_vals.size > 8 and odd_vals.size > 8:
+                    wake_checker_rel = float(abs(np.mean(even_vals) - np.mean(odd_vals)) / rms)
+                else:
+                    wake_checker_rel = float("nan")
+
+                pair_mask = reg_mask[:, :-1] & reg_mask[:, 1:]
+                a = centered[:, :-1][pair_mask]
+                b = centered[:, 1:][pair_mask]
+                if a.size > 16 and np.std(a) > 1.0e-12 and np.std(b) > 1.0e-12:
+                    corr = float(np.corrcoef(a, b)[0, 1])
+                    wake_checker_anti_corr_x = float(max(0.0, -corr))
+                else:
+                    wake_checker_anti_corr_x = float("nan")
+            else:
+                wake_checker_rel = float("nan")
+                wake_checker_anti_corr_x = float("nan")
+        else:
+            wake_checker_rel = float("nan")
+            wake_checker_anti_corr_x = float("nan")
+
+        return {
+            "mass": mass,
+            "inlet_var": inlet_var,
+            "inlet_line_rel_l2": inlet_rel_l2,
+            "inlet_line_rel_linf": inlet_rel_linf,
+            "inlet_wave_ux_rel": inlet_wave_ux_rel,
+            "inlet_wave_rho_rel": inlet_wave_rho_rel,
+            "wake_checker_rel": wake_checker_rel,
+            "wake_checker_anti_corr_x": wake_checker_anti_corr_x,
+        }
+
+    f = host_ddf.reshape(nq, nz, ny, nx)
+    rho = np.sum(f, axis=0)
+    ux = np.zeros_like(rho)
+    uy = np.zeros_like(rho)
+    uz = np.zeros_like(rho)
+    cx = np.array([0, 1,-1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0])
+    cy = np.array([0, 0, 0, 1,-1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1])
+    cz = np.array([0, 0, 0, 0, 0, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1])
+    for i in range(nq):
+        ux += cx[i] * f[i]
+        uy += cy[i] * f[i]
+        uz += cz[i] * f[i]
+    rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
+    ux /= rho_safe
+    uy /= rho_safe
+    uz /= rho_safe
+    vel = np.sqrt(ux * ux + uy * uy + uz * uz)
+
+    fluid = cfg["flag"].reshape(nz, ny, nx) == FLUID_FLAG
+    mass = float(np.nansum(rho[fluid]))
+
+    x0 = 1
+    x1 = min(nx - 1, 17)
+    inlet_window = np.zeros_like(fluid)
+    inlet_window[:, 1:ny - 1, x0:x1] = True
+    win_mask = fluid & inlet_window
+    inlet_var = float(np.nanvar(vel[win_mask])) if np.any(win_mask) else float("nan")
+
+    cx = float(cfg.get("cx", 0.25 * nx))
+    radius = float(cfg.get("radius", ny / 12.0))
+    x_pre0 = 1
+    x_pre1 = min(nx - 2, max(x_pre0 + 4, int(cx - 1.5 * radius)))
+    col_u = []
+    col_r = []
+    for xp in range(x_pre0, x_pre1):
+        col_mask = fluid[:, 1:ny - 1, xp]
+        if np.any(col_mask):
+            col_u.append(float(np.mean(ux[:, 1:ny - 1, xp][col_mask])))
+            col_r.append(float(np.mean(rho[:, 1:ny - 1, xp][col_mask])))
+
+    if int(cfg.get("inlet_profile", 0)) == 0:
+        u_target_mean = float(cfg["u0"])
+    else:
+        y_int = np.arange(1, ny - 1, dtype=np.float32)
+        yy_int = (y_int - 0.5 * (ny - 1)) / (ny - 2.0)
+        target_int = float(cfg["u0"]) * 1.5 * (1.0 - 4.0 * yy_int * yy_int)
+        u_target_mean = float(np.mean(target_int)) if target_int.size > 0 else float(cfg["u0"])
+
+    if len(col_u) >= 4:
+        col_u_arr = np.array(col_u, dtype=np.float64)
+        col_r_arr = np.array(col_r, dtype=np.float64)
+        inlet_wave_ux_rel = float(np.std(col_u_arr - u_target_mean) / max(abs(u_target_mean), 1.0e-8))
+        rho_ref = max(abs(float(np.mean(col_r_arr))), 1.0e-8)
+        inlet_wave_rho_rel = float(np.std(col_r_arr) / rho_ref)
+    else:
+        inlet_wave_ux_rel = float("nan")
+        inlet_wave_rho_rel = float("nan")
+
+    return {
+        "mass": mass,
+        "inlet_var": inlet_var,
+        "inlet_wave_ux_rel": inlet_wave_ux_rel,
+        "inlet_wave_rho_rel": inlet_wave_rho_rel,
+        "wake_checker_rel": float("nan"),
+        "wake_checker_anti_corr_x": float("nan"),
+        "inlet_line_rel_l2": float("nan"),
+        "inlet_line_rel_linf": float("nan"),
+    }
+
+
+def compute_trt_les_fields_2d(cfg, host_ddf):
+    if cfg["nq"] != 9:
+        return None
+
+    nx, ny = cfg["nx"], cfg["ny"]
+    f = host_ddf.reshape(9, ny, nx).astype(np.float64)
+
+    cx = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1], dtype=np.float64).reshape(9, 1, 1)
+    cy = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1], dtype=np.float64).reshape(9, 1, 1)
+    w = lattice_weights(9).astype(np.float64).reshape(9, 1, 1)
+
+    rho = np.sum(f, axis=0)
+    rho_safe = np.where(np.abs(rho) > 1.0e-12, rho, 1.0)
+    ux = np.sum(f * cx, axis=0) / rho_safe
+    uy = np.sum(f * cy, axis=0) / rho_safe
+
+    u2 = ux * ux + uy * uy
+    cu = 3.0 * (cx * ux[None, :, :] + cy * uy[None, :, :])
+    feq = w * rho[None, :, :] * (1.0 + cu + 0.5 * cu * cu - 1.5 * u2[None, :, :])
+    fneq = f - feq
+
+    pixx = np.sum(fneq * cx * cx, axis=0)
+    piyy = np.sum(fneq * cy * cy, axis=0)
+    pixy = np.sum(fneq * cx * cy, axis=0)
+
+    omega0 = float(cfg["omega"])
+    omega_min = float(cfg.get("omega_collision_min", OMEGA_COLLISION_MIN_DEFAULT))
+    omega_max = float(cfg["omega_collision_max"])
+    tau0 = max(1.0 / max(omega0, 1.0e-6), 0.500001)
+    nu0 = (tau0 - 0.5) * (1.0 / 3.0)
+    tau_max = 1.0 / max(omega_min, 1.0e-6)
+
+    tau_eff = np.full((ny, nx), tau0, dtype=np.float64)
+    nut = np.zeros((ny, nx), dtype=np.float64)
+    rho_ref = np.maximum(rho_safe, 1.0e-12)
+    cs2 = 1.0 / 3.0
+    nut_cap = max(0.0, LES_POST_NUT_MAX_RATIO * nu0)
+    cs = float(cfg.get("les_cs", 0.16))
+
+    for _ in range(LES_POST_FP_ITERS):
+        denom = 2.0 * rho_ref * cs2 * np.maximum(tau_eff, 0.500001)
+        sxx = -pixx / denom
+        syy = -piyy / denom
+        sxy = -pixy / denom
+
+        tr = 0.5 * (sxx + syy)
+        sxx_dev = sxx - tr
+        syy_dev = syy - tr
+        sxy_dev = sxy
+        s_mag = np.sqrt(np.maximum(0.0, 2.0 * (sxx_dev * sxx_dev + syy_dev * syy_dev + 2.0 * sxy_dev * sxy_dev)))
+
+        nut = np.clip((cs * cs) * s_mag, 0.0, nut_cap)
+        tau_new = np.clip(0.5 + 3.0 * (nu0 + nut), 0.500001, tau_max)
+        tau_eff = LES_POST_FP_RELAX * tau_new + (1.0 - LES_POST_FP_RELAX) * tau_eff
+
+    omega_plus = np.clip(1.0 / tau_eff, omega_min, omega_max)
+    denom_odd = np.maximum(1.0 / omega_plus - 0.5, 1.0e-9)
+    lam = float(cfg.get("trt_magic_param", 0.1875))
+    omega_minus = 1.0 / (lam / denom_odd + 0.5)
+
+    fluid = cfg["flag"].reshape(ny, nx) == FLUID_FLAG
+    rho_mean = float(np.mean(rho[fluid])) if np.any(fluid) else float(np.mean(rho))
+    rho_prime = rho - rho_mean
+
+    return {
+        "fluid": fluid,
+        "omega_plus": omega_plus,
+        "omega_minus": omega_minus,
+        "nut": nut,
+        "rho_prime": rho_prime,
+    }
+
+
+def plot_trt_les_maps(cfg, host_ddf, out_dir):
+    if cfg["nq"] != 9 or cfg["collision_model"] != 1 or not cfg["use_les"]:
+        return None
+
+    fields = compute_trt_les_fields_2d(cfg, host_ddf)
+    if fields is None:
+        return None
+
+    fluid = fields["fluid"]
+    tag = make_case_tag(cfg)
+    out_path = os.path.join(out_dir, f"{tag}_trt_les_fields.png")
+
+    fig, axes = plt.subplots(2, 2, figsize=(12, 9))
+    panels = [
+        ("omega_plus", r"$\omega_+$", "viridis", False),
+        ("omega_minus", r"$\omega_-$", "viridis", False),
+        ("nut", r"$\nu_t$", "magma", False),
+        ("rho_prime", r"$\rho-\overline{\rho}$", "RdBu_r", True),
+    ]
+
+    for ax, (key, title, cmap, signed) in zip(axes.ravel(), panels):
+        arr = fields[key]
+        arr_m = np.ma.array(arr, mask=~fluid)
+        finite_vals = arr[fluid]
+        if finite_vals.size == 0:
+            vmin, vmax = 0.0, 1.0
+        elif signed:
+            span = np.percentile(np.abs(finite_vals), 99)
+            span = max(float(span), 1.0e-12)
+            vmin, vmax = -span, span
+        else:
+            q1 = float(np.percentile(finite_vals, 1))
+            q99 = float(np.percentile(finite_vals, 99))
+            if not np.isfinite(q1) or not np.isfinite(q99) or q99 <= q1:
+                q1 = float(np.min(finite_vals))
+                q99 = float(np.max(finite_vals))
+                if q99 <= q1:
+                    q99 = q1 + 1.0e-12
+            vmin, vmax = q1, q99
+
+        im = ax.imshow(arr_m, origin="lower", aspect="auto", cmap=cmap, vmin=vmin, vmax=vmax)
+        plt.colorbar(im, ax=ax, fraction=0.046, pad=0.04)
+        ax.set_title(title)
+        ax.set_xlabel("x")
+        ax.set_ylabel("y")
+
+    fig.suptitle(f"{tag} TRT-LES diagnostics")
+    fig.tight_layout()
+    fig.savefig(out_path, dpi=160)
+    plt.close(fig)
+    return out_path
+
+
 def plot_case(cfg, host_ddf, out_dir):
     nq = cfg["nq"]
     nx, ny, nz = cfg["nx"], cfg["ny"], cfg["nz"]
@@ -176,7 +544,7 @@ def compute_vis_omega(reynolds, diameter, u0):
 
 def set_macros(nx, ny, nz, dim, nq, vis, u0, collision_model, use_les, les_cs,
                outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
-               omega_collision_max):
+               omega_collision_max, inlet_profile, trt_magic_param):
     lines = compiler.read_lines(compiler.kernel_path("macros.h"))
     defs = {
         "MULT_GPU": "False",
@@ -203,11 +571,12 @@ def set_macros(nx, ny, nz, dim, nq, vis, u0, collision_model, use_les, les_cs,
         "USE_DDF_SHIFTING": 0,
         "USE_LES": int(use_les),
         "LES_CS": f"{les_cs:.6f}f",
-        "INLET_PROFILE": 0,
+        "INLET_PROFILE": int(inlet_profile),
         "OUTLET_MODE": int(outlet_mode),
         "OUTLET_BACKFLOW_CLAMP": int(outlet_backflow_clamp),
         "OUTLET_BLEND_ALPHA": f"{float(outlet_blend_alpha):.3f}f",
         "OMEGA_COLLISION_MAX": f"{float(omega_collision_max):.3f}f",
+        "TRT_MAGIC_PARAM": f"{float(trt_magic_param):.6f}f",
     }
     for name, value in defs.items():
         lines = compiler.modify_macro(lines, name, value)
@@ -261,6 +630,8 @@ def run_case(device_id, cfg):
         outlet_backflow_clamp=cfg["outlet_backflow_clamp"],
         outlet_blend_alpha=cfg["outlet_blend_alpha"],
         omega_collision_max=cfg["omega_collision_max"],
+        inlet_profile=cfg["inlet_profile"],
+        trt_magic_param=cfg["trt_magic_param"],
     )
     compiler.compile_kernel_v2()
 
@@ -310,10 +681,18 @@ def run_case(device_id, cfg):
         grid = ((nx + 127) // 128, ny, nz)
 
         init_fn(d_flag, d_fi, block=block, grid=grid)
-        cuda.memcpy_dtod(d_fi2, d_fi, fsize)
-
         cuda.memcpy_htod(d_flag, cfg["flag"])
 
+        host0 = np.empty(n * nq, dtype=np.float32)
+        cuda.memcpy_dtoh(host0, d_fi)
+        host0 = impose_rest_state_on_nonfluid(cfg, host0)
+        cuda.memcpy_htod(d_fi, host0)
+        cuda.memcpy_htod(d_fi2, host0)
+
+        diag0 = compute_case_diagnostics(cfg, host0)
+        mass0 = diag0["mass"]
+        inlet_var0 = diag0["inlet_var"]
+
         t0 = time.time()
         for step in range(cfg["steps"]):
             step_fn(d_flag, d_fi, d_fi2, d_indx, d_delta, d_action, d_obs, block=block, grid=grid)
@@ -347,10 +726,26 @@ def run_case(device_id, cfg):
             center = float(rho[nz // 2, ny // 2, nx // 2])
 
         ok, reason = validate_case(rho)
+        diag_end = compute_case_diagnostics(cfg, host)
+
+        if np.isfinite(mass0) and mass0 != 0.0 and np.isfinite(diag_end["mass"]):
+            mass_drift = abs(diag_end["mass"] - mass0) / abs(mass0)
+        else:
+            mass_drift = float("nan")
+
+        if np.isfinite(inlet_var0) and inlet_var0 > 0.0 and np.isfinite(diag_end["inlet_var"]):
+            inlet_var_ratio_to_init = diag_end["inlet_var"] / inlet_var0
+        else:
+            inlet_var_ratio_to_init = float("nan")
+
         plot_path = None
         if cfg.get("save_plot", True):
             plot_path = plot_case(cfg, host, cfg["out_dir"])
 
+        trt_les_map_path = None
+        if cfg.get("save_trt_les_maps", True):
+            trt_les_map_path = plot_trt_les_maps(cfg, host, cfg["out_dir"])
+
         return {
             "case_tag": make_case_tag(cfg),
             "name": cfg["name"],
@@ -361,6 +756,18 @@ def run_case(device_id, cfg):
             "rho_center": center,
             "rho_min": float(np.nanmin(rho)),
             "rho_max": float(np.nanmax(rho)),
+            "mass0": float(mass0),
+            "mass_end": float(diag_end["mass"]),
+            "mass_drift": float(mass_drift),
+            "inlet_var0": float(inlet_var0),
+            "inlet_var_end": float(diag_end["inlet_var"]),
+            "inlet_var_ratio_to_init": float(inlet_var_ratio_to_init),
+            "inlet_line_rel_l2": float(diag_end.get("inlet_line_rel_l2", float("nan"))),
+            "inlet_line_rel_linf": float(diag_end.get("inlet_line_rel_linf", float("nan"))),
+            "inlet_wave_ux_rel": float(diag_end.get("inlet_wave_ux_rel", float("nan"))),
+            "inlet_wave_rho_rel": float(diag_end.get("inlet_wave_rho_rel", float("nan"))),
+            "wake_checker_rel": float(diag_end.get("wake_checker_rel", float("nan"))),
+            "wake_checker_anti_corr_x": float(diag_end.get("wake_checker_anti_corr_x", float("nan"))),
             "omega": cfg["omega"],
             "vis": cfg["vis"],
             "collision_model": cfg["collision_model"],
@@ -369,10 +776,13 @@ def run_case(device_id, cfg):
             "outlet_mode": int(cfg["outlet_mode"]),
             "outlet_backflow_clamp": int(cfg["outlet_backflow_clamp"]),
             "outlet_blend_alpha": float(cfg["outlet_blend_alpha"]),
+            "inlet_profile": int(cfg["inlet_profile"]),
             "omega_collision_max": float(cfg["omega_collision_max"]),
+            "trt_magic_param": float(cfg["trt_magic_param"]),
             "pass": bool(ok),
             "reason": reason,
             "plot_path": plot_path,
+            "trt_les_map_path": trt_les_map_path,
         }
     finally:
         ctx.pop()
@@ -380,7 +790,7 @@ def run_case(device_id, cfg):
 
 def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
                   outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
-                  omega_collision_max):
+                  omega_collision_max, inlet_profile, trt_magic_param):
     nx, ny, nz = 512, 256, 1
     cx, cy, radius = 128.0, 128.0, 24.0
     u0 = 0.03
@@ -392,6 +802,9 @@ def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
         "nx": nx,
         "ny": ny,
         "nz": nz,
+        "cx": cx,
+        "cy": cy,
+        "radius": radius,
         "flag": build_flags_2d(nx, ny, cx, cy, radius),
         "u0": u0,
         "vis": vis,
@@ -404,15 +817,18 @@ def build_case_2d(re2d, steps2d, collision_model, use_les, les_cs, out_dir,
         "outlet_mode": int(outlet_mode),
         "outlet_backflow_clamp": int(outlet_backflow_clamp),
         "outlet_blend_alpha": float(outlet_blend_alpha),
+        "inlet_profile": int(inlet_profile),
         "omega_collision_max": float(omega_collision_max),
+        "trt_magic_param": float(trt_magic_param),
         "target_re": re2d,
         "save_plot": True,
+        "save_trt_les_maps": True,
         "out_dir": out_dir,
     }
 
 def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
                   outlet_mode, outlet_backflow_clamp, outlet_blend_alpha,
-                  omega_collision_max):
+                  omega_collision_max, inlet_profile, trt_magic_param):
     nx, ny, nz = 256, 128, 32
     cx, cy, radius = 64.0, 64.0, 12.0
     u0 = 0.04
@@ -424,6 +840,9 @@ def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
         "nx": nx,
         "ny": ny,
         "nz": nz,
+        "cx": cx,
+        "cy": cy,
+        "radius": radius,
         "flag": build_flags_3d(nx, ny, nz, cx, cy, radius),
         "u0": u0,
         "vis": vis,
@@ -436,9 +855,12 @@ def build_case_3d(re3d, steps3d, collision_model, use_les, les_cs, out_dir,
         "outlet_mode": int(outlet_mode),
         "outlet_backflow_clamp": int(outlet_backflow_clamp),
         "outlet_blend_alpha": float(outlet_blend_alpha),
+        "inlet_profile": int(inlet_profile),
         "omega_collision_max": float(omega_collision_max),
+        "trt_magic_param": float(trt_magic_param),
         "target_re": re3d,
         "save_plot": True,
+        "save_trt_les_maps": True,
         "out_dir": out_dir,
     }
 
@@ -454,14 +876,18 @@ def build_comprehensive_cases(args, out_dir):
                                        outlet_mode=args.outlet_mode,
                                        outlet_backflow_clamp=1,
                                        outlet_blend_alpha=args.outlet_blend_alpha,
-                                       omega_collision_max=args.omega_collision_max))
+                                       omega_collision_max=args.omega_collision_max,
+                                       inlet_profile=args.inlet_profile,
+                                       trt_magic_param=args.trt_magic_param))
             cases.append(build_case_3d(re3d=200.0, steps3d=args.matrix_steps3d,
                                        collision_model=cm, use_les=les,
                                        les_cs=args.les_cs, out_dir=out_dir,
                                        outlet_mode=args.outlet_mode,
                                        outlet_backflow_clamp=1,
                                        outlet_blend_alpha=args.outlet_blend_alpha,
-                                       omega_collision_max=args.omega_collision_max))
+                                       omega_collision_max=args.omega_collision_max,
+                                       inlet_profile=args.inlet_profile,
+                                       trt_magic_param=args.trt_magic_param))
     return cases
 
 
@@ -480,10 +906,14 @@ def main():
     parser.add_argument("--les-cs", type=float, default=0.16)
     parser.add_argument("--outlet-mode", type=int, default=0, choices=[0, 1, 2],
                         help="0=non-equilibrium extrapolation, 1=zero-gradient copy, 2=damped blend")
+    parser.add_argument("--inlet-profile", type=int, default=1, choices=[0, 1],
+                        help="0=uniform inlet, 1=parabolic inlet")
     parser.add_argument("--outlet-blend-alpha", type=float, default=0.70,
                         help="Blend alpha for outlet-mode 2")
     parser.add_argument("--omega-collision-max", type=float, default=1.999,
                         help="Upper clamp for collision omega")
+    parser.add_argument("--trt-magic-param", type=float, default=0.002,
+                        help="TRT magic parameter Lambda used in omega- mapping")
     parser.add_argument("--only", choices=["2d", "3d", "both"], default="both")
     parser.add_argument("--comprehensive", action="store_true",
                         help="Run coverage matrix: SRT/TRT/MRT x LES on/off for 2D and 3D")
@@ -504,7 +934,8 @@ def main():
         if args.only in ("2d", "both"):
             c2 = build_case_2d(args.re2d, args.steps2d, args.collision, args.use_les,
                                args.les_cs, out_dir, args.outlet_mode, 1,
-                               args.outlet_blend_alpha, args.omega_collision_max)
+                               args.outlet_blend_alpha, args.omega_collision_max,
+                               args.inlet_profile, args.trt_magic_param)
             print("\n=== Running 2D high-Re case ===")
             print(f"  target Re={args.re2d:.1f}, vis={c2['vis']:.6e}, omega={c2['omega']:.6f}")
             results.append(run_case(args.device, c2))
@@ -512,7 +943,8 @@ def main():
         if args.only in ("3d", "both"):
             c3 = build_case_3d(args.re3d, args.steps3d, args.collision, args.use_les,
                                args.les_cs, out_dir, args.outlet_mode, 1,
-                               args.outlet_blend_alpha, args.omega_collision_max)
+                               args.outlet_blend_alpha, args.omega_collision_max,
+                               args.inlet_profile, args.trt_magic_param)
             print("\n=== Running 3D high-Re case ===")
             print(f"  target Re={args.re3d:.1f}, vis={c3['vis']:.6e}, omega={c3['omega']:.6f}")
             results.append(run_case(args.device, c3))
@@ -534,9 +966,15 @@ def main():
                 n_pass += 1
             print(f"{r['name']}: nan={r['nan_count']}, rho_center={r['rho_center']:.6f}, "
                   f"rho[min,max]=[{r['rho_min']:.6f}, {r['rho_max']:.6f}], "
+                  f"mass_drift={r['mass_drift']:.3e}, inlet_var_end={r['inlet_var_end']:.3e}, "
+                                    f"inlet_relL2={r['inlet_line_rel_l2']:.3e}, inlet_relLinf={r['inlet_line_rel_linf']:.3e}, "
+                                    f"waveUx={r['inlet_wave_ux_rel']:.3e}, waveRho={r['inlet_wave_rho_rel']:.3e}, "
+                                    f"chkRel={r['wake_checker_rel']:.3e}, chkAntiX={r['wake_checker_anti_corr_x']:.3e}, "
                   f"MLUPS={r['mlups']:.1f}, pass={r['pass']} ({r['reason']})")
             if r.get("plot_path"):
                 print(f"  plot: {r['plot_path']}")
+            if r.get("trt_les_map_path"):
+                print(f"  trt_les_map: {r['trt_les_map_path']}")
         print(f"Pass rate: {n_pass}/{len(results)}")
         print(f"Saved: {out_json}")
     finally: