DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
RANS_ke.hpp
Go to the documentation of this file.
1/** @file RANS_ke.hpp
2 * @brief RANS two-equation turbulence model implementations for the DNDSR CFD solver.
3 *
4 * Provides template functions for computing eddy viscosity, viscous fluxes, and
5 * source terms for several widely-used RANS turbulence closures:
6 *
7 * - **Realizable k-epsilon** (Shih et al.) with f_mu damping and realizability constraint.
8 * - **k-omega SST** (Menter) with F1/F2 blending, cross-diffusion, and DES/DDES/IDDES support.
9 * - **k-omega Wilcox** (1988 / 2006 variants) with stress limiter and f_beta correction.
10 * - **Spalart-Allmaras** (SA) one-equation model with optional rotation correction,
11 * negative-nuTilde extension, and DES/DDES/IDDES/WMLES length-scale modification.
12 *
13 * All functions are templated on spatial dimension @c dim and accept Eigen-compatible
14 * expression types for the conservative state vector @c UMeanXy and the gradient tensor
15 * @c DiffUxy (dim × nVars). The convention for variable layout is:
16 * - Indices 0..dim: density and momentum components.
17 * - Index @c I4 = dim+1: total energy.
18 * - Index @c I4+1: first RANS variable (k or nuTilde).
19 * - Index @c I4+2: second RANS variable (epsilon or omega), where applicable.
20 *
21 * Compile-time macros control model variants and limiters; see definitions below.
22 *
23 * @note This file is included by other Euler module headers and should not normally
24 * be included directly by application code.
25 */
26#pragma once
27
28#include "Euler.hpp"
29
30/** @brief When set to 1, cap turbulent viscosity at 1e5 * mu_laminar for k-epsilon. */
31#define KE_LIMIT_MUT 1
32
33/** @brief Select Wilcox k-omega model version: 0 = original 1988, 1 = 2006 with stress limiter, 2 = simplified. */
34#define KW_WILCOX_VER 1
35/** @brief When set to 1, enable production limiters for Wilcox k-omega (CFL3D-style capping). */
36#define KW_WILCOX_PROD_LIMITS 1
37/** @brief When set to 1, cap turbulent viscosity at 1e5 * mu_laminar for Wilcox k-omega. */
38#define KW_WILCOX_LIMIT_MUT 1
39
40/** @brief When set to 1, cap turbulent viscosity at 1e5 * mu_laminar for SST. */
41#define KW_SST_LIMIT_MUT 1
42/** @brief When set to 1, enable production limiters for SST (CFL3D-style capping). */
43#define KW_SST_PROD_LIMITS 1
44/** @brief Select omega production formulation for SST: 0 = classical, 1 = strain-rate based, 2 = vorticity based. */
45#define KW_SST_PROD_OMEGA_VERSION 1
46
47/** @brief When set to 1, enable the ft2 laminar suppression term in Spalart-Allmaras. Setting to 0 disables ft2. */
48#define SA_USE_FT2_TERM 1
49
50/** @brief RANS turbulence model functions (eddy viscosity, viscous flux, source terms). */
51namespace DNDS::Euler::RANS
52{
53 /** @brief Spalart-Allmaras model constants used across multiple files.
54 *
55 * Canonical definitions from Spalart & Allmaras (1994).
56 * These constants are shared between the SA source term computation
57 * and boundary condition routines in other translation units.
58 */
59 namespace SA
60 {
61 static constexpr real cnu1 = 7.1; ///< SA model constant c_{v1} controlling the viscous damping function f_{v1}.
62 static constexpr real cn1 = 16.0; ///< SA model constant c_{n1} for the negative-nuTilde extension function f_n.
63 static constexpr real sigma = 2.0 / 3.0; ///< SA model diffusion coefficient sigma.
64 }
65
66 /** @brief Wall-omega boundary condition coefficient for k-omega family models.
67 *
68 * Applied at wall boundaries as: rhoOmegaWall = mu / d^2 * kWallOmegaCoeff,
69 * where d is the wall-normal distance of the first cell center.
70 */
71 static constexpr real kWallOmegaCoeff = 800.0;
72
73 /**
74 * @brief Compute turbulent viscosity mu_t from the Shih et al. realizable k-epsilon model.
75 *
76 * Evaluates the eddy viscosity using the realizable k-epsilon formulation with:
77 * - f_mu damping function based on the turbulent Reynolds number Re_t = rho * k^2 / (mu * epsilon).
78 * - Realizability constraint: mu_t <= phi * rho * k / S, where S is the strain-rate magnitude.
79 * - Optional upper bound of 1e5 * mu_laminar when @c KE_LIMIT_MUT is enabled.
80 *
81 * @tparam dim Spatial dimension (2 or 3).
82 * @tparam TU Eigen-compatible type for the conservative state vector.
83 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
84 * @param UMeanXy Conservative state vector [rho, rhoU, ..., rhoE, rho*k, rho*epsilon].
85 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
86 * @param muf Laminar (molecular) dynamic viscosity.
87 * @param d Wall distance (unused in this model but kept for interface consistency).
88 * @return Turbulent dynamic viscosity mu_t (>= 0).
89 */
90 template <int dim, class TU, class TDiffU>
91 real GetMut_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d)
92 {
93 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
94 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
95 static const auto I4 = dim + 1;
96 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
97 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
98 real cmu = 0.09;
99 real phi = 2. / 3.;
100 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
101 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
102 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
103 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
104 Eigen::Matrix<real, dim, dim> SS = diffU + diffU.transpose() - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity();
105 real rho = UMeanXy(0);
106 real k = UMeanXy(I4 + 1) / rho + verySmallReal_4;
107 real epsilon = UMeanXy(I4 + 2) / rho + verySmallReal_4;
108 real Ret = rho * sqr(k) / (muf * epsilon) + smallReal;
109 real S = std::sqrt(SS.squaredNorm() / 2) + verySmallReal_4;
110 real fmu = (1 - std::exp(-0.01 * Ret)) / (1 - exp(-std::sqrt(Ret))) * std::max(1., std::sqrt(2 / Ret));
111 real mut = cmu * fmu * rho * sqr(k) / epsilon;
112 mut = std::min(mut, phi * rho * k / S);
113#if KE_LIMIT_MUT == 1
114 mut = std::min(mut, 1e5 * muf); // CFL3D
115#endif
116 if (std::isnan(mut) || !std::isfinite(mut))
117 {
118 std::cerr << k << " " << epsilon << " " << Ret << " " << S << "\n";
119 std::cerr << fmu << "\n";
120 std::cerr << SS << std::endl;
121 DNDS_assert(false);
122 }
123 return mut;
124 }
125
126 /**
127 * @brief Compute the RANS viscous flux for the realizable k-epsilon transport equations.
128 *
129 * Evaluates the diffusion (viscous) flux contributions for the k and epsilon equations
130 * projected onto the face normal direction. The effective diffusivities are:
131 * - k-equation: (mu_laminar + mu_t / sigma_k), sigma_k = 1.0
132 * - epsilon-equation: (mu_laminar + mu_t / sigma_e), sigma_e = 1.3
133 *
134 * @tparam dim Spatial dimension (2 or 3).
135 * @tparam TU Eigen-compatible type for the conservative state vector.
136 * @tparam TN Eigen-compatible type for the face unit normal vector.
137 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, primitive variables).
138 * @tparam TVFlux Eigen-compatible type for the output viscous flux vector.
139 * @param UMeanXy Conservative state vector.
140 * @param DiffUxyPrim Gradient tensor of *primitive* variables, size dim × nVars.
141 * @param uNorm Outward-pointing face unit normal vector (length dim).
142 * @param mut Turbulent dynamic viscosity (precomputed).
143 * @param d Wall distance (unused here, kept for interface consistency).
144 * @param muPhy Laminar (molecular) dynamic viscosity.
145 * @param[out] vFlux Output viscous flux vector; entries at I4+1 and I4+2 are set.
146 */
147 template <int dim, class TU, class TN, class TDiffU, class TVFlux>
148 void GetVisFlux_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxyPrim, TN &&uNorm, real mut, real d, real muPhy, TVFlux &vFlux)
149 {
150 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
151 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
152 static const auto I4 = dim + 1;
153 real sigK = 1.;
154 real sigE = 1.3;
155 vFlux(I4 + 1) = DiffUxyPrim(Seq012, {I4 + 1}).dot(uNorm) * (muPhy + mut / sigK);
156 vFlux(I4 + 2) = DiffUxyPrim(Seq012, {I4 + 2}).dot(uNorm) * (muPhy + mut / sigE);
157 }
158
159 /**
160 * @brief Compute source terms for the realizable k-epsilon transport equations.
161 *
162 * Evaluates production, destruction, and extra source contributions for the
163 * k and epsilon equations following the Shih et al. realizable k-epsilon model.
164 *
165 * Key terms computed:
166 * - Reynolds stress tensor tau_{ij} via Boussinesq hypothesis with realizability clipping.
167 * - Production of k: P_k = -tau_{ij} * dU_i/dx_j, capped by the Shih pphi limiter.
168 * - Turbulent time scale T_t with low-Re correction via zeta = sqrt(Re_t / 2).
169 * - Extra dissipation term E involving the gradient of the turbulent time scale.
170 *
171 * When @p mode == 0, the full source vector is returned:
172 * - source(I4+1) = P_k - rho*epsilon
173 * - source(I4+2) = (c_{e1} * P_k - c_{e2} * rho*epsilon + E) / T_t
174 *
175 * When @p mode == 1, the implicit diagonal contribution (positive) is returned for
176 * use in implicit time stepping:
177 * - source(I4+1) = 0
178 * - source(I4+2) = c_{e2} / T_t
179 *
180 * @tparam dim Spatial dimension (2 or 3).
181 * @tparam TU Eigen-compatible type for the conservative state vector.
182 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
183 * @tparam TSource Eigen-compatible type for the output source vector.
184 * @param UMeanXy Conservative state vector.
185 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
186 * @param muf Laminar (molecular) dynamic viscosity.
187 * @param d Wall distance (unused here, kept for interface consistency).
188 * @param[out] source Output source vector; entries at I4+1 and I4+2 are set.
189 * @param mode Source computation mode: 0 = full source, 1 = implicit diagonal (destruction only).
190 */
191 template <int dim, class TU, class TDiffU, class TSource>
192 void GetSource_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d, TSource &source, int mode)
193 {
194 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
195 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
196 static const auto I4 = dim + 1;
197 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
198 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
199
200 real cmu = 0.09;
201 real phi = 2. / 3.;
202 real ce1 = 1.44;
203 real ce2 = 1.92;
204 real AE = 0.3;
205 real pphi = 1.065;
206 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
207 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
208 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
209 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
210 Eigen::Matrix<real, dim, dim> SS = diffU + diffU.transpose() - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity();
211 real rho = UMeanXy(0);
212 real k = std::max(UMeanXy(I4 + 1) / rho, verySmallReal_4);
213 real epsilon = std::max(UMeanXy(I4 + 1) / rho, verySmallReal_4);
214 real Ret = rho * sqr(k) / (muf * epsilon) + verySmallReal_4;
215 real S = std::sqrt(SS.squaredNorm() / 2) + verySmallReal_4;
216 real fmu = (1 - std::exp(-0.01 * Ret)) / std::max(1 - exp(-std::sqrt(Ret)), verySmallReal_4) * std::max(1., std::sqrt(2 / Ret));
217 real mut = cmu * fmu * rho * sqr(k) / epsilon;
218 mut = std::min(mut, phi * rho * k / S);
219
220 Eigen::Matrix<real, dim, dim> rhoMuiuj = Eigen::Matrix<real, dim, dim>::Identity() * UMeanXy(I4 + 1) * (2. / 3.) - mut * SS;
221 real Pk = -(rhoMuiuj.array() * diffU.array()).sum();
222 Pk = std::min(Pk, pphi * cmu * sqr(UMeanXy(I4 + 1)) / mut);
223 real zeta = std::sqrt(Ret / 2);
224 real Tt = k / epsilon * std::max(1., 1. / zeta) + smallReal;
225
226 Eigen::Matrix<real, dim, 2> diffRhoKe = DiffUxy(Seq012, {I4 + 1, I4 + 2});
227 Eigen::Matrix<real, dim, 2> diffKe = (diffRhoKe - 1. / UMeanXy(0) * diffRho * UMeanXy({I4 + 1, I4 + 2}).transpose()) / UMeanXy(0);
228 Eigen::Matrix<real, dim, 1> diffTau = diffKe(Seq012, 0) / epsilon - k / sqr(epsilon) * diffKe(Seq012, 1);
229 real Psi = std::max(0., diffKe(Seq012, 0).dot(diffTau));
230 real E = AE * rho * std::sqrt(epsilon * Tt) * Psi * std::max(std::sqrt(k), std::pow(muf * epsilon / rho, 0.25));
231
232 if (mode == 0)
233 {
234 source(I4 + 1) = Pk - UMeanXy(I4 + 2);
235 source(I4 + 2) = (ce1 * Pk - ce2 * UMeanXy(I4 + 2) + E) / Tt;
236 }
237 else
238 {
239 source(I4 + 1) = 0;
240 source(I4 + 2) = (ce2) / Tt;
241 }
242 if (!source.allFinite() || source.hasNaN())
243 {
244 std::cerr << source.transpose() << "\n";
245 std::cerr << UMeanXy.transpose() << "\n";
246 std::cerr << DiffUxy << "\n";
247 std::cerr << S << "\n";
248 std::cerr << mut << "\n";
249
250 std::cout << std::endl;
251
252 DNDS_assert(false);
253 }
254 }
255
256 /**
257 * @brief Attempt to find equilibrium epsilon for zero-gradient (freestream) conditions via Newton iteration.
258 *
259 * Iteratively solves for the epsilon value that zeroes the k-equation source term
260 * when all spatial gradients are zero. Uses a simple Newton-Raphson loop with
261 * finite-difference Jacobian evaluation.
262 *
263 * @warning This function is noted as **not applicable** in practice — the zero-gradient
264 * equilibrium assumption does not hold for the realizable k-epsilon model.
265 * Retained for reference and experimentation only.
266 *
267 * @tparam dim Spatial dimension (2 or 3).
268 * @tparam TU Eigen-compatible type for the conservative state vector.
269 * @param[in,out] u Conservative state vector; u(I4+2) (rho*epsilon) is updated in place.
270 * @param muPhy Laminar (molecular) dynamic viscosity.
271 * @return A tuple (rhs_initial, rhs_final) giving the k-equation residual before and
272 * after the Newton iteration, for convergence diagnostics.
273 */
274 template <int dim, class TU>
275 std::tuple<real, real> SolveZeroGradEquilibrium(TU &u, real muPhy) // this is tested to be not applicable!
276 {
277 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
278 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
279 static const auto I4 = dim + 1;
280 Eigen::Matrix<real, dim, 7> diffU;
281 diffU.setZero();
282
283 Eigen::Vector<real, Eigen::Dynamic> src = u;
284 Eigen::Vector<real, Eigen::Dynamic> uc = u;
285 src.setZero();
286
287 real distEps = std::sqrt(smallReal);
288 real distV = 1 + distEps;
289
290 auto getDE = [&](real re) -> real
291 {
292 uc(I4 + 2) = re;
293 GetSource_RealizableKe<dim>(uc, diffU, muPhy, src);
294 return src(I4 + 1);
295 };
296 std::cout << "Mu " << muPhy << std::endl;
297
298 real re = u(I4 + 2);
299 real rhs0{0}, rhs{0};
300 for (int i = 0; i < 1000; i++)
301 {
302 rhs = getDE(re);
303 real dRhsDRe = (getDE(re * distV) - rhs) / (re * distEps);
304 real ren = re - rhs / (dRhsDRe + verySmallReal);
305 re = std::max(ren, verySmallReal);
306 std::cout << rhs << std::endl;
307 if (i == 0)
308 rhs0 = rhs;
309 }
310 u(I4 + 2) = re;
311 return std::make_tuple(rhs0, rhs);
312 }
313
314 /**
315 * @brief Compute the analytical diagonal of the source Jacobian for the realizable k-epsilon model.
316 *
317 * Evaluates -dS/dU diagonal entries for implicit time stepping of the k and epsilon
318 * transport equations. The diagonal entries approximate the linearised destruction terms:
319 * - source(I4+1) = -P_k / (rho*k) (destruction rate of k)
320 * - source(I4+2) = c_{e2} / T_t (destruction rate of epsilon)
321 *
322 * These positive values are used to augment the implicit operator diagonal,
323 * improving stability of the coupled RANS system.
324 *
325 * @tparam dim Spatial dimension (2 or 3).
326 * @tparam TU Eigen-compatible type for the conservative state vector.
327 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
328 * @tparam TSource Eigen-compatible type for the output Jacobian diagonal vector.
329 * @param UMeanXy Conservative state vector.
330 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
331 * @param muf Laminar (molecular) dynamic viscosity.
332 * @param[out] source Output Jacobian diagonal vector; entries at I4+1 and I4+2 are set.
333 */
334 template <int dim, class TU, class TDiffU, class TSource>
335 void GetSourceJacobianDiag_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, TSource &source)
336 {
337 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
338 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
339 static const auto I4 = dim + 1;
340
341 real cmu = 0.09;
342 real phi = 2. / 3.;
343 real ce1 = 1.44;
344 real ce2 = 1.92;
345 real AE = 0.3;
346 real pphi = 1.065;
347 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
348 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
349 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
350 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
351 Eigen::Matrix<real, dim, dim> SS = diffU + diffU.transpose() - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity();
352 real rho = UMeanXy(0);
353 real k = UMeanXy(I4 + 1) / rho + verySmallReal;
354 real epsilon = UMeanXy(I4 + 2) / rho + verySmallReal;
355 real Ret = rho * sqr(k) / (muf * epsilon) + verySmallReal;
356 real S = std::sqrt(SS.squaredNorm() / 2) + verySmallReal;
357 real fmu = (1 - std::exp(-0.01 * Ret)) / (1 - exp(-std::sqrt(Ret))) * std::max(1., std::sqrt(2 / Ret));
358 real mut = cmu * fmu * rho * sqr(k) / epsilon;
359 mut = std::min(mut, phi * rho * k / S);
360
361 Eigen::Matrix<real, dim, dim> rhoMuiuj = Eigen::Matrix<real, dim, dim>::Identity() * UMeanXy(I4 + 1) * (2. / 3.) - mut * SS;
362 real Pk = -(rhoMuiuj.array() * diffU.array()).sum();
363 Pk = std::min(Pk, pphi * cmu * sqr(UMeanXy(I4 + 1)) / mut);
364 real zeta = std::sqrt(Ret / 2);
365 real Tt = k / epsilon * std::max(1., 1. / zeta) + smallReal;
366
367 Eigen::Matrix<real, dim, 2> diffRhoKe = DiffUxy(Seq012, {I4 + 1, I4 + 2});
368 Eigen::Matrix<real, dim, 2> diffKe = (diffRhoKe - 1. / UMeanXy(0) * diffRho * UMeanXy({I4 + 1, I4 + 2}).transpose()) / UMeanXy(0);
369 Eigen::Matrix<real, dim, 1> diffTau = diffKe(Seq012, 0) / epsilon - k / sqr(epsilon) * diffKe(Seq012, 1);
370 real Psi = std::max(0., diffKe(Seq012, 0).dot(diffTau));
371 real E = AE * rho * std::sqrt(epsilon * Tt) * Psi * std::max(std::sqrt(k), std::pow(muf * epsilon / rho, 0.25));
372
373 real dPk = -((Eigen::Matrix<real, dim, dim>::Identity() * (2. / 3.)).array() * diffU.array()).sum();
374
375 source(I4 + 1) = -Pk / (UMeanXy(I4 + 1) + verySmallReal);
376 source(I4 + 2) = (ce2) / Tt;
377 }
378
379 /**
380 * @brief Compute the numerical (finite-difference) diagonal of the source Jacobian for the realizable k-epsilon model.
381 *
382 * Evaluates -dS_i/dU_i by forward finite-difference perturbation of rho*k and rho*epsilon
383 * independently. The perturbation step is proportional to sqrt(smallReal) times the variable
384 * magnitude. Results are clamped to non-negative values and amplified by a factor of 10
385 * for enhanced implicit stability.
386 *
387 * This is a fallback when the analytical Jacobian diagonal
388 * (GetSourceJacobianDiag_RealizableKe) is suspected of inaccuracy.
389 *
390 * @tparam dim Spatial dimension (2 or 3).
391 * @tparam TU Eigen-compatible type for the conservative state vector.
392 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
393 * @tparam TSource Eigen-compatible type for the output Jacobian diagonal vector.
394 * @param UMeanXy Conservative state vector.
395 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
396 * @param muf Laminar (molecular) dynamic viscosity.
397 * @param[out] source Output numerical Jacobian diagonal vector; entries at I4+1 and I4+2 are set.
398 */
399 template <int dim, class TU, class TDiffU, class TSource>
400 void GetSourceJacobianDiag_RealizableKe_ND(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, TSource &source)
401 {
402 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
403 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
404 static const auto I4 = dim + 1;
405 Eigen::VectorXd u0 = UMeanXy;
406 Eigen::VectorXd u1 = UMeanXy;
407 Eigen::VectorXd sb = source;
408 Eigen::VectorXd s0 = source;
409 Eigen::VectorXd s1 = source;
410 real epsRk = (u0(I4 + 1) + smallReal) * std::sqrt(smallReal);
411 real epsRe = (u0(I4 + 2) + smallReal) * std::sqrt(smallReal);
412 u0(I4 + 1) += epsRk;
413 u1(I4 + 2) += epsRe;
414 GetSource_RealizableKe<dim>(UMeanXy, DiffUxy, muf, sb);
415 GetSource_RealizableKe<dim>(u0, DiffUxy, muf, s0);
416 GetSource_RealizableKe<dim>(u1, DiffUxy, muf, s1);
417 source(I4 + 1) = -(s0(I4 + 1) - sb(I4 + 1)) / epsRk;
418 source(I4 + 2) = -(s1(I4 + 2) - sb(I4 + 2)) / epsRe;
419 source(I4 + 1) = std::max(source(I4 + 1), 0.) * 10;
420 source(I4 + 2) = std::max(source(I4 + 2), 0.) * 10;
421 }
422
423 /**
424 * @brief Compute turbulent viscosity mu_t from Menter's k-omega SST model.
425 *
426 * Evaluates the SST eddy viscosity with:
427 * - Vorticity magnitude Omega = ||(grad U)^T - grad U|| / sqrt(2).
428 * - Blending function F2 based on wall distance and turbulent Reynolds number.
429 * - Bradshaw's structural constraint: mu_t = a1 * k / max(Omega * F2, a1 * omega),
430 * where a1 = 0.31.
431 * - Optional upper bound of 1e5 * mu_laminar when @c KW_SST_LIMIT_MUT is enabled.
432 *
433 * @tparam dim Spatial dimension (2 or 3).
434 * @tparam TU Eigen-compatible type for the conservative state vector.
435 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
436 * @param UMeanXy Conservative state vector [rho, rhoU, ..., rhoE, rho*k, rho*omega].
437 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
438 * @param muf Laminar (molecular) dynamic viscosity.
439 * @param d Wall distance.
440 * @return Turbulent dynamic viscosity mu_t (>= 0).
441 */
442 template <int dim, class TU, class TDiffU>
443 real GetMut_SST(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d)
444 {
445 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
446 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
447 static const auto I4 = dim + 1;
448 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
449 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
450
451 real a1 = 0.31;
452 real betaStar = 0.09;
453 real sigOmega2 = 0.856;
454 real cmu = 0.09;
455 real phi = 2. / 3.;
456 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
457 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
458 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
459 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
460 Eigen::Matrix<real, dim, dim> SR2 = diffU + diffU.transpose(); // 2 times strain rate
461 Eigen::Matrix<real, dim, dim> SS = SR2 - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity(); // 2 times shear strain rate
462 Eigen::Matrix<real, dim, dim> OmegaM = (diffU.transpose() - diffU) * 0.5;
463 real OmegaMag = OmegaM.norm() * std::sqrt(2);
464 real rho = UMeanXy(0);
465 real k = std::max(UMeanXy(I4 + 1) / rho, verySmallReal_4);
466 real omegaaa = std::max(UMeanXy(I4 + 2) / rho, verySmallReal_4);
467 real S = std::sqrt(SS.squaredNorm() / 2) + verySmallReal_4;
468 real nuPhy = muf / rho;
469 real F2 = std::tanh(sqr(std::max(2 * std::sqrt(k) / (betaStar * omegaaa * d), 500 * nuPhy / (sqr(d) * omegaaa))));
470 // F2 = 0;
471 real mut = a1 * k / std::max(OmegaMag * F2, a1 * omegaaa) * rho;
472#if KW_SST_LIMIT_MUT == 1
473 mut = std::min(mut, 1e5 * muf); // CFL3D
474#endif
475
476 if (std::isnan(mut) || !std::isfinite(mut))
477 {
478 std::cerr << k << " " << omegaaa << " " << mut << " " << S << "\n";
479 std::cerr << SS << std::endl;
480 DNDS_assert(false);
481 }
482 return mut;
483 }
484
485 /**
486 * @brief Compute the RANS viscous flux for the k-omega SST transport equations.
487 *
488 * Evaluates the diffusion (viscous) flux contributions for the k and omega equations
489 * projected onto the face normal direction. The effective diffusion coefficients are
490 * blended between the inner-layer (set 1) and outer-layer (set 2) values using the
491 * SST blending function F1:
492 * - sigma_k = sigma_k1 * F1 + sigma_k2 * (1 - F1)
493 * - sigma_omega = sigma_omega1 * F1 + sigma_omega2 * (1 - F1)
494 *
495 * Effective diffusivity for each equation: mu_laminar + mu_t * sigma_{k,omega}.
496 *
497 * @tparam dim Spatial dimension (2 or 3).
498 * @tparam TU Eigen-compatible type for the conservative state vector.
499 * @tparam TN Eigen-compatible type for the face unit normal vector.
500 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, primitive variables).
501 * @tparam TVFlux Eigen-compatible type for the output viscous flux vector.
502 * @param UMeanXy Conservative state vector.
503 * @param DiffUxyPrim Gradient tensor of *primitive* variables, size dim × nVars.
504 * @param uNorm Outward-pointing face unit normal vector (length dim).
505 * @param mutIn Turbulent dynamic viscosity (precomputed, used for diffusion coefficients).
506 * @param d Wall distance.
507 * @param muf Laminar (molecular) dynamic viscosity.
508 * @param[out] vFlux Output viscous flux vector; entries at I4+1 and I4+2 are set.
509 */
510 template <int dim, class TU, class TN, class TDiffU, class TVFlux>
511 void GetVisFlux_SST(TU &&UMeanXy, TDiffU &&DiffUxyPrim, TN &&uNorm, real mutIn, real d, real muf, TVFlux &vFlux)
512 {
513 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
514 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
515 static const auto I4 = dim + 1;
516 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
517 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
518
519 real a1 = 0.31;
520 real betaStar = 0.09;
521 real sigK1 = 0.85;
522 real sigK2 = 1;
523 real sigO1 = 0.5;
524 real sigO2 = 0.856;
525 real beta1 = 0.075;
526 real beta2 = 0.0828;
527 real kap = 0.41;
528 real gamma1 = beta1 / betaStar - sigO1 * sqr(kap) / std::sqrt(betaStar);
529 real gamma2 = beta2 / betaStar - sigO2 * sqr(kap) / std::sqrt(betaStar);
530 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
531 Eigen::Matrix<real, dim, 1> diffRho = DiffUxyPrim(Seq012, {0});
532 Eigen::Matrix<real, dim, dim> diffU = DiffUxyPrim(Seq012, Seq123);
533 Eigen::Matrix<real, dim, dim> SR2 = diffU + diffU.transpose(); // 2 times strain rate
534 Eigen::Matrix<real, dim, dim> SS = SR2 - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity(); // 2 times shear strain rate
535 Eigen::Matrix<real, dim, 2> diffKO = DiffUxyPrim(Seq012, {I4 + 1, I4 + 2});
536
537 Eigen::Matrix<real, dim, dim> OmegaM = (diffU.transpose() - diffU) * 0.5;
538 real OmegaMag = OmegaM.norm() * std::sqrt(2);
539 real rho = UMeanXy(0);
540 real k = std::max(UMeanXy(I4 + 1) / rho, verySmallReal_4);
541 real omegaaa = std::max(UMeanXy(I4 + 2) / rho, verySmallReal_4);
542 real S = std::sqrt(SS.squaredNorm() / 2) + verySmallReal_4;
543 real nuPhy = muf / rho;
544 real CDKW = std::max(2 * rho * sigO2 / omegaaa * diffKO(EigenAll, 0).dot(diffKO(EigenAll, 1)), 1e-10);
545 real F1 = std::tanh(std::pow(
546 std::min(std::max(std::sqrt(k) / (betaStar * omegaaa * d), 500 * nuPhy / (sqr(d) * omegaaa)),
547 4 * rho * sigO2 * k / (CDKW * sqr(d))),
548 4));
549 real F2 = std::tanh(sqr(std::max(2 * std::sqrt(k) / (betaStar * omegaaa * d), 500 * nuPhy / (sqr(d) * omegaaa))));
550 // F2 = 0;
551 real mut = a1 * k / std::max(OmegaMag * F2, a1 * omegaaa) * rho;
552#if KW_SST_LIMIT_MUT == 1
553 mut = std::min(mut, 1e5 * muf); // CFL3D
554#endif
555
556 real sigK = sigK1 * F1 + sigK2 * (1 - F1);
557 real sigO = sigO1 * F1 + sigO2 * (1 - F1);
558
559 vFlux(I4 + 1) = diffKO(Seq012, 0).dot(uNorm) * (muf + mutIn * sigK);
560 vFlux(I4 + 2) = diffKO(Seq012, 1).dot(uNorm) * (muf + mutIn * sigO);
561
562 if (!vFlux.allFinite() || vFlux.hasNaN())
563 {
564 std::cerr << vFlux << "\n";
565 std::cerr << sigK << " " << sigO << "\n";
566 std::cerr << F1 << " " << F2 << "\n";
567 std::cerr << CDKW << "\n";
568 std::cerr << k << " " << omegaaa << "\n";
569 std::cerr << muf << " " << mut << "\n";
570 std::cerr << std::endl;
571 }
572 }
573
574 /**
575 * @brief Compute source terms for the k-omega SST transport equations.
576 *
577 * Evaluates production, destruction, and cross-diffusion source contributions for
578 * the k and omega equations following Menter's SST model. Supports DES/DDES
579 * capability via the RANS-to-LES length scale ratio.
580 *
581 * Key terms:
582 * - Production P_k from the Boussinesq Reynolds stress, optionally capped at
583 * 20 * beta* * rho * k * omega (when @c KW_SST_PROD_LIMITS is enabled).
584 * - Omega production P_omega with version-dependent formulation controlled by
585 * @c KW_SST_PROD_OMEGA_VERSION (0 = classical, 1 = strain-rate, 2 = vorticity).
586 * - Destruction terms: beta* * rho * k * omega for k, beta * rho * omega^2 for omega.
587 * - Cross-diffusion term: 2 * (1 - F1) * rho * sigma_omega2 / omega * grad(k) . grad(omega).
588 * - DES shielding: F_DES = max(1, l_RANS / l_LES * (1 - F2)), reducing the destruction
589 * of k outside the RANS region.
590 *
591 * When @p mode == 0, the full source vector is returned.
592 * When @p mode == 1, the implicit diagonal contribution (positive) is returned:
593 * - source(I4+1) = beta* * omega * F_DES
594 * - source(I4+2) = 2 * beta * omega
595 *
596 * @tparam dim Spatial dimension (2 or 3).
597 * @tparam TU Eigen-compatible type for the conservative state vector.
598 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
599 * @tparam TSource Eigen-compatible type for the output source vector.
600 * @param UMeanXy Conservative state vector.
601 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
602 * @param muf Laminar (molecular) dynamic viscosity.
603 * @param d Wall distance.
604 * @param lLES LES length scale for DES/DDES shielding (set to a large value to disable DES).
605 * @param[out] source Output source vector; entries at I4+1 and I4+2 are set.
606 * @param mode Source computation mode: 0 = full source, 1 = implicit diagonal (destruction only).
607 */
608 template <int dim, class TU, class TDiffU, class TSource>
609 void GetSource_SST(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d, real lLES, TSource &source, int mode)
610 {
611 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
612 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
613 static const auto I4 = dim + 1;
614 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
615 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
616
617 real a1 = 0.31;
618 real betaStar = 0.09;
619 real sigK1 = 0.85;
620 real sigK2 = 1;
621 real sigO1 = 0.5;
622 real sigO2 = 0.856;
623 real beta1 = 0.075;
624 real beta2 = 0.0828;
625 real kap = 0.41;
626 real gamma1 = beta1 / betaStar - sigO1 * sqr(kap) / std::sqrt(betaStar);
627 real gamma2 = beta2 / betaStar - sigO2 * sqr(kap) / std::sqrt(betaStar);
628 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
629 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
630 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
631 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
632 Eigen::Matrix<real, dim, dim> SR2 = diffU + diffU.transpose(); // 2 times strain rate
633 Eigen::Matrix<real, dim, dim> SS = SR2 - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity(); // 2 times shear strain rate
634 Eigen::Matrix<real, dim, 2> diffRhoKO = DiffUxy(Seq012, {I4 + 1, I4 + 2});
635 Eigen::Matrix<real, dim, 2> diffKO = (diffRhoKO - 1. / UMeanXy(0) * diffRho * UMeanXy({I4 + 1, I4 + 2}).transpose()) / UMeanXy(0);
636 Eigen::Matrix<real, dim, dim> OmegaM = (diffU.transpose() - diffU) * 0.5;
637 real OmegaMag = OmegaM.norm() * std::sqrt(2);
638 real rho = UMeanXy(0);
639 real k = std::max(UMeanXy(I4 + 1) / rho, verySmallReal_4);
640 real omegaaa = std::max(UMeanXy(I4 + 2) / rho, verySmallReal_4);
641 real S = std::sqrt(SS.squaredNorm() / 2) + verySmallReal_4;
642 real nuPhy = muf / rho;
643 real CDKW = std::max(2 * rho * sigO2 / omegaaa * diffKO(EigenAll, 0).dot(diffKO(EigenAll, 1)), 1e-10);
644 real F1 = std::tanh(std::pow(
645 std::min(std::max(std::sqrt(k) / (betaStar * omegaaa * d), 500 * nuPhy / (sqr(d) * omegaaa)),
646 4 * rho * sigO2 * k / (CDKW * sqr(d))),
647 4));
648 real F2 = std::tanh(sqr(std::max(2 * std::sqrt(k) / (betaStar * omegaaa * d), 500 * nuPhy / (sqr(d) * omegaaa))));
649 real lRANS = std::sqrt(k) / (betaStar * omegaaa);
650 real FDES = std::max(1.0, lRANS / (lLES + verySmallReal) * (1 - F2));
651 // F2 = 0;
652 real mut = a1 * k / std::max(OmegaMag * F2, a1 * omegaaa) * rho; // use S/OmegaMag for SST: S: CFD++, OmegaMag: Turbulence Modeling Validation, Testing, and Developmen
653#if KW_SST_LIMIT_MUT == 1
654 mut = std::min(mut, 1e5 * muf); // CFL3D
655#endif
656 real nutHat = std::max(mut / rho, 1e-8);
657
658 Eigen::Matrix<real, dim, dim> rhoMuiuj = Eigen::Matrix<real, dim, dim>::Identity() * UMeanXy(I4 + 1) * (2. / 3.) - mut * SS;
659 real Pk = -(rhoMuiuj.array() * diffU.array()).sum();
660 real PkTilde = Pk;
661#if KW_SST_PROD_LIMITS == 1
662 PkTilde = std::max(PkTilde, verySmallReal);
663 PkTilde = std::min(Pk, 20 * betaStar * rho * k * omegaaa); // CFD++'s limiting: 10 times
664#endif
665
666 real gammaC = gamma1 * F1 + gamma1 * (1 - F1);
667 real sigK = sigK1 * F1 + sigK2 * (1 - F1);
668 real sigO = sigO1 * F1 + sigO2 * (1 - F1);
669 real beta = beta1 * F1 + beta2 * (1 - F1);
670 real POmega = gammaC / nutHat * Pk;
671#if KW_SST_PROD_OMEGA_VERSION == 1
672 POmega =
673 0.5 * gammaC * rho * ((SR2 - SR2.trace() / 3. * Eigen::Matrix<real, dim, dim>::Identity()).array() * SR2.array()).sum();
674#elif KW_SST_PROD_OMEGA_VERSION == 2
675 POmega =
676 gammaC * rho * sqr(OmegaMag);
677#endif
678 if (mode == 0)
679 {
680 source(I4 + 1) = PkTilde - betaStar * rho * k * omegaaa * FDES;
681 source(I4 + 2) = POmega - beta * rho * sqr(omegaaa) +
682 2 * (1 - F1) * rho * sigO2 / omegaaa * diffKO(EigenAll, 0).dot(diffKO(EigenAll, 1));
683 // source(I4 + 2) = POmega - beta * rho * sqr(omegaaa) +
684 // 2 * (1 - F1) * rho * sigO2 / omegaaa * diffKO(EigenAll, 0).dot(diffKO(EigenAll, 1));
685 }
686 else
687 {
688 source(I4 + 1) = betaStar * omegaaa * FDES;
689 source(I4 + 2) = 2 * beta * omegaaa;
690 }
691 }
692
693 /**
694 * @brief Compute turbulent viscosity mu_t from the Wilcox k-omega model.
695 *
696 * Evaluates the eddy viscosity mu_t = rho * k / omega_tilde, where omega_tilde
697 * depends on the selected model version (@c KW_WILCOX_VER):
698 * - Version 0 (original 1988): omega_tilde = omega (no stress limiter).
699 * - Version 1 (2006): omega_tilde = max(omega, C_lim * sqrt(S_{ij}S_{ij} / beta*)),
700 * applying the stress limiter with C_lim = 7/8.
701 * - Version 2 (simplified): omega_tilde = omega (no stress limiter).
702 *
703 * Optional upper bound of 1e5 * mu_laminar when @c KW_WILCOX_LIMIT_MUT is enabled.
704 *
705 * @tparam dim Spatial dimension (2 or 3).
706 * @tparam TU Eigen-compatible type for the conservative state vector.
707 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
708 * @param UMeanXy Conservative state vector [rho, rhoU, ..., rhoE, rho*k, rho*omega].
709 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
710 * @param muf Laminar (molecular) dynamic viscosity.
711 * @param d Wall distance (unused in this model but kept for interface consistency).
712 * @return Turbulent dynamic viscosity mu_t (>= 0).
713 */
714 template <int dim, class TU, class TDiffU>
715 real GetMut_KOWilcox(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d)
716 {
717 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
718 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
719 static const auto I4 = dim + 1;
720 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
721 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
722
723 real CLim = 7. / 8.;
724 real betaS = 0.09;
725 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
726 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
727 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
728 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
729 Eigen::Matrix<real, dim, dim> SR2 = diffU + diffU.transpose();
730 real rho = UMeanXy(0);
731 real k = std::max(UMeanXy(I4 + 1) / rho, sqr(verySmallReal_4));
732 real omegaaa = std::max(UMeanXy(I4 + 2) / rho, verySmallReal_4);
733#if KW_WILCOX_VER == 0
734 real omegaaaTut = omegaaa;
735#else
736 real omegaaaTut = std::max(omegaaa, CLim * std::sqrt(0.5 * SR2.squaredNorm() / betaS));
737#endif
738 real mut = k / omegaaaTut * rho;
739#if KW_WILCOX_LIMIT_MUT == 1
740 mut = std::min(mut, 1e5 * muf); // CFL3D
741#endif
742
743 if (std::isnan(mut) || !std::isfinite(mut))
744 {
745 std::cerr << k << " " << omegaaa << " " << mut << "\n";
746 DNDS_assert(false);
747 }
748 return mut;
749 }
750
751 /**
752 * @brief Compute the RANS viscous flux for the Wilcox k-omega transport equations.
753 *
754 * Evaluates the diffusion (viscous) flux contributions for the k and omega equations
755 * projected onto the face normal direction. The Wilcox model uses constant diffusion
756 * coefficients:
757 * - sigma_k = 0.5
758 * - sigma_omega = 0.5
759 *
760 * Effective diffusivity: mu_laminar + mu_t * sigma_{k,omega}.
761 *
762 * @tparam dim Spatial dimension (2 or 3).
763 * @tparam TU Eigen-compatible type for the conservative state vector.
764 * @tparam TN Eigen-compatible type for the face unit normal vector.
765 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, primitive variables).
766 * @tparam TVFlux Eigen-compatible type for the output viscous flux vector.
767 * @param UMeanXy Conservative state vector.
768 * @param DiffUxyPrim Gradient tensor of *primitive* variables, size dim × nVars.
769 * @param uNorm Outward-pointing face unit normal vector (length dim).
770 * @param mutIn Turbulent dynamic viscosity (precomputed).
771 * @param d Wall distance (unused here, kept for interface consistency).
772 * @param muf Laminar (molecular) dynamic viscosity.
773 * @param[out] vFlux Output viscous flux vector; entries at I4+1 and I4+2 are set.
774 */
775 template <int dim, class TU, class TN, class TDiffU, class TVFlux>
776 void GetVisFlux_KOWilcox(TU &&UMeanXy, TDiffU &&DiffUxyPrim, TN &&uNorm, real mutIn, real d, real muf, TVFlux &vFlux)
777 {
778 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
779 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
780 static const auto I4 = dim + 1;
781 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
782 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
783
784 real alpha = 13. / 25.;
785 real betaS = 0.09;
786 real sigK = 0.5;
787 real sigO = 0.5;
788 real Prt = 0.9;
789 real CLim = 7. / 8.;
790 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
791 Eigen::Matrix<real, dim, 1> diffRho = DiffUxyPrim(Seq012, {0});
792 Eigen::Matrix<real, dim, dim> diffU = DiffUxyPrim(Seq012, Seq123);
793 Eigen::Matrix<real, dim, dim> SR2 = diffU + diffU.transpose(); // 2 times strain rate
794 Eigen::Matrix<real, dim, dim> SS = SR2 - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity(); // 2 times shear strain rate
795 Eigen::Matrix<real, dim, 2> diffKO = DiffUxyPrim(Seq012, {I4 + 1, I4 + 2});
796 real rho = UMeanXy(0);
797 real k = std::max(UMeanXy(I4 + 1) / rho, sqr(verySmallReal_4));
798 real omegaaa = std::max(UMeanXy(I4 + 2) / rho, verySmallReal_4);
799#if KW_WILCOX_VER == 0
800 real omegaaaTut = omegaaa;
801#else
802 real omegaaaTut = std::max(omegaaa, CLim * std::sqrt(0.5 * SR2.squaredNorm() / betaS));
803#endif
804 real mut = k / omegaaaTut * rho;
805#if KW_WILCOX_LIMIT_MUT == 1
806 mut = std::min(mut, 1e5 * muf); // CFL3D
807#endif
808
809 vFlux(I4 + 1) = diffKO(Seq012, 0).dot(uNorm) * (muf + mut * sigK);
810 vFlux(I4 + 2) = diffKO(Seq012, 1).dot(uNorm) * (muf + mut * sigO);
811 }
812
813 /**
814 * @brief Compute source terms for the Wilcox k-omega transport equations.
815 *
816 * Evaluates production, destruction, and cross-diffusion source contributions for
817 * the k and omega equations following the Wilcox k-omega model. The version-dependent
818 * behavior is controlled by @c KW_WILCOX_VER:
819 *
820 * - **Version 0 (1988):** alpha = 5/9, beta = 3/40, no stress limiter, no cross-diffusion.
821 * - **Version 1 (2006):** alpha = 13/25, beta = f_beta * beta_0, with the f_beta correction
822 * based on Chi_omega (mean rotation / strain invariant), stress limiter C_lim = 7/8,
823 * and cross-diffusion sigma_d term (sigma_d = 0.125 when grad(k) . grad(omega) > 0).
824 * - **Version 2 (simplified):** alpha = 13/25, beta = 0.075, vorticity-based production,
825 * no stress limiter, no cross-diffusion.
826 *
827 * Production of k is optionally capped at 20 * beta* * rho * k * omega (CFL3D-style)
828 * when @c KW_WILCOX_PROD_LIMITS is enabled.
829 *
830 * When @p mode == 0, the full source vector is returned:
831 * - source(I4+1) = P_k - beta* * rho * k * omega
832 * - source(I4+2) = P_omega - beta * rho * omega^2 + sigma_d / omega * grad(k) . grad(omega)
833 *
834 * When @p mode == 1, the implicit diagonal contribution (positive) is returned:
835 * - source(I4+1) = beta* * omega
836 * - source(I4+2) = 2 * beta * omega
837 *
838 * @tparam dim Spatial dimension (2 or 3).
839 * @tparam TU Eigen-compatible type for the conservative state vector.
840 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
841 * @tparam TSource Eigen-compatible type for the output source vector.
842 * @param UMeanXy Conservative state vector.
843 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
844 * @param muf Laminar (molecular) dynamic viscosity.
845 * @param d Wall distance (unused here, kept for interface consistency).
846 * @param[out] source Output source vector; entries at I4+1 and I4+2 are set.
847 * @param mode Source computation mode: 0 = full source, 1 = implicit diagonal (destruction only).
848 */
849 template <int dim, class TU, class TDiffU, class TSource>
850 void GetSource_KOWilcox(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d, TSource &source, int mode)
851 {
852 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
853 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
854 static const auto I4 = dim + 1;
855 static const auto verySmallReal_3 = std::pow(verySmallReal, 1. / 3);
856 static const auto verySmallReal_4 = std::pow(verySmallReal, 1. / 4);
857#if KW_WILCOX_VER == 0
858 real alpha = 5. / 9.;
859#else
860 real alpha = 13. / 25.;
861#endif
862 real betaS = 0.09;
863 real sigK = 0.5;
864 real sigO = 0.5;
865 real Prt = 0.9;
866 real CLim = 7. / 8.;
867 real betaO = 0.0708;
868 real kappa = 0.41;
869
870 Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
871 Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
872 Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
873 Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
874 Eigen::Matrix<real, dim, dim> SR2 = diffU + diffU.transpose(); // 2 times strain rate
875 Eigen::Matrix<real, dim, dim> OmegaM2 = diffU.transpose() - diffU; // 2 times rotation
876 Eigen::Matrix<real, dim, dim> SS = SR2 - (2. / 3.) * diffU.trace() * Eigen::Matrix<real, dim, dim>::Identity(); // 2 times shear strain rate
877 Eigen::Matrix<real, dim, 2> diffRhoKO = DiffUxy(Seq012, {I4 + 1, I4 + 2});
878 Eigen::Matrix<real, dim, 2> diffKO = (diffRhoKO - 1. / UMeanXy(0) * diffRho * UMeanXy({I4 + 1, I4 + 2}).transpose()) / UMeanXy(0);
879 real rho = UMeanXy(0);
880 real k = std::max(UMeanXy(I4 + 1) / rho, sqr(verySmallReal_4)); // make nu -> 0 when k,O->0
881 real omegaaa = std::max(UMeanXy(I4 + 2) / rho, verySmallReal_4);
882#if KW_WILCOX_VER == 0 || KW_WILCOX_VER == 2
883 real omegaaaTut = omegaaa;
884#else
885 real omegaaaTut = std::max(omegaaa, CLim * std::sqrt(0.5 * SR2.squaredNorm() / betaS));
886#endif
887 real mut = k / omegaaaTut * rho;
888#if KW_WILCOX_LIMIT_MUT == 1
889 mut = std::min(mut, 1e5 * muf); // CFL3D
890#endif
891
892 real ChiOmega = std::abs(((OmegaM2 * OmegaM2).array() * SR2.array()).sum() * 0.125 / cube(betaS * omegaaa));
893 real fBeta = (1 + 85 * ChiOmega) / (1 + 100 * ChiOmega);
894#if KW_WILCOX_VER == 0
895 real beta = 3. / 40.;
896#elif KW_WILCOX_VER == 1
897 real beta = fBeta * betaO;
898#else
899 real beta = 0.075;
900#endif
901 real crossDiff = diffKO(EigenAll, 0).dot(diffKO(EigenAll, 1));
902#if KW_WILCOX_VER == 0 || KW_WILCOX_VER == 2
903 real SigD = 0;
904#else
905 real SigD = crossDiff > 0 ? 0.125 : 0;
906#endif
907
908 Eigen::Matrix<real, dim, dim> rhoMuiuj = Eigen::Matrix<real, dim, dim>::Identity() * UMeanXy(I4 + 1) * (2. / 3.) - mut * SS;
909 real Pk = -(rhoMuiuj.array() * diffU.array()).sum();
910 real POmega = alpha * omegaaa / k * Pk;
911
912#if KW_WILCOX_VER == 2
913 Pk = OmegaM2.squaredNorm() * 0.5 * mut;
914 real gam = beta / betaS - sqr(kappa) / (sigO * std::sqrt(betaS));
915 POmega = gam * rho * OmegaM2.squaredNorm() * 0.5; // CFL3D ???
916#endif
917
918#if KW_WILCOX_PROD_LIMITS == 1
919 // Pk = std::min(Pk, OmegaM2.squaredNorm() * 0.5 * mut); // compare with CFL3D approx
920
921 Pk = std::max(Pk, verySmallReal);
922 Pk = std::min(Pk, betaS * rho * k * omegaaa * 20); // CFL3D
923#endif
924
925 // POmega = std::min(POmega, beta * rho * sqr(omegaaa) * 20); // CFL3D
926
927 if (mode == 0)
928 {
929 source(I4 + 1) = Pk - betaS * rho * k * omegaaa;
930 source(I4 + 2) = POmega - beta * rho * sqr(omegaaa) + SigD / omegaaa * crossDiff;
931 }
932 else
933 {
934 source(I4 + 1) = betaS * omegaaa;
935 source(I4 + 2) = 2 * beta * omegaaa;
936 }
937 }
938
939 /**
940 * @brief Compute source terms for the Spalart-Allmaras one-equation turbulence model.
941 *
942 * Evaluates the full SA source term including production, destruction, diffusion,
943 * and optional DES/DDES/IDDES/WMLES length-scale modification. The SA model solves
944 * a single transport equation for the modified kinematic viscosity nuTilde (stored
945 * as rho * nuTilde at index I4+1 in the conservative state vector).
946 *
947 * Key terms computed:
948 * - **Vorticity magnitude** S = ||Omega|| * sqrt(2), with optional rotation correction
949 * (SRotCor = c_rot * min(0, S - SS)) when @p rotCor is enabled.
950 * - **Modified vorticity** Sh = S + Sbar, where Sbar = nuTilde * f_{v2} / (kappa^2 * d^2).
951 * - **Production** P = c_{b1} * (1 - f_{t2}) * Sh * nuTilde.
952 * - **Destruction** D = (c_{w1} * f_w - c_{b1}/kappa^2 * f_{t2}) * (nuTilde / d)^2.
953 * - **Diffusion** term: c_{b2}/sigma * |grad(nuTilde)|^2 (conservative form contribution).
954 * - **f_{t2} laminar suppression** term controlled by @c SA_USE_FT2_TERM.
955 *
956 * DES/DDES/IDDES length-scale modification:
957 * - **DESMode 0 (RANS):** l_DES = d (pure RANS, wall distance).
958 * - **DESMode 1 (IDDES):** Blends l_RANS and l_LES using the IDDES shielding functions
959 * f_d, f_B, f_e, with psi correction for grid-induced separation avoidance.
960 * - **DESMode 2 (WMLES):** Uses the wall-modeled LES length scale
961 * l_WMLES = f_B * (1 + f_e) * l_RANS + (1 - f_B) * l_LES * psi.
962 *
963 * When @p mode == 0, the full source is returned at source(I4+1).
964 * When @p mode == 1, the implicit diagonal contribution -dS/dU (positive, suitable for
965 * augmenting the implicit operator) is returned at source(I4+1).
966 *
967 * @tparam dim Spatial dimension (2 or 3).
968 * @tparam TU Eigen-compatible type for the conservative state vector.
969 * @tparam TDiffU Eigen-compatible type for the gradient tensor (dim × nVars, conservative variables).
970 * @tparam TSource Eigen-compatible type for the output source vector.
971 * @param UMeanXy Conservative state vector [rho, rhoU, ..., rhoE, rho*nuTilde].
972 * @param DiffUxy Gradient tensor of conservative variables, size dim × nVars.
973 * @param muRef Reference viscosity scale used to non-dimensionalise nuTilde
974 * (nuTilde_physical = UMeanXy(I4+1) * muRef / rho).
975 * @param mufPhy Physical (dimensional) laminar dynamic viscosity.
976 * @param gamma Ratio of specific heats (used for pressure gradient in optional helicity correction).
977 * @param d Wall distance.
978 * @param lLES LES length scale for DES shielding (set to a large value to disable DES).
979 * @param hMax Maximum cell size (used in IDDES alpha parameter: alpha = 0.25 - d/hMax).
980 * @param DESMode DES operating mode: 0 = pure RANS, 1 = IDDES, 2 = WMLES.
981 * @param[out] source Output source vector; entry at I4+1 is set.
982 * @param rotCor Rotation correction flag: 0 = disabled, nonzero = enabled (c_rot = 2.0).
983 * @param mode Source computation mode: 0 = full source, 1 = implicit diagonal (positive).
984 */
985 template <int dim, class TU, class TDiffU, class TSource>
986 void GetSource_SA(TU &&UMeanXy, TDiffU &&DiffUxy, real muRef, real mufPhy, real gamma,
987 real d,
988 real lLES,
989 real hMax,
990 int DESMode,
991 TSource &source, int rotCor, int mode)
992 {
993 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
994 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
995 static const auto I4 = dim + 1;
996
997 static const real cb1 = 0.1355;
998 static const real cb2 = 0.622;
999 static const real sigma = SA::sigma;
1000 static const real cnu1 = SA::cnu1;
1001 static const real cnu2 = 0.7;
1002 static const real cnu3 = 0.9;
1003 static const real cw2 = 0.3;
1004 static const real cw3 = 2;
1005 static const real kappa = 0.41;
1006 static const real rlim = 10;
1007 static const real cw1 = cb1 / sqr(kappa) + (1 + cb2) / sigma;
1008
1009#if SA_USE_FT2_TERM
1010 static const real ct3 = 1.2;
1011 static const real ct4 = 0.5;
1012#else
1013 static const real ct3 = 0.0;
1014 static const real ct4 = 0.0;
1015#endif
1016
1017 const real nuh = UMeanXy(I4 + 1) * muRef / UMeanXy(0);
1018
1019 const real Chi = (UMeanXy(I4 + 1) * muRef / mufPhy);
1020 const real Chi3 = cube(Chi);
1021 const real fnu1 = Chi3 / (Chi3 + cube(cnu1));
1022 const real fnu2 = 1 - Chi / (1 + Chi * fnu1);
1023
1024 const Eigen::Matrix<real, dim, 1> velo = UMeanXy(Seq123) / UMeanXy(0);
1025 const Eigen::Matrix<real, dim, 1> diffRhoNu = DiffUxy(Seq012, {I4 + 1}) * muRef;
1026 const Eigen::Matrix<real, dim, 1> diffRho = DiffUxy(Seq012, {0});
1027 const Eigen::Matrix<real, dim, 1> diffNu = (diffRhoNu - nuh * diffRho) / UMeanXy(0);
1028 const Eigen::Matrix<real, dim, dim> diffRhoU = DiffUxy(Seq012, Seq123);
1029 const Eigen::Matrix<real, dim, dim> diffU = (diffRhoU - diffRho * velo.transpose()) / UMeanXy(0);
1030
1031 const Eigen::Matrix<real, dim, dim> Omega = 0.5 * (diffU.transpose() - diffU);
1032#ifndef USE_ABS_VELO_IN_ROTATION
1033 if (settings.frameConstRotation.enabled)
1034 Omega += Geom::CrossVecToMat(settings.frameConstRotation.vOmega())(Seq012, Seq012); // to static frame rotation
1035#endif
1036 const real S = Omega.norm() * std::sqrt(2.0); // is omega's magnitude
1037 const real SS = (diffU + diffU.transpose()).norm() * (1. / std::sqrt(2.0)); // is sqrt(2) * strainrate's norm
1038 real SRotCor = 0.0;
1039 if (rotCor)
1040 {
1041 // 2 is recommended but we use 1 to avoid negative production, see Diskin, Boris, Yi Liu, and Marshall C. Galbraith. "High-Fidelity CFD Verification Workshop 2024: Spalart-Allmaras QCR2000-R Turbulence Model." AIAA Scitech 2023 Forum. 2023.
1042 // we now use 2.0
1043 const real cRot = 2.0;
1044 SRotCor += cRot * std::min(0.0, S - SS);
1045 }
1046 real diffUNorm = diffU.norm();
1047
1048 const real ft2 = ct3 * std::exp(-ct4 * sqr(Chi));
1049 // {
1050 // Eigen::Matrix<real, dim, dim> sHat = 0.5 * (diffU.transpose() + diffU);
1051 // real sHatSqr = 2 * sHat.squaredNorm();
1052 // real rStar = std::sqrt(sHatSqr) / S;
1053 // real DD = 0.5 * (sHatSqr + sqr(S));
1054 // ! // need second derivatives for rotation term !(CFD++ user manual)
1055 // }
1056
1057 real lDES = d;
1058 // DDES shield
1059 {
1060 real lRANS = d;
1061 real fwStar = .424;
1062 real nuTur = mufPhy / UMeanXy(0) * std::max((Chi * fnu1), 0.0);
1063 real nufPhy = mufPhy / UMeanXy(0);
1064 real rd = (nuTur + nufPhy) / (sqr(kappa) * sqr(d) * std::max(smallReal, diffUNorm));
1065 real fd = 1. - std::tanh(cube(8 * rd));
1066 real psiSqr = 0.0;
1067 if (Chi <= 0)
1068 psiSqr = 100.0;
1069 else
1070 psiSqr = std::min(100.0, (1 - cb1 / (cw1 * sqr(kappa) * fwStar) * (ft2 + (1 - ft2) * fnu2)) /
1071 (fnu1 * std::max(smallReal, 1 - ft2)));
1072 // if (psiSqr < 1.01)
1073 // {
1074 // std::cout << psiSqr << "Chi " << Chi << " fnu1 " << fnu1 << " xx " << (fnu1 * std::max(smallReal, 1 - ft2)) << " xx " << cb1 / (cw1 * sqr(kappa) * fwStar) << std::endl;
1075 // }
1076 real psi = std::sqrt(std::max(psiSqr, 1.0));
1077 //! note that psi has lower bound of 1
1078 // psi = 1.0;
1079 lDES = lRANS - fd * std::max(0., lRANS - lLES * psi);
1080
1081 // IDDES switch
1082 if (DESMode == 1 || DESMode == 2)
1083 {
1084 const real alphaIDDES = 0.25 - d / hMax;
1085 const real fB = std::min(2 * std::exp(-9. * sqr(alphaIDDES)), 1.0);
1086
1087 const real cl = 3.55;
1088 const real ct = 1.63;
1089 const real rdt = nuTur / (sqr(kappa) * sqr(d) * std::max(smallReal, diffUNorm));
1090 const real rdl = nufPhy / (sqr(kappa) * sqr(d) * std::max(smallReal, diffUNorm));
1091 const real ft = std::tanh(cube(sqr(ct) * rdt));
1092 const real fl_tmp = sqr(sqr(cl) * rdl);
1093 const real fl = std::tanh(sqr(sqr(fl_tmp)) * fl_tmp); // power 5
1094
1095 const real fe1 = 2. * std::exp(-(alphaIDDES >= 0 ? 11.09 : 9.0) * sqr(alphaIDDES));
1096 const real fe2 = 1. - std::max(ft, fl);
1097 const real fe = std::max((fe1 - 1.), 0.) * psi * fe2;
1098 if (DESMode == 2)
1099 {
1100 const real lWMLES = fB * (1 + fe) * lRANS + (1 - fB) * lLES * psi;
1101 lDES = lWMLES;
1102 }
1103 else
1104 {
1105 const real fdTilde = std::max(fB, std::tanh(cube(8 * rdt)));
1106 real lIDDES = fdTilde * (1 + fe) * lRANS + (1 - fdTilde) * lLES * psi;
1107 lIDDES = std::max(lLES * smallReal, lIDDES);
1108 lDES = lIDDES;
1109 }
1110 }
1111 }
1112 d = lDES; //! super subs
1113
1114 const real Sbar = nuh / (sqr(kappa) * sqr(d)) * fnu2;
1115
1116 real Sh = 0.;
1117
1118 { // Lee, K., Wilson, M., and Vahdati, M. (April 16, 2018). "Validation of a Numerical Model for Predicting Stalled Flows in a Low-Speed Fan—Part I: Modification of Spalart–Allmaras Turbulence Model." ASME. J. Turbomach. May 2018; 140(5): 051008.
1119 // real betaSCor = 1;
1120 // real ch1 = 0.5;
1121 // real ch2 = 0.7;
1122 // real a1 = 3; //! is this good?
1123 // real a2 = 3;
1124 // Eigen::Vector<real, dim> diffP = (DiffUxy(Seq012, I4) - diffRhoU * velo - UMeanXy(0) * diffU * velo) * (gamma - 1);
1125 // real veloN = velo.norm();
1126 // Eigen::Vector<real, dim> uN = velo / (veloN + verySmallReal);
1127 // real pStar = diffP.dot(uN) / (sqr(UMeanXy(0)) * sqr(veloN) * veloN) * mufPhy;
1128 // Geom::tPoint omegaV = Geom::CrossMatToVec(Omega);
1129 // real HStar = omegaV.dot(velo) / (veloN * omegaV.norm() + verySmallReal);
1130 // real Cs = ch1 * std::tanh(a1 * sqr(pStar)) / std::tanh(1.0) + 1;
1131 // real Cvh = ch2 * std::tanh(a2 * sqr(HStar)) / std::tanh(1.0) + 1;
1132 // betaSCor = Cs * Cvh;
1133
1134 // S *= betaSCor;
1135 }
1136#ifdef USE_NS_SA_NEGATIVE_MODEL
1137 if (Sbar < -cnu2 * S)
1138 Sh = S + S * (sqr(cnu2) * S + cnu3 * Sbar) / ((cnu3 - 2 * cnu2) * S - Sbar);
1139 else //*negative fix
1140#endif
1141 Sh = S + Sbar;
1142 // here r is used for fw, we use real d instead of lDES
1143 const real r = std::min(nuh / (Sh * sqr(kappa * d) + verySmallReal), rlim);
1144 const real g = r + cw2 * (std::pow(r, 6) - r);
1145 const real fw = g * std::pow((1 + std::pow(cw3, 6)) / (std::pow(g, 6) + std::pow(cw3, 6)), 1. / 6.);
1146
1147 // if (d < 0.01)
1148 // std::cout << d << " " << lDES << " " << lLES << std::endl;
1149 // DDES
1150
1151#ifdef USE_NS_SA_NEGATIVE_MODEL
1152 real D = (cw1 * fw - cb1 / sqr(kappa) * ft2) * sqr(nuh / lDES); //! modified >>
1153 real P = cb1 * (1 - ft2) * (Sh + SRotCor) * nuh; //! modified >>
1154#else
1155 real D = (cw1 * fw - cb1 / sqr(kappa) * ft2) * sqr(nuh / lDES);
1156 real P = cb1 * (1 - ft2) * (Sh + SRotCor) * nuh;
1157#endif
1158 real fn = 1;
1159#ifdef USE_NS_SA_NEGATIVE_MODEL
1160 if (UMeanXy(I4 + 1) < 0)
1161 {
1162 real Chi = UMeanXy(I4 + 1) * muRef / mufPhy;
1163 fn = (SA::cn1 + std::pow(Chi, 3)) / (SA::cn1 - std::pow(Chi, 3));
1164 D = -cw1 * sqr(nuh / lDES);
1165 P = cb1 * (1 - ct3) * std::abs(S + SRotCor) * nuh;
1166 }
1167#endif
1168
1169 if (mode == 0)
1170 source(I4 + 1) = UMeanXy(0) * (P - D + diffNu.squaredNorm() * cb2 / sigma) / muRef -
1171 (UMeanXy(I4 + 1) * fn * muRef + mufPhy) / (UMeanXy(0) * sigma) * diffRho.dot(diffNu) / muRef;
1172 else
1173 source(I4 + 1) = -std::min(UMeanXy(0) * (std::min(P, 0.0) * 1 - D * 2) / muRef / (std::abs(UMeanXy(I4 + 1)) + verySmallReal), -verySmallReal);
1174
1175 if (!source.allFinite())
1176 {
1177 std::cout << P << std::endl;
1178 std::cout << D << std::endl;
1179 std::cout << UMeanXy(0) << std::endl;
1180 std::cout << Sh << std::endl;
1181 std::cout << nuh << std::endl;
1182 std::cout << g << std::endl;
1183 std::cout << r << std::endl;
1184 std::cout << S << std::endl;
1185 std::cout << d << std::endl;
1186 std::cout << fnu2 << std::endl;
1187 std::cout << mufPhy << std::endl;
1188 DNDS_assert(false);
1189 }
1190 }
1191}
DNDS_assert
#define DNDS_assert(expr)
Debug-only assertion (compiled out when DNDS_NDEBUG is defined). Prints the expression + file/line + ...
Definition Errors.hpp:108
Euler.hpp
Core type definitions, enumerations, and distributed array wrappers for compressible Navier-Stokes / ...
DNDS::Euler::RANS
RANS turbulence model functions (eddy viscosity, viscous flux, source terms).
Definition RANS_ke.hpp:52
DNDS::Euler::RANS::GetSourceJacobianDiag_RealizableKe
void GetSourceJacobianDiag_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, TSource &source)
Compute the analytical diagonal of the source Jacobian for the realizable k-epsilon model.
Definition RANS_ke.hpp:335
DNDS::Euler::RANS::GetSource_SA
void GetSource_SA(TU &&UMeanXy, TDiffU &&DiffUxy, real muRef, real mufPhy, real gamma, real d, real lLES, real hMax, int DESMode, TSource &source, int rotCor, int mode)
Compute source terms for the Spalart-Allmaras one-equation turbulence model.
Definition RANS_ke.hpp:986
DNDS::Euler::RANS::GetVisFlux_RealizableKe
void GetVisFlux_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxyPrim, TN &&uNorm, real mut, real d, real muPhy, TVFlux &vFlux)
Compute the RANS viscous flux for the realizable k-epsilon transport equations.
Definition RANS_ke.hpp:148
DNDS::Euler::RANS::GetSource_RealizableKe
void GetSource_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d, TSource &source, int mode)
Compute source terms for the realizable k-epsilon transport equations.
Definition RANS_ke.hpp:192
DNDS::Euler::RANS::GetMut_KOWilcox
real GetMut_KOWilcox(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d)
Compute turbulent viscosity mu_t from the Wilcox k-omega model.
Definition RANS_ke.hpp:715
DNDS::Euler::RANS::GetVisFlux_SST
void GetVisFlux_SST(TU &&UMeanXy, TDiffU &&DiffUxyPrim, TN &&uNorm, real mutIn, real d, real muf, TVFlux &vFlux)
Compute the RANS viscous flux for the k-omega SST transport equations.
Definition RANS_ke.hpp:511
DNDS::Euler::RANS::SolveZeroGradEquilibrium
std::tuple< real, real > SolveZeroGradEquilibrium(TU &u, real muPhy)
Attempt to find equilibrium epsilon for zero-gradient (freestream) conditions via Newton iteration.
Definition RANS_ke.hpp:275
DNDS::Euler::RANS::GetMut_RealizableKe
real GetMut_RealizableKe(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d)
Compute turbulent viscosity mu_t from the Shih et al. realizable k-epsilon model.
Definition RANS_ke.hpp:91
DNDS::Euler::RANS::GetSourceJacobianDiag_RealizableKe_ND
void GetSourceJacobianDiag_RealizableKe_ND(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, TSource &source)
Compute the numerical (finite-difference) diagonal of the source Jacobian for the realizable k-epsilo...
Definition RANS_ke.hpp:400
DNDS::Euler::RANS::GetSource_SST
void GetSource_SST(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d, real lLES, TSource &source, int mode)
Compute source terms for the k-omega SST transport equations.
Definition RANS_ke.hpp:609
DNDS::Euler::RANS::GetMut_SST
real GetMut_SST(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d)
Compute turbulent viscosity mu_t from Menter's k-omega SST model.
Definition RANS_ke.hpp:443
DNDS::Euler::RANS::GetSource_KOWilcox
void GetSource_KOWilcox(TU &&UMeanXy, TDiffU &&DiffUxy, real muf, real d, TSource &source, int mode)
Compute source terms for the Wilcox k-omega transport equations.
Definition RANS_ke.hpp:850
DNDS::Euler::RANS::GetVisFlux_KOWilcox
void GetVisFlux_KOWilcox(TU &&UMeanXy, TDiffU &&DiffUxyPrim, TN &&uNorm, real mutIn, real d, real muf, TVFlux &vFlux)
Compute the RANS viscous flux for the Wilcox k-omega transport equations.
Definition RANS_ke.hpp:776
DNDS::Geom::CrossVecToMat
tGPoint CrossVecToMat(const tPoint &axn)
Definition Geometric.hpp:195
DNDS::cube
constexpr T cube(const T &a)
a * a * a, constexpr.
Definition Defines.hpp:525
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::sqr
constexpr T sqr(const T &a)
a * a, constexpr. Works for all arithmetic types.
Definition Defines.hpp:517
DNDS::real
double real
Canonical floating-point scalar used throughout DNDSR (double precision).
Definition Defines.hpp:105
DNDS::smallReal
DNDS_CONSTANT const real smallReal
Loose lower bound (for iterative-solver tolerances etc.).
Definition Defines.hpp:196
r
tVec r(NCells)
alpha
real alpha
Definition test_GasThermo.cpp:454
velo
Eigen::Vector3d velo
Definition test_GasThermo.cpp:282