DNDSR 0.1.0.dev1+gcd065ad
Distributed Numeric Data Structure for CFV
Loading...
Searching...
No Matches
Gas.hpp
Go to the documentation of this file.
1/** @file Gas.hpp
2 * @brief Ideal-gas Riemann solvers, flux functions, and thermodynamic utilities
3 * for the compressible Euler / Navier-Stokes equations.
4 *
5 * Provides:
6 * - Right / left eigenvector matrices for the 1-D Euler system.
7 * - Conservative ↔ primitive variable conversions for an ideal gas.
8 * - Inviscid flux evaluation in Cartesian and face-normal directions (single
9 * and batched variants).
10 * - Roe, HLLC, and HLLEP approximate Riemann solvers with multiple entropy
11 * fix strategies selectable at compile time via the @c eigScheme template
12 * parameter, and at run time through RiemannSolverType dispatch.
13 * - Viscous (Navier-Stokes) flux with Sutherland viscosity, optional QCR
14 * correction, and adiabatic wall treatment.
15 *
16 * All functions are templated on the spatial dimension @c dim (2 or 3) and
17 * accept Eigen expression types so that fixed-size and dynamic vectors /
18 * matrices are handled without copies.
19 */
20#pragma once
21#include "DNDS/Defines.hpp"
22#include "DNDS/IdealGasPhysics.hpp"
23
24#include "DNDS/JsonUtil.hpp"
25#include "DNDS/ConfigEnum.hpp"
26#include <fmt/core.h>
27
28namespace DNDS::Euler::Gas
29{
30 using tVec = Eigen::Vector3d; ///< Convenience alias for 3-D real vector.
31 using tVec2 = Eigen::Vector2d; ///< Convenience alias for 2-D real vector.
32
33 /// Maximum column count for batched (SIMD-style) Riemann solver evaluation.
34 static const int MaxBatch = 16;
35
36 // Shared constants for Riemann solver entropy fixes and low-dissipation schemes.
37 // Delegated from DNDS::IdealGas.
38 static constexpr real kScaleHartenYee = IdealGas::kScaleHartenYee; ///< Base Harten-Yee entropy-fix threshold (0.05).
39 static constexpr real kScaleLD = IdealGas::kScaleLD; ///< Low-dissipation speed-of-sound scaling (0.2).
40 static constexpr real kScaleHFix = IdealGas::kScaleHFix; ///< H-correction transverse-wave scaling (0.25).
41
42 /**
43 * @brief Selects the approximate Riemann solver and its entropy-fix variant.
44 *
45 * | Enumerator | Solver / fix scheme |
46 * |------------|--------------------------------------------------------------|
47 * | Roe | Standard Roe with Harten-Yee entropy fix (eigScheme 0). |
48 * | HLLC | Harten-Lax-van Leer-Contact (known accuracy issues, see IV). |
49 * | HLLEP | HLL with Enhanced Pressure (type 0). |
50 * | HLLEP_V1 | HLLEP variant 1 (type 1). |
51 * | Roe_M1 | Roe + cLLF entropy fix (eigScheme 1). |
52 * | Roe_M2 | Roe + Lax-Friedrichs / vanilla Lax (eigScheme 2, early exit).|
53 * | Roe_M3 | Roe + LD (low-dissipation) Roe fix (eigScheme 3). |
54 * | Roe_M4 | Roe + ID (intermediate-dissipation) Roe fix (eigScheme 4). |
55 * | Roe_M5 | Roe + LD cLLF fix (eigScheme 5). |
56 * | Roe_M6 | Roe + H-correction only (eigScheme 6). |
57 * | Roe_M7 | Roe + Harten-Yee fix only, no H-correction (eigScheme 7). |
58 * | Roe_M8 | Roe + H-correction + Harten-Yee fix (eigScheme 8). |
59 * | Roe_M9 | Reserved (eigScheme 9, currently asserts false). |
60 */
61 enum RiemannSolverType
62 {
63 UnknownRS = 0,
64 Roe = 1,
65 HLLC = 2,
66 HLLEP = 3,
67 HLLEP_V1 = 21,
68 Roe_M1 = 11,
69 Roe_M2 = 12,
70 Roe_M3 = 13,
71 Roe_M4 = 14,
72 Roe_M5 = 15,
73 Roe_M6 = 16,
74 Roe_M7 = 17,
75 Roe_M8 = 18,
76 Roe_M9 = 19,
77 };
78
79 DNDS_DEFINE_ENUM_JSON(
80 RiemannSolverType,
81 {
82 {UnknownRS, "UnknownRS"},
83 {Roe, "Roe"},
84 {HLLC, "HLLC"},
85 {HLLEP, "HLLEP"},
86 {HLLEP_V1, "HLLEP_V1"},
87 {Roe_M1, "Roe_M1"},
88 {Roe_M2, "Roe_M2"},
89 {Roe_M3, "Roe_M3"},
90 {Roe_M4, "Roe_M4"},
91 {Roe_M5, "Roe_M5"},
92 {Roe_M6, "Roe_M6"},
93 {Roe_M7, "Roe_M7"},
94 {Roe_M8, "Roe_M8"},
95 {Roe_M9, "Roe_M9"},
96 })
97
98 /**
99 * @brief Fills the right eigenvector matrix for the 1-D Euler system in the
100 * x-direction.
101 *
102 * The matrix has (dim+2) columns, one per characteristic wave:
103 * - column 0: left-running acoustic wave (λ = u − a)
104 * - column 1: entropy wave (λ = u)
105 * - columns 2..dim: shear waves (λ = u)
106 * - column dim+1: right-running acoustic wave (λ = u + a)
107 *
108 * @warning @p ReV must be pre-allocated to size (dim+2)×(dim+2).
109 *
110 * @tparam dim Spatial dimension (2 or 3).
111 * @tparam TVec Velocity vector type (Eigen expression, length @c dim).
112 * @tparam TeV Eigenvector matrix type (Eigen expression, (dim+2)×(dim+2)).
113 * @param velo Roe-averaged (or cell-centred) velocity vector.
114 * @param Vsqr Squared velocity magnitude (= velo.squaredNorm()).
115 * @param H Total specific enthalpy.
116 * @param a Speed of sound.
117 * @param[out] ReV Right eigenvector matrix, overwritten on output.
118 */
119 template <int dim = 3, class TVec, class TeV>
120 inline void EulerGasRightEigenVector(const TVec &velo, real Vsqr, real H, real a, TeV &ReV)
121 {
122 ReV.setZero();
123 ReV(0, {0, 1, dim + 1}).setConstant(1);
124 ReV(Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>), Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>))
125 .diagonal()
126 .setConstant(1);
127
128 ReV(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), {0, 1, dim + 1}).colwise() = velo;
129 ReV(1, 0) -= a;
130 ReV(1, dim + 1) += a;
131
132 // Last Row
133 ReV(dim + 1, 0) = H - velo(0) * a;
134 ReV(dim + 1, dim + 1) = H + velo(0) * a;
135 ReV(dim + 1, 1) = 0.5 * Vsqr;
136
137 ReV(dim + 1, Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>)) =
138 velo(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim - 1>));
139 }
140
141 /**
142 * @brief Fills the left eigenvector matrix (inverse of the right eigenvector
143 * matrix) for the 1-D Euler system in the x-direction.
144 *
145 * Satisfies LeV * ReV = I for the (dim+2)×(dim+2) system.
146 *
147 * @warning @p LeV must be pre-allocated to size (dim+2)×(dim+2).
148 *
149 * @tparam dim Spatial dimension (2 or 3).
150 * @tparam TVec Velocity vector type.
151 * @tparam TeV Eigenvector matrix type.
152 * @param velo Roe-averaged (or cell-centred) velocity vector.
153 * @param Vsqr Squared velocity magnitude.
154 * @param H Total specific enthalpy.
155 * @param a Speed of sound.
156 * @param gamma Ratio of specific heats (γ).
157 * @param[out] LeV Left eigenvector matrix, overwritten on output.
158 */
159 template <int dim = 3, class TVec, class TeV>
160 inline void EulerGasLeftEigenVector(const TVec &velo, real Vsqr, real H, real a, real gamma, TeV &LeV)
161 {
162 LeV.setZero();
163 real gammaBar = gamma - 1;
164 LeV(0, 0) = H + a / gammaBar * (velo(0) - a);
165 LeV(0, 1) = -velo(0) - a / gammaBar;
166 LeV(0, Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>)) =
167 -velo(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim - 1>));
168 LeV(0, dim + 1) = 1;
169
170 LeV(1, 0) = -2 * H + 4 / gammaBar * (a * a);
171 LeV(1, Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) = velo.transpose() * 2;
172 LeV(1, dim + 1) = -2;
173
174 LeV(Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>), 0) =
175 velo(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim - 1>)) * (-2 * (a * a) / gammaBar);
176 LeV(Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>), Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>))
177 .diagonal()
178 .setConstant(2 * (a * a) / gammaBar);
179
180 LeV(dim + 1, 0) = H - a / gammaBar * (velo(0) + a);
181 LeV(dim + 1, 1) = -velo(0) + a / gammaBar;
182 LeV(dim + 1, Eigen::seq(Eigen::fix<2>, Eigen::fix<dim>)) =
183 -velo(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim - 1>));
184 LeV(dim + 1, dim + 1) = 1;
185
186 LeV *= gammaBar / (2 * a * a);
187 }
188
189 /// @brief Thin wrapper delegating to IdealGas::IdealGasThermal.
190 inline void IdealGasThermal(
191 real E, real rho, real vSqr, real gamma, real &p, real &asqr, real &H)
192 {
193 IdealGas::IdealGasThermal(E, rho, vSqr, gamma, p, asqr, H);
194 }
195
196 /**
197 * @brief Pre-computed Roe-averaged quantities shared by all Riemann solvers.
198 *
199 * Holds the Lm/Rm thermodynamic state and the Roe averages derived from them.
200 * Computed once by ComputeRoePreamble(), consumed by HLLEPFlux, HLLCFlux,
201 * and RoeFlux.
202 */
203 template <int dim>
204 struct RoePreamble
205 {
206 using TVec = Eigen::Vector<real, dim>;
207 // L/R mean-state thermodynamics
208 TVec veloLm, veloRm; ///< Left/right mean-state velocity vectors.
209 real asqrLm, asqrRm, pLm, pRm, HLm, HRm; ///< Speed-of-sound², pressure, and enthalpy for L/R mean states.
210 real vsqrLm, vsqrRm; ///< Squared velocity magnitudes for L/R mean states.
211 // Roe averages
212 TVec veloRoe; ///< Roe-averaged velocity vector.
213 real sqrtRhoLm, sqrtRhoRm; ///< √ρ for L/R states (used in Roe weighting).
214 real vsqrRoe, HRoe, asqrRoe, rhoRoe, aRoe; ///< Roe-averaged |v|², H, a², ρ, and speed of sound.
215 };
216
217 /**
218 * @brief Compute Roe-averaged quantities from mean-state L/R vectors.
219 *
220 * Extracts velocities, calls IdealGasThermal for both sides, then computes
221 * the Roe-averaged velocity, enthalpy, and speed of sound.
222 *
223 * @param ULm Left mean-state conservative vector
224 * @param URm Right mean-state conservative vector
225 * @param gamma Ratio of specific heats
226 * @param dumpInfo Callable invoked before assertion on invalid asqrRoe
227 * @return RoePreamble<dim> with all fields populated
228 */
229 template <int dim = 3, typename TULm, typename TURm, typename TFdumpInfo>
230 RoePreamble<dim> ComputeRoePreamble(const TULm &ULm, const TURm &URm,
231 real gamma, const TFdumpInfo &dumpInfo)
232 {
233 RoePreamble<dim> rp;
234
235 rp.veloLm = (ULm(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / ULm(0)).matrix();
236 rp.veloRm = (URm(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / URm(0)).matrix();
237 rp.vsqrLm = rp.veloLm.squaredNorm();
238 rp.vsqrRm = rp.veloRm.squaredNorm();
239 IdealGasThermal(ULm(dim + 1), ULm(0), rp.vsqrLm, gamma, rp.pLm, rp.asqrLm, rp.HLm);
240 IdealGasThermal(URm(dim + 1), URm(0), rp.vsqrRm, gamma, rp.pRm, rp.asqrRm, rp.HRm);
241
242 rp.sqrtRhoLm = std::sqrt(ULm(0));
243 rp.sqrtRhoRm = std::sqrt(URm(0));
244
245 rp.veloRoe = (rp.sqrtRhoLm * rp.veloLm + rp.sqrtRhoRm * rp.veloRm) / (rp.sqrtRhoLm + rp.sqrtRhoRm);
246 rp.vsqrRoe = rp.veloRoe.squaredNorm();
247 rp.HRoe = (rp.sqrtRhoLm * rp.HLm + rp.sqrtRhoRm * rp.HRm) / (rp.sqrtRhoLm + rp.sqrtRhoRm);
248 rp.asqrRoe = (gamma - 1) * (rp.HRoe - 0.5 * rp.vsqrRoe);
249 rp.rhoRoe = rp.sqrtRhoLm * rp.sqrtRhoRm;
250
251 if (!(rp.asqrRoe > 0))
252 {
253 dumpInfo();
254 }
255 DNDS_assert(rp.asqrRoe > 0);
256 rp.aRoe = std::sqrt(rp.asqrRoe);
257
258 return rp;
259 }
260
261 /**
262 * @brief Converts conservative variables to primitive variables for an
263 * ideal gas.
264 *
265 * Conservative layout: U = [ρ, ρu, ρv, (ρw,) E, ...]
266 * Primitive layout: prim = [ρ, u, v, (w,) p, ...]
267 *
268 * The primitive variable at index I4 = dim+1 stores **pressure** (Euler
269 * module convention), not temperature.
270 *
271 * @tparam dim Spatial dimension (2 or 3).
272 * @tparam TCons Conservative vector type (Eigen expression, length ≥ dim+2).
273 * @tparam TPrim Primitive vector type (Eigen expression, same size as TCons).
274 * @param U Conservative state vector (input).
275 * @param[out] prim Primitive state vector (output).
276 * @param gamma Ratio of specific heats (γ).
277 */
278 template <int dim = 3, class TCons, class TPrim>
279 inline void IdealGasThermalConservative2Primitive(
280 const TCons &U, TPrim &prim, real gamma)
281 {
282 prim = U / U(0);
283 real vSqr = (U(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) / U(0)).squaredNorm();
284 real rho = U(0);
285 real E = U(1 + dim);
286 real p = (gamma - 1) * (E - rho * 0.5 * vSqr);
287 prim(0) = rho;
288 prim(1 + dim) = p;
289 DNDS_assert(rho > 0);
290 }
291
292 /**
293 * @brief Converts primitive variables to conservative variables for an
294 * ideal gas.
295 *
296 * Inverse of IdealGasThermalConservative2Primitive().
297 *
298 * @tparam dim Spatial dimension (2 or 3).
299 * @tparam TCons Conservative vector type.
300 * @tparam TPrim Primitive vector type.
301 * @param prim Primitive state vector (input).
302 * @param[out] U Conservative state vector (output).
303 * @param gamma Ratio of specific heats (γ).
304 */
305 template <int dim = 3, class TCons, class TPrim>
306 inline void IdealGasThermalPrimitive2Conservative(
307 const TPrim &prim, TCons &U, real gamma)
308 {
309 U = prim * prim(0);
310 real vSqr = prim(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).squaredNorm();
311 real rho = prim(0);
312 real p = prim(dim + 1);
313 real E = p / (gamma - 1) + (rho * 0.5 * vSqr);
314 U(0) = rho;
315 U(dim + 1) = E;
316 DNDS_assert(rho > 0);
317 }
318
319 /**
320 * @brief Computes total (stagnation) pressure p0 and temperature T0 from
321 * a primitive state using isentropic relations.
322 *
323 * Uses p0 = p * (1 + (γ−1)/2 · M²)^(γ/(γ−1)) and T0 = T * (1 + (γ−1)/2 · M²).
324 *
325 * @tparam dim Spatial dimension (2 or 3).
326 * @tparam TPrim Primitive vector type.
327 * @param prim Primitive state [ρ, u, v, (w,) p, ...].
328 * @param gamma Ratio of specific heats (γ).
329 * @param rg Specific gas constant R_gas = Cp − Cv.
330 * @return std::tuple<real,real> (p0, T0).
331 */
332 template <int dim = 3, class TPrim>
333 std::tuple<real, real> IdealGasThermalPrimitiveGetP0T0(
334 const TPrim &prim, real gamma, real rg)
335 {
336 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
337 static const auto I4 = dim + 1;
338 real T = prim(I4) / (prim(0) * rg + verySmallReal);
339 real vsqr = prim(Seq123).squaredNorm();
340 real asqr = gamma * prim(I4) / prim(0);
341 real Msqr = vsqr / (asqr + verySmallReal);
342 real p0 = std::pow(1 + (gamma - 1) * 0.5 * Msqr, gamma / (gamma - 1)) * prim(I4);
343 real T0 = (1 + (gamma - 1) * 0.5 * Msqr) * T;
344 return std::make_tuple(p0, T0);
345 }
346
347 /**
348 * @brief Computes the inviscid (Euler) flux in the x-direction for a
349 * moving grid.
350 *
351 * The flux vector is written to the first (dim+2) entries of @p F:
352 * F = [ρ(u−ug), ρu(u−ug)+p, ρv(u−ug), (ρw(u−ug),) (E+p)u − E·ug]
353 * where ug = vg(0) is the grid velocity in the x-direction.
354 *
355 * @note Passive-scalar (RANS) entries beyond index dim+1 are **not** filled.
356 *
357 * @tparam dim Spatial dimension (2 or 3).
358 * @param U Conservative state vector.
359 * @param velo Fluid velocity vector.
360 * @param vg Grid velocity vector (for ALE / moving-mesh formulations).
361 * @param p Pressure.
362 * @param[out] F Flux vector (first dim+2 entries overwritten).
363 */
364 template <int dim = 3, typename TU, typename TF, class TVec, class TVecVG>
365 inline void GasInviscidFlux(const TU &U, const TVec &velo, const TVecVG &vg, real p, TF &F)
366 {
367 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = U(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) * (velo(0) - vg(0)); // note that additional flux are unattended!
368 F(1) += p;
369 F(dim + 1) += velo(0) * p;
370 // original form: F(dim + 1) += (velo(0) - vg(0)) * p + vg(0) * p;
371 }
372
373 /**
374 * @brief Computes the inviscid flux projected onto an arbitrary face normal
375 * @p n, accounting for grid motion @p vg.
376 *
377 * This is the main face-flux function used during spatial discretisation.
378 * The normal velocity is V_n = (velo − vg) · n.
379 *
380 * @note Passive-scalar entries beyond index dim+1 are **not** filled.
381 *
382 * @tparam dim Spatial dimension.
383 * @param U Conservative state vector.
384 * @param velo Fluid velocity vector.
385 * @param vg Grid velocity vector.
386 * @param n Outward face-normal vector (unit or area-weighted).
387 * @param p Pressure.
388 * @param[out] F Flux vector (first dim+2 entries overwritten).
389 */
390 template <int dim = 3, typename TU, typename TF, class TVec, class TVecN, class TVecVG>
391 inline void GasInviscidFlux_XY(const TU &U, const TVec &velo, const TVecVG &vg, const TVecN &n, real p, TF &F)
392 {
393 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = U(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) * (velo - vg).dot(n); // note that additional flux are unattended!
394 F(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) += p * n;
395 F(dim + 1) += velo.dot(n) * p;
396 // original form: F(dim + 1) += (velo(0) - vg(0)) * p + vg(0) * p;
397 }
398
399 /**
400 * @brief Batched x-direction inviscid flux for column-major state matrices.
401 *
402 * Each column of @p U, @p velo, @p vg represents one face in the batch;
403 * @p p is a row-vector of pressures. Output columns of @p F are filled
404 * independently.
405 *
406 * @tparam dim Spatial dimension.
407 * @param U Conservative states [(dim+2) × nBatch].
408 * @param velo Velocity vectors [dim × nBatch].
409 * @param vg Grid velocities [dim × nBatch].
410 * @param p Pressures [1 × nBatch].
411 * @param[out] F Flux matrix [(dim+2) × nBatch].
412 */
413 template <int dim = 3, typename TU, typename TF, class TVec, class TVecVG, class TP>
414 inline void GasInviscidFlux_Batch(const TU &U, const TVec &velo, const TVecVG &vg, TP &&p, TF &F)
415 {
416 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll) = U(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll).array().rowwise() * (velo(0, EigenAll) - vg(0, EigenAll)).array(); // note that additional flux are unattended!
417 F(1, EigenAll).array() += p.array();
418 F(dim + 1, EigenAll).array() += velo(0, EigenAll).array() * p.array();
419 // original form: F(dim + 1) += (velo(0) - vg(0)) * p + vg(0) * p;
420 }
421
422 /**
423 * @brief Batched face-normal inviscid flux for column-major state matrices.
424 *
425 * Each column represents one face in the batch. The flux at each column is
426 * projected onto the corresponding column of the normal matrix @p n.
427 *
428 * @tparam dim Spatial dimension.
429 * @param U Conservative states [(dim+2) × nBatch].
430 * @param velo Velocity vectors [dim × nBatch].
431 * @param vg Grid velocities [dim × nBatch].
432 * @param n Face normals [dim × nBatch].
433 * @param p Pressures [1 × nBatch].
434 * @param[out] F Flux matrix [(dim+2) × nBatch].
435 */
436 template <int dim = 3, typename TU, typename TF, class TVec, class TVecVG, class TVecN, class TP>
437 inline void GasInviscidFlux_XY_Batch(const TU &U, const TVec &velo, const TVecVG &vg, const TVecN &n, TP &&p, TF &F)
438 {
439 auto vn = (velo.array() * n.array()).colwise().sum();
440 auto vgn = (vg.array() * n.array()).colwise().sum();
441 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll) = U(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll).array().rowwise() * (vn - vgn).array(); // note that additional flux are unattended!
442 F(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll).array() += n.array().rowwise() * p.array();
443 F(dim + 1, EigenAll).array() += vn * p.array();
444 // original form: F(dim + 1) += (velo(0) - vg(0)) * p + vg(0) * p;
445 }
446
447 /**
448 * @brief Computes velocity and pressure increments from a conservative-state
449 * increment, used for the Lax-flux Jacobian computation.
450 *
451 * Given a state U and its increment dU, computes
452 * dVelo = d(momentum/ρ) and dp = (γ−1)(dE − ½(dMomentum·v + momentum·dv)).
453 *
454 * @tparam dim Spatial dimension.
455 * @param U Conservative state vector.
456 * @param dU Conservative state increment.
457 * @param velo Velocity vector (= U[1..dim] / U[0]).
458 * @param gamma Ratio of specific heats (γ).
459 * @param[out] dVelo Velocity increment (output).
460 * @param[out] dp Pressure increment (output).
461 */
462 template <int dim = 3, typename TU, class TVec>
463 inline void IdealGasUIncrement(const TU &U, const TU &dU, const TVec &velo, real gamma, TVec &dVelo, real &dp)
464 {
465 dVelo = (dU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) -
466 U(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) / U(0) * dU(0)) /
467 U(0);
468 dp = (gamma - 1) * (dU(dim + 1) -
469 0.5 * (dU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(velo) +
470 U(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(dVelo)));
471 } // For Lax-Flux jacobian
472
473 /**
474 * @brief Computes the increment of the facial inviscid flux from state,
475 * velocity, and pressure increments.
476 *
477 * Used in the implicit time-stepping Jacobian (e.g. LU-SGS).
478 * Writes the full flux vector @p F including passive-scalar entries
479 * beyond index dim+1 (valid for RANS and extra-equation models).
480 *
481 * \note now this function writes to the whole F, assuming SeqI52Last part is passive scalar
482 * this is valid for RANS and EX
483 *
484 * @tparam dim Spatial dimension.
485 * @param U Conservative state vector.
486 * @param dU Conservative state increment.
487 * @param unitNorm Face unit-normal vector.
488 * @param velo Velocity vector.
489 * @param dVelo Velocity increment (from IdealGasUIncrement()).
490 * @param vg Grid velocity vector.
491 * @param dp Pressure increment (from IdealGasUIncrement()).
492 * @param p Pressure.
493 * @param[out] F Incremental facial flux (all entries written).
494 */
495 template <int dim = 3, typename TU, typename TF, class TVec, class TVecVG>
496 inline void GasInviscidFluxFacialIncrement(const TU &U, const TU &dU,
497 const TVec &unitNorm,
498 const TVecVG &velo, const TVec &dVelo, const TVec &vg,
499 real dp, real p,
500 TF &F)
501 {
502 real vn = velo.dot(unitNorm);
503 real dvn = dVelo.dot(unitNorm);
504 F(0) = dU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(unitNorm);
505 F(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) =
506 dU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) * vn +
507 U(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) * dvn + unitNorm * dp;
508 F(dim + 1) = (dU(dim + 1) + dp) * vn + (U(dim + 1) + p) * dvn;
509 static const auto SeqI52Last = Eigen::seq(Eigen::fix<dim + 2>, EigenLast);
510 if constexpr (TU::RowsAtCompileTime == Eigen::Dynamic || TU::RowsAtCompileTime > dim + 2)
511 F(SeqI52Last) = dU(SeqI52Last) * vn + U(SeqI52Last) * dvn;
512 F -= dU * vg.dot(unitNorm);
513 }
514
515 /**
516 * @brief Convenience wrapper that computes the right eigenvector matrix
517 * directly from a conservative state vector U.
518 *
519 * Extracts velocity and thermodynamic quantities from U, then delegates to
520 * EulerGasRightEigenVector().
521 *
522 * @tparam dim Spatial dimension.
523 * @param U Conservative state vector.
524 * @param gamma Ratio of specific heats (γ).
525 * @return Eigen::Matrix<real, dim+2, dim+2> Right eigenvector matrix.
526 */
527 template <int dim = 3, typename TU>
528 inline auto IdealGas_EulerGasRightEigenVector(const TU &U, real gamma)
529 {
530 DNDS_assert(U(0) > 0);
531 Eigen::Matrix<real, dim + 2, dim + 2> ReV;
532 Eigen::Vector<real, dim> velo =
533 (U(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / U(0)).matrix();
534 real vsqr = velo.squaredNorm();
535 real asqr, p, H;
536 IdealGasThermal(U(dim + 1), U(0), vsqr, gamma, p, asqr, H);
537 DNDS_assert(asqr >= 0);
538 EulerGasRightEigenVector<dim>(velo, vsqr, H, std::sqrt(asqr), ReV);
539 return ReV;
540 }
541
542 /**
543 * @brief Convenience wrapper that computes the left eigenvector matrix
544 * directly from a conservative state vector U.
545 *
546 * Extracts velocity and thermodynamic quantities from U, then delegates to
547 * EulerGasLeftEigenVector().
548 *
549 * @tparam dim Spatial dimension.
550 * @param U Conservative state vector.
551 * @param gamma Ratio of specific heats (γ).
552 * @return Eigen::Matrix<real, dim+2, dim+2> Left eigenvector matrix.
553 */
554 template <int dim = 3, typename TU>
555 inline auto IdealGas_EulerGasLeftEigenVector(const TU &U, real gamma)
556 {
557 DNDS_assert(U(0) > 0);
558 Eigen::Matrix<real, dim + 2, dim + 2> LeV;
559 Eigen::Vector<real, dim> velo =
560 (U(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / U(0)).matrix();
561 real vsqr = velo.squaredNorm();
562 real asqr, p, H;
563 IdealGasThermal(U(dim + 1), U(0), vsqr, gamma, p, asqr, H);
564 DNDS_assert(asqr >= 0);
565 EulerGasLeftEigenVector<dim>(velo, vsqr, H, std::sqrt(asqr), gamma, LeV);
566 return LeV;
567 }
568
569 /**
570 * @brief HLLEP (HLL with Enhanced Pressure) approximate Riemann solver for
571 * an ideal gas on a moving grid.
572 *
573 * Uses Roe-decomposed wave strengths combined with HLL wave-speed estimates.
574 * When @p type = 0, computes the standard HLLEP flux; when @p type = 1,
575 * computes the HLLEP_V1 variant which uses modified eigenvalue scaling
576 * factors δ1, δ2, δ3 instead of the HLL blending formula.
577 *
578 * Left/right face states (@p UL, @p UR) are used for the actual flux and
579 * jump, while the mean states (@p ULm, @p URm) are used for Roe averaging
580 * (these may differ under MUSCL or WENO reconstruction).
581 *
582 * @tparam dim Spatial dimension (2 or 3).
583 * @tparam type 0 = HLLEP, 1 = HLLEP_V1.
584 * @param UL, UR Left/right conservative face states.
585 * @param ULm, URm Left/right mean-state conservative vectors for Roe averaging.
586 * @param vg Grid velocity vector.
587 * @param n Face-normal vector (unit or area-weighted).
588 * @param gamma Ratio of specific heats (γ).
589 * @param[out] F Numerical flux (first dim+2 entries written).
590 * @param dLambda H-correction transverse eigenvalue estimate from
591 * neighbouring faces.
592 * @param fixScale User-configurable scaling factor for entropy fix
593 * (from settings.rsFixScale).
594 * @param dumpInfo Callable invoked before assertion on invalid state
595 * (negative density, NaN, etc.).
596 * @param[out] lam0 Absolute eigenvalue |V_n − a| after entropy fixing.
597 * @param[out] lam123 Absolute eigenvalue |V_n| after entropy fixing.
598 * @param[out] lam4 Absolute eigenvalue |V_n + a| after entropy fixing.
599 */
600 // #define DNDS_GAS_HLLEP_USE_V1
601 template <int dim = 3, int type = 0,
602 typename TUL, typename TUR,
603 typename TULm, typename TURm, typename TVecVG, typename TVecN,
604 typename TF, typename TFdumpInfo>
605 void HLLEPFlux_IdealGas(const TUL &UL, const TUR &UR, const TULm &ULm, const TURm &URm,
606 const TVecVG &vg, const TVecN &n,
607 real gamma, TF &F, real dLambda, real fixScale,
608 const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
609 {
610 using TVec = Eigen::Vector<real, dim>;
611
612 if (!(UL(0) > 0 && UR(0) > 0))
613 {
614 dumpInfo();
615 }
616 DNDS_assert(UL(0) > 0 && UR(0) > 0);
617 TVec veloL = (UL(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / UL(0)).matrix();
618 TVec veloR = (UR(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / UR(0)).matrix();
619
620 real asqrL, asqrR, pL, pR, HL, HR;
621 real vsqrL = veloL.squaredNorm();
622 real vsqrR = veloR.squaredNorm();
623 IdealGasThermal(UL(dim + 1), UL(0), vsqrL, gamma, pL, asqrL, HL);
624 IdealGasThermal(UR(dim + 1), UR(0), vsqrR, gamma, pR, asqrR, HR);
625
626 auto rp = ComputeRoePreamble<dim>(ULm, URm, gamma, dumpInfo);
627
628 Eigen::Vector<real, dim + 2> FL, FR;
629 GasInviscidFlux_XY<dim>(UL, veloL, vg, n, pL, FL);
630 GasInviscidFlux_XY<dim>(UR, veloR, vg, n, pR, FR);
631
632 real veloRoeN = rp.veloRoe.dot(n);
633 real vgN = vg.dot(n);
634 real veloRoe0 = veloRoeN - vgN;
635 lam0 = std::abs(veloRoe0 - rp.aRoe);
636 lam123 = std::abs(veloRoe0);
637 lam4 = std::abs(veloRoe0 + rp.aRoe);
638 real veloLm0 = (rp.veloLm - vg).dot(n);
639 real veloRm0 = (rp.veloRm - vg).dot(n);
640
641 Eigen::Vector<real, dim + 2> incU =
642 UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) -
643 UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)); //! not using m, for this is accuracy-limited!
644 real incU123N = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(n);
645
646 TVec alpha23V = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) - incU(0) * rp.veloRoe;
647 TVec alpha23VT = alpha23V - n * alpha23V.dot(n);
648 real incU4b = incU(dim + 1) - alpha23VT.dot(rp.veloRoe);
649 real alpha1 = (gamma - 1) / rp.asqrRoe *
650 (incU(0) * (rp.HRoe - veloRoeN * veloRoeN) +
651 veloRoeN * incU123N - incU4b);
652 real alpha0 = (incU(0) * (veloRoeN + rp.aRoe) - incU123N - rp.aRoe * alpha1) / (2 * rp.aRoe);
653 real alpha4 = incU(0) - (alpha0 + alpha1);
654
655 // std::cout << alpha.transpose() << std::endl;
656 // std::cout << std::endl;
657
658 real SL = std::min(veloRoe0 - rp.aRoe, veloLm0 - std::sqrt(asqrL));
659 real SR = std::max(veloRoe0 + rp.aRoe, veloRm0 + std::sqrt(asqrR));
660 real UU = std::abs(veloRoe0);
661 real dfix = rp.aRoe / (rp.aRoe + std::max(UU, dLambda * fixScale));
662 real SP = std::max(SR, 0.0);
663 real SM = std::min(SL, 0.0);
664
665 if constexpr (type != 1) // ? HLLEP
666 {
667 Eigen::Vector<real, dim + 2> ReV1;
668 {
669 ReV1(0) = 1;
670 ReV1(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) = rp.veloRoe;
671 ReV1(dim + 1) = rp.vsqrRoe * 0.5;
672 }
673 real div = SP - SM;
674 div += signP(div) * verySmallReal;
675 // F = (SP * FL - SM * FR) / div + (SP * SM / div) * (UR - UL - dfix * ReVRoe(EigenAll, {1, 2, 3}) * alpha({1, 2, 3}));
676 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) =
677 (SP * FL - SM * FR) / div +
678 (SP * SM / div) *
679 (UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) -
680 UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) -
681 dfix * ReV1 * alpha1);
682 }
683 else
684 {
685 // ? HLLEP_V1
686 // real delta1 = aRoe / (UU + aRoe + verySmallReal);
687 // real delta1 = 0.5;
688 real delta1 = dfix;
689 real delta2 = 0.0;
690 real delta3 = 0.0;
691
692 lam123 = ((SP + SM) * (veloRoe0)-2.0 * (1.0 - delta1) * (SP * SM)) / (SP - SM);
693 lam4 = ((SP + SM) * (veloRoe0 + rp.aRoe) - 2.0 * (1.0 - delta2) * (SP * SM)) / (SP - SM);
694 lam0 = ((SP + SM) * (veloRoe0 - rp.aRoe) - 2.0 * (1.0 - delta3) * (SP * SM)) / (SP - SM);
695 lam123 = std::abs(lam123);
696 lam0 = std::abs(lam0);
697 lam4 = std::abs(lam4);
698
699 alpha0 *= lam0;
700 alpha1 *= lam123;
701 alpha23VT *= lam123;
702 alpha4 *= lam4; // here becomes alpha_i * lam_i
703
704 Eigen::Vector<real, dim + 2> incF;
705 incF(0) = alpha0 + alpha1 + alpha4;
706 incF(dim + 1) = (rp.HRoe - veloRoeN * rp.aRoe) * alpha0 + 0.5 * rp.vsqrRoe * alpha1 +
707 (rp.HRoe + veloRoeN * rp.aRoe) * alpha4 + alpha23VT.dot(rp.veloRoe);
708 incF(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) =
709 (rp.veloRoe - rp.aRoe * n) * alpha0 + (rp.veloRoe + rp.aRoe * n) * alpha4 +
710 rp.veloRoe * alpha1 + alpha23VT;
711
712 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = (FL + FR) * 0.5 - 0.5 * incF;
713 }
714 }
715
716 /**
717 * @brief HLLC (Harten-Lax-van Leer-Contact) approximate Riemann solver for
718 * an ideal gas on a moving grid.
719 *
720 * @warning This solver has known accuracy issues (observed in the IV test
721 * problem). Prefer Roe or HLLEP for production runs.
722 *
723 * @tparam dim Spatial dimension (2 or 3).
724 * @param UL, UR Left/right conservative face states.
725 * @param ULm, URm Left/right mean states for Roe averaging.
726 * @param vg Grid velocity vector.
727 * @param n Face-normal vector.
728 * @param gamma Ratio of specific heats (γ).
729 * @param[out] F Numerical flux (first dim+2 entries written).
730 * @param dLambda H-correction transverse eigenvalue estimate.
731 * @param fixScale Entropy-fix scaling factor.
732 * @param dumpInfo Debug dump callable for assertion failures.
733 * @param[out] lam0 |V_n − a| (Roe-averaged, before HLLC wave selection).
734 * @param[out] lam123 |V_n| (Roe-averaged).
735 * @param[out] lam4 |V_n + a| (Roe-averaged).
736 */
737 template <int dim = 3, typename TUL, typename TUR, typename TULm, typename TURm,
738 typename TVecVG, typename TVecN,
739 typename TF, typename TFdumpInfo>
740 void HLLCFlux_IdealGas_HartenYee(const TUL &UL, const TUR &UR, const TULm &ULm, const TURm &URm,
741 const TVecVG &vg, const TVecN &n,
742 real gamma, TF &F, real dLambda, real fixScale,
743 const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
744 {
745 //! warning: has accuracy issue (see IV test)
746 using TVec = Eigen::Vector<real, dim>;
747
748 if (!(UL(0) > 0 && UR(0) > 0))
749 {
750 dumpInfo();
751 }
752 DNDS_assert(UL(0) > 0 && UR(0) > 0);
753 TVec veloL = (UL(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / UL(0)).matrix();
754 TVec veloR = (UR(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / UR(0)).matrix();
755
756 real asqrL, asqrR, pL, pR, HL, HR;
757 real vsqrL = veloL.squaredNorm();
758 real vsqrR = veloR.squaredNorm();
759 IdealGasThermal(UL(dim + 1), UL(0), vsqrL, gamma, pL, asqrL, HL);
760 IdealGasThermal(UR(dim + 1), UR(0), vsqrR, gamma, pR, asqrR, HR);
761
762 auto rp = ComputeRoePreamble<dim>(ULm, URm, gamma, dumpInfo);
763
764 Eigen::Vector<real, dim + 2> FL, FR;
765 GasInviscidFlux_XY<dim>(UL, veloL, vg, n, pL, FL);
766 GasInviscidFlux_XY<dim>(UR, veloR, vg, n, pR, FR);
767
768 real veloRoeN = rp.veloRoe.dot(n);
769 real vgN = vg.dot(n);
770 real veloRoe0 = veloRoeN - vgN;
771 lam0 = std::abs(veloRoe0 - rp.aRoe);
772 lam123 = std::abs(veloRoe0);
773 lam4 = std::abs(veloRoe0 + rp.aRoe);
774 real veloLm0 = (rp.veloLm - vg).dot(n);
775 real veloRm0 = (rp.veloRm - vg).dot(n);
776
777 auto HLLCq = [&](real p, real pS)
778 {
779 real q = std::sqrt(1 + (gamma + 1) / 2 / gamma * (pS / p - 1));
780 if (pS <= p)
781 q = 1;
782 return q;
783 };
784 real pS = 0.5 * (rp.pLm + rp.pRm) - 0.5 * (veloLm0 - veloRm0) * rp.rhoRoe * rp.aRoe;
785 pS = std::max(0.0, pS);
786 real SL = veloLm0 - std::sqrt(rp.asqrLm) * HLLCq(rp.pLm, pS);
787 real SR = veloRm0 + std::sqrt(rp.asqrRm) * HLLCq(rp.pRm, pS);
788
789 dLambda += verySmallReal;
790 dLambda *= 2.0;
791
792 // * E-Fix
793 // SL += sign(SL) * std::exp(-std::abs(SL) / dLambda) * dLambda;
794 // SR += sign(SR) * std::exp(-std::abs(SR) / dLambda) * dLambda;
795
796 if (0 <= SL)
797 {
798 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = FL;
799 return;
800 }
801 if (SR <= 0)
802 {
803 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = FR;
804 return;
805 }
806 real SS = 0;
807 real div = (UL(0) * (SL - veloLm0) - UR(0) * (SR - veloRm0));
808 if (std::abs(div) > verySmallReal)
809 SS = (pR - pL +
810 (UL(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(n) - vgN * UL(0)) * (SL - veloLm0) -
811 (UR(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(n) - vgN * UR(0)) * (SR - veloRm0)) /
812 div;
813 //! is this right for moving mesh?
814 Eigen::Vector<real, dim + 2> DS;
815 DS.setZero();
816 DS(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) = n;
817 DS(dim + 1) = SS;
818 // SS += sign(SS) * std::exp(-std::abs(SS) / dLambda) * dLambda;
819 if (SS >= 0)
820 {
821 real div = SL - SS;
822 if (std::abs(div) < verySmallReal)
823 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = FL;
824 else
825 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) =
826 ((UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) * SL - FL) * SS +
827 DS * ((pL + UL(0) * (SL - veloLm0) * (SS - veloLm0)) * SL)) /
828 div;
829 }
830 else
831 {
832 real div = SR - SS;
833 if (std::abs(div) < verySmallReal)
834 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = FR;
835 else
836 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) =
837 ((UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) * SR - FR) * SS +
838 DS * ((pR + UR(0) * (SR - veloRm0) * (SS - veloRm0)) * SR)) /
839 div;
840 }
841 }
842
843 /**
844 * @brief Template-dispatched entropy-fix for Roe-type Riemann solvers.
845 *
846 * Modifies the absolute Roe eigenvalues (lam0, lam123, lam4) in-place to
847 * prevent the sonic glitch and control numerical dissipation. The fix
848 * strategy is selected at compile time via @p eigScheme:
849 *
850 * | eigScheme | Strategy |
851 * |-----------|-------------------------------------------------------------|
852 * | 0 | Harten-Yee (H2) — threshold = max(0, λ−λ_L, λ_R−λ) · s. |
853 * | 1 | cLLF — componentwise local Lax-Friedrichs. |
854 * | 3 | LD Roe — low-dissipation (Fleischmann et al. 2020). |
855 * | 4 | ID Roe — intermediate-dissipation (Fleischmann et al. 2020). |
856 * | 5 | LD cLLF — low-dissipation componentwise LLF. |
857 * | 6 | H-correction only (floor by dLambda · scaleHFix). |
858 * | 7 | Harten-Yee only, no H-correction. |
859 * | 8 | H-correction + Harten-Yee combined. |
860 *
861 * @tparam eigScheme Compile-time entropy-fix scheme selector (0–8).
862 * @param aL, aR Left/right mean-state speed of sound.
863 * @param aAve Roe-averaged speed of sound.
864 * @param uL, uR Left/right normal velocities relative to grid.
865 * @param uAve Roe-averaged normal velocity relative to grid.
866 * @param VL, VR Left/right total velocity magnitudes relative to grid.
867 * @param VAve Roe-averaged total velocity magnitude relative to grid.
868 * @param dLambda H-correction transverse eigenvalue estimate.
869 * @param fixScale User scaling factor applied to the base entropy-fix
870 * constants (kScaleHartenYee, kScaleLD, kScaleHFix).
871 * @param incFScale Incremental flux dissipation scale (settings.rsIncFScale),
872 * controls overall dissipation level in the contact/shear
873 * wave.
874 * @param[in,out] lam0 |V_n − a| — modified in place.
875 * @param[in,out] lam123 |V_n| — modified in place.
876 * @param[in,out] lam4 |V_n + a| — modified in place.
877 */
878 template <int eigScheme>
879 void Roe_EntropyFixer(const real aL, const real aR, const real aAve,
880 const real uL, const real uR, const real uAve,
881 const real VL, const real VR, const real VAve, // V is magnitude of velo
882 real dLambda, real fixScale, real incFScale,
883 real &lam0, real &lam123, real &lam4)
884 {
885 const real scaleHartenYee = kScaleHartenYee * fixScale;
886 const real scaleLD = kScaleLD * fixScale;
887 const real scaleHFix = kScaleHFix * fixScale;
888 static const real incFScaleA = 1.0;
889 if constexpr (eigScheme == 0)
890 {
891 //*H2
892 auto Flim = [=](real v, real lam, real lamL, real lamR)
893 {
894 const real thresH = std::max<real>({0., lam - lamL, lamR - lam}) * scaleHartenYee;
895 if (v < thresH)
896 return (sqr(v) / thresH + thresH) * 0.5;
897 else
898 return v;
899 };
900
901 lam0 = Flim(lam0 * incFScaleA, uAve - aAve, uL - aL, uR - aR);
902 lam123 = Flim(lam123 * incFScale, uAve, uL, uR);
903 lam4 = Flim(lam4 * incFScaleA, uAve + aAve, uL + aL, uR + aR);
904 //*H2
905
906 // lam0 = std::max(lam0, dLambda * scaleHFix);
907 // lam123 = std::max(lam123, dLambda * scaleHFix);
908 // lam4 = std::max(lam4, dLambda * scaleHFix);
909 }
910 else if constexpr (eigScheme == 1)
911 {
912 //* cLLF
913 const real aLm = aL * incFScaleA;
914 const real aRm = aR * incFScaleA;
915 const real veloLm0 = uL;
916 const real veloRm0 = uR;
917 const real uLm = signTol(veloLm0, aLm * smallReal) * std::max(std::abs(veloLm0) * incFScale, aLm * scaleLD);
918 const real uRm = signTol(veloRm0, aRm * smallReal) * std::max(std::abs(veloRm0) * incFScale, aRm * scaleLD);
919 lam0 = std::max(std::abs(uLm - aLm), std::abs(uRm - aRm));
920 lam123 = std::max(std::abs(uLm), std::abs(uRm));
921 lam4 = std::max(std::abs(uLm + aLm), std::abs(uRm + aRm));
922 }
923 else if constexpr (eigScheme == 3)
924 {
925 //*LD, Roe_M
926 /**
927 * Nico Fleischmann, Stefan Adami, Xiangyu Y. Hu, Nikolaus A. Adams, A low dissipation method to cure the grid-aligned shock instability, 2020
928 */
929 const real LDthreshold = std::abs(uAve) / scaleLD;
930 const real aRoeStar = std::min(LDthreshold, aAve * incFScaleA);
931 lam0 = std::abs(uAve * incFScale - aRoeStar);
932 lam4 = std::abs(uAve * incFScale + aRoeStar);
933 }
934 else if constexpr (eigScheme == 4)
935 {
936 //*ID, Roe_M
937 /**
938 * Nico Fleischmann, Stefan Adami, Xiangyu Y. Hu, Nikolaus A. Adams, A low dissipation method to cure the grid-aligned shock instability, 2020
939 */
940#ifdef USE_SIGN_MINUS_AT_ROE_M4_FLUX
941 const real uStar = signM(uAve) * std::max(aAve * scaleLD, std::abs(uAve));
942#else
943 const real uStar = signTol(uAve, aAve * smallReal) * std::max(aAve * scaleLD, std::abs(uAve) * incFScale); //! why signM here?
944#endif
945 lam0 = std::abs(uStar - aAve * incFScaleA);
946 lam123 = std::abs(uStar);
947 lam4 = std::abs(uStar + aAve * incFScaleA);
948 }
949 else if constexpr (eigScheme == 5)
950 {
951 //*LD, cLLF_M
952 /**
953 * Nico Fleischmann, Stefan Adami, Xiangyu Y. Hu, Nikolaus A. Adams, A low dissipation method to cure the grid-aligned shock instability, 2020
954 */
955 const real veloLm0 = uL * incFScale;
956 const real veloRm0 = uR * incFScale;
957 const real aLm = std::min(aL * incFScaleA, std::abs(veloLm0) / scaleLD);
958 const real aRm = std::min(aR * incFScaleA, std::abs(veloRm0) / scaleLD);
959 lam0 = std::max(std::abs(veloLm0 - aLm), std::abs(veloRm0 - aRm));
960 lam123 = std::max(std::abs(veloLm0), std::abs(veloRm0));
961 lam4 = std::max(std::abs(veloLm0 + aLm), std::abs(veloRm0 + aRm));
962 //*LD, cLLF_M
963 }
964 else if constexpr (eigScheme == 6) // simply H-correction
965 {
966 lam0 = std::max(std::abs(uAve * incFScale - aAve * incFScaleA), dLambda * scaleHFix);
967 lam4 = std::max(std::abs(uAve * incFScale + aAve * incFScaleA), dLambda * scaleHFix);
968 lam123 = std::max(std::abs(uAve * incFScale), dLambda * scaleHFix);
969 }
970 else if constexpr (eigScheme == 7) // simply Harten-Yee-fix
971 {
972 lam0 *= std::abs(uAve * incFScale - aAve * incFScaleA);
973 lam4 *= std::abs(uAve * incFScale + aAve * incFScaleA);
974 lam123 *= std::abs(uAve * incFScale);
975 //*HY
976 const real thresholdHartenYee = scaleHartenYee * (VAve + aAve);
977 const real thresholdHartenYeeS = sqr(thresholdHartenYee);
978 if (lam0 < thresholdHartenYee)
979 lam0 = (sqr(lam0) + thresholdHartenYeeS) / (2 * thresholdHartenYee);
980 if (lam4 < thresholdHartenYee)
981 lam4 = (sqr(lam4) + thresholdHartenYeeS) / (2 * thresholdHartenYee);
982 if (lam123 < thresholdHartenYee)
983 lam123 = (sqr(lam123) + thresholdHartenYeeS) / (2 * thresholdHartenYee);
984 //*HY
985 }
986 else if constexpr (eigScheme == 8) // H-cor + Harten-Yee-fix
987 {
988 lam0 = std::max(std::abs(uAve * incFScale - aAve * incFScaleA), dLambda * scaleHFix);
989 lam4 = std::max(std::abs(uAve * incFScale + aAve * incFScaleA), dLambda * scaleHFix);
990 lam123 = std::max(std::abs(uAve * incFScale), dLambda * scaleHFix);
991
992 //*HY
993 const real thresholdHartenYee = scaleHartenYee * (VAve + aAve);
994 const real thresholdHartenYeeS = sqr(thresholdHartenYee);
995 if (lam0 < thresholdHartenYee)
996 lam0 = (sqr(lam0) + thresholdHartenYeeS) / (2 * thresholdHartenYee);
997 if (lam4 < thresholdHartenYee)
998 lam4 = (sqr(lam4) + thresholdHartenYeeS) / (2 * thresholdHartenYee);
999 if (lam123 < thresholdHartenYee)
1000 lam123 = (sqr(lam123) + thresholdHartenYeeS) / (2 * thresholdHartenYee);
1001 //*HY
1002 }
1003 else
1004 {
1005 DNDS_assert(false);
1006 }
1007 }
1008
1009 /**
1010 * @brief Core Roe approximate Riemann solver with a selectable entropy-fix
1011 * scheme for an ideal gas on a moving grid.
1012 *
1013 * Computes the numerical inter-cell flux as
1014 * F = ½(F_L + F_R) − ½ |A_Roe| (U_R − U_L)
1015 * where |A_Roe| is the Roe matrix with eigenvalues modified by the chosen
1016 * entropy fix (see Roe_EntropyFixer()).
1017 *
1018 * When @p eigScheme = 2 (Lax-Friedrichs), the function takes an early-exit
1019 * path returning the vanilla Lax-Friedrichs flux without Roe decomposition.
1020 *
1021 * @tparam dim Spatial dimension (2 or 3).
1022 * @tparam eigScheme Compile-time entropy-fix selector (0–8); default 0
1023 * (Harten-Yee).
1024 * @param UL, UR Left/right conservative face states.
1025 * @param ULm, URm Left/right mean states for Roe averaging.
1026 * @param vg Grid velocity vector.
1027 * @param n Face-normal vector.
1028 * @param gamma Ratio of specific heats (γ).
1029 * @param[out] F Numerical flux (first dim+2 entries written).
1030 * @param dLambda H-correction transverse eigenvalue estimate.
1031 * @param fixScale Entropy-fix scaling factor (settings.rsFixScale).
1032 * @param incFScale Incremental flux dissipation scale (settings.rsIncFScale).
1033 * @param dumpInfo Debug dump callable for assertion failures.
1034 * @param[out] lam0 |V_n − a| after entropy fixing.
1035 * @param[out] lam123 |V_n| after entropy fixing.
1036 * @param[out] lam4 |V_n + a| after entropy fixing.
1037 */
1038 template <int dim = 3, int eigScheme = 0,
1039 typename TUL, typename TUR,
1040 typename TULm, typename TURm,
1041 typename TVecVG, typename TVecN,
1042 typename TF, typename TFdumpInfo>
1043 void RoeFlux_IdealGas_HartenYee(const TUL &UL, const TUR &UR,
1044 const TULm &ULm, const TURm &URm,
1045 const TVecVG &vg, const TVecN &n,
1046 real gamma, TF &F, real dLambda, real fixScale,
1047 real incFScale,
1048 const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
1049 {
1050 using TVec = Eigen::Vector<real, dim>;
1051
1052 if (!(UL(0) > 0 && UR(0) > 0))
1053 {
1054 dumpInfo();
1055 }
1056 DNDS_assert(UL(0) > 0 && UR(0) > 0);
1057 TVec veloL = (UL(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / UL(0)).matrix();
1058 TVec veloR = (UR(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / UR(0)).matrix();
1059
1060 real asqrL, asqrR, pL, pR, HL, HR;
1061 real vsqrL = veloL.squaredNorm();
1062 real vsqrR = veloR.squaredNorm();
1063 IdealGasThermal(UL(dim + 1), UL(0), vsqrL, gamma, pL, asqrL, HL);
1064 IdealGasThermal(UR(dim + 1), UR(0), vsqrR, gamma, pR, asqrR, HR);
1065
1066 auto rp = ComputeRoePreamble<dim>(ULm, URm, gamma, dumpInfo);
1067
1068 Eigen::Vector<real, dim + 2> FL, FR;
1069 GasInviscidFlux_XY<dim>(UL, veloL, vg, n, pL, FL);
1070 GasInviscidFlux_XY<dim>(UR, veloR, vg, n, pR, FR);
1071
1072 real veloRoeN = rp.veloRoe.dot(n);
1073 real vgN = vg.dot(n);
1074 real veloRoe0 = veloRoeN - vgN;
1075 lam0 = std::abs(veloRoe0 - rp.aRoe);
1076 lam123 = std::abs(veloRoe0);
1077 lam4 = std::abs(veloRoe0 + rp.aRoe);
1078 real veloLm0 = (rp.veloLm - vg).dot(n);
1079 real veloRm0 = (rp.veloRm - vg).dot(n);
1080
1081 if constexpr (eigScheme == 2)
1082 {
1083 // *vanilla Lax
1084 // lam0 = lam123 = lam4 = std::max({lam0, lam123, lam4});
1085 lam0 = lam123 = lam4 = std::max(std::abs(veloLm0) + std::sqrt(rp.asqrLm), std::abs(veloRm0) + std::sqrt(rp.asqrRm));
1086 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) =
1087 (FL + FR) * 0.5 -
1088 0.5 * lam0 * (UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) - UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)));
1089 return; //* early exit
1090 }
1091 else
1092 Roe_EntropyFixer<eigScheme>(std::sqrt(rp.asqrLm), std::sqrt(rp.asqrRm), rp.aRoe,
1093 veloLm0, veloRm0, veloRoe0,
1094 (rp.veloLm - vg).norm(), (rp.veloRm - vg).norm(), (rp.veloRoe - vg).norm(),
1095 dLambda, fixScale, incFScale,
1096 lam0, lam123, lam4);
1097
1098 Eigen::Vector<real, dim + 2> incU =
1099 UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) -
1100 UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)); //! not using m, for this is accuracy-limited!
1101 real incU123N = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(n);
1102
1103 TVec alpha23V = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) - incU(0) * rp.veloRoe;
1104 TVec alpha23VT = alpha23V - n * alpha23V.dot(n);
1105 real incU4b = incU(dim + 1) - alpha23VT.dot(rp.veloRoe);
1106 real alpha1 = (gamma - 1) / rp.asqrRoe *
1107 (incU(0) * (rp.HRoe - veloRoeN * veloRoeN) +
1108 veloRoeN * incU123N - incU4b);
1109 real alpha0 = (incU(0) * (veloRoeN + rp.aRoe) - incU123N - rp.aRoe * alpha1) / (2 * rp.aRoe);
1110 real alpha4 = incU(0) - (alpha0 + alpha1);
1111
1112 alpha0 *= lam0;
1113 alpha1 *= lam123;
1114 alpha23VT *= lam123;
1115 alpha4 *= lam4; // here becomes alpha_i * lam_i
1116
1117 Eigen::Vector<real, dim + 2> incF;
1118 incF(0) = alpha0 + alpha1 + alpha4;
1119 incF(dim + 1) = (rp.HRoe - veloRoeN * rp.aRoe) * alpha0 + 0.5 * rp.vsqrRoe * alpha1 +
1120 (rp.HRoe + veloRoeN * rp.aRoe) * alpha4 + alpha23VT.dot(rp.veloRoe);
1121 incF(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) =
1122 (rp.veloRoe - rp.aRoe * n) * alpha0 + (rp.veloRoe + rp.aRoe * n) * alpha4 +
1123 rp.veloRoe * alpha1 + alpha23VT;
1124
1125 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>)) = (FL + FR) * 0.5 - 0.5 * incF;
1126 }
1127
1128 /**
1129 * @brief Computes the Roe-averaged state vector including passive scalars
1130 * (e.g. RANS turbulence variables).
1131 *
1132 * The Roe average for the Euler components (density, momentum, energy) is
1133 * the standard √ρ-weighted mean. Passive scalars beyond index dim+1 are
1134 * also averaged using the same √ρ weighting.
1135 *
1136 * @tparam dim Spatial dimension (2 or 3).
1137 * @param UL, UR Left/right conservative state vectors (full length
1138 * including any passive scalars).
1139 * @param gamma Ratio of specific heats (γ).
1140 * @param[out] veloRoe Roe-averaged velocity vector.
1141 * @param[out] vsqrRoe Roe-averaged squared velocity magnitude.
1142 * @param[out] aRoe Roe-averaged speed of sound.
1143 * @param[out] asqrRoe Roe-averaged speed of sound squared.
1144 * @param[out] HRoe Roe-averaged total specific enthalpy.
1145 * @param[out] UOut Roe-averaged state vector (same layout as UL/UR,
1146 * including passive scalars).
1147 */
1148 template <int dim = 3, typename TUL, typename TUR, typename TVecV, typename TUOut>
1149 void GetRoeAverage(const TUL &UL, const TUR &UR, real gamma,
1150 TVecV &veloRoe, real &vsqrRoe, real &aRoe, real &asqrRoe, real &HRoe, TUOut &UOut)
1151 {
1152 using TVec = Eigen::Vector<real, dim>;
1153 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
1154 static const auto I4 = dim + 1;
1155 TVec veloLm = (UL(Seq123).array() / UL(0)).matrix();
1156 TVec veloRm = (UR(Seq123).array() / UR(0)).matrix();
1157 real asqrLm, asqrRm, pLm, pRm, HLm, HRm;
1158 real vsqrLm = veloLm.squaredNorm();
1159 real vsqrRm = veloRm.squaredNorm();
1160 IdealGasThermal(UL(dim + 1), UL(0), vsqrLm, gamma, pLm, asqrLm, HLm);
1161 IdealGasThermal(UR(dim + 1), UR(0), vsqrRm, gamma, pRm, asqrRm, HRm);
1162 DNDS_assert(UL(0) >= 0 && UR(0) >= 0);
1163 real sqrtRhoLm = std::sqrt(UL(0));
1164 real sqrtRhoRm = std::sqrt(UR(0));
1165
1166 veloRoe = (sqrtRhoLm * veloLm + sqrtRhoRm * veloRm) / (sqrtRhoLm + sqrtRhoRm);
1167 vsqrRoe = veloRoe.squaredNorm();
1168 HRoe = (sqrtRhoLm * HLm + sqrtRhoRm * HRm) / (sqrtRhoLm + sqrtRhoRm);
1169 asqrRoe = (gamma - 1) * (HRoe - 0.5 * vsqrRoe);
1170 DNDS_assert(asqrRoe >= 0);
1171 aRoe = std::sqrt(asqrRoe);
1172
1173 UOut(0) = sqrtRhoLm * sqrtRhoRm;
1174 UOut(Seq123) = veloRoe * UOut(0);
1175 real pRoe = asqrRoe * UOut(0) / gamma;
1176 UOut(I4) = pRoe / (gamma - 1) + 0.5 * vsqrRoe * UOut(0);
1177 UOut(Eigen::seq(I4 + 1, EigenLast)) =
1178 UOut(0) / (sqrtRhoLm + sqrtRhoRm) *
1179 (UL(Eigen::seq(I4 + 1, EigenLast)) / sqrtRhoLm +
1180 UR(Eigen::seq(I4 + 1, EigenLast)) / sqrtRhoRm);
1181 }
1182
1183 /**
1184 * @brief Computes the dissipation part of the Roe flux |A|·dU, given
1185 * pre-computed Roe averages and entropy-fixed eigenvalues.
1186 *
1187 * This function **accumulates** into @p incF (does not zero it first).
1188 * It handles the Euler components (indices 0..dim+1) via the standard Roe
1189 * eigenvector decomposition, and passive-scalar components (indices
1190 * dim+2..end) with the contact/shear eigenvalue lam123.
1191 *
1192 * Used by the LUSGS implicit solver where the dissipation term is needed
1193 * separately from the central flux.
1194 *
1195 * @tparam dim Spatial dimension (2 or 3).
1196 * @param incU State jump U_R − U_L (full length, including scalars).
1197 * @param n Face-normal vector.
1198 * @param veloRoe Roe-averaged velocity.
1199 * @param vsqrRoe Roe-averaged |v|².
1200 * @param aRoe Roe-averaged speed of sound.
1201 * @param asqrRoe Roe-averaged a².
1202 * @param HRoe Roe-averaged total enthalpy.
1203 * @param lam0 Entropy-fixed eigenvalue |V_n − a|.
1204 * @param lam123 Entropy-fixed eigenvalue |V_n|.
1205 * @param lam4 Entropy-fixed eigenvalue |V_n + a|.
1206 * @param gamma Ratio of specific heats (γ).
1207 * @param[in,out] incF Dissipation flux; accumulated (not zeroed).
1208 */
1209 template <int dim = 3, typename TDU, typename TDF, typename TVecV, typename TVecN>
1210 void RoeFluxIncFDiff(const TDU &incU, const TVecN &n, const TVecV &veloRoe,
1211 real vsqrRoe, real aRoe, real asqrRoe, real HRoe,
1212 real lam0, real lam123, real lam4, real gamma,
1213 TDF &incF)
1214 {
1215 static const auto SeqI52Last = Eigen::seq(Eigen::fix<dim + 2>, EigenLast);
1216 using TVec = Eigen::Vector<real, dim>;
1217 real veloRoeN = veloRoe.dot(n);
1218
1219 real incU123N = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).dot(n);
1220
1221 TVec alpha23V = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) - incU(0) * veloRoe;
1222 TVec alpha23VT = alpha23V - n * alpha23V.dot(n);
1223 real incU4b = incU(dim + 1) - alpha23VT.dot(veloRoe);
1224 real alpha1 = (gamma - 1) / asqrRoe *
1225 (incU(0) * (HRoe - veloRoeN * veloRoeN) +
1226 veloRoeN * incU123N - incU4b);
1227 real alpha0 = (incU(0) * (veloRoeN + aRoe) - incU123N - aRoe * alpha1) / (2 * aRoe);
1228 real alpha4 = incU(0) - (alpha0 + alpha1);
1229
1230 alpha0 *= lam0;
1231 alpha1 *= lam123;
1232 alpha23VT *= lam123;
1233 alpha4 *= lam4; // here becomes alpha_i * lam_i
1234
1235 incF(0) += alpha0 + alpha1 + alpha4;
1236 incF(dim + 1) += (HRoe - veloRoeN * aRoe) * alpha0 + 0.5 * vsqrRoe * alpha1 +
1237 (HRoe + veloRoeN * aRoe) * alpha4 + alpha23VT.dot(veloRoe);
1238 incF(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)) +=
1239 (veloRoe - aRoe * n) * alpha0 + (veloRoe + aRoe * n) * alpha4 +
1240 veloRoe * alpha1 + alpha23VT;
1241 // ! handling the passive part of U!
1242 incF(SeqI52Last) += incU(SeqI52Last) * lam123;
1243 }
1244
1245 /**
1246 * @brief Batched Roe flux with selectable entropy fix for column-major
1247 * state matrices.
1248 *
1249 * Evaluates the Roe flux for @c nB = UL.cols() faces simultaneously. The
1250 * Roe averaging is performed from the single-column mean states @p ULm,
1251 * @p URm, sharing one set of Roe eigenvalues across the batch. Per-face
1252 * normal vectors @p n and grid velocities @p vg allow each column to have
1253 * a distinct orientation and motion; @p nm and @p vgm are the mean-state
1254 * normal/grid velocity used for the shared Roe average.
1255 *
1256 * @tparam dim Spatial dimension (2 or 3).
1257 * @tparam eigScheme Compile-time entropy-fix selector (0–8); default 0.
1258 * @param UL, UR Left/right face states [(dim+2) × nBatch].
1259 * @param ULm, URm Mean-state vectors for Roe averaging (single column).
1260 * @param vg, vgm Per-face and mean-state grid velocities.
1261 * @param n, nm Per-face and mean-state face normals.
1262 * @param gamma Ratio of specific heats (γ).
1263 * @param[out] F Numerical flux matrix [(dim+2) × nBatch].
1264 * @param dLambda H-correction transverse eigenvalue estimate.
1265 * @param fixScale Entropy-fix scaling factor.
1266 * @param incFScale Incremental flux dissipation scale.
1267 * @param dumpInfo Debug dump callable.
1268 * @param[out] lam0 |V_n − a| after entropy fixing.
1269 * @param[out] lam123 |V_n| after entropy fixing.
1270 * @param[out] lam4 |V_n + a| after entropy fixing.
1271 */
1272 template <int dim = 3, int eigScheme = 0,
1273 typename TUL, typename TUR,
1274 typename TULm, typename TURm,
1275 typename TVecVG, typename TVecVGm,
1276 typename TVecN, typename TVecNm,
1277 typename TF, typename TFdumpInfo>
1278 void RoeFlux_IdealGas_HartenYee_Batch(const TUL &UL, const TUR &UR,
1279 const TULm &ULm, const TURm &URm,
1280 const TVecVG &vg, const TVecVGm &vgm,
1281 const TVecN &n, const TVecNm &nm,
1282 real gamma, TF &F,
1283 real dLambda, real fixScale,
1284 real incFScale,
1285 const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
1286 {
1287 using TVec = Eigen::Vector<real, dim>;
1288 using TVec_Batch = Eigen::Matrix<real, dim, -1, Eigen::ColMajor, dim, MaxBatch>;
1289 using TReal_Batch = Eigen::Matrix<real, 1, -1, Eigen::RowMajor, 1, MaxBatch>;
1290 using TU5_Batch = Eigen::Matrix<real, dim + 2, -1, Eigen::ColMajor, dim + 2, MaxBatch>;
1291
1292 int nB = UL.cols();
1293 for (int iB = 0; iB < nB; iB++)
1294 if (!(UL(0, iB) > 0 && UR(0, iB) > 0))
1295 {
1296 dumpInfo();
1297 DNDS_assert(false);
1298 }
1299 TVec_Batch veloL = (UL(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll).array().rowwise() / UL(0, EigenAll).array()).matrix();
1300 TVec_Batch veloR = (UR(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll).array().rowwise() / UR(0, EigenAll).array()).matrix();
1301 TVec veloLm = (ULm(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / ULm(0)).matrix();
1302 TVec veloRm = (URm(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>)).array() / URm(0)).matrix();
1303
1304 TReal_Batch pL, pR;
1305 pL.resize(nB), pR.resize(nB);
1306 for (int iB = 0; iB < nB; iB++)
1307 {
1308 real asqrL, asqrR, HL, HR;
1309 real vsqrL = veloL(EigenAll, iB).squaredNorm();
1310 real vsqrR = veloR(EigenAll, iB).squaredNorm();
1311 IdealGasThermal(UL(dim + 1, iB), UL(0, iB), vsqrL, gamma, pL(iB), asqrL, HL);
1312 IdealGasThermal(UR(dim + 1, iB), UR(0, iB), vsqrR, gamma, pR(iB), asqrR, HR);
1313 }
1314
1315 real asqrLm, asqrRm, pLm, pRm, HLm, HRm;
1316 real vsqrLm = veloLm.squaredNorm();
1317 real vsqrRm = veloRm.squaredNorm();
1318 IdealGasThermal(ULm(dim + 1), ULm(0), vsqrLm, gamma, pLm, asqrLm, HLm);
1319 IdealGasThermal(URm(dim + 1), URm(0), vsqrRm, gamma, pRm, asqrRm, HRm);
1320
1321 real sqrtRhoLm = std::sqrt(ULm(0));
1322 real sqrtRhoRm = std::sqrt(URm(0));
1323
1324 TVec veloRoe = (sqrtRhoLm * veloLm + sqrtRhoRm * veloRm) / (sqrtRhoLm + sqrtRhoRm);
1325 real vsqrRoe = veloRoe.squaredNorm();
1326 real HRoe = (sqrtRhoLm * HLm + sqrtRhoRm * HRm) / (sqrtRhoLm + sqrtRhoRm);
1327 real asqrRoe = (gamma - 1) * (HRoe - 0.5 * vsqrRoe);
1328 real rhoRoe = sqrtRhoLm * sqrtRhoRm;
1329 real veloRoeN = veloRoe.dot(nm);
1330 real vgmN = vgm.dot(nm);
1331 real veloRoeRN = veloRoeN - vgmN;
1332 real veloLm0 = (veloLm - vgm).dot(nm);
1333 real veloRm0 = (veloRm - vgm).dot(nm);
1334
1335 TU5_Batch FL, FR;
1336 FL.resize(Eigen::NoChange, UL.cols());
1337 FR.resize(Eigen::NoChange, UL.cols());
1338 GasInviscidFlux_XY_Batch<dim>(UL, veloL, vg, n, pL, FL);
1339 GasInviscidFlux_XY_Batch<dim>(UR, veloR, vg, n, pR, FR);
1340
1341 if (!(asqrRoe > 0))
1342 {
1343 dumpInfo();
1344 }
1345 DNDS_assert(asqrRoe > 0);
1346 real aRoe = std::sqrt(asqrRoe);
1347
1348 lam0 = std::abs(veloRoeRN - aRoe);
1349 lam123 = std::abs(veloRoeRN);
1350 lam4 = std::abs(veloRoeRN + aRoe);
1351
1352 if constexpr (eigScheme == 2)
1353 {
1354 // *vanilla Lax
1355 // lam0 = lam123 = lam4 = std::max({lam0, lam123, lam4});
1356 lam0 = lam123 = lam4 = std::max(std::abs(veloLm0) + std::sqrt(asqrLm), std::abs(veloRm0) + std::sqrt(asqrRm));
1357 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll) =
1358 (FL + FR) * 0.5 -
1359 0.5 * lam0 * (UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll) - UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll));
1360 return; //* early exit
1361 }
1362 else
1363 Roe_EntropyFixer<eigScheme>(std::sqrt(asqrLm), std::sqrt(asqrRm), aRoe,
1364 veloLm0, veloRm0, veloRoeRN,
1365 (veloLm - vgm).norm(), (veloRm - vgm).norm(), (veloRoe - vgm).norm(),
1366 dLambda, fixScale, incFScale,
1367 lam0, lam123, lam4);
1368
1369 TU5_Batch incU =
1370 UR(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll) -
1371 UL(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll); //! not using m, for this is accuracy-limited!
1372 TReal_Batch incU123N =
1373 (incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll).array() * n.array()).colwise().sum();
1374 TVec_Batch alpha23V = incU(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll) - veloRoe * incU(0, EigenAll);
1375 TVec_Batch alpha23VT = alpha23V.array() - n.array().rowwise() * (alpha23V.array() * n.array()).colwise().sum();
1376 TReal_Batch incU4b =
1377 incU(dim + 1, EigenAll) -
1378 veloRoe.transpose() * alpha23VT;
1379 TReal_Batch alpha1 =
1380 (gamma - 1) / asqrRoe *
1381 (incU(0, EigenAll) * (HRoe - veloRoeN * veloRoeN) +
1382 veloRoeN * incU123N - incU4b);
1383 TReal_Batch alpha0 =
1384 (incU(0, EigenAll) * (veloRoeN + aRoe) - incU123N - aRoe * alpha1) / (2 * aRoe);
1385 TReal_Batch alpha4 =
1386 incU(0, EigenAll) - (alpha0 + alpha1);
1387
1388 alpha0 *= lam0;
1389 alpha1 *= lam123;
1390 alpha23VT *= lam123;
1391 alpha4 *= lam4; // here becomes alpha_i * lam_i
1392
1393 TU5_Batch incF;
1394 incF.resize(Eigen::NoChange, UL.cols());
1395 incF(0, EigenAll) = alpha0 + alpha1 + alpha4;
1396 incF(dim + 1, EigenAll) = (HRoe - veloRoeN * aRoe) * alpha0 + 0.5 * vsqrRoe * alpha1 +
1397 (HRoe + veloRoeN * aRoe) * alpha4 + veloRoe.transpose() * alpha23VT;
1398 incF(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll) =
1399 ((-aRoe * n).array().colwise() + veloRoe.array()).rowwise() * alpha0.array() +
1400 ((aRoe * n).array().colwise() + veloRoe.array()).rowwise() * alpha4.array() +
1401 alpha23VT.array() +
1402 (veloRoe * alpha1).array();
1403 // incF(Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>), EigenAll) =
1404 // (veloRoe.array() - (aRoe * n).array().colwise()) * alpha0 * (veloRoe.array() + (aRoe * n).array().colwise()) * alpha4;
1405
1406 F(Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>), EigenAll) = (FL + FR) * 0.5 - 0.5 * incF;
1407 }
1408
1409 /**
1410 * @brief Runtime dispatcher from RiemannSolverType to the correct Roe /
1411 * HLLC / HLLEP template instantiation.
1412 *
1413 * Maps RiemannSolverType enumerators to the corresponding compile-time
1414 * eigScheme / type template arguments. For Roe variants, calls
1415 * RoeFlux_IdealGas_HartenYee; for HLLEP, calls HLLEPFlux_IdealGas; for
1416 * HLLC, calls HLLCFlux_IdealGas_HartenYee.
1417 *
1418 * @tparam dim Spatial dimension (2 or 3).
1419 * @param type Riemann solver selector (runtime).
1420 * @param UL, UR Left/right conservative face states.
1421 * @param ULm, URm Mean states for Roe averaging.
1422 * @param vg Grid velocity.
1423 * @param n Face-normal.
1424 * @param gamma Ratio of specific heats (γ).
1425 * @param[out] F Numerical flux.
1426 * @param dLambda H-correction transverse eigenvalue estimate.
1427 * @param fixScale Entropy-fix scaling factor.
1428 * @param incFScale Incremental flux dissipation scale.
1429 * @param dumpInfo Debug dump callable.
1430 * @param[out] lam0 |V_n − a| after entropy fixing.
1431 * @param[out] lam123 |V_n| after entropy fixing.
1432 * @param[out] lam4 |V_n + a| after entropy fixing.
1433 */
1434 template <int dim = 3,
1435 typename TUL, typename TUR,
1436 typename TULm, typename TURm,
1437 typename TVecVG, typename TVecN,
1438 typename TF, typename TFdumpInfo>
1439 void InviscidFlux_IdealGas_Dispatcher(
1440 RiemannSolverType type,
1441 TUL &&UL, TUR &&UR, TULm &&ULm, TURm &&URm,
1442 TVecVG &&vg, TVecN &&n, real gamma, TF &&F,
1443 real dLambda, real fixScale,
1444 real incFScale,
1445 TFdumpInfo &&dumpInfo, real &lam0, real &lam123, real &lam4)
1446 {
1447#define DNDS_GAS_CALL_ROE(type) \
1448 RoeFlux_IdealGas_HartenYee<dim, type>( \
1449 UL, UR, ULm, URm, vg, n, gamma, F, dLambda, fixScale, incFScale, \
1450 dumpInfo, lam0, lam123, lam4)
1451
1452 if (type == Roe)
1453 RoeFlux_IdealGas_HartenYee<dim>(
1454 UL, UR, ULm, URm, vg, n, gamma, F, dLambda, fixScale, incFScale,
1455 dumpInfo, lam0, lam123, lam4);
1456 else if (type == Roe_M1)
1457 DNDS_GAS_CALL_ROE(1);
1458 else if (type == Roe_M2)
1459 DNDS_GAS_CALL_ROE(2);
1460 else if (type == Roe_M3)
1461 DNDS_GAS_CALL_ROE(3);
1462 else if (type == Roe_M4)
1463 DNDS_GAS_CALL_ROE(4);
1464 else if (type == Roe_M5)
1465 DNDS_GAS_CALL_ROE(5);
1466 else if (type == Roe_M6)
1467 DNDS_GAS_CALL_ROE(6);
1468 else if (type == Roe_M7)
1469 DNDS_GAS_CALL_ROE(7);
1470 else if (type == Roe_M8)
1471 DNDS_GAS_CALL_ROE(8);
1472 else if (type == Roe_M9)
1473 DNDS_GAS_CALL_ROE(9);
1474 else if (type == HLLEP)
1475 HLLEPFlux_IdealGas<dim, 0>(
1476 UL, UR, ULm, URm, vg, n, gamma, F, dLambda, fixScale,
1477 dumpInfo, lam0, lam123, lam4);
1478 else if (type == HLLEP_V1)
1479 HLLEPFlux_IdealGas<dim, 1>(
1480 UL, UR, ULm, URm, vg, n, gamma, F, dLambda, fixScale,
1481 dumpInfo, lam0, lam123, lam4);
1482 else if (type == HLLC)
1483 HLLCFlux_IdealGas_HartenYee<dim>(
1484 UL, UR, ULm, URm, vg, n, gamma, F, dLambda, fixScale,
1485 dumpInfo, lam0, lam123, lam4);
1486 else
1487 DNDS_assert_info(false, "the rs type is invalid");
1488#undef DNDS_GAS_CALL_ROE
1489 }
1490
1491 /**
1492 * @brief Runtime dispatcher for batched Roe flux from RiemannSolverType to
1493 * the correct RoeFlux_IdealGas_HartenYee_Batch template instantiation.
1494 *
1495 * Only Roe variants (Roe, Roe_M1..M9) are supported in batch mode; HLLC
1496 * and HLLEP are not available and will trigger an assertion failure.
1497 *
1498 * @tparam dim Spatial dimension (2 or 3).
1499 * @param type Riemann solver selector (runtime).
1500 * @param UL, UR Left/right face states [(dim+2) × nBatch].
1501 * @param ULm, URm Mean states for Roe averaging (single column).
1502 * @param vg, vgm Per-face and mean-state grid velocities.
1503 * @param n, nm Per-face and mean-state face normals.
1504 * @param gamma Ratio of specific heats (γ).
1505 * @param[out] F Numerical flux matrix [(dim+2) × nBatch].
1506 * @param dLambda H-correction transverse eigenvalue estimate.
1507 * @param fixScale Entropy-fix scaling factor.
1508 * @param incFScale Incremental flux dissipation scale.
1509 * @param dumpInfo Debug dump callable.
1510 * @param[out] lam0 |V_n − a| after entropy fixing.
1511 * @param[out] lam123 |V_n| after entropy fixing.
1512 * @param[out] lam4 |V_n + a| after entropy fixing.
1513 */
1514 template <int dim = 3,
1515 typename TUL, typename TUR,
1516 typename TULm, typename TURm,
1517 typename TVecVG, typename TVecVGm,
1518 typename TVecN, typename TVecNm,
1519 typename TF, typename TFdumpInfo>
1520 void InviscidFlux_IdealGas_Batch_Dispatcher(
1521 RiemannSolverType type,
1522 TUL &&UL, TUR &&UR,
1523 TULm &&ULm, TURm &&URm,
1524 TVecVG &&vg, TVecVGm &&vgm,
1525 TVecN &&n, TVecNm &&nm,
1526 real gamma, TF &&F, real dLambda, real fixScale,
1527 real incFScale,
1528 TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
1529 {
1530#define DNDS_GAS_CALL_ROE(type) \
1531 RoeFlux_IdealGas_HartenYee_Batch<dim, type>( \
1532 UL, UR, ULm, URm, vg, vgm, n, nm, gamma, F, dLambda, fixScale, incFScale, \
1533 dumpInfo, lam0, lam123, lam4)
1534
1535 if (type == Roe)
1536 RoeFlux_IdealGas_HartenYee_Batch<dim>(
1537 UL, UR, ULm, URm, vg, vgm, n, nm, gamma, F, dLambda, fixScale, incFScale,
1538 dumpInfo, lam0, lam123, lam4);
1539 else if (type == Roe_M1)
1540 DNDS_GAS_CALL_ROE(1);
1541 else if (type == Roe_M2)
1542 DNDS_GAS_CALL_ROE(2);
1543 else if (type == Roe_M3)
1544 DNDS_GAS_CALL_ROE(3);
1545 else if (type == Roe_M4)
1546 DNDS_GAS_CALL_ROE(4);
1547 else if (type == Roe_M5)
1548 DNDS_GAS_CALL_ROE(5);
1549 else if (type == Roe_M6)
1550 DNDS_GAS_CALL_ROE(6);
1551 else if (type == Roe_M7)
1552 DNDS_GAS_CALL_ROE(7);
1553 else if (type == Roe_M8)
1554 DNDS_GAS_CALL_ROE(8);
1555 else if (type == Roe_M9)
1556 DNDS_GAS_CALL_ROE(9);
1557 else
1558 DNDS_assert_info(false, "the rs type is invalid (for batch version)");
1559#undef DNDS_GAS_CALL_ROE
1560 }
1561
1562 /**
1563 * @brief Computes the viscous (Navier-Stokes) flux projected onto a face
1564 * normal for an ideal gas.
1565 *
1566 * Evaluates the viscous stress tensor τ (including Stokes' hypothesis
1567 * λ = −2/3 μ) and the heat flux vector q = −k ∇T, then projects onto
1568 * @p norm. Supports:
1569 * - Sutherland-law viscosity (via the externally provided @p mu).
1570 * - Turbulent eddy viscosity through @p mutRatio (μ_t / μ).
1571 * - QCR (Quadratic Constitutive Relation) correction when @p mutQCRFix
1572 * is true, using ccr1 = 0.3.
1573 * - Adiabatic wall treatment: when @p adiabatic is true, the wall-normal
1574 * component of ∇T is removed.
1575 *
1576 * @note @p GradUPrim has size [dim × nVars]; it is the gradient of the
1577 * **primitive** variables (ρ, u, v, [w,] p, ...) in physical space.
1578 * "3×5 TGradU" in the original comment refers to dim=3 with 5 Euler
1579 * variables.
1580 *
1581 * @tparam dim Spatial dimension (2 or 3).
1582 * @param U Conservative state vector.
1583 * @param GradUPrim Gradient of primitive variables [dim × nVars].
1584 * @param norm Face-normal vector (unit or area-weighted).
1585 * @param adiabatic If true, zero the wall-normal heat flux component.
1586 * @param gamma Ratio of specific heats (γ).
1587 * @param mu Dynamic (molecular + turbulent) viscosity μ + μ_t.
1588 * @param mutRatio Turbulent-to-molecular viscosity ratio μ_t / μ.
1589 * @param mutQCRFix Enable QCR stress correction.
1590 * @param k Thermal conductivity (molecular + turbulent).
1591 * @param Cp Specific heat at constant pressure.
1592 * @param[out] Flux Viscous flux vector (dim+2 entries, Flux[0] = 0).
1593 */
1594 template <int dim = 3, typename TU, typename TGradU, typename TFlux, typename TNorm>
1595 void ViscousFlux_IdealGas(const TU &U, const TGradU &GradUPrim, TNorm norm, bool adiabatic, real gamma,
1596 real mu, real mutRatio, bool mutQCRFix, real k, real Cp, TFlux &Flux)
1597 {
1598 static const auto Seq01234 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>);
1599 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
1600 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
1601
1602 Eigen::Vector<real, dim> velo = U(Seq123) / U(0);
1603 static const real lambda = -2. / 3.;
1604 Eigen::Matrix<real, dim, dim> diffVelo = GradUPrim(Seq012, Seq123); // dU_j/dx_i
1605 Eigen::Vector<real, dim> GradP = GradUPrim(Seq012, dim + 1);
1606 real vSqr = velo.squaredNorm();
1607 real p = (gamma - 1) * (U(dim + 1) - U(0) * 0.5 * vSqr);
1608 Eigen::Vector<real, dim> GradT = (gamma / ((gamma - 1) * Cp * U(0) * U(0))) *
1609 (U(0) * GradP - p * GradUPrim(Seq012, 0)); // GradU(:,0) is grad rho no matter prim or not
1610
1611 if (adiabatic) //! is this fix reasonable?
1612 GradT -= GradT.dot(norm) * norm;
1613
1614 Eigen::Matrix<real, dim, dim> vStress = (diffVelo + diffVelo.transpose()) * mu +
1615 Eigen::Matrix<real, dim, dim>::Identity() * (lambda * mu * diffVelo.trace());
1616 if (mutQCRFix)
1617 {
1618 real b = std::sqrt((diffVelo.array() * diffVelo.array()).sum());
1619 Eigen::Matrix<real, dim, dim> O = (diffVelo.transpose() - diffVelo) / (b + verySmallReal); // dU_i/dx_j-dU_j/dx_i
1620 real ccr1 = 0.3;
1621 Eigen::Matrix<real, dim, dim> vStressQCRFix;
1622 vStressQCRFix.setZero();
1623 vStressQCRFix.diagonal() = (vStress.array() * O.array()).rowwise().sum();
1624 vStressQCRFix(0, 1) = O(0, 1) * (vStress(1, 1) - vStress(0, 0));
1625 if (dim == 3)
1626 {
1627 vStressQCRFix(0, 1) += O(1, 2) * vStress(0, 2) + O(0, 2) * vStress(1, 2);
1628 vStressQCRFix(0, 2) = O(0, 2) * (vStress(2, 2) - vStress(0, 0)) + O(0, 1) * vStress(2, 1) + O(2, 1) * vStress(0, 1);
1629 vStressQCRFix(1, 2) = O(1, 2) * (vStress(2, 2) - vStress(1, 1)) + O(1, 0) * vStress(2, 0) + O(2, 0) * vStress(1, 0);
1630 }
1631 Eigen::Matrix<real, dim, dim> vStressQCRFixFull = vStressQCRFix + vStressQCRFix.transpose();
1632 vStress -= ccr1 * mutRatio * vStressQCRFixFull;
1633 }
1634 Flux(0) = 0;
1635 Flux(Seq123) = vStress * norm;
1636 Flux(dim + 1) = (vStress * velo + k * GradT).dot(norm);
1637 if (!Flux.allFinite())
1638 {
1639 std::cout << "U\n"
1640 << U.transpose() << "\n";
1641 std::cout << "GradUPrim\n"
1642 << GradUPrim << "\n";
1643 std::cout << "norm\n"
1644 << norm << "\n";
1645 std::cout << "gamma\n"
1646 << gamma << "\n";
1647 std::cout << "mu\n"
1648 << mu << "\n";
1649 std::cout << "k\n"
1650 << k << "\n";
1651 std::cout << "Cp\n"
1652 << Cp << "\n";
1653 DNDS_assert(false);
1654 }
1655 }
1656
1657 /**
1658 * @brief Converts the gradient of conservative variables to the gradient of
1659 * primitive variables for an ideal gas.
1660 *
1661 * Given GradU = ∂U/∂x (the spatial gradient of the conservative vector),
1662 * computes GradUPrim = ∂(ρ,u,v,[w,]p,...)/∂x in-place. The density
1663 * gradient column (column 0) is unchanged, velocity gradient columns are
1664 * transformed via the quotient rule, and the pressure gradient column uses
1665 * the chain rule through the equation of state.
1666 *
1667 * @note The input @p GradU has size [dim × nVars], where nVars ≥ dim+2.
1668 * Passive-scalar gradient columns beyond index dim+1 are also
1669 * converted (divided by ρ after removing the density contribution).
1670 *
1671 * @tparam dim Spatial dimension (2 or 3).
1672 * @param U Conservative state vector.
1673 * @param GradU Gradient of conservative variables [dim × nVars].
1674 * @param[out] GradUPrim Gradient of primitive variables [dim × nVars].
1675 * @param gamma Ratio of specific heats (γ).
1676 */
1677 template <int dim = 3, typename TU, typename TGradU, typename TGradUPrim>
1678 void GradientCons2Prim_IdealGas(const TU &U, const TGradU &GradU, TGradUPrim &GradUPrim, real gamma)
1679 {
1680 static const auto Seq01234 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>);
1681 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
1682 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
1683 static const auto I4 = dim + 1;
1684
1685 Eigen::Vector<real, dim> velo = U(Seq123) / U(0);
1686 GradUPrim = GradU;
1687
1688 GradUPrim(Seq012, Seq123) = (1.0 / sqr(U(0))) *
1689 (U(0) * GradU(Seq012, Seq123) -
1690 GradU(Seq012, 0) * Eigen::RowVector<real, dim>(U(Seq123))); // dU_j/dx_i
1691 GradUPrim(Seq012, I4) = (gamma - 1) *
1692 (GradU(Seq012, dim + 1) -
1693 0.5 *
1694 (GradU(Seq012, Seq123) * velo +
1695 GradUPrim(Seq012, Seq123) * Eigen::Vector<real, dim>(U(Seq123))));
1696 GradUPrim(Seq012, Eigen::seq(Eigen::fix<I4 + 1>, EigenLast)) -= GradU(Seq012, 0) * U(Eigen::seq(Eigen::fix<I4 + 1>, EigenLast)).transpose() / U(0);
1697 GradUPrim(Seq012, Eigen::seq(Eigen::fix<I4 + 1>, EigenLast)) /= U(0);
1698 }
1699
1700 /**
1701 * @brief Extracts the velocity gradient tensor from the conservative-variable
1702 * gradient using the quotient rule.
1703 *
1704 * Returns diffVelo(i,j) = ∂u_j/∂x_i, a [dim × dim] matrix, computed as
1705 * (ρ · ∂(ρu_j)/∂x_i − (ρu_j) · ∂ρ/∂x_i) / ρ².
1706 *
1707 * @tparam dim Spatial dimension (2 or 3).
1708 * @param U Conservative state vector.
1709 * @param GradU Gradient of conservative variables [dim × nVars].
1710 * @return Eigen::Matrix<real, dim, dim> Velocity gradient tensor (dU_j/dx_i).
1711 */
1712 template <int dim, typename TU, typename TGradU>
1713 auto GetGradVelo(const TU &U, const TGradU &GradU)
1714 {
1715 static const auto Seq01234 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>);
1716 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
1717 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
1718 Eigen::Matrix<real, dim, dim> diffVelo = (1.0 / sqr(U(0))) *
1719 (U(0) * GradU(Seq012, Seq123) -
1720 GradU(Seq012, 0) * Eigen::RowVector<real, dim>(U(Seq123))); // dU_j/dx_i
1721 return diffVelo;
1722 }
1723
1724 /**
1725 * @brief Computes the maximum safe compression ratio (scaling factor α ∈ [0,1])
1726 * for a conservative-variable increment that keeps internal energy
1727 * above a prescribed positive floor.
1728 *
1729 * Given a state @p u and an increment @p uInc, finds the largest α such that
1730 * the internal energy e = E − ½ρ|v|² of (u + α·uInc) remains above
1731 * @p newrhoEinteralNew. Two schemes are available:
1732 * - scheme 0: analytic quadratic solve with iterative safety fallback.
1733 * - scheme 1: convex (linear) estimation only.
1734 *
1735 * Used during implicit time stepping and limiting to prevent negative
1736 * pressures.
1737 *
1738 * TODO: vectorize
1739 *
1740 * @tparam dim Spatial dimension (2 or 3).
1741 * @tparam scheme 0 = quadratic solve, 1 = convex estimate.
1742 * @tparam nVarsFixed Compile-time size of the state vector.
1743 * @param u Current conservative state.
1744 * @param uInc Proposed conservative state increment.
1745 * @param newrhoEinteralNew Desired minimum internal energy floor ρ·e_min
1746 * (= p_min / (γ−1)).
1747 * @return The safe scaling factor α in [0, 1].
1748 */
1749 template <int dim = 3, int scheme = 0, int nVarsFixed, typename TU, typename TUInc>
1750 real IdealGasGetCompressionRatioPressure(const TU &u, const TUInc &uInc, real newrhoEinteralNew)
1751 {
1752 static const real safetyRatio = 1 - 1e-5;
1753 static const auto Seq01234 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim + 1>);
1754 static const auto Seq012 = Eigen::seq(Eigen::fix<0>, Eigen::fix<dim - 1>);
1755 static const auto Seq123 = Eigen::seq(Eigen::fix<1>, Eigen::fix<dim>);
1756 static const auto I4 = dim + 1;
1757
1758 Eigen::Vector<real, nVarsFixed> ret = uInc;
1759 Eigen::Vector<real, nVarsFixed> uNew = u + uInc;
1760 real rhoEOld = u(I4) - u(Seq123).squaredNorm() / (u(0) + verySmallReal) * 0.5;
1761 newrhoEinteralNew = std::max(smallReal * rhoEOld, newrhoEinteralNew);
1762 real rhoENew = uNew(I4) - uNew(Seq123).squaredNorm() / (uNew(0) + verySmallReal) * 0.5;
1763 real alphaEst1 = (rhoEOld - newrhoEinteralNew) / std::max(-rhoENew + rhoEOld, verySmallReal);
1764 if (rhoENew > rhoEOld)
1765 alphaEst1 = 1;
1766 alphaEst1 = std::min(alphaEst1, 1.);
1767 alphaEst1 = std::max(alphaEst1, 0.);
1768 real alpha = alphaEst1; //! using convex estimation
1769
1770 real alphaL, alphaR, c0, c1, c2;
1771 alphaL = alphaR = c0 = c1 = c2 = 0;
1772 if constexpr (scheme == 0)
1773 {
1774 c0 = 2 * u(I4) * u(0) - u(Seq123).squaredNorm() - 2 * u(0) * newrhoEinteralNew;
1775 c1 = 2 * u(I4) * ret(0) + 2 * u(0) * ret(I4) - 2 * u(Seq123).dot(ret(Seq123)) - 2 * ret(0) * newrhoEinteralNew;
1776 c2 = 2 * ret(I4) * ret(0) - ret(Seq123).squaredNorm();
1777 c2 += signP(c2) * verySmallReal;
1778 real deltaC = sqr(c1) - 4 * c0 * c2;
1779 if (deltaC <= -sqr(c0) * smallReal)
1780 {
1781 std::cout << std::scientific << std::setprecision(5);
1782 std::cout << u.transpose() << std::endl;
1783 std::cout << uInc.transpose() << std::endl;
1784 std::cout << newrhoEinteralNew << std::endl;
1785 std::cout << fmt::format("{} {} {}", c0, c1, c2) << std::endl;
1786
1787 DNDS_assert(false);
1788 }
1789 deltaC = std::max(0., deltaC);
1790 real alphaL = (-std::sqrt(deltaC) - c1) / (2 * c2);
1791 real alphaR = (std::sqrt(deltaC) - c1) / (2 * c2);
1792 // if (c2 > 0)
1793 // DNDS_assert(alphaL > 0);
1794 // DNDS_assert(alphaR > 0);
1795 // DNDS_assert(alphaL < 1);
1796 // if (c2 < 0)
1797 // DNDS_assert(alphaR < 1);
1798
1799 if (std::abs(c2) < 1e-10 * c0)
1800 {
1801 if (std::abs(c1) < 1e-10 * c0)
1802 {
1803 alpha = 0;
1804 }
1805 else
1806 {
1807 alpha = std::min(-c0 / c1, 1.);
1808 }
1809 }
1810 else
1811 {
1812 alpha = std::min((c2 > 0 ? alphaL : alphaL), 1.);
1813 }
1814 alpha = std::max(0., alpha);
1815 alpha *= safetyRatio;
1816 if (alpha < smallReal)
1817 alpha = 0;
1818 }
1819 else if constexpr (scheme == 1)
1820 {
1821 // has used convex
1822 alpha *= safetyRatio;
1823 if (alpha < smallReal)
1824 alpha = 0;
1825 }
1826
1827 ret *= alpha;
1828
1829 // Eigen::Vector<real, nVarsFixed> uNew = u + ret;
1830 // real eNew = uNew(I4) - 0.5 * uNew(Seq123).squaredNorm() / uNew(0);
1831
1832 real decay = 1 - 1e-2;
1833 int iter;
1834 for (iter = 0; iter < 1000; iter++)
1835 {
1836 real ek = 0.5 * (u(Seq123) + ret(Seq123)).squaredNorm() / (u(0) + ret(0) + verySmallReal);
1837 if (ret(I4) + u(I4) - ek < newrhoEinteralNew)
1838 {
1839
1840 ret *= decay, alpha *= decay;
1841 }
1842 else
1843 break;
1844 }
1845 if (iter >= 1000)
1846 {
1847 real ek = 0.5 * (u(Seq123) + ret(Seq123)).squaredNorm() / (u(0) + ret(0) + verySmallReal);
1848 {
1849 std::cout << std::scientific << std::setprecision(5);
1850 std::cout << fmt::format("alphas: {}, {}, {}\n", alpha, alphaL, alphaR);
1851 std::cout << fmt::format("ABC: {}, {}, {}\n", c2, c1, c0);
1852 std::cout << u.transpose() << std::endl;
1853 std::cout << uInc.transpose() << std::endl;
1854 std::cout << fmt::format("eks: {} {}\n", ret(I4) + u(I4) - ek, newrhoEinteralNew);
1855 DNDS_assert(false);
1856 }
1857 }
1858 // std::cout << fmt::format("{} {} {} {} {}", c0, c1, c2, alphaL, alphaR) << std::endl;
1859
1860 return alpha;
1861 }
1862}
ConfigEnum.hpp
Extended enum-to-JSON serialization macro that also exposes allowed string values for JSON Schema gen...
DNDS_DEFINE_ENUM_JSON
#define DNDS_DEFINE_ENUM_JSON(EnumType_,...)
Define JSON serialization for an enum AND expose its allowed string values.
Definition ConfigEnum.hpp:126
Defines.hpp
Core type aliases, constants, and metaprogramming utilities for the DNDS framework.
DNDS_assert_info
#define DNDS_assert_info(expr, info)
Debug-only assertion with an extra std::string info message.
Definition Errors.hpp:113
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
DNDS_GAS_CALL_ROE
#define DNDS_GAS_CALL_ROE(type)
IdealGasPhysics.hpp
Shared ideal-gas thermodynamics and Roe-flux primitives.
JsonUtil.hpp
JSON-to-Eigen conversion utilities and nlohmann_json helper macros.
DNDS::Euler::Gas::HLLEPFlux_IdealGas
void HLLEPFlux_IdealGas(const TUL &UL, const TUR &UR, const TULm &ULm, const TURm &URm, const TVecVG &vg, const TVecN &n, real gamma, TF &F, real dLambda, real fixScale, const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
HLLEP (HLL with Enhanced Pressure) approximate Riemann solver for an ideal gas on a moving grid.
Definition Gas.hpp:605
DNDS::Euler::Gas::GradientCons2Prim_IdealGas
void GradientCons2Prim_IdealGas(const TU &U, const TGradU &GradU, TGradUPrim &GradUPrim, real gamma)
Converts the gradient of conservative variables to the gradient of primitive variables for an ideal g...
Definition Gas.hpp:1678
DNDS::Euler::Gas::RoeFluxIncFDiff
void RoeFluxIncFDiff(const TDU &incU, const TVecN &n, const TVecV &veloRoe, real vsqrRoe, real aRoe, real asqrRoe, real HRoe, real lam0, real lam123, real lam4, real gamma, TDF &incF)
Computes the dissipation part of the Roe flux |A|·dU, given pre-computed Roe averages and entropy-fix...
Definition Gas.hpp:1210
DNDS::Euler::Gas::ComputeRoePreamble
RoePreamble< dim > ComputeRoePreamble(const TULm &ULm, const TURm &URm, real gamma, const TFdumpInfo &dumpInfo)
Compute Roe-averaged quantities from mean-state L/R vectors.
Definition Gas.hpp:230
DNDS::Euler::Gas::RiemannSolverType
RiemannSolverType
Selects the approximate Riemann solver and its entropy-fix variant.
Definition Gas.hpp:62
DNDS::Euler::Gas::GasInviscidFlux_XY_Batch
void GasInviscidFlux_XY_Batch(const TU &U, const TVec &velo, const TVecVG &vg, const TVecN &n, TP &&p, TF &F)
Batched face-normal inviscid flux for column-major state matrices.
Definition Gas.hpp:437
DNDS::Euler::Gas::tVec2
Eigen::Vector2d tVec2
Convenience alias for 2-D real vector.
Definition Gas.hpp:31
DNDS::Euler::Gas::RoeFlux_IdealGas_HartenYee_Batch
void RoeFlux_IdealGas_HartenYee_Batch(const TUL &UL, const TUR &UR, const TULm &ULm, const TURm &URm, const TVecVG &vg, const TVecVGm &vgm, const TVecN &n, const TVecNm &nm, real gamma, TF &F, real dLambda, real fixScale, real incFScale, const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
Batched Roe flux with selectable entropy fix for column-major state matrices.
Definition Gas.hpp:1278
DNDS::Euler::Gas::GasInviscidFlux
void GasInviscidFlux(const TU &U, const TVec &velo, const TVecVG &vg, real p, TF &F)
Computes the inviscid (Euler) flux in the x-direction for a moving grid.
Definition Gas.hpp:365
DNDS::Euler::Gas::ViscousFlux_IdealGas
void ViscousFlux_IdealGas(const TU &U, const TGradU &GradUPrim, TNorm norm, bool adiabatic, real gamma, real mu, real mutRatio, bool mutQCRFix, real k, real Cp, TFlux &Flux)
Computes the viscous (Navier-Stokes) flux projected onto a face normal for an ideal gas.
Definition Gas.hpp:1595
DNDS::Euler::Gas::HLLCFlux_IdealGas_HartenYee
void HLLCFlux_IdealGas_HartenYee(const TUL &UL, const TUR &UR, const TULm &ULm, const TURm &URm, const TVecVG &vg, const TVecN &n, real gamma, TF &F, real dLambda, real fixScale, const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
HLLC (Harten-Lax-van Leer-Contact) approximate Riemann solver for an ideal gas on a moving grid.
Definition Gas.hpp:740
DNDS::Euler::Gas::IdealGasThermal
void IdealGasThermal(real E, real rho, real vSqr, real gamma, real &p, real &asqr, real &H)
Thin wrapper delegating to IdealGas::IdealGasThermal.
Definition Gas.hpp:190
DNDS::Euler::Gas::GasInviscidFluxFacialIncrement
void GasInviscidFluxFacialIncrement(const TU &U, const TU &dU, const TVec &unitNorm, const TVecVG &velo, const TVec &dVelo, const TVec &vg, real dp, real p, TF &F)
Computes the increment of the facial inviscid flux from state, velocity, and pressure increments.
Definition Gas.hpp:496
DNDS::Euler::Gas::RoeFlux_IdealGas_HartenYee
void RoeFlux_IdealGas_HartenYee(const TUL &UL, const TUR &UR, const TULm &ULm, const TURm &URm, const TVecVG &vg, const TVecN &n, real gamma, TF &F, real dLambda, real fixScale, real incFScale, const TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
Core Roe approximate Riemann solver with a selectable entropy-fix scheme for an ideal gas on a moving...
Definition Gas.hpp:1043
DNDS::Euler::Gas::InviscidFlux_IdealGas_Dispatcher
void InviscidFlux_IdealGas_Dispatcher(RiemannSolverType type, TUL &&UL, TUR &&UR, TULm &&ULm, TURm &&URm, TVecVG &&vg, TVecN &&n, real gamma, TF &&F, real dLambda, real fixScale, real incFScale, TFdumpInfo &&dumpInfo, real &lam0, real &lam123, real &lam4)
Runtime dispatcher from RiemannSolverType to the correct Roe / HLLC / HLLEP template instantiation.
Definition Gas.hpp:1439
DNDS::Euler::Gas::GasInviscidFlux_XY
void GasInviscidFlux_XY(const TU &U, const TVec &velo, const TVecVG &vg, const TVecN &n, real p, TF &F)
Computes the inviscid flux projected onto an arbitrary face normal n, accounting for grid motion vg.
Definition Gas.hpp:391
DNDS::Euler::Gas::InviscidFlux_IdealGas_Batch_Dispatcher
void InviscidFlux_IdealGas_Batch_Dispatcher(RiemannSolverType type, TUL &&UL, TUR &&UR, TULm &&ULm, TURm &&URm, TVecVG &&vg, TVecVGm &&vgm, TVecN &&n, TVecNm &&nm, real gamma, TF &&F, real dLambda, real fixScale, real incFScale, TFdumpInfo &dumpInfo, real &lam0, real &lam123, real &lam4)
Runtime dispatcher for batched Roe flux from RiemannSolverType to the correct RoeFlux_IdealGas_Harten...
Definition Gas.hpp:1520
DNDS::Euler::Gas::EulerGasLeftEigenVector
void EulerGasLeftEigenVector(const TVec &velo, real Vsqr, real H, real a, real gamma, TeV &LeV)
Fills the left eigenvector matrix (inverse of the right eigenvector matrix) for the 1-D Euler system ...
Definition Gas.hpp:160
DNDS::Euler::Gas::IdealGasThermalPrimitive2Conservative
void IdealGasThermalPrimitive2Conservative(const TPrim &prim, TCons &U, real gamma)
Converts primitive variables to conservative variables for an ideal gas.
Definition Gas.hpp:306
DNDS::Euler::Gas::IdealGasGetCompressionRatioPressure
real IdealGasGetCompressionRatioPressure(const TU &u, const TUInc &uInc, real newrhoEinteralNew)
Computes the maximum safe compression ratio (scaling factor α ∈ [0,1]) for a conservative-variable in...
Definition Gas.hpp:1750
DNDS::Euler::Gas::IdealGas_EulerGasRightEigenVector
auto IdealGas_EulerGasRightEigenVector(const TU &U, real gamma)
Convenience wrapper that computes the right eigenvector matrix directly from a conservative state vec...
Definition Gas.hpp:528
DNDS::Euler::Gas::GasInviscidFlux_Batch
void GasInviscidFlux_Batch(const TU &U, const TVec &velo, const TVecVG &vg, TP &&p, TF &F)
Batched x-direction inviscid flux for column-major state matrices.
Definition Gas.hpp:414
DNDS::Euler::Gas::IdealGasThermalPrimitiveGetP0T0
std::tuple< real, real > IdealGasThermalPrimitiveGetP0T0(const TPrim &prim, real gamma, real rg)
Computes total (stagnation) pressure p0 and temperature T0 from a primitive state using isentropic re...
Definition Gas.hpp:333
DNDS::Euler::Gas::tVec
Eigen::Vector3d tVec
Convenience alias for 3-D real vector.
Definition Gas.hpp:30
DNDS::Euler::Gas::IdealGasThermalConservative2Primitive
void IdealGasThermalConservative2Primitive(const TCons &U, TPrim &prim, real gamma)
Converts conservative variables to primitive variables for an ideal gas.
Definition Gas.hpp:279
DNDS::Euler::Gas::EulerGasRightEigenVector
class TeV inline void EulerGasRightEigenVector(const TVec &velo, real Vsqr, real H, real a, TeV &ReV)
Definition Gas.hpp:120
DNDS::Euler::Gas::IdealGasUIncrement
void IdealGasUIncrement(const TU &U, const TU &dU, const TVec &velo, real gamma, TVec &dVelo, real &dp)
Computes velocity and pressure increments from a conservative-state increment, used for the Lax-flux ...
Definition Gas.hpp:463
DNDS::Euler::Gas::IdealGas_EulerGasLeftEigenVector
auto IdealGas_EulerGasLeftEigenVector(const TU &U, real gamma)
Convenience wrapper that computes the left eigenvector matrix directly from a conservative state vect...
Definition Gas.hpp:555
DNDS::Euler::Gas::GetGradVelo
auto GetGradVelo(const TU &U, const TGradU &GradU)
Extracts the velocity gradient tensor from the conservative-variable gradient using the quotient rule...
Definition Gas.hpp:1713
DNDS::Euler::Gas::Roe_EntropyFixer
void Roe_EntropyFixer(const real aL, const real aR, const real aAve, const real uL, const real uR, const real uAve, const real VL, const real VR, const real VAve, real dLambda, real fixScale, real incFScale, real &lam0, real &lam123, real &lam4)
Template-dispatched entropy-fix for Roe-type Riemann solvers.
Definition Gas.hpp:879
DNDS::Euler::Gas::GetRoeAverage
void GetRoeAverage(const TUL &UL, const TUR &UR, real gamma, TVecV &veloRoe, real &vsqrRoe, real &aRoe, real &asqrRoe, real &HRoe, TUOut &UOut)
Computes the Roe-averaged state vector including passive scalars (e.g. RANS turbulence variables).
Definition Gas.hpp:1149
DNDS::IdealGas::IdealGasThermal
DNDS_DEVICE_CALLABLE void IdealGasThermal(real E, real rho, real vSqr, real gamma, real &p, real &asqr, real &H)
Compute pressure, speed-of-sound squared, and specific enthalpy from total energy,...
Definition IdealGasPhysics.hpp:38
DNDS::signM
constexpr real signM(real a)
"Signum, biased toward -1": treats 0 as negative.
Definition Defines.hpp:550
DNDS::verySmallReal
DNDS_CONSTANT const real verySmallReal
Catch-all lower bound ("effectively zero").
Definition Defines.hpp:194
DNDS::signTol
constexpr real signTol(real a, real tol)
Tolerant signum: returns 0 inside [-tol, tol].
Definition Defines.hpp:538
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::signP
constexpr real signP(real a)
"Signum, biased toward +1": treats 0 as positive.
Definition Defines.hpp:544
DNDS::smallReal
DNDS_CONSTANT const real smallReal
Loose lower bound (for iterative-solver tolerances etc.).
Definition Defines.hpp:196
DNDS::Euler::Gas::RoePreamble
Pre-computed Roe-averaged quantities shared by all Riemann solvers.
Definition Gas.hpp:205
DNDS::Euler::Gas::RoePreamble::veloRoe
TVec veloRoe
Roe-averaged velocity vector.
Definition Gas.hpp:212
DNDS::Euler::Gas::RoePreamble::TVec
Eigen::Vector< real, dim > TVec
Definition Gas.hpp:206
DNDS::Euler::Gas::RoePreamble::sqrtRhoRm
real sqrtRhoRm
√ρ for L/R states (used in Roe weighting).
Definition Gas.hpp:213
DNDS::Euler::Gas::RoePreamble::veloRm
TVec veloRm
Left/right mean-state velocity vectors.
Definition Gas.hpp:208
DNDS::Euler::Gas::RoePreamble::aRoe
real aRoe
Roe-averaged |v|², H, a², ρ, and speed of sound.
Definition Gas.hpp:214
DNDS::Euler::Gas::RoePreamble::HRm
real HRm
Speed-of-sound², pressure, and enthalpy for L/R mean states.
Definition Gas.hpp:209
DNDS::Euler::Gas::RoePreamble::vsqrRm
real vsqrRm
Squared velocity magnitudes for L/R mean states.
Definition Gas.hpp:210
setConstant
dof1 setConstant(0.0)
v
Eigen::Matrix< real, 5, 1 > v
Definition test_ArrayDOF.cpp:468
b
tVec b(NCells)
alpha
real alpha
Definition test_GasThermo.cpp:454
vg
Eigen::Vector3d vg(0, 0, 0)
velo
Eigen::Vector3d velo
Definition test_GasThermo.cpp:282
dU
Eigen::Vector< real, 5 > dU
Definition test_GasThermo.cpp:452
n
Eigen::Vector3d n(1.0, 0.0, 0.0)