38 lines
1.9 KiB
Diff
38 lines
1.9 KiB
Diff
--- 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
|