- Fixed tke integration in thermodynamic variables

Common/include/CConfig.hpp

@@ -1961,7 +1961,7 @@ class CConfig {
    * \brief Get the value of the non-dimensionalized freestream energy.
    * \return Non-dimensionalized freestream energy.
    */
-  su2double GetEnergy_FreeStreamND(void) const { cout << "La chiedo non-dimensionale " << endl; return Energy_FreeStreamND; }
+  su2double GetEnergy_FreeStreamND(void) const { return Energy_FreeStreamND; }
 
   /*!
    * \brief Get the value of the non-dimensionalized freestream viscosity.

SU2_CFD/include/numerics_simd/flow/convection/common.hpp

@@ -103,7 +103,6 @@ FORCEINLINE void musclEdgeLimited(Int iPoint,
 template<class ReconVarType, class PrimVarType, size_t nDim, class VariableType>
 FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int jPoint,
                                                       bool muscl, LIMITER limiterType,
-                                                      bool musclTurb, LIMITER limiterTypeTurb,
                                                       const CPair<PrimVarType>& V1st,
                                                       const VectorDbl<nDim>& vector_ij,
                                                       const VariableType& solution,
@@ -113,9 +112,6 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
   const auto& gradients = solution.GetGradient_Reconstruction();
   const auto& limiters = solution.GetLimiter_Primitive();
 
-  // const auto& gradientsTurb = turbSolution.GetGradient_Reconstruction();
-  // const auto& limitersTurb = turbSolution.GetLimiter_Primitive();
-
   CPair<ReconVarType> V;
 
   for (size_t iVar = 0; iVar < ReconVarType::nVar; ++iVar) {
@@ -153,10 +149,9 @@ FORCEINLINE CPair<ReconVarType> reconstructPrimitives(Int iEdge, Int iPoint, Int
     for (size_t iDim = 0; iDim < nDim; ++iDim) {
       v_squared += pow(R*V.j.velocity(iDim) + V.i.velocity(iDim), 2);
     }
-    Double tke = R*V1st.j.allTurb(0) + V1st.i.allTurb(0);
     /*--- Multiply enthalpy by R+1 since v^2 was not divided by (R+1)^2.
      * Note: a = sqrt((gamma-1) * (H - 0.5 * v^2)) ---*/
-    const Double neg_sound_speed = enthalpy * (R+1) < (0.5 * v_squared - tke);
+    const Double neg_sound_speed = enthalpy * (R+1) < 0.5 * v_squared;
 
     /*--- Revert to first order if the state is non-physical. ---*/
     Double bad_recon = fmax(neg_p_or_rho, neg_sound_speed);

SU2_CFD/include/numerics_simd/flow/convection/roe.hpp

@@ -60,9 +60,7 @@ class CRoeBase : public Base {
   const bool finestGrid;
   const bool dynamicGrid;
   const bool muscl;
-  const bool musclTurb;
   const LIMITER typeLimiter;
-  const LIMITER typeLimiterTurb;
 
   using Base::turbVars;
 
@@ -77,9 +75,7 @@ class CRoeBase : public Base {
     finestGrid(iMesh == MESH_0),
     dynamicGrid(config.GetDynamic_Grid()),
     muscl(finestGrid && config.GetMUSCL_Flow()),
-    musclTurb(finestGrid && config.GetMUSCL_Turb()),
-    typeLimiter(config.GetKind_SlopeLimit_Flow()),
-    typeLimiterTurb(config.GetKind_SlopeLimit_Turb()) {
+    typeLimiter(config.GetKind_SlopeLimit_Flow()) {
   }
 
 public:
@@ -132,7 +128,7 @@ class CRoeBase : public Base {
     }
 
     auto V = reconstructPrimitives<CCompressiblePrimitives<nDim,nPrimVarGrad> >(
-                 iEdge, iPoint, jPoint, muscl, typeLimiter, musclTurb, typeLimiterTurb, V1st, vector_ij, solution, tkeNeeded);
+                 iEdge, iPoint, jPoint, muscl, typeLimiter, V1st, vector_ij, solution, tkeNeeded);
 
     /*--- Compute conservative variables. ---*/
 

SU2_CFD/include/numerics_simd/flow/variables.hpp

@@ -136,7 +136,7 @@ FORCEINLINE CRoeVariables<nDim> roeAveragedVariables(Double gamma,
   }
   roeAvg.enthalpy = (R*V.j.enthalpy() + V.i.enthalpy()) * D;
   roeAvg.tke = (R*V.j.tke() + V.i.tke()) * D;
-  roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity) - roeAvg.tke));
+  roeAvg.speedSound = sqrt((gamma-1) * (roeAvg.enthalpy - 0.5*squaredNorm(roeAvg.velocity)));
   roeAvg.projVel = dot(roeAvg.velocity, normal);
   return roeAvg;
 }

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -5494,7 +5494,9 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain
   unsigned short nSpanWiseSections = geometry->GetnSpanWiseSections(config->GetMarker_All_TurbomachineryFlag(val_marker));
   bool viscous = config->GetViscous();
   bool gravity = (config->GetGravityForce());
-  bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
+  bool tkeNeeded = ((config->GetKind_Turb_Model() == TURB_MODEL::SST) && !(config->GetSSTParsedOptions().modified));
+  CVariable* turbNodes = nullptr;
+  if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   su2double *Normal, *turboNormal, *UnitNormal, *FlowDirMix, FlowDirMixMag, *turboVelocity;
   Normal = new su2double[nDim];
@@ -5568,6 +5570,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);
 
@@ -5613,11 +5616,12 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain
               turboVelocity[iDim] = sqrt(Velocity2_e)*Flow_Dir[iDim];
             ComputeBackVelocity(turboVelocity,turboNormal, Velocity_e, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone));
             StaticEnthalpy_e = Enthalpy_e - 0.5 * Velocity2_e;
+            if(tkeNeeded) StaticEnthalpy_e -= turbNodes->GetSolution(iPoint, 0);
             GetFluidModel()->SetTDState_hs(StaticEnthalpy_e, Entropy_e);
             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:
@@ -5648,10 +5652,12 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain
             ComputeBackVelocity(turboVelocity,turboNormal, Velocity_e, config->GetMarker_All_TurbomachineryFlag(val_marker),config->GetKind_TurboMachinery(iZone));
 
             StaticEnthalpy_e = Enthalpy_e - 0.5 * Velocity2_e;
+            if(tkeNeeded) StaticEnthalpy_e -= turbNodes->GetSolution(iPoint, 0);
             GetFluidModel()->SetTDState_hs(StaticEnthalpy_e, Entropy_e);
             Density_e = GetFluidModel()->GetDensity();
             StaticEnergy_e = GetFluidModel()->GetStaticEnergy();
             Energy_e = StaticEnergy_e + 0.5 * Velocity2_e;
+            if(tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint, 0);
             // if (tkeNeeded) Energy_e += GetTke_Inf();
             break;
 
@@ -5670,6 +5676,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:
@@ -5687,6 +5694,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;
 
 
@@ -5704,6 +5712,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:
@@ -6915,7 +6924,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
   const su2double Gas_Constant = config->GetGas_ConstantND();
   const auto Kind_Inlet = config->GetKind_Inlet();
   const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker);
-  const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
+  const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified));
 
   /*--- Loop over all the vertices on this boundary marker ---*/
 

SU2_CFD/src/solvers/CIncNSSolver.cpp

@@ -287,7 +287,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, c
   const TURB_MODEL turb_model = config->GetKind_Turb_Model();
   const SPECIES_MODEL species_model = config->GetKind_Species_Model();
 
-  bool tkeNeeded = (turb_model == TURB_MODEL::SST);
+  bool tkeNeeded = (turb_model == TURB_MODEL::SST && !(config->GetSSTParsedOptions().modified));
 
   AD::StartNoSharedReading();