From c0dc0cee396ba32632edde575713d7629c8c652c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 13:52:07 +0200 Subject: [PATCH 1/7] Include the SST k of the total energy in the convective numerics With SST the total energy contains k (CNSVariable subtracts it from the static energy, the free-stream and boundary energies add it), but the convective schemes computed the speed of sound, the eigenvector matrices and the Jacobians as if it did not. Add a k argument (held fixed, default 0) to GetPMatrix, GetPMatrix_inv and GetInviscidProjJac, and to the SIMD pMatrix, pMatrixInv, inviscidProjJac and Roe average. Remove k from the speed of sound of Roe, HLLC, AUSM-type and FVS schemes, use the Roe-averaged k in the Roe matrices, advect the part of Delta(rho E) due to Delta k with the contact wave (Roe, JST-mat), and include k in the accurate AUSM and HLLC Jacobians and in the preaccumulation inputs. Upwind_Residual passes the k of the cells to the numerics, adds it to the total enthalpy of the consistent and low-Mach extrapolations, and checks the speed of sound of each MUSCL face state as well as of the Roe average. Unit tests for the Roe and AUSM contacts with a k jump and the HLLC Jacobians with k. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/include/numerics/CNumerics.hpp | 24 ++- .../include/numerics/flow/convection/roe.hpp | 1 + .../flow/convection/centered.hpp | 33 +++- .../numerics_simd/flow/convection/common.hpp | 29 ++- .../numerics_simd/flow/convection/upwind.hpp | 44 +++-- .../include/numerics_simd/flow/variables.hpp | 9 +- SU2_CFD/include/solvers/CEulerSolver.hpp | 10 +- SU2_CFD/src/numerics/CNumerics.cpp | 25 +-- .../numerics/flow/convection/ausm_slau.cpp | 56 +++--- SU2_CFD/src/numerics/flow/convection/fvs.cpp | 17 +- SU2_CFD/src/numerics/flow/convection/hllc.cpp | 44 +++-- SU2_CFD/src/numerics/flow/convection/roe.cpp | 61 ++++-- SU2_CFD/src/solvers/CEulerSolver.cpp | 39 +++- UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp | 173 +++++++++++++++++ .../SU2_CFD/numerics/hllc_jacobian_tests.cpp | 11 +- UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp | 176 ++++++++++++++++++ UnitTests/meson.build | 2 + 17 files changed, 625 insertions(+), 129 deletions(-) create mode 100644 UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp create mode 100644 UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp diff --git a/SU2_CFD/include/numerics/CNumerics.hpp b/SU2_CFD/include/numerics/CNumerics.hpp index ec511a3fa1db..b26e911c3542 100644 --- a/SU2_CFD/include/numerics/CNumerics.hpp +++ b/SU2_CFD/include/numerics/CNumerics.hpp @@ -1146,10 +1146,11 @@ class CNumerics { * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. * \param[in] val_scale - Scale of the projection. * \param[out] val_Proj_Jac_tensor - Pointer to the projected inviscid Jacobian. + * \param[in] val_tke - Turbulent kinetic energy contained in the total energy (SST), held fixed. */ void GetInviscidProjJac(const su2double *val_velocity, const su2double *val_energy, const su2double *val_normal, su2double val_scale, - su2double **val_Proj_Jac_tensor) const; + su2double **val_Proj_Jac_tensor, su2double val_tke = 0.0) const; /*! * \brief Compute the projection of the inviscid Jacobian matrices (incompressible). @@ -1222,11 +1223,12 @@ class CNumerics { * \param[in] val_normal - Normal vector, the norm of the vector is the area of the face. * \param[in] val_scale - Scale of the projection. * \param[out] val_Proj_Jac_tensor - Pointer to the projected inviscid Jacobian. + * \param[in] val_tke - Turbulent kinetic energy contained in the total energy (SST), held fixed. */ void GetInviscidProjJac(const su2double *val_velocity, const su2double *val_enthalphy, const su2double *val_chi, const su2double *val_kappa, const su2double *val_normal, su2double val_scale, - su2double **val_Proj_Jac_tensor) const; + su2double **val_Proj_Jac_tensor, su2double val_tke = 0.0) const; /*! * \brief Mapping between primitives variables P and conservatives variables C. @@ -1253,7 +1255,7 @@ class CNumerics { void GetPMatrix(const su2double *val_density, const su2double *val_velocity, const su2double *val_soundspeed, const su2double *val_enthalpy, const su2double *val_chi, const su2double *val_kappa, - const su2double *val_normal, su2double **val_p_tensor) const; + const su2double *val_normal, su2double **val_p_tensor, su2double val_tke = 0.0) const; /*! * \brief Computation of the matrix P, this matrix diagonalize the conservative Jacobians in @@ -1267,12 +1269,13 @@ class CNumerics { template void GetPMatrix(const su2double& density, const su2double* velocity, const su2double& soundspeed, const su2double* normal, - Matrix& p_tensor) const { + Matrix& p_tensor, su2double tke = 0.0) const { + /*--- With SST the total energy contains k (tke): it is added to the kinetic energy of the energy row. ---*/ const su2double rhooc = density / soundspeed; const su2double rhoxc = density * soundspeed; if (nDim == 2) { - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity) + tke; const su2double projvel = GeometryToolbox::DotProduct(2, velocity, normal); p_tensor[0][0] = 1.0; @@ -1295,7 +1298,7 @@ class CNumerics { p_tensor[3][2] = 0.5 * (ke * rhooc + density * projvel + rhoxc / Gamma_Minus_One); p_tensor[3][3] = 0.5 * (ke * rhooc - density * projvel + rhoxc / Gamma_Minus_One); } else { - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity) + tke; const su2double projvel = GeometryToolbox::DotProduct(3, velocity, normal); p_tensor[0][0] = normal[0]; @@ -1430,7 +1433,7 @@ class CNumerics { void GetPMatrix_inv(su2double **val_invp_tensor, const su2double *val_density, const su2double *val_velocity, const su2double *val_soundspeed, const su2double *val_chi, const su2double *val_kappa, - const su2double *val_normal) const; + const su2double *val_normal, su2double val_tke = 0.0) const; /*! * \brief Computation of the matrix P^{-1}, this matrix diagonalize the conservative Jacobians @@ -1444,7 +1447,8 @@ class CNumerics { template void GetPMatrix_inv(const su2double& density, const su2double* velocity, const su2double& soundspeed, const su2double* normal, - Matrix& inv_p_tensor) const { + Matrix& inv_p_tensor, su2double tke = 0.0) const { + /*--- With SST the total energy contains k (tke), dp/drho = (gamma-1) (|u|^2/2 - k). ---*/ const su2double rhoxc = density * soundspeed; const su2double c2 = pow(soundspeed, 2); const su2double gm1 = Gamma_Minus_One; @@ -1455,7 +1459,7 @@ class CNumerics { if (nDim == 3) { const su2double k2orho = normal[2] / density; - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(3, velocity) - tke; const su2double projvel_o_rho = GeometryToolbox::DotProduct(3, velocity, normal) / density; inv_p_tensor[0][0] = normal[0] + k1orho * velocity[2] - k2orho * velocity[1] - normal[0] * gm1_o_c2 * ke; @@ -1488,7 +1492,7 @@ class CNumerics { inv_p_tensor[4][3] = -k2orho - gm1_o_rhoxc * velocity[2]; inv_p_tensor[4][4] = gm1_o_rhoxc; } else { - const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity); + const su2double ke = 0.5 * GeometryToolbox::SquaredNorm(2, velocity) - tke; const su2double projvel_o_rho = GeometryToolbox::DotProduct(2, velocity, normal) / density; inv_p_tensor[0][0] = 1 - gm1_o_c2 * ke; diff --git a/SU2_CFD/include/numerics/flow/convection/roe.hpp b/SU2_CFD/include/numerics/flow/convection/roe.hpp index 01c0874f8e98..8664c60e49e1 100644 --- a/SU2_CFD/include/numerics/flow/convection/roe.hpp +++ b/SU2_CFD/include/numerics/flow/convection/roe.hpp @@ -45,6 +45,7 @@ class CUpwRoeBase_Flow : public CNumerics { su2double *ProjFlux_j = nullptr, *Conservatives_j = nullptr; su2double **P_Tensor = nullptr, **invP_Tensor = nullptr; su2double RoeDensity, RoeEnthalpy, RoeSoundSpeed, ProjVelocity, RoeSoundSpeed2, kappa; + su2double RoeTke = 0.0; /*!< \brief Roe-averaged k (SST, contained in the total energy). */ su2double* Flux = nullptr; /*!< \brief The flux accross the face. */ su2double** Jacobian_i = nullptr; /*!< \brief The Jacobian w.r.t. point i after computation. */ diff --git a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp index faaf56912351..193b5ec7624f 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/centered.hpp @@ -52,6 +52,7 @@ class CCenteredBase : public Base { const su2double fixFactor; const bool dynamicGrid; const su2double stretchParam = 0.3; + const CVariable* const tkeVars; /*!< \brief SST turbulence variables (k is part of the total energy), else nullptr. */ /*! * \brief Constructor, store some constants and forward args to base. @@ -60,7 +61,8 @@ class CCenteredBase : public Base { CCenteredBase(const CConfig& config, Ts&... args) : Base(config, args...), gamma(config.GetGamma()), fixFactor(config.GetCent_Jac_Fix_Factor()), - dynamicGrid(config.GetDynamic_Grid()) { + dynamicGrid(config.GetDynamic_Grid()), + tkeVars(config.GetKind_Turb_Model() == TURB_MODEL::SST ? findTurbVars(args...) : nullptr) { } /*! @@ -140,8 +142,14 @@ class CCenteredBase : public Base { MatrixDbl jac_i, jac_j; if (implicit) { - jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, 0.5); - jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, 0.5); + /*--- With SST the total energy contains k, held fixed in the Jacobian. ---*/ + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } + jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, 0.5, tke_i); + jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, 0.5, tke_j); } /*--- Grid motion. ---*/ @@ -283,6 +291,7 @@ class CJSTmatScheme : public CCenteredBase,Decorator> { using Base::nVar; using Base::gamma; using Base::fixFactor; + using Base::tkeVars; const su2double kappa2; const su2double kappa4; const su2double entropyFix; @@ -350,11 +359,24 @@ class CJSTmatScheme : public CCenteredBase,Decorator> { const auto unitProjVel = dot(avgV.velocity(), unitNormal); + /*--- With SST the total energy contains k, averaged between the cells. ---*/ + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } + const Double avgTke = 0.5 * (tke_i + tke_j); + + /*--- Delta(rho k) = avg(rho) Delta k + avg(k) Delta rho. The part with Delta k is not a pressure jump: remove it + from the 2nd order dissipation before the projection and advect it with the contact wave. The 4th order term + keeps it (the Laplacian of k is not available). ---*/ + const Double rhoDeltaTke = eps2 * avgV.density() * (tke_i - tke_j); + scalarDissip(nVar-1) -= rhoDeltaTke; auto pMat = pMatrix(gamma, avgV.density(), avgV.velocity(), - unitProjVel, avgV.speedSound(), unitNormal); + unitProjVel, avgV.speedSound(), unitNormal, avgTke); auto pMatInv = pMatrixInv(gamma, avgV.density(), avgV.velocity(), - unitProjVel, avgV.speedSound(), unitNormal); + unitProjVel, avgV.speedSound(), unitNormal, avgTke); /*--- Compute limited absolute eigenvalues (times area). ---*/ @@ -394,6 +416,7 @@ class CJSTmatScheme : public CCenteredBase,Decorator> { } } } + flux(nVar-1) += lambda(0) * rhoDeltaTke; } }; diff --git a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp index 800167d4b0f5..2d09e62d6b3b 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/common.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/common.hpp @@ -134,9 +134,11 @@ FORCEINLINE CPair reconstructPrimitives(const Int& iEdge, */ template FORCEINLINE MatrixDbl pMatrix(const Double& gamma, const Double& density, const RandomAccessIterator& velocity, - const Double& projVel, const Double& speedSound, const VectorDbl& normal) { + const Double& projVel, const Double& speedSound, const VectorDbl& normal, + const Double& tke = Double(0.0)) { MatrixDbl pMat; - const Double vel2 = 0.5*squaredNorm(velocity); + /*--- With SST the total energy contains k (tke). ---*/ + const Double vel2 = 0.5*squaredNorm(velocity) + tke; if (nDim == 2) { pMat(0,0) = 1.0; @@ -197,11 +199,13 @@ FORCEINLINE MatrixDbl pMatrix(const Double& gamma, const Double& density template FORCEINLINE MatrixDbl pMatrixInv(const Double& gamma, const Double& density, const RandomAccessIterator& velocity, const Double& projVel, - const Double& speedSound, const VectorDbl& normal) { + const Double& speedSound, const VectorDbl& normal, + const Double& tke = Double(0.0)) { MatrixDbl pMatInv; const Double c2 = pow(speedSound,2); - const Double vel2 = 0.5*squaredNorm(velocity); + /*--- With SST the total energy contains k (tke), dp/drho = (gamma-1) (|u|^2/2 - k). ---*/ + const Double vel2 = 0.5*squaredNorm(velocity) - tke; const Double oneOnRho = 1 / density; if (nDim == 2) { @@ -273,19 +277,30 @@ FORCEINLINE VectorDbl inviscidProjFlux(const PrimVarType& V, return flux; } +/*! + * \brief Turbulence variables among the constructor arguments of a scheme (nullptr if there are none). + */ +inline const CVariable* findTurbVars() { return nullptr; } +template +const CVariable* findTurbVars(T& first, Ts&... rest) { + if constexpr (std::is_convertible::value) return first; + else return findTurbVars(rest...); +} + /*! * \brief Jacobian of the convective flux (compressible flow, ideal gas). + * \note With SST the total energy contains k (tke), which is held fixed. */ template FORCEINLINE MatrixDbl inviscidProjJac(const Double& gamma, RandomAccessIterator velocity, const Double& energy, const VectorDbl& normal, - const Double& scale) { + const Double& scale, const Double& tke = Double(0.0)) { MatrixDbl jac; Double projVel = dot(velocity, normal); Double gamma_m_1 = gamma-1; - Double phi = 0.5*gamma_m_1*squaredNorm(velocity); - Double a1 = gamma*energy - phi; + Double phi = gamma_m_1*(0.5*squaredNorm(velocity) - tke); // dp/drho, k held fixed + Double a1 = gamma*energy - gamma_m_1*(0.5*squaredNorm(velocity) + tke); // total enthalpy jac(0,0) = 0.0; for (size_t iDim = 0; iDim < nDim; ++iDim) { diff --git a/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp b/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp index 38ec02c902b5..f8f0e43e983a 100644 --- a/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp +++ b/SU2_CFD/include/numerics_simd/flow/convection/upwind.hpp @@ -61,6 +61,7 @@ class CUpwindBase : public Base { const bool muscl; const su2double umusclKappa; const LIMITER typeLimiter; + const CVariable* const tkeVars; /*!< \brief SST turbulence variables (k is part of the total energy), else nullptr. */ /*! * \brief Constructor, store some constants and forward args to base. @@ -73,7 +74,8 @@ class CUpwindBase : public Base { dynamicGrid(config.GetDynamic_Grid()), muscl(finestGrid && config.GetMUSCL_Flow()), umusclKappa(config.GetMUSCL_Kappa_Flow()), - typeLimiter(config.GetKind_SlopeLimit_Flow()) { + typeLimiter(config.GetKind_SlopeLimit_Flow()), + tkeVars(config.GetKind_Turb_Model() == TURB_MODEL::SST ? findTurbVars(args...) : nullptr) { } public: @@ -178,6 +180,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { using Base::gamma; using Base::gasConst; using Base::dynamicGrid; + using Base::tkeVars; public: /*! @@ -190,6 +193,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { typeDissip(static_cast(config.GetKind_RoeLowDiss())) { } + /*! * \brief Updates flux and Jacobians with standard Roe dissipation. * \note "Ts" is here just in case other schemes in the family need extra args. @@ -210,9 +214,14 @@ class CRoeScheme : public CUpwindBase, Decorator> { const CEulerVariable& solution, const CGeometry& geometry, Ts&...) const { - /*--- Roe averaged variables. ---*/ + /*--- Roe averaged variables. With SST the total enthalpy of the cells contains k. ---*/ - auto roeAvg = roeAveragedVariables(gamma, V, unitNormal); + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } + auto roeAvg = roeAveragedVariables(gamma, V, unitNormal, tke_i, tke_j); /*--- Grid motion. ---*/ @@ -249,8 +258,8 @@ class CRoeScheme : public CUpwindBase, Decorator> { flux(iVar) = 0.5 * (flux_i(iVar) + flux_j(iVar)); } if (implicit) { - jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, kappa); - jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, kappa); + jac_i = inviscidProjJac(gamma, V.i.velocity(), U.i.energy(), normal, kappa, tke_i); + jac_j = inviscidProjJac(gamma, V.j.velocity(), U.j.energy(), normal, kappa, tke_j); } /*--- Correct for grid motion. ---*/ @@ -270,12 +279,12 @@ class CRoeScheme : public CUpwindBase, Decorator> { /*--- P tensor. ---*/ auto pMat = pMatrix(gamma, roeAvg.density, roeAvg.velocity, - roeAvg.projVel, roeAvg.speedSound, unitNormal); + roeAvg.projVel, roeAvg.speedSound, unitNormal, roeAvg.tke); /*--- Inverse P tensor. ---*/ auto pMatInv = pMatrixInv(gamma, roeAvg.density, roeAvg.velocity, - roeAvg.projVel, roeAvg.speedSound, unitNormal); + roeAvg.projVel, roeAvg.speedSound, unitNormal, roeAvg.tke); /*--- Diference between conservative variables at jPoint and iPoint. ---*/ @@ -284,6 +293,11 @@ class CRoeScheme : public CUpwindBase, Decorator> { deltaU(iVar) = U.j.all(iVar) - U.i.all(iVar); } + /*--- With SST rho*E contains rho*k, and Delta(rho k) = rho_roe Delta k + k_roe Delta rho. The part with Delta k + is not a pressure jump: remove it before the projection and advect it with the contact wave. ---*/ + const Double rhoDeltaTke = roeAvg.density * (tke_j - tke_i); + deltaU(nVar-1) -= rhoDeltaTke; + /*--- Dissipation terms. ---*/ Double dissipation = roeDissipation(iPoint, jPoint, typeDissip, solution); @@ -309,6 +323,7 @@ class CRoeScheme : public CUpwindBase, Decorator> { } } } + flux(nVar-1) -= (1-kappa) * area * dissipation * lambda(0) * rhoDeltaTke; } }; @@ -330,6 +345,7 @@ class CMSWScheme : public CUpwindBase, Decorator> { using Base::gamma; using Base::gasConst; using Base::dynamicGrid; + using Base::tkeVars; public: /*! @@ -397,10 +413,16 @@ class CMSWScheme : public CUpwindBase, Decorator> { lambda(nDim) = fmax(projVel_i + soundSpeed_i, 0); lambda(nDim+1) = fmax(projVel_i - soundSpeed_i, 0); + /*--- With SST the total energy contains k, taken from the cells. ---*/ + Double tke_i = 0.0, tke_j = 0.0; + if (tkeVars) { + tke_i = gatherVariables(iPoint, tkeVars->GetSolution()); + tke_j = gatherVariables(jPoint, tkeVars->GetSolution()); + } auto pMat = pMatrix(gamma, Vweighted.i.density(), Vweighted.i.velocity(), - projVel_i, soundSpeed_i, unitNormal); + projVel_i, soundSpeed_i, unitNormal, tke_i); auto pMatInv = pMatrixInv(gamma, Vweighted.i.density(), Vweighted.i.velocity(), - projVel_i, soundSpeed_i, unitNormal); + projVel_i, soundSpeed_i, unitNormal, tke_i); auto updateFlux = [&](const auto& u, auto& jac) { for (size_t iVar = 0; iVar < nVar; ++iVar) { @@ -429,9 +451,9 @@ class CMSWScheme : public CUpwindBase, Decorator> { lambda(nDim+1) = fmin(projVel_j - soundSpeed_j, 0); pMat = pMatrix(gamma, Vweighted.j.density(), Vweighted.j.velocity(), - projVel_j, soundSpeed_j, unitNormal); + projVel_j, soundSpeed_j, unitNormal, tke_j); pMatInv = pMatrixInv(gamma, Vweighted.j.density(), Vweighted.j.velocity(), - projVel_j, soundSpeed_j, unitNormal); + projVel_j, soundSpeed_j, unitNormal, tke_j); updateFlux(U.j.all, jac_j); } }; diff --git a/SU2_CFD/include/numerics_simd/flow/variables.hpp b/SU2_CFD/include/numerics_simd/flow/variables.hpp index 4406a06abf7c..abeabf049ed1 100644 --- a/SU2_CFD/include/numerics_simd/flow/variables.hpp +++ b/SU2_CFD/include/numerics_simd/flow/variables.hpp @@ -107,15 +107,19 @@ struct CRoeVariables { Double enthalpy; Double speedSound; Double projVel; + Double tke; /*!< \brief Roe-averaged k (SST, part of the total enthalpy). */ }; /*! * \brief Compute Roe-averaged variables from pair of primitive variables. + * \note With SST the total enthalpy contains k (tke_i, tke_j), which is not part of the speed of sound. */ template FORCEINLINE CRoeVariables roeAveragedVariables(const Double& gamma, const CPair& V, - const VectorDbl& normal) { + const VectorDbl& normal, + const Double& tke_i = Double(0.0), + const Double& tke_j = Double(0.0)) { CRoeVariables roeAvg; Double R = sqrt(V.j.density() / V.i.density()); Double D = 1 / (R+1); @@ -124,7 +128,8 @@ FORCEINLINE CRoeVariables roeAveragedVariables(const Double& gamma, roeAvg.velocity(iDim) = (R*V.j.velocity(iDim) + V.i.velocity(iDim)) * D; } roeAvg.enthalpy = (R*V.j.enthalpy() + V.i.enthalpy()) * D; - roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity))); + roeAvg.tke = (R*tke_j + tke_i) * D; + roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity) - roeAvg.tke)); roeAvg.projVel = dot(roeAvg.velocity, normal); return roeAvg; } diff --git a/SU2_CFD/include/solvers/CEulerSolver.hpp b/SU2_CFD/include/solvers/CEulerSolver.hpp index 345be6d3d3a4..ed4b5db05748 100644 --- a/SU2_CFD/include/solvers/CEulerSolver.hpp +++ b/SU2_CFD/include/solvers/CEulerSolver.hpp @@ -419,11 +419,13 @@ class CEulerSolver : public CFVMFlowSolverBase CUpwAUSMPLUS_SLAU_Base_Flow::ComputeResidual(const CCo AD::SetPreaccIn(Normal, nDim); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy /*--- Variables for the general form and primitives for mass flux and pressure calculation. ---*/ /*--- F_{1/2} = ||A|| ( 0.5 * mdot * (psi_i+psi_j) - 0.5 * |mdot| * (psi_i-psi_j) + N * pf ) ---*/ @@ -413,8 +416,9 @@ void CUpwAUSMPLUSUP_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su /*--- Compute interface speed of sound (aF) ---*/ - su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_i); - su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_j); + /*--- With SST the total enthalpy contains k, which is not part of the critical speed of sound. ---*/ + su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_i-turb_ke_i)); + su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_j-turb_ke_j)); su2double ahatL = astarL*astarL/max(astarL, ProjVelocity_i); su2double ahatR = astarR*astarR/max(astarR,-ProjVelocity_j); @@ -584,7 +588,7 @@ void CUpwAUSMPLUSUP_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su astar_b *= 2.0*tmp; Vn_i_b -= tmp*tmp * aF_b; } - H_i_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*Enthalpy_i)) * astar_b; + H_i_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*(Enthalpy_i-turb_ke_i))) * astar_b; H_j_b = 0.0; } else { @@ -593,7 +597,7 @@ void CUpwAUSMPLUSUP_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su astar_b *= 2.0*tmp; Vn_j_b += tmp*tmp * aF_b; } - H_j_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*Enthalpy_j)) * astar_b; + H_j_b = sqrt(0.5*(Gamma-1.0)/((Gamma+1.0)*(Enthalpy_j-turb_ke_j))) * astar_b; H_i_b = 0.0; } @@ -639,8 +643,9 @@ void CUpwAUSMPLUSUP2_Flow::ComputeMassAndPressureFluxes(const CConfig* config, s /*--- Compute interface speed of sound (aF) ---*/ - su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_i); - su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*Enthalpy_j); + /*--- With SST the total enthalpy contains k, which is not part of the critical speed of sound. ---*/ + su2double astarL = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_i-turb_ke_i)); + su2double astarR = sqrt(2.0*(Gamma-1.0)/(Gamma+1.0)*(Enthalpy_j-turb_ke_j)); su2double ahatL = astarL*astarL/max(astarL, ProjVelocity_i); su2double ahatR = astarR*astarR/max(astarR,-ProjVelocity_j); @@ -731,11 +736,12 @@ void CUpwSLAU_Flow::ComputeMassAndPressureFluxes(const CConfig* config, su2doubl sq_velj += Velocity_j[iDim]*Velocity_j[iDim]; } + /*--- With SST the total energy contains k, which is not part of the speed of sound. ---*/ su2double Energy_i = Enthalpy_i - Pressure_i/Density_i; - SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_veli))); + SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_veli-turb_ke_i))); su2double Energy_j = Enthalpy_j - Pressure_j/Density_j; - SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_velj))); + SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_velj-turb_ke_j))); /*--- Compute interface speed of sound (aF), and left/right Mach number ---*/ @@ -842,6 +848,7 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) AD::SetPreaccIn(Normal, nDim); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy /*--- Face area (norm or the normal vector) ---*/ Area = GeometryToolbox::Norm(nDim, Normal); @@ -860,7 +867,7 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) Density_i = V_i[nDim+2]; Enthalpy_i = V_i[nDim+3]; Energy_i = Enthalpy_i - Pressure_i/Density_i; - SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_vel))); + SoundSpeed_i = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_i-0.5*sq_vel-turb_ke_i))); // SST: E contains k /*--- Primitive variables at point j ---*/ sq_vel = 0.0; @@ -872,7 +879,7 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) Density_j = V_j[nDim+2]; Enthalpy_j = V_j[nDim+3]; Energy_j = Enthalpy_j - Pressure_j/Density_j; - SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_vel))); + SoundSpeed_j = sqrt(fabs(Gamma*Gamma_Minus_One*(Energy_j-0.5*sq_vel-turb_ke_j))); // SST: E contains k /*--- Projected velocities ---*/ ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; @@ -925,10 +932,11 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim]; } RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1); - RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel))); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel-(R*turb_ke_j+turb_ke_i)/(R+1)))); /*--- Compute P and Lambda (do it with the Normal) ---*/ - GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor); + GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor, (R*turb_ke_j+turb_ke_i)/(R+1)); ProjVelocity = 0.0; ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -944,11 +952,11 @@ CNumerics::ResidualType<> CUpwAUSM_Flow::ComputeResidual(const CConfig* config) Lambda[nVar-1] = ProjVelocity - RoeSoundSpeed; /*--- Compute inverse P ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, (R*turb_ke_j+turb_ke_i)/(R+1)); /*--- Jacobias of the inviscid flux, scale = 0.5 because val_residual ~ 0.5*(fc_i+fc_j)*Normal ---*/ - GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i); - GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i, turb_ke_i); + GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j, turb_ke_j); /*--- Roe's Flux approximation ---*/ for (iVar = 0; iVar < nVar; iVar++) { diff --git a/SU2_CFD/src/numerics/flow/convection/fvs.cpp b/SU2_CFD/src/numerics/flow/convection/fvs.cpp index 7655b793c851..83684cb34b7a 100644 --- a/SU2_CFD/src/numerics/flow/convection/fvs.cpp +++ b/SU2_CFD/src/numerics/flow/convection/fvs.cpp @@ -45,6 +45,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim + 4); AD::SetPreaccIn(V_j, nDim + 4); + AD::SetPreaccIn(turb_ke_i); // SST: k in the total energy + AD::SetPreaccIn(turb_ke_j); AD::SetPreaccIn(Sensor_i, Sensor_j); AD::SetPreaccIn(Normal, nDim); if (dynamic_grid) { @@ -71,8 +73,9 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { /*--- Recompute the speed of sound because it is not MUSCL-reconstructed. ---*/ const su2double sqvel_i = GeometryToolbox::SquaredNorm(nDim, V_i + 1); const su2double sqvel_j = GeometryToolbox::SquaredNorm(nDim, V_j + 1); - const su2double c_i = sqrt(fmax((Gamma - 1) * (H_i - 0.5 * sqvel_i), EPS)); - const su2double c_j = sqrt(fmax((Gamma - 1) * (H_j - 0.5 * sqvel_j), EPS)); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + const su2double c_i = sqrt(fmax((Gamma - 1) * (H_i - 0.5 * sqvel_i - turb_ke_i), EPS)); + const su2double c_j = sqrt(fmax((Gamma - 1) * (H_j - 0.5 * sqvel_j - turb_ke_j), EPS)); /*--- Recompute conservatives ---*/ @@ -109,6 +112,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { } Vst_i[nDim + 4] = onemw * c_i + w * c_j; Vst_j[nDim + 4] = onemw * c_j + w * c_i; + const su2double tkest_i = onemw * turb_ke_i + w * turb_ke_j; + const su2double tkest_j = onemw * turb_ke_j + w * turb_ke_i; su2double Velst_i[MAXNDIM] = {}, Velst_j[MAXNDIM] = {}; su2double ProjVelst_i{}, ProjVelst_j{}; @@ -132,8 +137,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { /*--- Compute projected P, invP, and Lambda ---*/ su2double P_Tensor[MAXNVAR][MAXNVAR], invP_Tensor[MAXNVAR][MAXNVAR]; - GetPMatrix(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, P_Tensor); - GetPMatrix_inv(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, invP_Tensor); + GetPMatrix(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, P_Tensor, tkest_i); + GetPMatrix_inv(Vst_i[nDim + 2], Velst_i, Vst_i[nDim + 4], UnitNormal, invP_Tensor, tkest_i); /*--- Projected flux (f+) at i ---*/ @@ -167,8 +172,8 @@ CNumerics::ResidualType<> CUpwMSW_Flow::ComputeResidual(const CConfig* config) { /*--- Compute projected P, invP, and Lambda ---*/ - GetPMatrix(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, P_Tensor); - GetPMatrix_inv(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, invP_Tensor); + GetPMatrix(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, P_Tensor, tkest_j); + GetPMatrix_inv(Vst_j[nDim + 2], Velst_j, Vst_j[nDim + 4], UnitNormal, invP_Tensor, tkest_j); /*--- Projected flux (f-) ---*/ diff --git a/SU2_CFD/src/numerics/flow/convection/hllc.cpp b/SU2_CFD/src/numerics/flow/convection/hllc.cpp index 18b210e72381..4fa644208e74 100644 --- a/SU2_CFD/src/numerics/flow/convection/hllc.cpp +++ b/SU2_CFD/src/numerics/flow/convection/hllc.cpp @@ -177,8 +177,9 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) Energy_i = Enthalpy_i - Pressure_i / Density_i; Energy_j = Enthalpy_j - Pressure_j / Density_j; - SoundSpeed_i = sqrt((Enthalpy_i - 0.5 * sq_vel_i) * Gamma_Minus_One); - SoundSpeed_j = sqrt((Enthalpy_j - 0.5 * sq_vel_j) * Gamma_Minus_One); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + SoundSpeed_i = sqrt((Enthalpy_i - 0.5 * sq_vel_i - turb_ke_i) * Gamma_Minus_One); + SoundSpeed_j = sqrt((Enthalpy_j - 0.5 * sq_vel_j - turb_ke_j) * Gamma_Minus_One); /*--- Projected velocities ---*/ @@ -202,7 +203,8 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Roe-averaged speed of sound ---*/ - RoeSoundSpeed = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe ) ); + const su2double RoeTke = ( sqrt(Density_j) * turb_ke_j + sqrt(Density_i) * turb_ke_i ) / Rrho; + RoeSoundSpeed = sqrt( Gamma_Minus_One * ( RoeEnthalpy - 0.5 * sq_velRoe - RoeTke ) ); /*--- Speed of sound at L and R ---*/ @@ -310,7 +312,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) for (jVar = 0; jVar < nVar; jVar++) Jacobian_j[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_i, &Energy_i, UnitNormal, 1.0, Jacobian_i); + GetInviscidProjJac(Velocity_i, &Energy_i, UnitNormal, 1.0, Jacobian_i, turb_ke_i); } else { @@ -326,7 +328,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Computing pressure derivatives d/dU_L (PI) ---*/ - dPI_dU[0] = 0.5 * Gamma_Minus_One * sq_vel_i; + dPI_dU[0] = Gamma_Minus_One * (0.5 * sq_vel_i - turb_ke_i); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Gamma_Minus_One * Velocity_i[iDim]; dPI_dU[nVar-1] = Gamma_Minus_One; @@ -399,7 +401,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Computing d/dU_R (Sm) ---*/ - dSm_dU[0] = ( ProjVelocity_j * ProjVelocity_j - sM * sR - 0.5 * Gamma_Minus_One * sq_vel_j ) / RHO; + dSm_dU[0] = ( ProjVelocity_j * ProjVelocity_j - sM * sR - Gamma_Minus_One * (0.5 * sq_vel_j - turb_ke_j) ) / RHO; for (iDim = 0; iDim < nDim; iDim++) dSm_dU[iDim+1] = - ( UnitNormal[iDim] * ( 2 * ProjVelocity_j - sR - sM) - Gamma_Minus_One * Velocity_j[iDim] ) / RHO; dSm_dU[nVar-1] = - Gamma_Minus_One / RHO; @@ -445,7 +447,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) for (jVar = 0; jVar < nVar; jVar++) Jacobian_i[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_j, &Energy_j, UnitNormal, 1.0, Jacobian_j); + GetInviscidProjJac(Velocity_j, &Energy_j, UnitNormal, 1.0, Jacobian_j, turb_ke_j); } else { @@ -461,7 +463,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Computing d/dU_L (Sm) ---*/ - dSm_dU[0] = ( - ProjVelocity_i * ProjVelocity_i + sM * sL + 0.5 * Gamma_Minus_One * sq_vel_i ) / RHO; + dSm_dU[0] = ( - ProjVelocity_i * ProjVelocity_i + sM * sL + Gamma_Minus_One * (0.5 * sq_vel_i - turb_ke_i) ) / RHO; for (iDim = 0; iDim < nDim; iDim++) dSm_dU[iDim+1] = ( UnitNormal[iDim] * ( 2 * ProjVelocity_i - sL - sM ) - Gamma_Minus_One * Velocity_i[iDim] ) / RHO; dSm_dU[nVar-1] = Gamma_Minus_One / RHO; @@ -504,7 +506,7 @@ CNumerics::ResidualType<> CUpwHLLC_Flow::ComputeResidual(const CConfig* config) /*--- Computing pressure derivatives d/dU_R (PI) ---*/ - dPI_dU[0] = 0.5 * Gamma_Minus_One * sq_vel_j; + dPI_dU[0] = Gamma_Minus_One * (0.5 * sq_vel_j - turb_ke_j); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Gamma_Minus_One * Velocity_j[iDim]; dPI_dU[nVar-1] = Gamma_Minus_One; @@ -695,8 +697,9 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c } Energy_i = Enthalpy_i - Pressure_i / Density_i; - StaticEnthalpy_i = Enthalpy_i - 0.5 * sq_vel_i; - StaticEnergy_i = Energy_i - 0.5 * sq_vel_i; + /*--- With SST the total energy and enthalpy contain k. ---*/ + StaticEnthalpy_i = Enthalpy_i - 0.5 * sq_vel_i - turb_ke_i; + StaticEnergy_i = Energy_i - 0.5 * sq_vel_i - turb_ke_i; Kappa_i = S_i[1] / Density_i; Chi_i = S_i[0] - Kappa_i * StaticEnergy_i; @@ -704,8 +707,8 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c Energy_j = Enthalpy_j - Pressure_j / Density_j; - StaticEnthalpy_j = Enthalpy_j - 0.5 * sq_vel_j; - StaticEnergy_j = Energy_j - 0.5 * sq_vel_j; + StaticEnthalpy_j = Enthalpy_j - 0.5 * sq_vel_j - turb_ke_j; + StaticEnergy_j = Energy_j - 0.5 * sq_vel_j - turb_ke_j; Kappa_j = S_j[1] / Density_j; Chi_j = S_j[0] - Kappa_j * StaticEnergy_j; @@ -747,7 +750,8 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Roe-averaged speed of sound ---*/ //RoeSoundSpeed2 = RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe ); - RoeSoundSpeed = sqrt( RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe ) ); + RoeSoundSpeed = sqrt( RoeChi + RoeKappa * ( RoeEnthalpy - 0.5 * sq_velRoe + - ( sqrt(Density_j) * turb_ke_j + sqrt(Density_i) * turb_ke_i ) / Rrho ) ); /*--- Speed of sound at L and R ---*/ @@ -855,7 +859,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c Jacobian_j[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, 1.0, Jacobian_i); + GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, 1.0, Jacobian_i, turb_ke_i); } else { @@ -871,7 +875,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_L (PI) ---*/ - dPI_dU[0] = Chi_i + 0.5 * Kappa_i * sq_vel_i; + dPI_dU[0] = Chi_i + Kappa_i * (0.5 * sq_vel_i - turb_ke_i); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_i * Velocity_i[iDim]; dPI_dU[nVar-1] = Kappa_i; @@ -945,7 +949,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_R (PI) ---*/ - dPI_dU[0] = Chi_j + 0.5 * Kappa_j * sq_vel_j; + dPI_dU[0] = Chi_j + Kappa_j * (0.5 * sq_vel_j - turb_ke_j); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_j * Velocity_j[iDim]; dPI_dU[nVar-1] = Kappa_j; @@ -999,7 +1003,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c for (jVar = 0; jVar < nVar; jVar++) Jacobian_i[iVar][jVar] = 0; - GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, UnitNormal, 1.0, Jacobian_j); + GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, UnitNormal, 1.0, Jacobian_j, turb_ke_j); } else { @@ -1015,7 +1019,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_L (PI) ---*/ - dPI_dU[0] = Chi_i + 0.5 * Kappa_i * sq_vel_i; + dPI_dU[0] = Chi_i + Kappa_i * (0.5 * sq_vel_i - turb_ke_i); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_i * Velocity_i[iDim]; dPI_dU[nVar-1] = Kappa_i; @@ -1066,7 +1070,7 @@ CNumerics::ResidualType<> CUpwGeneralHLLC_Flow::ComputeResidual(const CConfig* c /*--- Computing pressure derivatives d/dU_R (PI) ---*/ - dPI_dU[0] = Chi_j + 0.5 * Kappa_j * sq_vel_j; + dPI_dU[0] = Chi_j + Kappa_j * (0.5 * sq_vel_j - turb_ke_j); // k (SST) held fixed for (iDim = 0; iDim < nDim; iDim++) dPI_dU[iDim+1] = - Kappa_j * Velocity_j[iDim]; dPI_dU[nVar-1] = Kappa_j; diff --git a/SU2_CFD/src/numerics/flow/convection/roe.cpp b/SU2_CFD/src/numerics/flow/convection/roe.cpp index 37c697515cfb..ad965e01621d 100644 --- a/SU2_CFD/src/numerics/flow/convection/roe.cpp +++ b/SU2_CFD/src/numerics/flow/convection/roe.cpp @@ -101,6 +101,7 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy if (dynamic_grid) { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); } @@ -146,7 +147,9 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim]; } RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1); - RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + RoeTke = (R*turb_ke_j+turb_ke_i)/(R+1); + RoeSoundSpeed2 = (Gamma-1)*(RoeEnthalpy-0.5*sq_vel-RoeTke); /*--- Negative RoeSoundSpeed^2, the jump variables is too large, clear fluxes and exit. ---*/ @@ -170,7 +173,7 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi /*--- P tensor ---*/ - GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor); + GetPMatrix(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, P_Tensor, RoeTke); /*--- Projected velocity adjusted for mesh motion ---*/ @@ -222,8 +225,8 @@ CNumerics::ResidualType<> CUpwRoeBase_Flow::ComputeResidual(const CConfig* confi Flux[iVar] = 0.5*(ProjFlux_i[iVar]+ProjFlux_j[iVar]); if (implicit) { - GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i); - GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i, turb_ke_i); + GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j, turb_ke_j); } /*--- Finalize in children class ---*/ @@ -261,12 +264,17 @@ void CUpwRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_Jac unsigned short iVar, jVar, kVar; /*--- Compute inverse P tensor ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, RoeTke); /*--- Diference between conservative variables at jPoint and iPoint ---*/ for (iVar = 0; iVar < nVar; iVar++) Diff_U[iVar] = Conservatives_j[iVar]-Conservatives_i[iVar]; + /*--- With SST rho*E contains rho*k, and Delta(rho k) = RoeDensity Delta k + RoeTke Delta rho. The part with + Delta k is not a pressure jump: remove it before the projection and advect it with the contact wave. ---*/ + const su2double rhoDeltaTke = RoeDensity*(turb_ke_j-turb_ke_i); + Diff_U[nVar-1] -= rhoDeltaTke; + /*--- Low dissipation formulation ---*/ if (roe_low_dissipation) Dissipation_ij = GetRoe_Dissipation(Dissipation_i, Dissipation_j, Sensor_i, Sensor_j, config); @@ -291,6 +299,7 @@ void CUpwRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_Jac } } } + val_residual[nVar-1] -= (1.0-kappa)*Lambda[0]*rhoDeltaTke*Area*Dissipation_ij; } @@ -350,11 +359,14 @@ void CUpwL2Roe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J for (kVar = 0; kVar < nVar; kVar++) val_residual[iVar] -= (1.0-kappa)*Lambda[kVar]*delta_wave[kVar]*P_Tensor[iVar][kVar]*Area; + /*--- With SST, the jump of k in rho*E is advected with the contact wave. ---*/ + val_residual[nVar-1] -= (1.0-kappa)*Lambda[0]*RoeDensity*(turb_ke_j-turb_ke_i)*Area; + if (!implicit) return; /*--- If implicit use the Jacobians of the standard Roe scheme as an approximation ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, RoeTke); for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { @@ -426,11 +438,14 @@ void CUpwLMRoe_Flow::FinalizeResidual(su2double *val_residual, su2double **val_J for (kVar = 0; kVar < nVar; kVar++) val_residual[iVar] -= (1.0-kappa)*Lambda[kVar]*delta_wave[kVar]*P_Tensor[iVar][kVar]*Area; + /*--- With SST, the jump of k in rho*E is advected with the contact wave. ---*/ + val_residual[nVar-1] -= (1.0-kappa)*Lambda[0]*RoeDensity*(turb_ke_j-turb_ke_i)*Area; + if (!implicit) return; /*--- If implicit use the Jacobians of the standard Roe scheme as an approximation ---*/ - GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor); + GetPMatrix_inv(RoeDensity, RoeVelocity, RoeSoundSpeed, UnitNormal, invP_Tensor, RoeTke); for (iVar = 0; iVar < nVar; iVar++) { for (jVar = 0; jVar < nVar; jVar++) { @@ -568,7 +583,8 @@ CNumerics::ResidualType<> CUpwTurkel_Flow::ComputeResidual(const CConfig* config sq_vel += RoeVelocity[iDim]*RoeVelocity[iDim]; } RoeEnthalpy = (R*Enthalpy_j+Enthalpy_i)/(R+1); - RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel))); + /*--- With SST the total enthalpy contains k, which is not part of the speed of sound. ---*/ + RoeSoundSpeed = sqrt(fabs((Gamma-1)*(RoeEnthalpy-0.5*sq_vel-(R*turb_ke_j+turb_ke_i)/(R+1)))); RoePressure = RoeDensity/Gamma*RoeSoundSpeed*RoeSoundSpeed; /*--- Compute ProjFlux_i ---*/ @@ -635,8 +651,8 @@ CNumerics::ResidualType<> CUpwTurkel_Flow::ComputeResidual(const CConfig* config if (implicit) { /*--- Jacobians of the inviscid flux, scaled by 0.5 because Flux ~ 0.5*(fc_i+fc_j)*Normal ---*/ - GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i); - GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_i, &Energy_i, Normal, 0.5, Jacobian_i, turb_ke_i); + GetInviscidProjJac(Velocity_j, &Energy_j, Normal, 0.5, Jacobian_j, turb_ke_j); } for (iVar = 0; iVar < nVar; iVar ++) { @@ -747,6 +763,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co AD::StartPreacc(); AD::SetPreaccIn(V_i, nDim+4); AD::SetPreaccIn(V_j, nDim+4); AD::SetPreaccIn(Normal, nDim); + AD::SetPreaccIn(turb_ke_i); AD::SetPreaccIn(turb_ke_j); // SST: k in the total energy AD::SetPreaccIn(S_i, 2); AD::SetPreaccIn(S_j, 2); if (dynamic_grid) { AD::SetPreaccIn(GridVel_i, nDim); AD::SetPreaccIn(GridVel_j, nDim); @@ -774,7 +791,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co Density_i = V_i[nDim+2]; Enthalpy_i = V_i[nDim+3]; Energy_i = Enthalpy_i - Pressure_i/Density_i; - StaticEnthalpy_i = Enthalpy_i - 0.5*Velocity2_i; + StaticEnthalpy_i = Enthalpy_i - 0.5*Velocity2_i - turb_ke_i; // SST: H contains k StaticEnergy_i = StaticEnthalpy_i - Pressure_i/Density_i; Kappa_i = S_i[1]/Density_i; @@ -794,7 +811,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co Enthalpy_j = V_j[nDim+3]; Energy_j = Enthalpy_j - Pressure_j/Density_j; - StaticEnthalpy_j = Enthalpy_j - 0.5*Velocity2_j; + StaticEnthalpy_j = Enthalpy_j - 0.5*Velocity2_j - turb_ke_j; // SST: H contains k StaticEnergy_j = StaticEnthalpy_j - Pressure_j/Density_j; Kappa_j = S_j[1]/Density_j; @@ -837,7 +854,8 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co /*--- Compute P and Lambda (do it with the Normal) ---*/ - GetPMatrix(&RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeEnthalpy, &RoeChi, &RoeKappa, UnitNormal, P_Tensor); + GetPMatrix(&RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeEnthalpy, &RoeChi, &RoeKappa, UnitNormal, P_Tensor, + (R*turb_ke_j+turb_ke_i)/(R+1)); ProjVelocity = 0.0; ProjVelocity_i = 0.0; ProjVelocity_j = 0.0; for (iDim = 0; iDim < nDim; iDim++) { @@ -923,6 +941,8 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co for (jVar = 0; jVar < nVar; jVar++) Flux[iVar] -= 0.5*Lambda[jVar]*delta_wave[jVar]*P_Tensor[iVar][jVar]*Area; } + /*--- With SST, the jump of k in rho*E is advected with the contact wave. ---*/ + Flux[nVar-1] -= 0.5*Lambda[0]*RoeDensity*(turb_ke_j-turb_ke_i)*Area; /*--- Flux contribution due to grid motion ---*/ if (dynamic_grid) { @@ -938,20 +958,26 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co /*--- Compute inverse P ---*/ - GetPMatrix_inv(invP_Tensor, &RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeChi , &RoeKappa, UnitNormal); + GetPMatrix_inv(invP_Tensor, &RoeDensity, RoeVelocity, &RoeSoundSpeed, &RoeChi , &RoeKappa, UnitNormal, + (R*turb_ke_j+turb_ke_i)/(R+1)); /*--- Jacobians of the inviscid flux, scaled by 0.5 because val_resconv ~ 0.5*(fc_i+fc_j)*Normal ---*/ - GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, Normal, 0.5, Jacobian_i); + GetInviscidProjJac(Velocity_i, &Enthalpy_i, &Chi_i, &Kappa_i, Normal, 0.5, Jacobian_i, turb_ke_i); - GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, Normal, 0.5, Jacobian_j); + GetInviscidProjJac(Velocity_j, &Enthalpy_j, &Chi_j, &Kappa_j, Normal, 0.5, Jacobian_j, turb_ke_j); /*--- Diference variables iPoint and jPoint ---*/ for (iVar = 0; iVar < nVar; iVar++) Diff_U[iVar] = U_j[iVar]-U_i[iVar]; + /*--- With SST, the part of the jump of rho*E due to the jump of k is not a pressure jump: remove it before the + projection and advect it with the contact wave. ---*/ + const su2double rhoDeltaTke = RoeDensity*(turb_ke_j-turb_ke_i); + Diff_U[nVar-1] -= rhoDeltaTke; + /*--- Roe's Flux approximation ---*/ for (iVar = 0; iVar < nVar; iVar++) { Flux[iVar] = 0.5*(ProjFlux_i[iVar]+ProjFlux_j[iVar]); @@ -968,6 +994,7 @@ CNumerics::ResidualType<> CUpwGeneralRoe_Flow::ComputeResidual(const CConfig* co Jacobian_j[iVar][jVar] -= (1.0-kappa)*Proj_ModJac_Tensor_ij*Area; } } + Flux[nVar-1] -= (1.0-kappa)*Lambda[0]*rhoDeltaTke*Area; /*--- Jacobian contributions due to grid motion ---*/ if (dynamic_grid) { @@ -1026,6 +1053,6 @@ void CUpwGeneralRoe_Flow::ComputeRoeAverage() { // // } - RoeSoundSpeed2 = RoeChi + RoeKappa*(RoeEnthalpy-0.5*sq_vel); + RoeSoundSpeed2 = RoeChi + RoeKappa*(RoeEnthalpy-0.5*sq_vel-(R*turb_ke_j+turb_ke_i)/(R+1)); // SST: H contains k } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 6b8784f228aa..e556d0e4aecc 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -1843,6 +1843,9 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain const bool muscl = (config->GetMUSCL_Flow() && (iMesh == MESH_0)); const bool limiter = (config->GetKind_SlopeLimit_Flow() != LIMITER::NONE); + /*--- SST: the total enthalpy of the cells contains k. ---*/ + const bool tkeInEnergy = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + const CVariable* turbNodes = tkeInEnergy ? solver_container[TURB_SOL]->GetNodes() : nullptr; const bool van_albada = (config->GetKind_SlopeLimit_Flow() == LIMITER::VAN_ALBADA_EDGE); const su2double kappa = config->GetMUSCL_Kappa_Flow(); @@ -1883,6 +1886,13 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain numerics->SetNormal(geometry->edges->GetNormal(iEdge)); + su2double tke_i = 0.0, tke_j = 0.0; + if (tkeInEnergy) { + tke_i = turbNodes->GetSolution(iPoint, 0); + tke_j = turbNodes->GetSolution(jPoint, 0); + numerics->SetTurbKineticEnergy(tke_i, tke_j); + } + auto Coord_i = geometry->nodes->GetCoord(iPoint); auto Coord_j = geometry->nodes->GetCoord(jPoint); @@ -1958,17 +1968,18 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain if (!neg_pres_or_rho_i && !neg_pres_or_rho_j) { - /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. ---*/ + /*--- Recompute the reconstructed quantities in a thermodynamically consistent way. + With SST the total enthalpy contains k, taken from the cells (as the vectorized reconstruction does). ---*/ if (!ideal_gas || low_mach_corr) { - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i); - ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_i, Secondary_i, tke_i); + ComputeConsistentExtrapolation(GetFluidModel(), nDim, Primitive_j, Secondary_j, tke_j); } /*--- Low-Mach number correction. ---*/ if (low_mach_corr) { - LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j); + LowMachPrimitiveCorrection(GetFluidModel(), nDim, Primitive_i, Primitive_j, tke_i, tke_j); } /*--- The correction changes the face states, check them again. ---*/ @@ -1988,7 +1999,14 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } su2double RoeEnthalpy = (R * Primitive_j[prim_idx.Enthalpy()] + Primitive_i[prim_idx.Enthalpy()]) / (R+1); - const bool neg_sound_speed = ((Gamma-1)*(RoeEnthalpy-0.5*sq_vel) < 0.0); + const su2double RoeTke = (R * tke_j + tke_i) / (R+1); // SST: the total enthalpy contains k + /*--- Speed of sound of the Roe average and of each side (schemes such as HLLC use both). ---*/ + auto negSoundSpeed2 = [&](const su2double* prim, su2double tke) { + const su2double vel2 = GeometryToolbox::SquaredNorm(nDim, &prim[prim_idx.Velocity()]); + return (Gamma-1)*(prim[prim_idx.Enthalpy()] - 0.5*vel2 - tke) < 0.0; + }; + const bool neg_sound_speed = ((Gamma-1)*(RoeEnthalpy-0.5*sq_vel-RoeTke) < 0.0) || + negSoundSpeed2(Primitive_i, tke_i) || negSoundSpeed2(Primitive_j, tke_j); bool bad_recon = neg_sound_speed || neg_pres_or_rho_i || neg_pres_or_rho_j; bad_recon = nodes->UpdateNonPhysicalEdgeCounter(iEdge, bad_recon); counter_local += bad_recon; @@ -2048,7 +2066,7 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain } void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsigned short nDim, - su2double *primitive, su2double *secondary) { + su2double *primitive, su2double *secondary, su2double tke) { SU2_ZONE_SCOPED const CEulerVariable::CIndices prim_idx(nDim, 0); const su2double density = primitive[prim_idx.Density()]; @@ -2058,7 +2076,7 @@ void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsig fluidModel->SetTDState_Prho(pressure, density); primitive[prim_idx.Temperature()] = fluidModel->GetTemperature(); - primitive[prim_idx.Enthalpy()] = fluidModel->GetStaticEnergy() + pressure / density + 0.5*velocity2; + primitive[prim_idx.Enthalpy()] = fluidModel->GetStaticEnergy() + pressure / density + 0.5*velocity2 + tke; primitive[prim_idx.SoundSpeed()] = fluidModel->GetSoundSpeed(); secondary[0] = fluidModel->GetdPdrho_e(); secondary[1] = fluidModel->GetdPde_rho(); @@ -2066,7 +2084,8 @@ void CEulerSolver::ComputeConsistentExtrapolation(CFluidModel *fluidModel, unsig } void CEulerSolver::LowMachPrimitiveCorrection(CFluidModel *fluidModel, unsigned short nDim, - su2double *primitive_i, su2double *primitive_j) { + su2double *primitive_i, su2double *primitive_j, + su2double tke_i, su2double tke_j) { SU2_ZONE_SCOPED unsigned short iDim; @@ -2097,10 +2116,10 @@ void CEulerSolver::LowMachPrimitiveCorrection(CFluidModel *fluidModel, unsigned } fluidModel->SetEnergy_Prho(primitive_i[nDim+1], primitive_i[nDim+2]); - primitive_i[nDim+3]= fluidModel->GetStaticEnergy() + primitive_i[nDim+1]/primitive_i[nDim+2] + 0.5*velocity2_i; + primitive_i[nDim+3]= fluidModel->GetStaticEnergy() + primitive_i[nDim+1]/primitive_i[nDim+2] + 0.5*velocity2_i + tke_i; fluidModel->SetEnergy_Prho(primitive_j[nDim+1], primitive_j[nDim+2]); - primitive_j[nDim+3]= fluidModel->GetStaticEnergy() + primitive_j[nDim+1]/primitive_j[nDim+2] + 0.5*velocity2_j; + primitive_j[nDim+3]= fluidModel->GetStaticEnergy() + primitive_j[nDim+1]/primitive_j[nDim+2] + 0.5*velocity2_j + tke_j; } diff --git a/UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp b/UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp new file mode 100644 index 000000000000..1a8decc8a0a0 --- /dev/null +++ b/UnitTests/SU2_CFD/numerics/ausm_tke_tests.cpp @@ -0,0 +1,173 @@ +/*! + * \file ausm_tke_tests.cpp + * \brief SLAU, SLAU2 and AUSM with the turbulent kinetic energy of SST in the total energy. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include +#include +#include +#include +#include "../../../SU2_CFD/include/numerics/flow/convection/ausm_slau.hpp" +#include "../../../SU2_CFD/include/numerics/flow/convection/fvs.hpp" + +namespace { + +constexpr su2double gamma = 1.4; + +std::unique_ptr MakeConfig(bool implicit = false) { + std::stringstream options; + options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= SLAU\nMACH_NUMBER= 0.5\n"; + options << "TIME_DISCRE_FLOW= " + << (implicit ? "EULER_IMPLICIT\nUSE_ACCURATE_FLUX_JACOBIANS= YES\n" : "EULER_EXPLICIT\n"); + return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); +} + +/*--- Primitive variables (T, u, v, [w], p, rho, h, c) of an ideal gas whose total energy contains k. ---*/ +void Primitives(unsigned short nDim, su2double rho, const su2double* vel, su2double p, su2double k, su2double* V) { + su2double q2 = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + V[iDim + 1] = vel[iDim]; + q2 += vel[iDim] * vel[iDim]; + } + V[0] = p / (rho * 287.058); + V[nDim + 1] = p; + V[nDim + 2] = rho; + V[nDim + 3] = gamma / (gamma - 1) * p / rho + 0.5 * q2 + k; + V[nDim + 4] = sqrt(gamma * p / rho); +} + +/*--- k changes neither the pressure nor the speed of sound: with the same density, velocity and pressure, the mass + and momentum fluxes must not depend on k, and the energy flux changes only by the advected k. ---*/ +void CheckTkeIndependence(CNumerics& numerics, const CConfig* config, unsigned short nDim, bool uniformTke = false) { + const su2double vel_i[3] = {0.8, 0.1, -0.05}, vel_j[3] = {0.6, -0.1, 0.05}, n[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double n2[2] = {0.6, 0.8}; + su2double Vi[10] = {0.0}, Vj[10] = {0.0}; + numerics.SetNormal(nDim == 2 ? n2 : n); + + su2double flux[2][5] = {{0.0}}; + /*--- MSW mixes the two states, so k must be uniform for the fluxes to be independent of it. ---*/ + const su2double k_i[2] = {0.0, 0.3}, k_j[2] = {0.0, uniformTke ? 0.3 : 0.2}; + for (int c = 0; c < 2; c++) { + Primitives(nDim, 1.2, vel_i, 1.1, k_i[c], Vi); + Primitives(nDim, 0.9, vel_j, 0.8, k_j[c], Vj); + numerics.SetPrimitive(Vi, Vj); + numerics.SetTurbKineticEnergy(k_i[c], k_j[c]); + const auto res = numerics.ComputeResidual(config); + for (unsigned short iVar = 0; iVar < nDim + 2; iVar++) flux[c][iVar] = res.residual[iVar]; + } + CAPTURE(nDim); + for (unsigned short iVar = 0; iVar < nDim + 1; iVar++) CHECK(flux[1][iVar] == Approx(flux[0][iVar]).epsilon(1e-12)); + /*--- Energy: the extra flux is the mass flux times the upwind k (the flow is from i to j here). ---*/ + CHECK(flux[1][nDim + 1] - flux[0][nDim + 1] == Approx(flux[0][0] * k_i[1]).epsilon(1e-10)); +} + +/*--- Largest difference between the (accurate) Jacobians and central differences of the flux with respect to the + conservative variables, k held fixed, relative to the largest entry. ---*/ +su2double JacobianError(CNumerics& numerics, const CConfig* config, su2double k) { + constexpr unsigned short nDim = 2, nVar = 4; + using State = std::array; + const su2double normal[2] = {0.6, 0.8}; + numerics.SetNormal(normal); + auto evaluate = [&](const State& Ui, const State& Uj, su2double* flux, su2double(*jac_i)[nVar], + su2double(*jac_j)[nVar]) { + su2double Vi[10] = {0.0}, Vj[10] = {0.0}; + for (int side = 0; side < 2; side++) { + const State& U = side ? Uj : Ui; + const su2double vel[2] = {U[1] / U[0], U[2] / U[0]}; + const su2double p = (gamma - 1) * (U[3] - 0.5 * (U[1] * U[1] + U[2] * U[2]) / U[0] - U[0] * k); + Primitives(nDim, U[0], vel, p, k, side ? Vj : Vi); + } + numerics.SetPrimitive(Vi, Vj); + numerics.SetTurbKineticEnergy(k, k); + const auto res = numerics.ComputeResidual(config); + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + flux[iVar] = res.residual[iVar]; + if (jac_i) + for (unsigned short jVar = 0; jVar < nVar; jVar++) { + jac_i[iVar][jVar] = res.jacobian_i[iVar][jVar]; + jac_j[iVar][jVar] = res.jacobian_j[iVar][jVar]; + } + } + }; + const State Ui = {1.2, 0.72, 0.24, 1.1 / (gamma - 1) + 0.5 * (0.72 * 0.72 + 0.24 * 0.24) / 1.2 + 1.2 * k}; + const State Uj = {0.9, 0.36, -0.18, 0.8 / (gamma - 1) + 0.5 * (0.36 * 0.36 + 0.18 * 0.18) / 0.9 + 0.9 * k}; + su2double flux[nVar], jac_i[nVar][nVar], jac_j[nVar][nVar]; + evaluate(Ui, Uj, flux, jac_i, jac_j); + su2double maxErr = 0.0, maxJac = 0.0; + for (int side = 0; side < 2; side++) { + for (unsigned short jVar = 0; jVar < nVar; jVar++) { + State Up = side ? Uj : Ui, Um = Up; + const su2double h = 1e-6; + Up[jVar] += h; + Um[jVar] -= h; + su2double fp[nVar], fm[nVar]; + if (side) { + evaluate(Ui, Up, fp, nullptr, nullptr); + evaluate(Ui, Um, fm, nullptr, nullptr); + } else { + evaluate(Up, Uj, fp, nullptr, nullptr); + evaluate(Um, Uj, fm, nullptr, nullptr); + } + for (unsigned short iVar = 0; iVar < nVar; iVar++) { + const su2double jac = side ? jac_j[iVar][jVar] : jac_i[iVar][jVar]; + maxErr = std::max(maxErr, std::abs(jac - (fp[iVar] - fm[iVar]) / (2 * h))); + maxJac = std::max(maxJac, std::abs(jac)); + } + } + } + return maxErr / maxJac; +} + +} // namespace + +TEST_CASE("Accurate SLAU and AUSM+up Jacobians hold k fixed", "[AUSM][SST]") { + auto config = MakeConfig(true); + for (const su2double k : {0.0, 0.2}) { + CAPTURE(k); + CUpwSLAU_Flow slau(2, 4, config.get(), false); + CHECK(JacobianError(slau, config.get(), k) < 1e-4); + CUpwAUSMPLUSUP_Flow ausmup(2, 4, config.get()); + CHECK(JacobianError(ausmup, config.get(), k) < 1e-4); + } +} + +TEST_CASE("SLAU, SLAU2, AUSM, AUSM+up, AUSM+up2 and MSW do not use k in the speed of sound", "[AUSM][SST]") { + auto config = MakeConfig(); + for (const unsigned short nDim : {2, 3}) { + CUpwSLAU_Flow slau(nDim, nDim + 2, config.get(), false); + CheckTkeIndependence(slau, config.get(), nDim); + CUpwSLAU2_Flow slau2(nDim, nDim + 2, config.get(), false); + CheckTkeIndependence(slau2, config.get(), nDim); + CUpwAUSM_Flow ausm(nDim, nDim + 2, config.get()); + CheckTkeIndependence(ausm, config.get(), nDim); + CUpwAUSMPLUSUP_Flow ausmup(nDim, nDim + 2, config.get()); + CheckTkeIndependence(ausmup, config.get(), nDim); + CUpwAUSMPLUSUP2_Flow ausmup2(nDim, nDim + 2, config.get()); + CheckTkeIndependence(ausmup2, config.get(), nDim); + CUpwMSW_Flow msw(nDim, nDim + 2, config.get()); + CheckTkeIndependence(msw, config.get(), nDim, true); + } +} diff --git a/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp index 0a14682612ed..050de8503571 100644 --- a/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp +++ b/UnitTests/SU2_CFD/numerics/hllc_jacobian_tests.cpp @@ -1,7 +1,8 @@ /*! * \file hllc_jacobian_tests.cpp * \brief Finite-difference verification of the HLLC Jacobians (ideal and general gas, 2D and 3D), on fixed and moving - * faces, exact (USE_ACCURATE_FLUX_JACOBIANS) or default. + * faces, with the turbulent kinetic energy of SST in the total energy, exact (USE_ACCURATE_FLUX_JACOBIANS) or + * default. * \version 8.5.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -211,9 +212,11 @@ void CheckAllBranches(unsigned short nDim, bool moving, su2double k, bool accura TEST_CASE("HLLC Jacobians match finite differences", "[HLLC]") { for (const unsigned short nDim : {2, 3}) { for (const bool moving : {false, true}) { - for (const bool accurate : {false, true}) { - CheckAllBranches(nDim, moving, 0.0, accurate); - CheckAllBranches(nDim, moving, 0.0, accurate); + for (const su2double k : {0.0, 0.2}) { + for (const bool accurate : {false, true}) { + CheckAllBranches(nDim, moving, k, accurate); + CheckAllBranches(nDim, moving, k, accurate); + } } } } diff --git a/UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp b/UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp new file mode 100644 index 000000000000..c33e4ae38735 --- /dev/null +++ b/UnitTests/SU2_CFD/numerics/roe_tke_tests.cpp @@ -0,0 +1,176 @@ +/*! + * \file roe_tke_tests.cpp + * \brief Roe eigenvector matrices with the turbulent kinetic energy of SST in the total energy. + * \version 8.5.0 "Harrier" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2026, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include +#include +#include "../../../SU2_CFD/include/numerics/flow/convection/roe.hpp" + +namespace { + +constexpr su2double gamma = 1.4; + +std::unique_ptr MakeConfig(bool implicit = true) { + std::stringstream options; + options << "SOLVER= EULER\nCONV_NUM_METHOD_FLOW= ROE\nENTROPY_FIX_COEFF= 0.0\nROE_KAPPA= 0.5\n"; + options << "TIME_DISCRE_FLOW= " << (implicit ? "EULER_IMPLICIT" : "EULER_EXPLICIT") << "\n"; + return std::make_unique(options, SU2_COMPONENT::SU2_CFD, false); +} + +/*--- Primitive variables (T, u, v, [w], p, rho, h, c) of an ideal gas whose total energy contains k. ---*/ +void Primitives(unsigned short nDim, su2double rho, const su2double* vel, su2double p, su2double k, su2double* V) { + su2double q2 = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + V[iDim + 1] = vel[iDim]; + q2 += vel[iDim] * vel[iDim]; + } + V[0] = p / (rho * 287.058); + V[nDim + 1] = p; + V[nDim + 2] = rho; + V[nDim + 3] = gamma / (gamma - 1) * p / rho + 0.5 * q2 + k; + V[nDim + 4] = sqrt(gamma * p / rho); +} + +/*--- Secondary variables (dp/drho_e, dp/de_rho) of the ideal gas, for the general-gas schemes. ---*/ +void Secondary(su2double rho, su2double p, su2double* S) { + S[0] = p / rho; + S[1] = (gamma - 1) * rho; +} + +/*--- Common constructor for the schemes of the family. ---*/ +template +std::unique_ptr MakeScheme(unsigned short nDim, unsigned short nVar, const CConfig* config) { + return std::make_unique(nDim, nVar, config); +} +template <> +std::unique_ptr MakeScheme(unsigned short nDim, unsigned short nVar, const CConfig* config) { + return std::make_unique(nDim, nVar, config, false); +} + +/*--- Roe flux between two states with the same pressure and normal velocity and different density, k and (with + shear) tangential velocity. With the entropy fix off this is a contact discontinuity, whose exact (upwind) flux is + the flux of the upstream state. ---*/ +template +void CheckContactWithTkeJump(const char* name, unsigned short nDim, const su2double* vel, bool implicit, bool shear) { + auto config = MakeConfig(implicit); + const unsigned short nVar = nDim + 2; + auto numerics = MakeScheme(nDim, nVar, config.get()); + + const su2double p = 1.0, rho_i = 1.4, rho_j = 0.7, k_i = 0.1, k_j = 0.3; + const su2double n2[2] = {0.6, 0.8}, n3[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double* n = (nDim == 2) ? n2 : n3; + const su2double t2[2] = {0.16, -0.12}, t3[3] = {0.15, -0.1, 0.0}; // tangential velocity jumps + const su2double* t = (nDim == 2) ? t2 : t3; + su2double vel_j[3] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) vel_j[iDim] = vel[iDim] + (shear ? t[iDim] : 0.0); + su2double Vi[10] = {0.0}, Vj[10] = {0.0}, Si[2], Sj[2]; + Primitives(nDim, rho_i, vel, p, k_i, Vi); + Primitives(nDim, rho_j, vel_j, p, k_j, Vj); + Secondary(rho_i, p, Si); + Secondary(rho_j, p, Sj); + numerics->SetPrimitive(Vi, Vj); + numerics->SetSecondary(Si, Sj); + numerics->SetNormal(n); + numerics->SetTurbKineticEnergy(k_i, k_j); + const auto res = numerics->ComputeResidual(config.get()); + + su2double un = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) un += vel[iDim] * n[iDim]; + const su2double rho_up = (un >= 0.0) ? rho_i : rho_j; + const su2double* V_up = (un >= 0.0) ? Vi : Vj; + + CAPTURE(name, nDim, un, implicit, shear); + CHECK(res.residual[0] == Approx(rho_up * un).margin(1e-12)); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + CHECK(res.residual[iDim + 1] == Approx(rho_up * V_up[iDim + 1] * un + p * n[iDim]).margin(1e-12)); + CHECK(res.residual[nVar - 1] == Approx(rho_up * V_up[nDim + 3] * un).margin(1e-12)); +} + +template +void CheckContactsWithTkeJump(const char* name, bool implicit = true, bool shear = true) { + const su2double zero[3] = {0.0, 0.0, 0.0}, moving[3] = {0.3, -0.1, 0.2}, reverse[3] = {-0.3, 0.1, -0.2}; + for (const unsigned short nDim : {2, 3}) { + for (const auto* vel : {zero, moving, reverse}) { + CheckContactWithTkeJump(name, nDim, vel, implicit, false); + if (shear) CheckContactWithTkeJump(name, nDim, vel, implicit, true); + } + } +} + +} // namespace + +TEST_CASE("Roe schemes resolve a contact with a jump of k exactly", "[Roe][SST]") { + CheckContactsWithTkeJump("Roe"); + CheckContactsWithTkeJump("general Roe", true); + CheckContactsWithTkeJump("general Roe", false); + CheckContactsWithTkeJump("LMRoe"); + /*--- L2Roe scales the tangential velocity jump at low Mach number (not upwind by design). ---*/ + CheckContactsWithTkeJump("L2Roe", true, false); +} + +TEST_CASE("Roe eigenvector matrices are inverse with k in the total energy", "[Roe][SST]") { + auto config = MakeConfig(); + for (const unsigned short nDim : {2, 3}) { + const unsigned short nVar = nDim + 2; + CUpwRoe_Flow numerics(nDim, nVar, config.get(), false); + const su2double velocity[3] = {0.4, -0.3, 0.2}, normal[3] = {2.0 / 7, 3.0 / 7, 6.0 / 7}; + const su2double n2[3] = {0.6, 0.8, 0.0}; + su2double P[5][5], Pinv[5][5]; + numerics.GetPMatrix(1.1, velocity, 1.2, nDim == 3 ? normal : n2, P, 0.3); + numerics.GetPMatrix_inv(1.1, velocity, 1.2, nDim == 3 ? normal : n2, Pinv, 0.3); + for (unsigned short i = 0; i < nVar; i++) { + for (unsigned short j = 0; j < nVar; j++) { + su2double prod = 0.0; + for (unsigned short l = 0; l < nVar; l++) prod += P[i][l] * Pinv[l][j]; + CAPTURE(nDim, i, j); + CHECK(prod == Approx(i == j ? 1.0 : 0.0).margin(1e-12)); + } + } + } +} + +TEST_CASE("Roe does not dissipate a stationary contact with k in the total energy", "[Roe][SST]") { + /*--- Zero velocity, same pressure and k on both sides, different density: a contact discontinuity at rest. + The Roe flux must be the pressure only (no mass or energy flux) when the entropy fix is off. ---*/ + auto config = MakeConfig(); + constexpr unsigned short nDim = 2, nVar = nDim + 2; + CUpwRoe_Flow numerics(nDim, nVar, config.get(), false); + + const su2double zero[2] = {0.0, 0.0}, p = 1.0, k = 0.2; + su2double Vi[10] = {0.0}, Vj[10] = {0.0}, normal[2] = {0.6, 0.8}; + Primitives(nDim, 1.4, zero, p, k, Vi); + Primitives(nDim, 0.7, zero, p, k, Vj); + numerics.SetPrimitive(Vi, Vj); + numerics.SetNormal(normal); + numerics.SetTurbKineticEnergy(k, k); + const auto res = numerics.ComputeResidual(config.get()); + + CHECK(res.residual[0] == Approx(0.0).margin(1e-12)); + CHECK(res.residual[1] == Approx(p * normal[0]).epsilon(1e-12)); + CHECK(res.residual[2] == Approx(p * normal[1]).epsilon(1e-12)); + CHECK(res.residual[3] == Approx(0.0).margin(1e-12)); +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 8c2e61614e74..7bc1de9e4b5e 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -17,6 +17,8 @@ su2_cfd_tests = files(['Common/CConfig_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', 'SU2_CFD/numerics/roe_contact_tests.cpp', 'SU2_CFD/numerics/hllc_jacobian_tests.cpp', + 'SU2_CFD/numerics/roe_tke_tests.cpp', + 'SU2_CFD/numerics/ausm_tke_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/solvers/CNEMOEulerSolver_tests.cpp', From 3eb9dc89be1cc001d2ef714d4df5149f7539b76b Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 13:54:02 +0200 Subject: [PATCH 2/7] Use the k of the boundary node in the boundary states of the compressible solver With SST the boundary states added the free-stream k to the total energy, while the domain next to the boundary has its own k (the free-stream value decays right away, #1851): the Riemann problem at the boundary then sees a jump of total energy that is not a jump of pressure. Use the k of the boundary node in the far-field, Riemann, turbo Riemann, inlet, outlet, supersonic inlet, engine and actuator-disk states, and pass it to the convective numerics of every boundary. Where the interior state is used, remove k from it: static energy of the interior state in BC_Riemann, BC_TurboRiemann and BC_Giles, energy of the DENSITY_VELOCITY Riemann state, total speed of sound at the inlet, engine exhaust and actuator disk. The P, P^-1 and Jacobian matrices of the Riemann boundaries take k, and the boundary total enthalpy of BC_Giles includes it. The initial check for non-physical points also subtracts the free-stream k from the energy. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/src/solvers/CEulerSolver.cpp | 181 ++++++++++++++++++++++----- 1 file changed, 148 insertions(+), 33 deletions(-) diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index e556d0e4aecc..76e26dafe04c 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -65,6 +65,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND); const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING); const bool adjoint = config->GetContinuous_Adjoint() || config->GetDiscrete_Adjoint(); + const bool tkeNeeded = (rans && config->GetKind_Turb_Model() == TURB_MODEL::SST); int Unst_RestartIter = 0; unsigned long iPoint, iMarker, counter_local = 0, counter_global = 0; @@ -246,6 +247,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Density_Inf = config->GetDensity_FreeStreamND(); Energy_Inf = config->GetEnergy_FreeStreamND(); Mach_Inf = config->GetMach(); + const su2double TKE_Inf = config->GetTke_FreeStreamND(); /*--- Initialize the secondary values for direct derivative approximations ---*/ @@ -315,6 +317,7 @@ CEulerSolver::CEulerSolver(CGeometry *geometry, CConfig *config, Velocity2 += pow(nodes->GetSolution(iPoint,iDim+1)/Density,2); StaticEnergy= nodes->GetEnergy(iPoint) - 0.5*Velocity2; + if (tkeNeeded) StaticEnergy -= TKE_Inf; GetFluidModel()->SetTDState_rhoe(Density, StaticEnergy); Pressure= GetFluidModel()->GetPressure(); @@ -5050,7 +5053,9 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, bool implicit = config->GetKind_TimeIntScheme() == EULER_IMPLICIT; bool viscous = config->GetViscous(); - bool tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST; + bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); auto *Normal = new su2double[nDim]; @@ -5200,7 +5205,8 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, } Pressure = Density*SoundSpeed*SoundSpeed/Gamma; Energy = Pressure/(Gamma_Minus_One*Density) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Store new primitive state for computing the flux. ---*/ @@ -5217,6 +5223,16 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container, conv_numerics->SetPrimitive(V_domain, V_infty); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + + conv_numerics->SetTurbKineticEnergy(tke, tke); + + } + if (dynamic_grid) { conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -5297,6 +5313,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, tkeNeeded = config->GetKind_Turb_Model() == TURB_MODEL::SST, ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || config->GetKind_FluidModel() == IDEAL_GAS; + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double **P_Tensor = new su2double*[nVar], **invP_Tensor = new su2double*[nVar]; @@ -5342,7 +5360,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Density_i = nodes->GetDensity(iPoint); const auto Energy_i = nodes->GetEnergy(iPoint); - const su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + su2double StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + if (tkeNeeded) StaticEnergy_i -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5396,7 +5415,6 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); break; case STATIC_SUPERSONIC_INFLOW_PT: @@ -5422,7 +5440,6 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); break; case STATIC_SUPERSONIC_INFLOW_PD: @@ -5449,7 +5466,6 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); break; case DENSITY_VELOCITY: @@ -5465,7 +5481,9 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, for (auto iDim = 0u; iDim < nDim; iDim++) Velocity_e[iDim] = VelMag_e*Flow_Dir[iDim]; + /*--- The total energy of the node already contains k, which is added below to every case. ---*/ Energy_e = Energy_i; + if (tkeNeeded) Energy_e -= turbNodes->GetSolution(iPoint, 0); break; case STATIC_PRESSURE: @@ -5489,11 +5507,18 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, break; } + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). + It is also removed from the static energy of the boundary state. ---*/ + const su2double Tke_e = tkeNeeded ? turbNodes->GetSolution(iPoint, 0) : 0.0; + Energy_e += Tke_e; + /*--- Compute P (matrix of right eigenvectors) ---*/ - conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor, + Tke_e); /*--- Compute inverse P (matrix of left eigenvectors)---*/ - conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal, + Tke_e); /*--- eigenvalues contribution due to grid motion ---*/ su2double ProjVelocity_i = GeometryToolbox::DotProduct(nDim, Velocity_i, UnitNormal); @@ -5542,7 +5567,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, const auto Velocity2_b = GeometryToolbox::SquaredNorm(nDim, Velocity_b); const su2double Energy_b = u_b[nVar-1]/Density_b; - const su2double StaticEnergy_b = Energy_b - 0.5*Velocity2_b; + const su2double StaticEnergy_b = Energy_b - 0.5*Velocity2_b - Tke_e; GetFluidModel()->SetTDState_rhoe(Density_b, StaticEnergy_b); /*--- Store number of Newton iterations at BC ---*/ @@ -5596,7 +5621,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, } /*--- Compute flux Jacobian in state b ---*/ - conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b); + conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b, Tke_e); /*--- Jacobian contribution due to grid motion ---*/ if (dynamic_grid){ @@ -5686,7 +5711,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container, /*--- Turbulent kinetic energy ---*/ - if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + if (tkeNeeded) visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0), solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0)); @@ -5742,6 +5767,8 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); const bool ideal_gas = config->GetKind_FluidModel() == STANDARD_AIR || config->GetKind_FluidModel() == IDEAL_GAS; + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); su2double *Normal, *turboNormal, *UnitNormal, *FlowDirMix, FlowDirMixMag, *turboVelocity; Normal = new su2double[nDim]; @@ -5815,6 +5842,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Energy_i = nodes->GetEnergy(iPoint); StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + if (tkeNeeded) StaticEnergy_i -= turbNodes->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -5878,7 +5906,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; case MIXING_IN: @@ -5913,7 +5941,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Density_e = GetFluidModel()->GetDensity(); StaticEnergy_e = GetFluidModel()->GetStaticEnergy(); Energy_e = StaticEnergy_e + 0.5 * Velocity2_e; - // if (tkeNeeded) Energy_e += GetTke_Inf(); + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; @@ -5931,6 +5959,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; case STATIC_PRESSURE: @@ -5948,6 +5977,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; @@ -5965,6 +5995,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_e += Velocity_e[iDim]*Velocity_e[iDim]; } Energy_e = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_e; + if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0); break; default: @@ -5973,10 +6004,13 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain } /*--- Compute P (matrix of right eigenvectors) ---*/ - conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor); + const su2double Tke_i = tkeNeeded ? turbNodes->GetSolution(iPoint, 0) : 0.0; // SST: k in the total energy + conv_numerics->GetPMatrix(&Density_i, Velocity_i, &SoundSpeed_i, &Enthalpy_i, &Chi_i, &Kappa_i, UnitNormal, P_Tensor, + Tke_i); /*--- Compute inverse P (matrix of left eigenvectors)---*/ - conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal); + conv_numerics->GetPMatrix_inv(invP_Tensor, &Density_i, Velocity_i, &SoundSpeed_i, &Chi_i, &Kappa_i, UnitNormal, + Tke_i); /*--- eigenvalues contribution due to grid motion ---*/ if (dynamic_grid){ @@ -6040,7 +6074,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain Velocity2_b += Velocity_b[iDim]*Velocity_b[iDim]; } Energy_b = u_b[nVar-1]/Density_b; - StaticEnergy_b = Energy_b - 0.5*Velocity2_b; + StaticEnergy_b = Energy_b - 0.5*Velocity2_b - Tke_i; // SST: k in the total energy GetFluidModel()->SetTDState_rhoe(Density_b, StaticEnergy_b); Pressure_b = GetFluidModel()->GetPressure(); Temperature_b = GetFluidModel()->GetTemperature(); @@ -6096,7 +6130,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain } /*--- Compute flux Jacobian in state b ---*/ - conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b); + conv_numerics->GetInviscidProjJac(Velocity_b, &Enthalpy_b, &Chi_b, &Kappa_b, Normal, 1.0, Jacobian_b, Tke_i); /*--- Jacobian contribution due to grid motion ---*/ if (dynamic_grid) @@ -6798,6 +6832,9 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Energy_i = nodes->GetEnergy(iPoint); StaticEnergy_i = Energy_i - 0.5*Velocity2_i; + /*--- SST: the total energy contains k. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + StaticEnergy_i -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); GetFluidModel()->SetTDState_rhoe(Density_i, StaticEnergy_i); @@ -7054,6 +7091,9 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu Energy_b = GetFluidModel()->GetStaticEnergy() + 0.5*Velocity2_b; Temperature_b= GetFluidModel()->GetTemperature(); Enthalpy_b = Energy_b + Pressure_b/Density_b; + /*--- SST: the boundary state carries the k of the node in its total enthalpy, as the interior one. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) + Enthalpy_b += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); /*--- Primitive variables, using the derived quantities ---*/ V_boundary[0] = Temperature_b; @@ -7071,6 +7111,16 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container, CNu /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_boundary); + + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + + conv_numerics->SetTurbKineticEnergy(tke, tke); + + } conv_numerics->SetSecondary(S_domain, S_boundary); @@ -7294,6 +7344,8 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + /*--- SST: the total enthalpy of the domain contains k, the imposed total enthalpy does not. ---*/ + if (tkeNeeded) Energy -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); Pressure = V_domain[nDim+1]; H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; SoundSpeed2 = Gamma*Pressure/Density; @@ -7364,7 +7416,6 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Using pressure, density, & velocity, compute the energy ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ @@ -7375,6 +7426,9 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[nDim+2] = Density; V_inlet[nDim+3] = Energy + Pressure/Density; + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) V_inlet[nDim+3] += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + break; } /*--- Mass flow has been specified at the inlet. ---*/ @@ -7436,7 +7490,6 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, /*--- Energy for the fictitious inlet state ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Vel_Mag*Vel_Mag; - if (tkeNeeded) Energy += GetTke_Inf(); /*--- Primitive variables, using the derived quantities ---*/ Temperature = Pressure / ( Gas_Constant * Density); @@ -7447,6 +7500,9 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, V_inlet[nDim+2] = Density; V_inlet[nDim+3] = Energy + Pressure/Density; + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) V_inlet[nDim+3] += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + break; } default: @@ -7458,6 +7514,16 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + + conv_numerics->SetTurbKineticEnergy(tke, tke); + + } + if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -7495,6 +7561,9 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, bool gravity = (config->GetGravityForce()); bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST); + CVariable* turbNodes = nullptr; + if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes(); + auto *Normal = new su2double[nDim]; /*--- Loop over all the vertices on this boundary marker ---*/ @@ -7576,7 +7645,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = P_Exit/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0); /*--- Conservative variables, using the derived quantities ---*/ V_outlet[0] = Pressure / ( Gas_Constant * Density); @@ -7590,6 +7659,11 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container, /*--- Set various quantities in the solver class ---*/ conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), geometry->nodes->GetGridVel(iPoint)); @@ -7687,13 +7761,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con /*--- Compute the energy from the specified state. ---*/ const su2double Velocity2 = GeometryToolbox::SquaredNorm(int(MAXNDIM), Velocity); - su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; - if (tkeNeeded) { - const su2double* Turb_Properties = config->GetInlet_TurbVal(Marker_Tag); - const su2double Intensity = Turb_Properties[0]; - const su2double Tke = 3.0 / 2.0 * (Velocity2 * pow(Intensity, 2)); - Energy += Tke; - } + const su2double Energy = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2; /*--- Primitive variables, using the derived quantities. ---*/ @@ -7705,6 +7773,9 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con for (unsigned short iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+prim_idx.Velocity()] = Velocity[iDim]; + /*--- k of the boundary node, not the free-stream value, which decays right away in the domain (#1851). ---*/ + if (tkeNeeded) V_inlet[prim_idx.Enthalpy()] += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); + /*--- Current solution at this boundary node. ---*/ const auto* V_domain = nodes->GetPrimitive(iPoint); @@ -7719,6 +7790,11 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -7800,6 +7876,11 @@ void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_co conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } if (dynamic_grid) conv_numerics->SetGridVel(geometry->nodes->GetGridVel(iPoint), @@ -8016,7 +8097,7 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai } Energy = Inflow_Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Conservative variables, using the derived quantities ---*/ @@ -8032,6 +8113,11 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_inflow); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } /*--- Set grid movement ---*/ @@ -8181,6 +8267,8 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + /*--- SST: the total enthalpy of the domain contains k, the imposed total enthalpy does not. ---*/ + if (tkeNeeded) Energy -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); Pressure = V_domain[nDim+1]; H_Exhaust = (Gamma*Gas_Constant/Gamma_Minus_One)*Exhaust_Temperature; SoundSpeed2 = Gamma*Pressure/Density; @@ -8253,7 +8341,7 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta /*--- Using pressure, density, & velocity, compute the energy ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Primitive variables, using the derived quantities ---*/ @@ -8284,6 +8372,11 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta conv_numerics->SetNormal(Normal); conv_numerics->SetPrimitive(V_domain, V_exhaust); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } /*--- Set grid movement ---*/ @@ -8545,7 +8638,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Conservative variables, using the derived quantities ---*/ @@ -8557,6 +8650,11 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C V_inlet[nDim+3] = Energy + Pressure/Density; V_inlet[nDim+4] = SoundSpeed; conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } @@ -8632,6 +8730,8 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = V_domain[nDim+3] - V_domain[nDim+1]/V_domain[nDim+2]; + /*--- SST: the total enthalpy of the domain contains k, the imposed total enthalpy does not. ---*/ + if (tkeNeeded) Energy -= solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); Pressure = V_domain[nDim+1]; H_Total = (Gamma*Gas_Constant/Gamma_Minus_One)*T_Total; SoundSpeed2 = Gamma*Pressure/Density; @@ -8700,7 +8800,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C /*--- Using pressure, density, & velocity, compute the energy ---*/ Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Primitive variables, using the derived quantities ---*/ @@ -8712,6 +8812,11 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C V_outlet[nDim+3] = Energy + Pressure/Density; V_outlet[nDim+4] = sqrt(SoundSpeed2); conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } @@ -8948,7 +9053,7 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver Velocity2 += Velocity[iDim]*Velocity[iDim]; } Energy = Pressure/(Density*Gamma_Minus_One) + 0.5*Velocity2; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 /*--- Conservative variables, using the derived quantities ---*/ @@ -8959,6 +9064,11 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver V_inlet[nDim+3] = Energy + Pressure/Density; V_inlet[nDim+4] = SoundSpeed; conv_numerics->SetPrimitive(V_domain, V_inlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } else { /*--- Acoustic Riemann invariant extrapolation form the interior domain. ---*/ @@ -8996,7 +9106,7 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver /*--- Computation of the enthalpy, total energy, temperature and speed of sound. ---*/ H_out = H_in/Density_in + Fa/Density_out; Energy = H_out - Pressure_out/Density_out; - if (tkeNeeded) Energy += GetTke_Inf(); + if (tkeNeeded) Energy += solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0); // node k, see #1851 Temperature_out = (Energy-0.5*Velocity2/(pow(Density_out,2)))*(Gamma_Minus_One/Gas_Constant); SoS_out = sqrt(Gamma*Gas_Constant*Temperature_out); @@ -9010,6 +9120,11 @@ void CEulerSolver::BC_ActDisk_VariableLoad(CGeometry *geometry, CSolver **solver V_outlet[nDim+3] = H_out; V_outlet[nDim+4] = SoS_out; conv_numerics->SetPrimitive(V_domain, V_outlet); + /*--- SST: both states contain the k of the node in their total enthalpy. ---*/ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST) { + const su2double tke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint, 0); + conv_numerics->SetTurbKineticEnergy(tke, tke); + } } /*--- Grid Movement (NOT TESTED!)---*/ From 8536636377cbf7237192a84a5eef5b87c87a0266 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 13:54:51 +0200 Subject: [PATCH 3/7] Include the SST k in the pressure of the axisymmetric source, remove a wrong viscous term The axisymmetric source computed the pressure from the total energy without removing k (SST). Subtract it in the pressure and in the energy row of the Jacobian (k held fixed). The viscous part added mu_tot*v*2/3*rho*k/y to the energy equation for every RANS model (with the SA variable in place of k): this term does not come from the stress tensor of the flow solver, remove it. k is set only for SST. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- SU2_CFD/include/numerics/flow/flow_sources.hpp | 3 ++- SU2_CFD/src/numerics/flow/flow_sources.cpp | 11 ++++++----- SU2_CFD/src/solvers/CEulerSolver.cpp | 2 +- 3 files changed, 9 insertions(+), 7 deletions(-) diff --git a/SU2_CFD/include/numerics/flow/flow_sources.hpp b/SU2_CFD/include/numerics/flow/flow_sources.hpp index 83c3811326cd..ca0613a0ebfc 100644 --- a/SU2_CFD/include/numerics/flow/flow_sources.hpp +++ b/SU2_CFD/include/numerics/flow/flow_sources.hpp @@ -72,7 +72,8 @@ class CSourceBase_Flow : public CNumerics { */ class CSourceAxisymmetric_Flow : public CSourceBase_Flow { protected: - bool implicit, viscous, rans; + bool implicit, viscous; + bool tkeInEnergy; /*!< \brief SST: k is part of the total energy (turb_ke_i is set by the solver). */ su2double yinv{0.0}; /*! diff --git a/SU2_CFD/src/numerics/flow/flow_sources.cpp b/SU2_CFD/src/numerics/flow/flow_sources.cpp index 87784e73f39f..44e5cf579dcf 100644 --- a/SU2_CFD/src/numerics/flow/flow_sources.cpp +++ b/SU2_CFD/src/numerics/flow/flow_sources.cpp @@ -54,7 +54,7 @@ CSourceAxisymmetric_Flow::CSourceAxisymmetric_Flow(unsigned short val_nDim, unsi implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT); viscous = config->GetViscous(); - rans = (config->GetKind_Turb_Model() != TURB_MODEL::NONE); + tkeInEnergy = (config->GetKind_Turb_Model() == TURB_MODEL::SST); } @@ -74,7 +74,9 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi Velocity_i = U_i[iDim+1] / U_i[0]; sq_vel += Velocity_i * Velocity_i; } - Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel); + /*--- With SST the total energy contains k. ---*/ + const su2double tke = tkeInEnergy ? turb_ke_i : 0.0; + Pressure_i = Gamma_Minus_One*U_i[0]*(U_i[nDim+1]/U_i[0]-0.5*sq_vel-tke); Enthalpy_i = (U_i[nDim+1] + Pressure_i) / U_i[0]; /*--- Smooth blending between gradient formulation and standard formulation ---*/ @@ -140,6 +142,7 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi jacobian[3][1] = -(Gamma-1)*U_i[2]*U_i[1]/(U_i[0]*U_i[0]); jacobian[3][2] = Gamma*U_i[3]/U_i[0] - 0.5*(Gamma-1)*( (U_i[1]*U_i[1]+U_i[2]*U_i[2])/(U_i[0]*U_i[0]) + 2*U_i[2]*U_i[2]/(U_i[0]*U_i[0]) ); jacobian[3][3] = Gamma*U_i[2]/U_i[0]; + jacobian[3][2] -= Gamma_Minus_One*tke; // k is a variable of the turbulence solver for (iVar=0; iVar < nVar; iVar++) for (jVar=0; jVar < nVar; jVar++) @@ -161,8 +164,6 @@ CNumerics::ResidualType<> CSourceAxisymmetric_Flow::ComputeResidual(const CConfi void CSourceAxisymmetric_Flow::ResidualDiffusion(){ - if (!rans){ turb_ke_i = 0.0; } - su2double laminar_viscosity_i = V_i[nDim+5]; su2double eddy_viscosity_i = V_i[nDim+6]; su2double thermal_conductivity_i = V_i[nDim+7]; @@ -181,7 +182,7 @@ void CSourceAxisymmetric_Flow::ResidualDiffusion(){ -TWO3*AuxVar_Grad_i[0][1]); residual[3] -= Volume*(yinv*(total_viscosity_i*(u*(PrimVar_Grad_i[2][0]+PrimVar_Grad_i[1][1]) +v*TWO3*(2*PrimVar_Grad_i[2][1]-PrimVar_Grad_i[1][0] - -v*yinv+U_i[0]*turb_ke_i)) + -v*yinv)) +total_conductivity_i*PrimVar_Grad_i[0][1]) -TWO3*(AuxVar_Grad_i[1][1]+AuxVar_Grad_i[2][0])); } diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 76e26dafe04c..8d490cb2c9d1 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2239,7 +2239,7 @@ void CEulerSolver::Source_Residual(CGeometry *geometry, CSolver **solver_contain numerics->SetAuxVarGrad(nodes->GetAuxVarGradient(iPoint), nullptr); /*--- Set turbulence kinetic energy ---*/ - if (rans){ + if (config->GetKind_Turb_Model() == TURB_MODEL::SST){ CVariable* turbNodes = solver_container[TURB_SOL]->GetNodes(); numerics->SetTurbKineticEnergy(turbNodes->GetSolution(iPoint,0), turbNodes->GetSolution(iPoint,0)); } From 580800a8ec34eafd88fb43d57eca214b56dd1d3a Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 13:56:04 +0200 Subject: [PATCH 4/7] Apply the SST compressibility corrections together with the production modifiers COMPRESSIBILITY-WILCOX and COMPRESSIBILITY-SARKAR were stored as the production option, overwriting VORTICITY, KATO-LAUNDER or UQ when combined with them (and their switch cases used the strain production). They are already stored in compWilcox and compSarkar: use these flags for the correction, independently of the production modifier. Use beta_star in place of the literal 0.09 in the corrected dissipation of omega. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- Common/include/option_structure.hpp | 7 ++--- Common/src/CConfig.cpp | 8 ++--- .../numerics/turbulent/turb_sources.hpp | 30 ++++++++----------- 3 files changed, 16 insertions(+), 29 deletions(-) diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bcd2e9c06328..b8222da53784 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -1200,13 +1200,10 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne SSTParsedOptions.production = SST_OPTIONS::UQ; } - // Parse compressibility options + // Parse compressibility options, stored in compWilcox and compSarkar so that they can be combined with a + // production modifier if (sst_compWilcox && sst_compSarkar) { SU2_MPI::Error("Please select only one compressibility correction (COMPRESSIBILITY-WILCOX or COMPRESSIBILITY-SARKAR).", CURRENT_FUNCTION); - } else if (sst_compWilcox) { - SSTParsedOptions.production = SST_OPTIONS::COMP_Wilcox; - } else if (sst_compSarkar) { - SSTParsedOptions.production = SST_OPTIONS::COMP_Sarkar; } SSTParsedOptions.sust = sst_sust; diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 3a5bfee8f713..a4c17da1b4d9 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -6713,16 +6713,12 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) { cout << "\nperturbing the Reynold's Stress Matrix towards " << eig_val_comp << " component turbulence"; if (uq_permute) cout << " (permuting eigenvectors)"; break; - case SST_OPTIONS::COMP_Wilcox: - cout << " with compressibility correction of Wilcox"; - break; - case SST_OPTIONS::COMP_Sarkar: - cout << " with compressibility correction of Sarkar"; - break; default: cout << " with no production modification"; break; } + if (sstParsedOptions.compWilcox) cout << ", with compressibility correction of Wilcox"; + if (sstParsedOptions.compSarkar) cout << ", with compressibility correction of Sarkar"; if (sstParsedOptions.dll){ cout << "\nusing non dimensional lower limits relative to infinity values clipping by Coefficients:" ; diff --git a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp index 9f72b6e37cc4..d1cfd2f09054 100644 --- a/SU2_CFD/include/numerics/turbulent/turb_sources.hpp +++ b/SU2_CFD/include/numerics/turbulent/turb_sources.hpp @@ -938,33 +938,27 @@ class CSourcePieceWise_TurbSST final : public CNumerics { P_Base = sqrt(StrainMag_i*VorticityMag); break; - case SST_OPTIONS::COMP_Wilcox: - P_Base = StrainMag_i; - if (Mt >= 0.25) { - zetaFMt = 2.0 * (Mt * Mt - 0.25 * 0.25); - } - break; - - case SST_OPTIONS::COMP_Sarkar: - P_Base = StrainMag_i; - if (Mt >= 0.25) { - zetaFMt = 0.5 * (Mt * Mt); - } - break; - default: /*--- Base production term for SST-1994 and SST-2003 ---*/ P_Base = StrainMag_i; break; } + /*--- Compressibility corrections, independent of the production modifier. ---*/ + if (sstParsedOptions.compWilcox && Mt >= 0.25) { + zetaFMt = 2.0 * (Mt * Mt - 0.25 * 0.25); + } + if (sstParsedOptions.compSarkar && Mt >= 0.25) { + zetaFMt = 0.5 * (Mt * Mt); + } + /*--- Production limiter. ---*/ const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0]; su2double P = Eddy_Viscosity_i * pow(P_Base, 2); su2double pk = max(0.0, min(P, prod_limit)); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; + const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i; const su2double zeta = max(ScalarVar_i[1], eddy_visc_var * F2_i / a1); /*--- Production limiter only for V2003, recompute for V1994. ---*/ @@ -988,7 +982,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { pw = max(pw, sust_w); } - if (sstParsedOptions.production == SST_OPTIONS::COMP_Sarkar) { + if (sstParsedOptions.compSarkar) { const su2double Dilatation_Sarkar = -0.15 * pk * Mt + 0.2 * beta_star * (1.0 +zetaFMt) * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * Mt * Mt; pk += Dilatation_Sarkar; } @@ -996,7 +990,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { /*--- Dissipation ---*/ su2double dk = beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0] * (1.0 + zetaFMt); - su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - 0.09/beta_blended * zetaFMt); + su2double dw = beta_blended * Density_i * ScalarVar_i[1] * ScalarVar_i[1] * (1.0 - beta_star/beta_blended * zetaFMt); /*--- LM model coupling with production and dissipation term for k transport equation---*/ if (config->GetKind_Trans_Model() == TURB_TRANS_MODEL::LM) { @@ -1025,7 +1019,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics { Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume * (1.0 + zetaFMt); Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume * (1.0 + zetaFMt); Jacobian_i[1][0] = 0.0; - Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - 0.09/beta_blended * zetaFMt); + Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume * (1.0 - beta_star/beta_blended * zetaFMt); } AD::SetPreaccOut(Residual, nVar); From a002e777053715f19f9e00e41f49db05dd527999 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 13:56:23 +0200 Subject: [PATCH 5/7] Check the SST free-stream values and markers, fix the omega of the turbo inlet - Stop with an error if the free-stream k or omega of SST is not positive and finite (zero turbulence intensity, viscosity ratio or velocity), which gave mu_t = 0/0 in the whole field. - Reject SST with engine and actuator-disk markers: the turbulence solver has no boundary condition for them (checked after the turbulence model of Euler zones is cleared). - BC_Inlet_Turbo computed omega from the averaged k of the flow solver instead of the imposed k. - Remove a dead write of the wall node solution in the SST wall function. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- Common/src/CConfig.cpp | 8 ++++++++ SU2_CFD/src/solvers/CTurbSSTSolver.cpp | 16 ++++++++++++---- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a4c17da1b4d9..852909cff588 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -5011,6 +5011,14 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i Kind_Solver == MAIN_SOLVER::FEM_EULER) Kind_Turb_Model = TURB_MODEL::NONE; + /*--- SST has no engine or actuator-disk boundary conditions: the faces would get no turbulence flux at all. + Checked after the turbulence model of Euler zones is cleared (multizone). ---*/ + if (Kind_Turb_Model == TURB_MODEL::SST && + (nMarker_EngineInflow + nMarker_EngineExhaust + nMarker_ActDiskInlet + nMarker_ActDiskOutlet) > 0) { + SU2_MPI::Error("MARKER_ENGINE_INFLOW, MARKER_ENGINE_EXHAUST and MARKER_ACTDISK are not supported with the SST model.", + CURRENT_FUNCTION); + } + Kappa_2nd_Flow = jst_coeff[0]; Kappa_4th_Flow = jst_coeff[1]; Kappa_2nd_AdjFlow = jst_adj_coeff[0]; diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index d6f98ecdee6b..1377defbaefc 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -124,6 +124,15 @@ CTurbSSTSolver::CTurbSSTSolver(CGeometry *geometry, CConfig *config, const CSolv su2double kine_Inf = 3.0/2.0*(VelMag2*Intensity*Intensity); su2double omega_Inf = rhoInf*kine_Inf/(muLamInf*viscRatio); + /*--- A zero turbulence intensity, viscosity ratio or free-stream velocity gives k = 0 or omega = 0 (or infinite), + and then mu_t = 0/0 in the whole field. ---*/ + if (!(kine_Inf > 0.0 && omega_Inf > 0.0 && std::isfinite(SU2_TYPE::GetValue(kine_Inf)) && + std::isfinite(SU2_TYPE::GetValue(omega_Inf)))) { + SU2_MPI::Error("SST: the free-stream k and omega must be positive and finite.\n" + "Set FREESTREAM_TURBULENCEINTENSITY > 0 and FREESTREAM_TURB2LAMVISCRATIO > 0 with a nonzero\n" + "free-stream velocity.", CURRENT_FUNCTION); + } + Solution_Inf[0] = kine_Inf; Solution_Inf[1] = omega_Inf; @@ -251,7 +260,7 @@ void CTurbSSTSolver::Postprocessing(CGeometry *geometry, CSolver **solver_contai const su2double kine = nodes->GetSolution(iPoint,0); const su2double omega = nodes->GetSolution(iPoint,1); - const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag; + const su2double eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag; const su2double muT = max(0.0, rho * a1 * kine / max(a1 * omega, eddy_visc_var * F2)); nodes->SetmuT(iPoint, muT); @@ -580,7 +589,6 @@ void CTurbSSTSolver::SetTurbVars_WF(CGeometry *geometry, CSolver **solver_contai su2double solution[MAXNVAR] = {k, omega}; nodes->SetSolution_Old(iPoint_Neighbor,solution); - nodes->SetSolution(iPoint,solution); LinSysRes.SetBlock_Zero(iPoint_Neighbor); @@ -777,15 +785,15 @@ void CTurbSSTSolver::BC_Inlet_Turbo(CGeometry *geometry, CSolver **solver_contai su2double rho = flowSolver->GetAverageDensity(val_marker, iSpan); su2double pressure = flowSolver->GetAveragePressure(val_marker, iSpan); - su2double kine = flowSolver->GetAverageKine(val_marker, iSpan); FluidModel->SetTDState_Prho(pressure, rho); su2double muLam = FluidModel->GetLaminarViscosity(); su2double VelMag2 = GeometryToolbox::SquaredNorm(nDim, flowSolver->GetAverageTurboVelocity(val_marker, iSpan)); + /*--- Imposed k from the intensity, omega from the imposed k and the viscosity ratio. ---*/ su2double kine_b = 3.0/2.0*(VelMag2*Intensity*Intensity); - su2double omega_b = rho*kine/(muLam*viscRatio); + su2double omega_b = rho*kine_b/(muLam*viscRatio); const su2double solution_j[] = {kine_b, omega_b}; From 37dcfedd0e07553af04066473fb386fb6af1516c Mon Sep 17 00:00:00 2001 From: rois1995 Date: Wed, 30 Sep 2026 14:07:43 +0200 Subject: [PATCH 6/7] Add a regression test of the Riemann DENSITY_VELOCITY inlet with SST Axisymmetric air nozzle with the inlet as a Riemann DENSITY_VELOCITY boundary, whose boundary state takes the energy of the node (containing k) and must remove k before adding it back with the other cases. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01LEL91DW5WPbPwgFtCvHga6 --- .../air_nozzle_density_velocity.cfg | 100 ++++++++++++++++++ TestCases/parallel_regression.py | 8 ++ TestCases/serial_regression.py | 8 ++ 3 files changed, 116 insertions(+) create mode 100644 TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg diff --git a/TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg b/TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg new file mode 100644 index 000000000000..208ad6b63159 --- /dev/null +++ b/TestCases/axisymmetric_rans/air_nozzle/air_nozzle_density_velocity.cfg @@ -0,0 +1,100 @@ +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% % +% SU2 configuration file % +% Case description: Axisymmetric supersonic converging-diverging air nozzle, % +% Riemann DENSITY_VELOCITY inlet (k in the boundary energy) % +% Author: Florian Dittmann % +% Date: 2021.12.02 % +% File Version 8.5.0 "Harrier" % +% % +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------% +% +SOLVER= RANS +KIND_TURB_MODEL= SST +RESTART_SOL= YES +AXISYMMETRIC= YES + +% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------% +% +MACH_NUMBER= 1E-9 +INIT_OPTION= TD_CONDITIONS +FREESTREAM_OPTION= TEMPERATURE_FS +FREESTREAM_PRESSURE= 1400000 +FREESTREAM_TEMPERATURE= 373.15 +REF_DIMENSIONALIZATION= DIMENSIONAL + +% ---- IDEAL GAS, POLYTROPIC, VAN DER WAALS AND PENG ROBINSON CONSTANTS -------% +% +FLUID_MODEL= STANDARD_AIR + +% --------------------------- VISCOSITY MODEL ---------------------------------% +% +VISCOSITY_MODEL= CONSTANT_VISCOSITY +MU_CONSTANT= 1.716E-5 + +% --------------------------- THERMAL CONDUCTIVITY MODEL ----------------------% +% +CONDUCTIVITY_MODEL= CONSTANT_PRANDTL +PRANDTL_LAM= 0.72 +PRANDTL_TURB= 0.90 + +% -------------------- BOUNDARY CONDITION DEFINITION --------------------------% +% +MARKER_HEATFLUX= ( WALL, 0.0 ) +MARKER_SYM= ( SYMMETRY ) +MARKER_RIEMANN= ( INFLOW, DENSITY_VELOCITY, 13.07, 100.0, 1.0, 0.0, 0.0, \ + OUTFLOW, STATIC_PRESSURE, 100000.0, 0.0, 0.0, 0.0, 0.0 ) +MARKER_MONITORING = (WALL) +MARKER_ANALYZE= (INFLOW) +% +% High free-stream turbulence, so that k is not negligible in the energy of the inlet state +FREESTREAM_TURBULENCEINTENSITY= 0.2 + +% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------% +% +NUM_METHOD_GRAD= GREEN_GAUSS +CFL_NUMBER= 1000.0 +CFL_ADAPT= NO +MAX_DELTA_TIME= 1E6 +OBJECTIVE_FUNCTION= DRAG + +% ----------- SLOPE LIMITER AND DISSIPATION SENSOR DEFINITION -----------------% +% +MUSCL_FLOW= YES +SLOPE_LIMITER_FLOW= NONE + +% ------------------------ LINEAR SOLVER DEFINITION ---------------------------% +% +LINEAR_SOLVER= FGMRES +LINEAR_SOLVER_PREC= ILU +LINEAR_SOLVER_ILU_FILL_IN= 0 +LINEAR_SOLVER_ERROR= 0.01 +LINEAR_SOLVER_ITER= 10 + +% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------% +% +CONV_NUM_METHOD_FLOW= ROE +ENTROPY_FIX_COEFF= 0.1 +TIME_DISCRE_FLOW= EULER_IMPLICIT + +% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------% +% +CONV_NUM_METHOD_TURB= SCALAR_UPWIND +TIME_DISCRE_TURB= EULER_IMPLICIT +CFL_REDUCTION_TURB= 1.0 + +% --------------------------- CONVERGENCE PARAMETERS --------------------------% +% +ITER= 3000 +CONV_RESIDUAL_MINVAL= -12 +CONV_STARTITER= 10 + +% ------------------------- INPUT/OUTPUT INFORMATION --------------------------% +% +MESH_FILENAME= nozzle.su2 +SOLUTION_FILENAME= solution_flow +RESTART_FILENAME= restart_flow +OUTPUT_WRT_FREQ= 1000 +SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_ENERGY, RMS_TKE, RMS_DISSIPATION, SURFACE_STATIC_TEMPERATURE) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index e85b5b7be660..8fff340267e9 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -624,6 +624,14 @@ def main(): axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) + # Riemann DENSITY_VELOCITY inlet with SST: k of the node in the energy of the boundary state + axi_rans_air_nozzle_density_velocity = TestCase('axi_rans_air_nozzle_density_velocity') + axi_rans_air_nozzle_density_velocity.cfg_dir = "axisymmetric_rans/air_nozzle" + axi_rans_air_nozzle_density_velocity.cfg_file = "air_nozzle_density_velocity.cfg" + axi_rans_air_nozzle_density_velocity.test_iter = 10 + axi_rans_air_nozzle_density_velocity.test_vals = [-1.879323, 4.022709, 0.514084, 4.414617, 702.55, 702.55] + test_list.append(axi_rans_air_nozzle_density_velocity) + ################################# ## Compressible RANS Restart ### ################################# diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index c34fe6f47225..77bf40b28fa1 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -372,6 +372,14 @@ def main(): axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) + # Riemann DENSITY_VELOCITY inlet with SST: k of the node in the energy of the boundary state + axi_rans_air_nozzle_density_velocity = TestCase('axi_rans_air_nozzle_density_velocity') + axi_rans_air_nozzle_density_velocity.cfg_dir = "axisymmetric_rans/air_nozzle" + axi_rans_air_nozzle_density_velocity.cfg_file = "air_nozzle_density_velocity.cfg" + axi_rans_air_nozzle_density_velocity.test_iter = 10 + axi_rans_air_nozzle_density_velocity.test_vals = [-1.879849, 4.023872, 0.501157, 4.419530, 704.80, 704.80] + test_list.append(axi_rans_air_nozzle_density_velocity) + # Axisymmetric air nozzle species axi_rans_air_nozzle_species = TestCase('axi_rans_air_nozzle_species') axi_rans_air_nozzle_species.cfg_dir = "axisymmetric_rans/air_nozzle" From bcdd766c017eff4e1bbef7ca86cbaecc2f65a6a2 Mon Sep 17 00:00:00 2001 From: rois1995 Date: Thu, 22 Oct 2026 09:44:45 +0200 Subject: [PATCH 7/7] Update the reference values from the CI Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_013YFdgQwrRqCQ2ihG446aNe --- TestCases/hybrid_regression.py | 46 +++++++++++------------ TestCases/hybrid_regression_AD.py | 2 +- TestCases/parallel_regression.py | 58 ++++++++++++++--------------- TestCases/parallel_regression_AD.py | 10 ++--- TestCases/serial_regression.py | 54 +++++++++++++-------------- TestCases/serial_regression_AD.py | 6 +-- TestCases/tutorials.py | 6 +-- TestCases/vandv.py | 10 ++--- 8 files changed, 96 insertions(+), 96 deletions(-) diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py index 60e5c01fbc7f..fa01e5b51118 100644 --- a/TestCases/hybrid_regression.py +++ b/TestCases/hybrid_regression.py @@ -149,7 +149,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.004080, -7.638116, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.004086, -7.637842, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.004276, -7.636719, -0.000000, 2.089953] test_list.append(poiseuille_profile) @@ -178,7 +178,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.461774, 5.883885, 0.630267, 0.020404, 100.000000] + rae2822_sst.test_vals = [-1.436649, 5.889182, 0.659764, 0.021350, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -194,7 +194,7 @@ def main(): turb_flatplate.cfg_dir = "rans/flatplate" turb_flatplate.cfg_file = "turb_SA_flatplate.cfg" turb_flatplate.test_iter = 20 - turb_flatplate.test_vals = [-0.187373, 0.003723, 10.000000, -1.443902] + turb_flatplate.test_vals = [-0.187373, 0.003723, 10.000000, -1.444398] test_list.append(turb_flatplate) # ONERA M6 Wing @@ -210,7 +210,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.038075, -16.332088, 1.080346, 0.018385, 20.000000, -2.873507, 0.000000, -14.250270, 0.000000] + turb_naca0012_sa.test_vals = [-12.038098, -16.332088, 1.080346, 0.018385, 20.000000, -2.873012, 0.000000, -14.250270, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.038091, -16.332090, 1.080346, 0.018385, 20.000000, -2.873236, 0.000000, -14.250271, 0.000000] test_list.append(turb_naca0012_sa) @@ -219,7 +219,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.094026, -15.250728, -5.906323, 1.070413, 0.015775, -2.854334, 0.000000] + turb_naca0012_sst.test_vals = [-5.948963, -10.295280, -3.740016, 1.069295, 0.015905, -3.049395, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.075928, -15.246732, -5.861249, 1.070036, 0.015841, -2.835263, 0] test_list.append(turb_naca0012_sst) @@ -228,7 +228,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.080919, -14.837175, -5.732906, 1.000893, 0.019109, -2.119961] + turb_naca0012_sst_sust.test_vals = [-9.951719, -11.674782, -5.092201, 1.000888, 0.019126, -1.639020] turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073210, -14.836724, -5.732627, 1.000050, 0.019144, -2.629689] test_list.append(turb_naca0012_sst_sust) @@ -245,7 +245,7 @@ def main(): turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012" turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg" turb_naca0012_sst_expliciteuler.test_iter = 10 - turb_naca0012_sst_expliciteuler.test_vals = [-3.532365, -3.157224, 3.743381, 1.124798, 0.501715, -16.000000] + turb_naca0012_sst_expliciteuler.test_vals = [-3.408465, -3.157224, 3.743383, 1.124790, 0.501715, -16.000000] test_list.append(turb_naca0012_sst_expliciteuler) # PROPELLER @@ -265,7 +265,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-5.069624, 0.505288, -2.816368, 1.728267, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-5.066047, 0.510609, -1.717316, 2.920473, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.140441, -9.154674, -10.886121, -5.806594, 0.000000] test_list.append(axi_rans_air_nozzle_restart) @@ -279,7 +279,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.551665, -5.057150, 0.830240, -0.008811, 0.078187] + turb_naca0012_sst_restart_mg.test_vals = [-3.689602, -5.056964, 0.830559, -0.008732, 0.078126] test_list.append(turb_naca0012_sst_restart_mg) ############################# @@ -291,7 +291,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.978722, 1.343530, 0.443909, -0.029125] + turb_naca0012_1c.test_vals = [-4.964592, 1.345541, 0.446118, -0.028439] turb_naca0012_1c.test_vals_aarch64 = [-4.976620, 1.345983, 0.433171, -0.033685] test_list.append(turb_naca0012_1c) @@ -300,7 +300,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.482844, 1.260870, 0.404376, -0.040361] + turb_naca0012_2c.test_vals = [-5.482846, 1.260731, 0.404913, -0.040192] turb_naca0012_2c.test_vals_aarch64 = [-5.485484, 1.263406, 0.411442, -0.040859] test_list.append(turb_naca0012_2c) @@ -309,7 +309,7 @@ def main(): turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" turb_naca0012_3c.test_iter = 10 - turb_naca0012_3c.test_vals = [-5.583730, 1.228732, 0.381967, -0.046233] + turb_naca0012_3c.test_vals = [-5.583734, 1.228655, 0.379717, -0.047003] turb_naca0012_3c.test_vals_aarch64 = [-5.583737, 1.232005, 0.390258, -0.046305] test_list.append(turb_naca0012_3c) @@ -318,7 +318,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.134181, 1.283464, 0.553971, 0.011916] + turb_naca0012_p1c1.test_vals = [-5.165277, 1.279520, 0.549254, 0.010996] turb_naca0012_p1c1.test_vals_aarch64 = [-5.114189, 1.285037, 0.406851, -0.043003] test_list.append(turb_naca0012_p1c1) @@ -327,7 +327,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.553988, 1.234031, 0.424168, -0.033501] + turb_naca0012_p1c2.test_vals = [-5.553200, 1.234059, 0.424695, -0.033351] turb_naca0012_p1c2.test_vals_aarch64 = [-5.548245, 1.236384, 0.381821, -0.050337] test_list.append(turb_naca0012_p1c2) @@ -501,7 +501,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-2.560678, -1.175981, 0.062204, 1.399351, 2.219223, 1.399297, 2.217470, 0.000000] + square_cylinder.test_vals = [-2.560698, -1.175955, 0.062131, 1.399399, 2.220356, 1.399347, 2.218595, 0.000000] square_cylinder.unsteady = True test_list.append(square_cylinder) @@ -602,7 +602,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.959252, -12.213981, -19.314222, -13.600248, -19.076384, -13.444335, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-0.170804, -0.162041, -8.236409, -2.533080, -7.364785, -1.876220, 73290.000000, 73290.000000, 0.020065, 82.291000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -610,7 +610,7 @@ def main(): axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D" axial_stage2D.cfg_file = "Axial_stage2D.cfg" axial_stage2D.test_iter = 20 - axial_stage2D.test_vals = [1.167512, 1.598494, -2.928579, 2.573642, -2.527390, 3.016171, 106370.000000, 106370.000000, 5.726800, 64.383000] + axial_stage2D.test_vals = [1.166038, 1.598534, -2.928564, 2.573650, -2.527328, 3.016255, 106370.000000, 106370.000000, 5.726800, 64.395000] test_list.append(axial_stage2D) # 2D transonic stator restart @@ -618,7 +618,7 @@ def main(): transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg" transonic_stator_restart.test_iter = 20 - transonic_stator_restart.test_vals = [-4.367144, -2.487733, -2.079075, 1.728173, -1.464952, 3.225025, -471620.000000, 94.839000, -0.051124] + transonic_stator_restart.test_vals = [-4.363280, -2.484853, -2.078097, 1.730856, -1.464954, 3.223391, -471620.000000, 94.839000, -0.098969] transonic_stator_restart.test_vals_aarch64 = [-4.442510, -2.561369, -2.165778, 1.652750, -1.355494, 3.172712, -471620.000000, 94.843000, -0.043825] test_list.append(transonic_stator_restart) @@ -640,7 +640,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 5 - uniform_flow.test_vals = [5.000000, 0.000000, -0.193990, -10.624443] + uniform_flow.test_vals = [5.000000, 0.000000, -0.193990, -10.624456] uniform_flow.unsteady = True uniform_flow.multizone = True test_list.append(uniform_flow) @@ -659,7 +659,7 @@ def main(): channel_2D.cfg_dir = "sliding_interface/channel_2D" channel_2D.cfg_file = "channel_2D_WA.cfg" channel_2D.test_iter = 2 - channel_2D.test_vals = [2.000000, 0.000000, 0.399109, 0.350302, 0.405582] + channel_2D.test_vals = [2.000000, 0.000000, 0.399108, 0.350302, 0.405584] channel_2D.unsteady = True channel_2D.multizone = True test_list.append(channel_2D) @@ -701,7 +701,7 @@ def main(): supersonic_vortex_shedding.cfg_dir = "sliding_interface/supersonic_vortex_shedding" supersonic_vortex_shedding.cfg_file = "sup_vor_shed_WA.cfg" supersonic_vortex_shedding.test_iter = 5 - supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.128138, 1.064252] + supersonic_vortex_shedding.test_vals = [5.000000, 0.000000, 1.128138, 1.064251] supersonic_vortex_shedding.unsteady = True supersonic_vortex_shedding.multizone = True test_list.append(supersonic_vortex_shedding) @@ -763,7 +763,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.328791, -4.152991, 0.000000, 103.000000] + dyn_fsi.test_vals = [-4.328791, -4.152991, 0.000000, 102.000000] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) @@ -773,7 +773,7 @@ def main(): fsi_cht_restart.cfg_dir = "fea_fsi/stat_fsi" fsi_cht_restart.cfg_file = "config_restart.cfg" fsi_cht_restart.test_iter = 0 - fsi_cht_restart.test_vals = [5.000000, 0.006352, -1.960362, -9.327033, -9.644570, -9.319618, 608.380000, -0.012974, 0.000000, 20.000000] + fsi_cht_restart.test_vals = [5.000000, 0.006352, -1.960362, -9.327033, -9.619476, -9.319378, 608.380000, -0.012974, 0.000000, 20.000000] fsi_cht_restart.multizone = True test_list.append(fsi_cht_restart) diff --git a/TestCases/hybrid_regression_AD.py b/TestCases/hybrid_regression_AD.py index 6a693812d903..5f168cba1b67 100644 --- a/TestCases/hybrid_regression_AD.py +++ b/TestCases/hybrid_regression_AD.py @@ -87,7 +87,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.201555, -0.175213, 3.045200, -0.041846] + discadj_rans_naca0012_sst.test_vals = [-2.201606, -0.172546, 3.045700, -0.041853] discadj_rans_naca0012_sst.test_vals_aarch64 = [-2.201855, -0.172443, 3.043400, -0.041820] test_list.append(discadj_rans_naca0012_sst) diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py index 8180634dddf9..05eb8db8ec33 100755 --- a/TestCases/parallel_regression.py +++ b/TestCases/parallel_regression.py @@ -336,7 +336,7 @@ def main(): ramp_msw.cfg_dir = "euler/ramp" ramp_msw.cfg_file = "inv_ramp_msw.cfg" ramp_msw.test_iter = 100 - ramp_msw.test_vals = [-7.242600, -1.471141, -0.077507, 0.054419] + ramp_msw.test_vals = [-6.972874, -1.216000, -0.077507, 0.054419] ramp_msw.tol = [0.2, 0.2, 0.00001, 0.00001] test_list.append(ramp_msw) @@ -471,7 +471,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.485957, 5.891491, 0.675200, 0.023427, 100.000000] + rae2822_sst.test_vals = [-1.493359, 5.883418, 0.657627, 0.022989, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -519,7 +519,7 @@ def main(): turb_flatplate_sst_roughBCKnopp.cfg_dir = "rans/flatplate/roughness/bc_knopp" turb_flatplate_sst_roughBCKnopp.cfg_file = "turb_SST_flatplate_roughBCKnopp.cfg" turb_flatplate_sst_roughBCKnopp.test_iter = 10 - turb_flatplate_sst_roughBCKnopp.test_vals = [-5.058634, -2.460850, -2.847064, 0.447200, -2.595042, 1.497149, -0.188079, 0.004571] + turb_flatplate_sst_roughBCKnopp.test_vals = [-4.619122, -2.026481, -2.392910, 0.885845, -2.501638, 1.532115, -0.188182, 0.004571] test_list.append(turb_flatplate_sst_roughBCKnopp) # FLAT PLATE, ROUGHNESS BC AUPOIX SST @@ -527,7 +527,7 @@ def main(): turb_flatplate_sst_roughBCAupoix.cfg_dir = "rans/flatplate/roughness/bc_aupoix" turb_flatplate_sst_roughBCAupoix.cfg_file = "turb_SST_flatplate_roughBCAupoix.cfg" turb_flatplate_sst_roughBCAupoix.test_iter = 10 - turb_flatplate_sst_roughBCAupoix.test_vals = [-5.278728, -2.302982, -2.890419, 0.227585, -1.393786, 3.192601, -0.188729, 0.0071821] + turb_flatplate_sst_roughBCAupoix.test_vals = [-4.844108, -2.112938, -2.523195, 0.663960, -1.393641, 3.192483, -0.188503, 0.007182] test_list.append(turb_flatplate_sst_roughBCAupoix) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SA, restart from a solution without wall functions @@ -571,7 +571,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.037552, -16.376949, 1.080346, 0.018385, 20.000000, -1.564082, 20.000000, -4.180917, 0.000000] + turb_naca0012_sa.test_vals = [-12.037491, -16.376949, 1.080346, 0.018385, 20.000000, -1.564168, 20.000000, -4.180902, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.037489, -16.376949, 1.080346, 0.018385, 20.000000, -1.564143, 20.000000, -4.180945, 0.000000] turb_naca0012_sa.timeout = 3200 test_list.append(turb_naca0012_sa) @@ -581,7 +581,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.094692, -15.251093, -5.906365, 1.070413, 0.015775, -2.376043, 0.000000] + turb_naca0012_sst.test_vals = [-5.949151, -10.294577, -3.763019, 1.069291, 0.015902, -3.459014, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000] turb_naca0012_sst.timeout = 3200 test_list.append(turb_naca0012_sst) @@ -591,7 +591,7 @@ def main(): turb_naca0012_sst_sust.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust.test_iter = 10 - turb_naca0012_sst_sust.test_vals = [-12.081966, -14.837177, -5.733436, 1.000893, 0.019109, -2.241123] + turb_naca0012_sst_sust.test_vals = [-9.963808, -11.675704, -5.093478, 1.000889, 0.019126, -1.603990] turb_naca0012_sst_sust.test_vals_aarch64 = [-12.073964, -14.836726, -5.732390, 1.000050, 0.019144, -2.229074] turb_naca0012_sst_sust.timeout = 3200 test_list.append(turb_naca0012_sst_sust) @@ -601,7 +601,7 @@ def main(): turb_naca0012_sst_2003_Vm.cfg_dir = "rans/naca0012" turb_naca0012_sst_2003_Vm.cfg_file = "turb_NACA0012_sst_2003-Vm.cfg" turb_naca0012_sst_2003_Vm.test_iter = 10 - turb_naca0012_sst_2003_Vm.test_vals = [-10.168739, -3.668709, 1.060563, 0.019147, -2.282800] + turb_naca0012_sst_2003_Vm.test_vals = [-10.076544, -3.586031, 1.059474, 0.019231, -3.088893] turb_naca0012_sst_2003_Vm.timeout = 3200 test_list.append(turb_naca0012_sst_2003_Vm) @@ -610,7 +610,7 @@ def main(): turb_naca0012_sst_1994_KLm.cfg_dir = "rans/naca0012" turb_naca0012_sst_1994_KLm.cfg_file = "turb_NACA0012_sst_1994-KLm.cfg" turb_naca0012_sst_1994_KLm.test_iter = 10 - turb_naca0012_sst_1994_KLm.test_vals = [-10.390871, -3.701584, 1.062164, 0.019076, -2.426123] + turb_naca0012_sst_1994_KLm.test_vals = [-10.196762, -3.596952, 1.061057, 0.019163, -3.186595] turb_naca0012_sst_1994_KLm.timeout = 3200 test_list.append(turb_naca0012_sst_1994_KLm) @@ -619,7 +619,7 @@ def main(): turb_naca0012_sst_fixedvalues.cfg_dir = "rans/naca0012" turb_naca0012_sst_fixedvalues.cfg_file = "turb_NACA0012_sst_fixedvalues.cfg" turb_naca0012_sst_fixedvalues.test_iter = 10 - turb_naca0012_sst_fixedvalues.test_vals = [-10.440016, 0.774146, 1.022363, 0.040546, -3.736437] + turb_naca0012_sst_fixedvalues.test_vals = [-10.440029, 0.774146, 1.022363, 0.040546, -3.736439] turb_naca0012_sst_fixedvalues.timeout = 3200 test_list.append(turb_naca0012_sst_fixedvalues) @@ -628,7 +628,7 @@ def main(): turb_naca0012_sst_expliciteuler.cfg_dir = "rans/naca0012" turb_naca0012_sst_expliciteuler.cfg_file = "turb_NACA0012_sst_expliciteuler.cfg" turb_naca0012_sst_expliciteuler.test_iter = 10 - turb_naca0012_sst_expliciteuler.test_vals = [-3.532365, -3.157224, 3.743381, 1.124798, 0.501715, -16.000000] + turb_naca0012_sst_expliciteuler.test_vals = [-3.408465, -3.157224, 3.743383, 1.124790, 0.501715, -16.000000] turb_naca0012_sst_expliciteuler.timeout = 3200 test_list.append(turb_naca0012_sst_expliciteuler) @@ -660,7 +660,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-5.130654, 0.443770, -2.934019, 1.603840, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-5.126464, 0.450358, -1.631056, 2.883065, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143310, -9.163287, -10.858232, -5.787715, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) @@ -683,7 +683,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 20 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.560190, -5.057149, 0.830238, -0.008812, 0.078167] + turb_naca0012_sst_restart_mg.test_vals = [-3.692233, -5.056998, 0.830505, -0.008795, 0.078083] turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 test_list.append(turb_naca0012_sst_restart_mg) @@ -965,7 +965,7 @@ def main(): turbmod_sa_neg_rae2822.cfg_dir = "turbulence_models/sa/rae2822" turbmod_sa_neg_rae2822.cfg_file = "turb_SA_NEG_RAE2822.cfg" turbmod_sa_neg_rae2822.test_iter = 10 - turbmod_sa_neg_rae2822.test_vals = [1.275994, 1.062588, -2.067822, 1.300662, 0.480881, 0.000000] + turbmod_sa_neg_rae2822.test_vals = [1.275049, 1.063047, -2.067311, 1.280691, 0.469057, 0.000000] turbmod_sa_neg_rae2822.test_vals_aarch64 = [-1.345593, 1.448310, 1.208721, -0.846597, 1.248410, 0.489117, 0.000000] test_list.append(turbmod_sa_neg_rae2822) @@ -1115,7 +1115,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.979554, 1.344291, 0.666805, 0.010247] + turb_naca0012_1c.test_vals = [-4.965424, 1.346117, 0.672687, 0.011931] turb_naca0012_1c.test_vals_aarch64 = [-4.981036, 1.345868, 0.673232, 0.010091] test_list.append(turb_naca0012_1c) @@ -1124,7 +1124,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.482691, 1.262918, 0.529442, -0.023584] + turb_naca0012_2c.test_vals = [-5.482676, 1.262851, 0.535634, -0.021752] turb_naca0012_2c.test_vals_aarch64 = [-5.484365, 1.264701, 0.501741, -0.033109] test_list.append(turb_naca0012_2c) @@ -1133,7 +1133,7 @@ def main(): turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" turb_naca0012_3c.test_iter = 10 - turb_naca0012_3c.test_vals = [-5.583617, 1.229594, 0.463384, -0.035444] + turb_naca0012_3c.test_vals = [-5.583621, 1.229520, 0.465410, -0.035039] test_list.append(turb_naca0012_3c) # NACA0012 p1c1 @@ -1141,7 +1141,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.129555, 1.283892, 0.808074, 0.047784] + turb_naca0012_p1c1.test_vals = [-5.164519, 1.279293, 0.810485, 0.048509] turb_naca0012_p1c1.test_vals_aarch64 = [-5.122100, 1.284478, 0.608744, -0.008593] test_list.append(turb_naca0012_p1c1) @@ -1150,7 +1150,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.554045, 1.234446, 0.603731, -0.008717] + turb_naca0012_p1c2.test_vals = [-5.553255, 1.234324, 0.608768, -0.007437] test_list.append(turb_naca0012_p1c2) ###################################### @@ -1211,7 +1211,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-1.175979, 0.062203, 1.399352, 2.219226, 1.399298, 2.217472, 0.000000] + square_cylinder.test_vals = [-1.175953, 0.062131, 1.399400, 2.220360, 1.399348, 2.218598, 0.000000] square_cylinder.unsteady = True test_list.append(square_cylinder) @@ -1318,7 +1318,7 @@ def main(): datadriven_fluidModel.cfg_dir = "nicf/datadriven" datadriven_fluidModel.cfg_file = "datadriven_nozzle.cfg" datadriven_fluidModel.test_iter = 50 - datadriven_fluidModel.test_vals = [-6.283594, -3.253961, -4.172935, -0.949415, -2.046039, 1.546933] + datadriven_fluidModel.test_vals = [-5.779669, -2.794058, -3.760455, -0.618981, -1.475343, 1.677566] test_list.append(datadriven_fluidModel) ###################################### @@ -1339,7 +1339,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.905382, -12.212100, -19.212143, -13.508268, -19.078783, -13.446389, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-0.170809, -0.162035, -8.236415, -2.533086, -7.364784, -1.876217, 73290.000000, 73290.000000, 0.020065, 82.291000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -1347,7 +1347,7 @@ def main(): axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D" axial_stage2D.cfg_file = "Axial_stage2D.cfg" axial_stage2D.test_iter = 20 - axial_stage2D.test_vals = [1.167500, 1.598546, -2.928573, 2.573647, -2.527331, 3.016254, 106370.000000, 106370.000000, 5.726800, 64.383000] + axial_stage2D.test_vals = [1.166016, 1.598547, -2.928563, 2.573651, -2.527330, 3.016254, 106370.000000, 106370.000000, 5.726800, 64.395000] test_list.append(axial_stage2D) # 2D transonic stator restart @@ -1355,7 +1355,7 @@ def main(): transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg" transonic_stator_restart.test_iter = 20 - transonic_stator_restart.test_vals = [-4.365636, -2.482248, -2.081252, 1.729763, -1.467305, 3.231073, -471620.000000, 94.837000, -0.046834] + transonic_stator_restart.test_vals = [-4.361754, -2.479177, -2.080217, 1.732474, -1.467446, 3.229306, -471620.000000, 94.837000, -0.093140] transonic_stator_restart.test_vals_aarch64 = [-4.437809, -2.553049, -2.164729, 1.657542, -1.356823, 3.178788, -471620.000000, 94.842000, -0.040365] test_list.append(transonic_stator_restart) @@ -1447,7 +1447,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.457912, -1.541047] + bars_SST_2D.test_vals = [13.000000, -0.457935, -1.540951] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1647,7 +1647,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-12.094692, -15.251093, -5.906365, 1.070413, 0.015775, -2.376043, 0.000000] + pywrapper_turb_naca0012_sst.test_vals = [-5.949151, -10.294577, -3.763019, 1.069291, 0.015902, -3.459014, 0.000000] pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.075620, -15.246688, -5.861276, 1.070036, 0.015841, -1.991001, 0.000000] pywrapper_turb_naca0012_sst.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f") pywrapper_turb_naca0012_sst.timeout = 3200 @@ -1658,7 +1658,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 10 - pywrapper_square_cylinder.test_vals = [-1.178521, -0.349772, 1.401402, 2.359105, 1.401729, 2.300974, 0.000000] + pywrapper_square_cylinder.test_vals = [-1.178744, -0.349499, 1.401128, 2.358038, 1.401474, 2.301063, 0.000000] pywrapper_square_cylinder.command = TestCase.Command("mpirun -np 2", "SU2_CFD.py", "--parallel -f") pywrapper_square_cylinder.unsteady = True test_list.append(pywrapper_square_cylinder) @@ -1709,7 +1709,7 @@ def main(): pywrapper_unsteadyCHT.cfg_dir = "py_wrapper/flatPlate_unsteady_CHT" pywrapper_unsteadyCHT.cfg_file = "unsteady_CHT_FlatPlate_Conf.cfg" pywrapper_unsteadyCHT.test_iter = 5 - pywrapper_unsteadyCHT.test_vals = [-1.614168, 2.260203, -0.008065, 0.201068] + pywrapper_unsteadyCHT.test_vals = [-1.614169, 2.260208, -0.008127, 0.201931] pywrapper_unsteadyCHT.command = TestCase.Command("mpirun -np 2", "python", "launch_unsteady_CHT_FlatPlate.py --parallel -f") pywrapper_unsteadyCHT.unsteady = True test_list.append(pywrapper_unsteadyCHT) @@ -1977,7 +1977,7 @@ def main(): cgns_writer.cfg_dir = "cgns_writer" cgns_writer.cfg_file = "config.cfg" cgns_writer.test_iter = 1 - cgns_writer.test_vals = [-2.974473, 0.640256, 5.371028, -6.732060] + cgns_writer.test_vals = [-2.974219, 0.639253, 5.371471, -6.721279] cgns_writer.new_output = True test_list.append(cgns_writer) diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py index ebd1193f4b0c..2f3d24727c4a 100644 --- a/TestCases/parallel_regression_AD.py +++ b/TestCases/parallel_regression_AD.py @@ -91,7 +91,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.226820, -0.256327, -2.299900, -0.003171] + discadj_rans_naca0012_sst.test_vals = [-2.226827, -0.253089, -2.299100, -0.003174] discadj_rans_naca0012_sst.test_vals_aarch64 = [-2.226852, -0.253100, -2.299300, -0.003164] test_list.append(discadj_rans_naca0012_sst) @@ -150,7 +150,7 @@ def main(): discadj_axisymmetric_rans_nozzle.cfg_dir = "axisymmetric_rans/air_nozzle" discadj_axisymmetric_rans_nozzle.cfg_file = "air_nozzle_restart.cfg" discadj_axisymmetric_rans_nozzle.test_iter = 10 - discadj_axisymmetric_rans_nozzle.test_vals = [9.738854, 5.135195, 7.108430, 2.487299] + discadj_axisymmetric_rans_nozzle.test_vals = [9.738700, 5.135180, 7.101795, 2.487548] discadj_axisymmetric_rans_nozzle.no_restart = True test_list.append(discadj_axisymmetric_rans_nozzle) @@ -251,7 +251,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, -7.555647, -10.335486, -10.356919, -13.629543] + discadj_trans_stator.test_vals = [79.000000, -7.554772, -10.335151, -10.357769, -13.628908] discadj_trans_stator.test_vals_aarch64 = [79.000000, -7.555647, -10.335486, -10.356919, -13.629543] test_list.append(discadj_trans_stator) @@ -260,7 +260,7 @@ def main(): discadj_axial_stage.cfg_dir = "disc_adj_turbomachinery/axial_stage_2D" discadj_axial_stage.cfg_file = "Axial_stage2D.cfg" discadj_axial_stage.test_iter = 79 - discadj_axial_stage.test_vals = [79.000000, -6.606356, -7.139726] + discadj_axial_stage.test_vals = [79.000000, -6.605179, -7.180296] discadj_axial_stage.test_vals_aarch64 = [79.000000, -6.606356, -7.139726] test_list.append(discadj_axial_stage) @@ -315,7 +315,7 @@ def main(): discadj_fsi2.cfg_dir = "disc_adj_fsi/Airfoil_2d" discadj_fsi2.cfg_file = "config.cfg" discadj_fsi2.test_iter = 8 - discadj_fsi2.test_vals = [-3.824591, 1.979595, -3.863368, 0.295450, 3.839800] + discadj_fsi2.test_vals = [-3.824635, 1.979542, -3.863368, 0.295450, 3.839800] discadj_fsi2.test_vals_aarch64 = [-3.824870, 1.979160, -3.863368, 0.295450, 3.839800] discadj_fsi2.tol = 0.00001 test_list.append(discadj_fsi2) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index 16c57318510a..cb017999420a 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -216,7 +216,7 @@ def main(): poiseuille_profile.cfg_dir = "navierstokes/poiseuille" poiseuille_profile.cfg_file = "profile_poiseuille.cfg" poiseuille_profile.test_iter = 10 - poiseuille_profile.test_vals = [-12.003115, -7.626023, -0.000000, 2.089953] + poiseuille_profile.test_vals = [-12.003110, -7.625808, -0.000000, 2.089953] poiseuille_profile.test_vals_aarch64 = [-12.009012, -7.262299, -0.000000, 2.089953] #last 4 columns test_list.append(poiseuille_profile) @@ -244,7 +244,7 @@ def main(): rae2822_sst.cfg_dir = "rans/rae2822" rae2822_sst.cfg_file = "turb_SST_RAE2822.cfg" rae2822_sst.test_iter = 20 - rae2822_sst.test_vals = [-1.916371, -1.461911, 5.894621, 0.611793, 0.016873, 100.000000] + rae2822_sst.test_vals = [-1.849587, -1.466490, 5.887578, 0.592053, 0.015314, 100.000000] test_list.append(rae2822_sst) # RAE2822 SST_SUST @@ -268,7 +268,7 @@ def main(): turb_wallfunction_flatplate_sst.cfg_dir = "wallfunctions/flatplate/compressible_SST" turb_wallfunction_flatplate_sst.cfg_file = "turb_SST_flatplate.cfg" turb_wallfunction_flatplate_sst.test_iter = 10 - turb_wallfunction_flatplate_sst.test_vals = [-4.036300, -1.916264, -1.821790, 1.443878, -1.579554, 1.521535, 10.000000, -2.351460, 0.030011, 0.002409] + turb_wallfunction_flatplate_sst.test_vals = [-4.021540, -1.941019, -1.835888, 1.458341, -1.581319, 1.513862, 10.000000, -2.562969, 0.029270, 0.002412] test_list.append(turb_wallfunction_flatplate_sst) # FLAT PLATE, ROUGHNESS BC WILCOX2006 SST @@ -276,7 +276,7 @@ def main(): turb_flatplate_sst_roughBCWilcox2006.cfg_dir = "rans/flatplate/roughness/bc_wilcox2006" turb_flatplate_sst_roughBCWilcox2006.cfg_file = "turb_SST_flatplate_roughBCWilcox2006.cfg" turb_flatplate_sst_roughBCWilcox2006.test_iter = 10 - turb_flatplate_sst_roughBCWilcox2006.test_vals = [-5.117900, -2.534224, -2.904279, 0.381710, -3.100346, 1.180161, -0.188797, 0.004029] + turb_flatplate_sst_roughBCWilcox2006.test_vals = [-4.679459, -2.104822, -2.446835, 0.818054, -2.775001, 1.266518, -0.188523, 0.004029] test_list.append(turb_flatplate_sst_roughBCWilcox2006) # FLAT PLATE, WALL FUNCTIONS, COMPRESSIBLE SA @@ -301,7 +301,7 @@ def main(): turb_naca0012_sa.cfg_dir = "rans/naca0012" turb_naca0012_sa.cfg_file = "turb_NACA0012_sa.cfg" turb_naca0012_sa.test_iter = 5 - turb_naca0012_sa.test_vals = [-12.037319, -16.384159, 1.080346, 0.018385, 20.000000, -3.456939, 20.000000, -4.641254, 0.000000] + turb_naca0012_sa.test_vals = [-12.037313, -16.384159, 1.080346, 0.018385, 20.000000, -3.457122, 20.000000, -4.641254, 0.000000] turb_naca0012_sa.test_vals_aarch64 = [-12.037297, -16.384158, 1.080346, 0.018385, 20.000000, -3.455886, 20.000000, -4.641247, 0.000000] turb_naca0012_sa.timeout = 3200 test_list.append(turb_naca0012_sa) @@ -311,7 +311,7 @@ def main(): turb_naca0012_sst.cfg_dir = "rans/naca0012" turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" turb_naca0012_sst.test_iter = 10 - turb_naca0012_sst.test_vals = [-12.094431, -15.251082, -5.906366, 1.070413, 0.015775, -3.178469, 0.000000] + turb_naca0012_sst.test_vals = [-5.948789, -10.293776, -3.762067, 1.069269, 0.015902, -3.345367, 0.000000] turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000] turb_naca0012_sst.timeout = 3200 test_list.append(turb_naca0012_sst) @@ -334,7 +334,7 @@ def main(): turb_naca0012_sst_2003m.cfg_dir = "rans/naca0012" turb_naca0012_sst_2003m.cfg_file = "turb_NACA0012_sst_2003m.cfg" turb_naca0012_sst_2003m.test_iter = 10 - turb_naca0012_sst_2003m.test_vals = [-7.129622, -10.410728, -3.699440, 1.061671, 0.019047, -2.435494, 0.000000] + turb_naca0012_sst_2003m.test_vals = [-5.953767, -10.205843, -3.595874, 1.060514, 0.019139, -3.188840, 0.000000] turb_naca0012_sst_2003m.timeout = 3200 test_list.append(turb_naca0012_sst_2003m) @@ -343,7 +343,7 @@ def main(): turb_naca0012_sst_sust_restart.cfg_dir = "rans/naca0012" turb_naca0012_sst_sust_restart.cfg_file = "turb_NACA0012_sst_sust.cfg" turb_naca0012_sst_sust_restart.test_iter = 10 - turb_naca0012_sst_sust_restart.test_vals = [-12.080496, -14.837169, -5.733461, 1.000893, 0.019109, -2.634008] + turb_naca0012_sst_sust_restart.test_vals = [-9.962094, -11.675499, -5.093415, 1.000888, 0.019126, -1.923274] turb_naca0012_sst_sust_restart.test_vals_aarch64 = [-12.074189, -14.836725, -5.732398, 1.000050, 0.019144, -3.315560] turb_naca0012_sst_sust_restart.timeout = 3200 test_list.append(turb_naca0012_sst_sust_restart) @@ -353,7 +353,7 @@ def main(): turb_naca0012_sst_fixedvalues.cfg_dir = "rans/naca0012" turb_naca0012_sst_fixedvalues.cfg_file = "turb_NACA0012_sst_fixedvalues.cfg" turb_naca0012_sst_fixedvalues.test_iter = 10 - turb_naca0012_sst_fixedvalues.test_vals = [-5.206618, -10.436758, 0.774095, 1.021995, 0.040553, -3.477597] + turb_naca0012_sst_fixedvalues.test_vals = [-5.206619, -10.436774, 0.774095, 1.021995, 0.040553, -3.477596] turb_naca0012_sst_fixedvalues.timeout = 3200 test_list.append(turb_naca0012_sst_fixedvalues) @@ -375,7 +375,7 @@ def main(): axi_rans_air_nozzle_restart.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_restart.cfg_file = "air_nozzle_restart.cfg" axi_rans_air_nozzle_restart.test_iter = 10 - axi_rans_air_nozzle_restart.test_vals = [-5.121517, 0.452804, -2.910659, 1.662143, 0.000000] + axi_rans_air_nozzle_restart.test_vals = [-5.118766, 0.457805, -1.643100, 2.889958, 0.000000] axi_rans_air_nozzle_restart.test_vals_aarch64 = [-14.143715, -9.170705, -10.848554, -5.776746, 0.000000] axi_rans_air_nozzle_restart.tol = 0.0001 test_list.append(axi_rans_air_nozzle_restart) @@ -393,7 +393,7 @@ def main(): axi_rans_air_nozzle_species.cfg_dir = "axisymmetric_rans/air_nozzle" axi_rans_air_nozzle_species.cfg_file = "air_nozzle_species.cfg" axi_rans_air_nozzle_species.test_iter = 10 - axi_rans_air_nozzle_species.test_vals = [-2.148980, 3.342062, -1.939182, 5.883411, -2.443141, 0.000000] + axi_rans_air_nozzle_species.test_vals = [-2.148943, 3.342084, -1.939303, 5.883410, -2.443272, 0.000000] axi_rans_air_nozzle_species.tol = 0.0001 test_list.append(axi_rans_air_nozzle_species) @@ -407,7 +407,7 @@ def main(): turb_naca0012_sst_restart_mg.cfg_file = "turb_NACA0012_sst_multigrid_restart.cfg" turb_naca0012_sst_restart_mg.test_iter = 50 turb_naca0012_sst_restart_mg.ntest_vals = 5 - turb_naca0012_sst_restart_mg.test_vals = [-6.570974, -5.081421, 0.810883, -0.008830, 0.077967] + turb_naca0012_sst_restart_mg.test_vals = [-3.649967, -5.081630, 0.810585, -0.009208, 0.078506] turb_naca0012_sst_restart_mg.timeout = 3200 turb_naca0012_sst_restart_mg.tol = 0.000001 test_list.append(turb_naca0012_sst_restart_mg) @@ -760,7 +760,7 @@ def main(): turb_naca0012_1c.cfg_dir = "rans_uq/naca0012" turb_naca0012_1c.cfg_file = "turb_NACA0012_uq_1c.cfg" turb_naca0012_1c.test_iter = 10 - turb_naca0012_1c.test_vals = [-4.981181, 1.343583, 0.597111, 0.017244] + turb_naca0012_1c.test_vals = [-4.964306, 1.346025, 0.604256, 0.019699] turb_naca0012_1c.test_vals_aarch64 = [-4.992791, 1.342873, 0.557941, 0.003269] test_list.append(turb_naca0012_1c) @@ -769,7 +769,7 @@ def main(): turb_naca0012_2c.cfg_dir = "rans_uq/naca0012" turb_naca0012_2c.cfg_file = "turb_NACA0012_uq_2c.cfg" turb_naca0012_2c.test_iter = 10 - turb_naca0012_2c.test_vals = [-5.482907, 1.262063, 0.459671, -0.026410] + turb_naca0012_2c.test_vals = [-5.482904, 1.261984, 0.464231, -0.024984] test_list.append(turb_naca0012_2c) # NACA0012 3c @@ -777,7 +777,7 @@ def main(): turb_naca0012_3c.cfg_dir = "rans_uq/naca0012" turb_naca0012_3c.cfg_file = "turb_NACA0012_uq_3c.cfg" turb_naca0012_3c.test_iter = 10 - turb_naca0012_3c.test_vals = [-5.583768, 1.229824, 0.426258, -0.033150] + turb_naca0012_3c.test_vals = [-5.583772, 1.229781, 0.426469, -0.033207] test_list.append(turb_naca0012_3c) # NACA0012 p1c1 @@ -785,7 +785,7 @@ def main(): turb_naca0012_p1c1.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c1.cfg_file = "turb_NACA0012_uq_p1c1.cfg" turb_naca0012_p1c1.test_iter = 10 - turb_naca0012_p1c1.test_vals = [-5.127701, 1.284875, 0.779301, 0.086111] + turb_naca0012_p1c1.test_vals = [-5.163575, 1.280570, 0.794472, 0.091341] turb_naca0012_p1c1.test_vals_aarch64 = [-5.119942, 1.283920, 0.486264, -0.021518] test_list.append(turb_naca0012_p1c1) @@ -794,7 +794,7 @@ def main(): turb_naca0012_p1c2.cfg_dir = "rans_uq/naca0012" turb_naca0012_p1c2.cfg_file = "turb_NACA0012_uq_p1c2.cfg" turb_naca0012_p1c2.test_iter = 10 - turb_naca0012_p1c2.test_vals = [-5.553987, 1.235071, 0.521197, -0.004487] + turb_naca0012_p1c2.test_vals = [-5.553130, 1.235035, 0.577646, 0.018241] test_list.append(turb_naca0012_p1c2) ###################################### @@ -854,7 +854,7 @@ def main(): square_cylinder.cfg_dir = "unsteady/square_cylinder" square_cylinder.cfg_file = "turb_square.cfg" square_cylinder.test_iter = 3 - square_cylinder.test_vals = [-2.560678, -1.175979, 0.062202, 1.399351, 2.219234, 1.399297, 2.217479, 0.000000] + square_cylinder.test_vals = [-2.560698, -1.175954, 0.062130, 1.399399, 2.220367, 1.399347, 2.218605, 0.000000] square_cylinder.unsteady = True test_list.append(square_cylinder) @@ -974,7 +974,7 @@ def main(): Jones_tc_restart.cfg_dir = "turbomachinery/APU_turbocharger" Jones_tc_restart.cfg_file = "Jones_restart.cfg" Jones_tc_restart.test_iter = 5 - Jones_tc_restart.test_vals = [-11.905642, -12.212181, -19.212554, -13.508053, -19.069991, -13.439494, 73286.000000, 73286.000000, 0.020056, 82.286000] + Jones_tc_restart.test_vals = [-0.170799, -0.162039, -8.236404, -2.533076, -7.364785, -1.876219, 73290.000000, 73290.000000, 0.020065, 82.291000] test_list.append(Jones_tc_restart) # 2D axial stage @@ -982,7 +982,7 @@ def main(): axial_stage2D.cfg_dir = "turbomachinery/axial_stage_2D" axial_stage2D.cfg_file = "Axial_stage2D.cfg" axial_stage2D.test_iter = 20 - axial_stage2D.test_vals = [1.167521, 1.598536, -2.928575, 2.573646, -2.527330, 3.016253, 106370.000000, 106370.000000, 5.726800, 64.383000] + axial_stage2D.test_vals = [1.166038, 1.598537, -2.928565, 2.573650, -2.527330, 3.016254, 106370.000000, 106370.000000, 5.726800, 64.395000] test_list.append(axial_stage2D) # 2D transonic stator restart @@ -990,7 +990,7 @@ def main(): transonic_stator_restart.cfg_dir = "turbomachinery/transonic_stator_2D" transonic_stator_restart.cfg_file = "transonic_stator_restart.cfg" transonic_stator_restart.test_iter = 20 - transonic_stator_restart.test_vals = [-4.367784, -2.492912, -2.082414, 1.727491, -1.466974, 3.224730, -471620.000000, 94.839000, -0.052082] + transonic_stator_restart.test_vals = [-4.363858, -2.489999, -2.081340, 1.730232, -1.467013, 3.223220, -471620.000000, 94.839000, -0.100140] transonic_stator_restart.test_vals_aarch64 = [-4.367784, -2.492912, -2.082414, 1.727491, -1.466974, 3.224730, -471620.000000, 94.839000, -0.052082] test_list.append(transonic_stator_restart) @@ -1020,7 +1020,7 @@ def main(): uniform_flow.cfg_dir = "sliding_interface/uniform_flow" uniform_flow.cfg_file = "uniform_NN.cfg" uniform_flow.test_iter = 2 - uniform_flow.test_vals = [2.000000, 0.000000, -0.230680, -13.251592] + uniform_flow.test_vals = [2.000000, 0.000000, -0.230680, -13.248947] uniform_flow.test_vals_aarch64 = [2.000000, 0.000000, -0.230641, -13.249000] uniform_flow.tol = 0.000001 uniform_flow.unsteady = True @@ -1094,7 +1094,7 @@ def main(): bars_SST_2D.cfg_dir = "sliding_interface/bars_SST_2D" bars_SST_2D.cfg_file = "bars.cfg" bars_SST_2D.test_iter = 13 - bars_SST_2D.test_vals = [13.000000, -0.456143, -1.541051] + bars_SST_2D.test_vals = [13.000000, -0.456166, -1.540955] bars_SST_2D.multizone = True test_list.append(bars_SST_2D) @@ -1172,7 +1172,7 @@ def main(): fsi_cht.cfg_dir = "fea_fsi/stat_fsi" fsi_cht.cfg_file = "config.cfg" fsi_cht.test_iter = 20 - fsi_cht.test_vals = [5.000000, -5.077002, -5.379450, -9.247804, -9.320014, -9.185034, 608.350000, -0.012973, 0.000000, 30.000000] + fsi_cht.test_vals = [5.000000, -5.077011, -5.379457, -9.247803, -9.320382, -9.185192, 608.350000, -0.012973, 0.000000, 29.000000] fsi_cht.multizone = True test_list.append(fsi_cht) @@ -1181,7 +1181,7 @@ def main(): dyn_fsi.cfg_dir = "fea_fsi/dyn_fsi" dyn_fsi.cfg_file = "config.cfg" dyn_fsi.test_iter = 4 - dyn_fsi.test_vals = [-4.328792, -4.153004, 0.000000, 86.000000] + dyn_fsi.test_vals = [-4.328792, -4.153004, 0.000000, 85.000000] dyn_fsi.multizone = True dyn_fsi.unsteady = True test_list.append(dyn_fsi) @@ -1696,7 +1696,7 @@ def main(): pywrapper_turb_naca0012_sst.cfg_dir = "rans/naca0012" pywrapper_turb_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" pywrapper_turb_naca0012_sst.test_iter = 10 - pywrapper_turb_naca0012_sst.test_vals = [-12.094431, -15.251082, -5.906366, 1.070413, 0.015775, -3.178469, 0.000000] + pywrapper_turb_naca0012_sst.test_vals = [-5.948789, -10.293776, -3.762067, 1.069269, 0.015902, -3.345367, 0.000000] pywrapper_turb_naca0012_sst.test_vals_aarch64 = [-12.076068, -15.246740, -5.861280, 1.070036, 0.015841, -3.297854, 0.000000] pywrapper_turb_naca0012_sst.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_turb_naca0012_sst.timeout = 3200 @@ -1710,7 +1710,7 @@ def main(): pywrapper_square_cylinder.cfg_dir = "unsteady/square_cylinder" pywrapper_square_cylinder.cfg_file = "turb_square.cfg" pywrapper_square_cylinder.test_iter = 3 - pywrapper_square_cylinder.test_vals = [-2.560678, -1.175979, 0.062202, 1.399351, 2.219234, 1.399297, 2.217479, 0.000000] + pywrapper_square_cylinder.test_vals = [-2.560698, -1.175954, 0.062130, 1.399399, 2.220367, 1.399347, 2.218605, 0.000000] pywrapper_square_cylinder.command = TestCase.Command(exec = "SU2_CFD.py", param = "-f") pywrapper_square_cylinder.timeout = 1600 pywrapper_square_cylinder.tol = 0.00001 diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py index 83e768e954e8..9f9206919255 100644 --- a/TestCases/serial_regression_AD.py +++ b/TestCases/serial_regression_AD.py @@ -96,7 +96,7 @@ def main(): discadj_rans_naca0012_sst.cfg_dir = "disc_adj_rans/naca0012" discadj_rans_naca0012_sst.cfg_file = "turb_NACA0012_sst.cfg" discadj_rans_naca0012_sst.test_iter = 10 - discadj_rans_naca0012_sst.test_vals = [-2.096075, -0.181237, 0.350380, -0.023435] + discadj_rans_naca0012_sst.test_vals = [-2.096038, -0.178557, 0.351620, -0.023443] test_list.append(discadj_rans_naca0012_sst) ####################################### @@ -200,7 +200,7 @@ def main(): discadj_trans_stator.cfg_dir = "disc_adj_turbomachinery/transonic_stator_2D" discadj_trans_stator.cfg_file = "transonic_stator.cfg" discadj_trans_stator.test_iter = 79 - discadj_trans_stator.test_vals = [79.000000, -7.308168, -9.891121, -10.038669, -13.368502] + discadj_trans_stator.test_vals = [79.000000, -7.306692, -9.889618, -10.037399, -13.367132] discadj_trans_stator.test_vals_aarch64 = [79.000000, -7.308167, -9.891117, -10.038669, -13.368501] test_list.append(discadj_trans_stator) @@ -209,7 +209,7 @@ def main(): discadj_axial_stage.cfg_dir = "disc_adj_turbomachinery/axial_stage_2D" discadj_axial_stage.cfg_file = "Axial_stage2D.cfg" discadj_axial_stage.test_iter = 79 - discadj_axial_stage.test_vals = [79.000000, -6.605593, -7.138207] + discadj_axial_stage.test_vals = [79.000000, -6.604508, -7.178602] discadj_axial_stage.test_vals_aarch64 = [79.000000, -6.605593, -7.138207] test_list.append(discadj_axial_stage) diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py index 22e9e55fdb03..40b9366fc0bf 100644 --- a/TestCases/tutorials.py +++ b/TestCases/tutorials.py @@ -249,7 +249,7 @@ def main(): tutorial_trans_flatplate_T3A.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A" tutorial_trans_flatplate_T3A.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3A.test_iter = 20 - tutorial_trans_flatplate_T3A.test_vals = [-5.790137, -2.054834, -3.894659, -0.255074, -1.747220, 5.119341, -3.493237, 0.393262] + tutorial_trans_flatplate_T3A.test_vals = [-5.790225, -2.054843, -3.895225, -0.255157, -1.746780, 5.119343, -3.493256, 0.392487] tutorial_trans_flatplate_T3A.test_vals_aarch64 = [-5.808996, -2.070606, -3.969765, -0.277943, -1.953289, 1.708472, -3.514943, 0.357411] tutorial_trans_flatplate_T3A.no_restart = True test_list.append(tutorial_trans_flatplate_T3A) @@ -259,7 +259,7 @@ def main(): tutorial_trans_flatplate_T3Am.cfg_dir = "../Tutorials/compressible_flow/Transitional_Flat_Plate/Langtry_and_Menter/T3A-" tutorial_trans_flatplate_T3Am.cfg_file = "transitional_LM_model_ConfigFile.cfg" tutorial_trans_flatplate_T3Am.test_iter = 20 - tutorial_trans_flatplate_T3Am.test_vals = [-5.587389, -1.700868, -3.093935, -0.102834, -3.750523, 3.287643, -2.394575, 1.119623] + tutorial_trans_flatplate_T3Am.test_vals = [-5.588634, -1.700867, -3.096990, -0.103933, -3.750523, 3.287644, -2.394576, 1.119624] tutorial_trans_flatplate_T3Am.test_vals_aarch64 = [-5.540938, -1.681627, -2.878831, -0.058224, -3.695533, 3.413628, -2.385345, 1.103633] tutorial_trans_flatplate_T3Am.no_restart = True test_list.append(tutorial_trans_flatplate_T3Am) @@ -314,7 +314,7 @@ def main(): tutorial_unst_naca0012.cfg_dir = "../Tutorials/compressible_flow/Unsteady_NACA0012" tutorial_unst_naca0012.cfg_file = "unsteady_naca0012.cfg" tutorial_unst_naca0012.test_iter = 520 - tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.292186, 0.000000, 0.304255, 0.792893, 0.001907, 0.010888] + tutorial_unst_naca0012.test_vals = [520.000000, 0.000000, -5.299254, 0.000000, 0.306206, 0.779691, 0.003296, 0.016165] tutorial_unst_naca0012.test_vals_aarch64 = [520, 0, -5.292359, 0, 0.284720, 0.766329, 0.000954, 0.007565] tutorial_unst_naca0012.unsteady = True test_list.append(tutorial_unst_naca0012) diff --git a/TestCases/vandv.py b/TestCases/vandv.py index 88602b2d2275..a67d0a574636 100644 --- a/TestCases/vandv.py +++ b/TestCases/vandv.py @@ -63,7 +63,7 @@ def main(): flatplate_sst1994m.cfg_dir = "vandv/rans/flatplate" flatplate_sst1994m.cfg_file = "turb_flatplate_sst.cfg" flatplate_sst1994m.test_iter = 5 - flatplate_sst1994m.test_vals = [-13.040943, -10.136971, -10.942764, -7.985181, -10.323857, -4.732487, 0.002801] + flatplate_sst1994m.test_vals = [-5.418449, -3.351150, -2.912790, 0.065870, -2.985271, 1.227315, 0.002802] flatplate_sst1994m.test_vals_aarch64 = [-13.021715, -9.534786, -10.401912, -7.501836, -9.750800, -4.850665, 0.002807] test_list.append(flatplate_sst1994m) @@ -72,7 +72,7 @@ def main(): bump_sst1994m.cfg_dir = "vandv/rans/bump_in_channel" bump_sst1994m.cfg_file = "turb_bump_sst.cfg" bump_sst1994m.test_iter = 5 - bump_sst1994m.test_vals = [-11.928292, -10.095796, -9.512953, -6.445652, -11.774088, -6.988752, 0.004931] + bump_sst1994m.test_vals = [-6.497322, -4.234388, -3.813607, -1.011693, -4.012643, 0.181789, 0.004933] bump_sst1994m.test_vals_aarch64 = [-13.042689, -10.812982, -10.604523, -7.655547, -10.816257, -5.308083, 0.004911] test_list.append(bump_sst1994m) @@ -91,7 +91,7 @@ def main(): swbli_sst.cfg_dir = "vandv/rans/swbli" swbli_sst.cfg_file = "config_sst.cfg" swbli_sst.test_iter = 5 - swbli_sst.test_vals = [-11.569218, -10.909086, -11.607984, -10.431163, -11.407588, -2.637660, 0.001816, -1.305818, -3.514590, 13.399000] + swbli_sst.test_vals = [-6.628941, -5.974370, -6.441345, -5.452139, -7.663134, 0.071743, 0.001815, -1.199316, -3.785242, 13.399000] test_list.append(swbli_sst) # DSMA661 - SA @@ -99,7 +99,7 @@ def main(): dsma661_sa.cfg_dir = "vandv/rans/dsma661" dsma661_sa.cfg_file = "dsma661_sa_config.cfg" dsma661_sa.test_iter = 5 - dsma661_sa.test_vals = [-11.255379, -8.243275, -9.037314, -5.927070, -10.737675, 0.155687, 0.024232] + dsma661_sa.test_vals = [-11.247085, -8.242815, -8.952043, -5.904549, -10.737668, 0.155687, 0.024232] dsma661_sa.test_vals_aarch64 = [-11.293183, -8.241775, -9.083761, -6.011398, -10.737680, 0.155687, 0.024232] test_list.append(dsma661_sa) @@ -108,7 +108,7 @@ def main(): dsma661_sst.cfg_dir = "vandv/rans/dsma661" dsma661_sst.cfg_file = "dsma661_sst_config.cfg" dsma661_sst.test_iter = 5 - dsma661_sst.test_vals = [-11.025153, -8.156995, -9.057021, -5.947228, -10.650874, -7.884423, 0.155882, 0.023344] + dsma661_sst.test_vals = [-2.616671, 0.009492, -0.244740, 2.840296, -2.180395, 1.413352, 0.156029, 0.023339] dsma661_sst.test_vals_aarch64 = [-10.977195, -8.403731, -8.747068, -5.808899, -10.522786, -7.369851, 0.155875, 0.023353] test_list.append(dsma661_sst)