--- a/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
+++ b/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh
@@ -5,7 +5,9 @@
 //   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.
+// |S| = sqrt(2 S_ij S_ij), with Delta = 1 in lattice units.
+// The stress estimate uses self-consistent tau_eff; tau0 is molecular.
+// Closed form uses full stress. Iteration uses deviatoric stress and a cap.
 // ============================================================================
 
 #ifndef CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
@@ -65,12 +67,12 @@
 
 #if LES_CLOSED_FORM
     // Closed-form Smagorinsky (no iteration):
-    //   ω = 2 / (τ₀ + √(τ₀² + 18·Cs²·√Q / ρ))
+    //   ω = 2 / (τ₀ + √(τ₀² + 18·√2·Cs²·√Q / ρ))
     // where Q = Πneq_ij · Πneq_ij (Frobenius norm squared)
     const float Q = pixx*pixx + piyy*piyy + 2.0f*pixy*pixy;
     const float sqrtQ = sqrtf(fmaxf(Q, 0.0f));
     const float smag_const_sq = LES_CS * LES_CS;
-    const float discriminant = tau0*tau0 + 18.0f * smag_const_sq * sqrtQ / rho_safe;
+    const float discriminant = tau0*tau0 + 18.0f * sqrtf(2.0f) * smag_const_sq * sqrtQ / rho_safe;
     const float tau_eff = 0.5f * (tau0 + sqrtf(fmaxf(discriminant, tau0*tau0)));
     return clamp_omega(1.0f / fmaxf(tau_eff, 0.500001f));
 #else
@@ -140,7 +142,7 @@
                   + 2.0f*(pixy*pixy + pixz*pixz + piyz*piyz);
     const float sqrtQ = sqrtf(fmaxf(Q, 0.0f));
     const float smag_const_sq = LES_CS * LES_CS;
-    const float discriminant = tau0*tau0 + 18.0f * smag_const_sq * sqrtQ / rho_safe;
+    const float discriminant = tau0*tau0 + 18.0f * sqrtf(2.0f) * smag_const_sq * sqrtQ / rho_safe;
     const float tau_eff = 0.5f * (tau0 + sqrtf(fmaxf(discriminant, tau0*tau0)));
     return clamp_omega(1.0f / fmaxf(tau_eff, 0.500001f));
 #else