[WIP] Inconsistencies and improvements to SST model

Common/include/CConfig.hpp

@@ -856,7 +856,8 @@ class CConfig {
   nMu_Temperature_Ref,           /*!< \brief Number of species reference temperature for Sutherland model. */
   nMu_S,                         /*!< \brief Number of species reference S for Sutherland model. */
   nThermal_Conductivity_Constant,/*!< \brief Number of species constant thermal conductivity. */
-  nPrandtl_Lam,                  /*!< \brief Number of species laminar Prandtl number. */
+  nPrandtl_Lam,                  /*!< \brief Number of species 
+  addDoubleOption("FREESTREAM_TURB2LAMVISCRATIO", TurbIntensityAndViscRatioFreeStream[1], 10.0); Prandtl number. */
   nPrandtl_Turb,                 /*!< \brief Number of species turbulent Prandtl number. */
   nConstant_Lewis_Number;       /*!< \brief Number of species Lewis Number. */
   su2double Diffusivity_Constant;   /*!< \brief Constant mass diffusivity for scalar transport.  */
@@ -868,6 +869,7 @@ class CConfig {
   array<su2double, N_POLY_COEFFS> MuPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for viscosity. */
   array<su2double, N_POLY_COEFFS> KtPolyCoefficientsND{{0.0}};  /*!< \brief Definition of the non-dimensional temperature polynomial coefficients for thermal conductivity. */
   su2double TurbIntensityAndViscRatioFreeStream[2]; /*!< \brief Freestream turbulent intensity and viscosity ratio for turbulence and transition models. */
+
   su2double Energy_FreeStream,     /*!< \brief Free-stream total energy of the fluid.  */
   ModVel_FreeStream,               /*!< \brief Magnitude of the free-stream velocity of the fluid.  */
   ModVel_FreeStreamND,             /*!< \brief Non-dimensional magnitude of the free-stream velocity of the fluid.  */
@@ -1165,6 +1167,10 @@ class CConfig {
   nHistoryOutput, nVolumeOutput;  /*!< \brief Number of variables printed to the history file. */
   bool Multizone_Residual;        /*!< \brief Determines if memory should be allocated for the multizone residual. */
   SST_ParsedOptions sstParsedOptions; /*!< \brief Additional parameters for the SST turbulence model. */
+  su2double lowerLimitTKE,
+            lowerLimitDissipation;
+  su2double prodLimConst;
+  su2double LDomain;
   SA_ParsedOptions saParsedOptions;   /*!< \brief Additional parameters for the SA turbulence model. */
   LM_ParsedOptions lmParsedOptions;   /*!< \brief Additional parameters for the LM transition model. */
   su2double uq_delta_b;         /*!< \brief Parameter used to perturb eigenvalues of Reynolds Stress Matrix */
@@ -9847,6 +9853,11 @@ class CConfig {
    */
   SST_ParsedOptions GetSSTParsedOptions() const { return sstParsedOptions; }
 
+  su2double GetLowerLimitTKE() const { return lowerLimitTKE; }
+  su2double GetLowerLimitDissipation() const { return lowerLimitDissipation; }
+  su2double GetProdLimConst() const { return prodLimConst; }
+  su2double GetLDomain() const { return LDomain; }
+
   /*!
    * \brief Get parsed SA option data structure.
    * \return SA option data structure.

Common/include/option_structure.hpp

@@ -990,18 +990,26 @@ enum class SST_OPTIONS {
   V,           /*!< \brief Menter k-w SST model with vorticity production terms. */
   KL,          /*!< \brief Menter k-w SST model with Kato-Launder production terms. */
   UQ,          /*!< \brief Menter k-w SST model with uncertainty quantification modifications. */
+  FULLPROD,    /*!< \brief Menter k-w SST model with full production term. */
+  LLT,         /*!< \brief Menter k-w SST model with Lower Limits Treatments. */
+  PRODLIM,     /*!< \brief Menter k-w SST model with user-defined production limiter constant. */
+  NEWBC,       /*!< \brief Menter k-w SST model with new boundary conditions. */
 };
 static const MapType<std::string, SST_OPTIONS> SST_Options_Map = {
   MakePair("NONE", SST_OPTIONS::NONE)
   MakePair("V1994m", SST_OPTIONS::V1994m)
   MakePair("V2003m", SST_OPTIONS::V2003m)
   /// TODO: For now we do not support "unmodified" versions of SST.
-  //MakePair("V1994", SST_OPTIONS::V1994)
-  //MakePair("V2003", SST_OPTIONS::V2003)
+  MakePair("V1994", SST_OPTIONS::V1994)
+  MakePair("V2003", SST_OPTIONS::V2003)
   MakePair("SUSTAINING", SST_OPTIONS::SUST)
   MakePair("VORTICITY", SST_OPTIONS::V)
   MakePair("KATO-LAUNDER", SST_OPTIONS::KL)
   MakePair("UQ", SST_OPTIONS::UQ)
+  MakePair("FULLPROD", SST_OPTIONS::FULLPROD)
+  MakePair("LLT", SST_OPTIONS::LLT)
+  MakePair("PRODLIM", SST_OPTIONS::PRODLIM)
+  MakePair("NEWBC", SST_OPTIONS::NEWBC)
 };
 
 /*!
@@ -1013,6 +1021,10 @@ struct SST_ParsedOptions {
   bool sust = false;                          /*!< \brief Bool for SST model with sustaining terms. */
   bool uq = false;                            /*!< \brief Bool for using uncertainty quantification. */
   bool modified = false;                      /*!< \brief Bool for modified (m) SST model. */
+  bool fullProd = false;                      /*!< \brief Bool for full production term. */
+  bool llt = false;                           /*!< \brief Bool for Lower Limits Treatment. */
+  bool prodLim = false;                       /*!< \brief Bool for user-defined production limiter constant. */
+  bool newBC = false;                         /*!< \brief Bool for new boundary conditions. */
 };
 
 /*!
@@ -1030,6 +1042,11 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
     return std::find(SST_Options, sst_options_end, option) != sst_options_end;
   };
 
+  const bool found_fullProd = IsPresent(SST_OPTIONS::FULLPROD);
+  const bool found_llt = IsPresent(SST_OPTIONS::LLT);
+  const bool found_prodLim = IsPresent(SST_OPTIONS::PRODLIM);
+  const bool found_newBC = IsPresent(SST_OPTIONS::NEWBC);
+
   const bool found_1994 = IsPresent(SST_OPTIONS::V1994);
   const bool found_2003 = IsPresent(SST_OPTIONS::V2003);
   const bool found_1994m = IsPresent(SST_OPTIONS::V1994m);
@@ -1071,6 +1088,10 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
   SSTParsedOptions.sust = sst_sust;
   SSTParsedOptions.modified = sst_m;
   SSTParsedOptions.uq = sst_uq;
+  SSTParsedOptions.fullProd = found_fullProd;
+  SSTParsedOptions.llt = found_llt;
+  SSTParsedOptions.prodLim = found_prodLim;
+  SSTParsedOptions.newBC = found_newBC;
   return SSTParsedOptions;
 }
 

Common/src/CConfig.cpp

@@ -1118,6 +1118,11 @@ void CConfig::SetConfig_Options() {
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify SA turbulence model options/corrections. \n Options: see \link SA_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumListOption("SA_OPTIONS", nSA_Options, SA_Options, SA_Options_Map);
 
+  addDoubleOption("LOWER_LIMIT_TKE", lowerLimitTKE, 1e-20);
+  addDoubleOption("LOWER_LIMIT_DISSIPATION", lowerLimitDissipation, 1e-6);
+  addDoubleOption("PROD_LIM_CONST", prodLimConst, 20.0);
+  addDoubleOption("L_DOMAIN", LDomain, 1.0);
+
   /*!\brief KIND_TRANS_MODEL \n DESCRIPTION: Specify transition model OPTIONS: see \link Trans_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumOption("KIND_TRANS_MODEL", Kind_Trans_Model, Trans_Model_Map, TURB_TRANS_MODEL::NONE);
   /*!\brief SST_OPTIONS \n DESCRIPTION: Specify LM transition model options/correlations. \n Options: see \link LM_Options_Map \endlink \n DEFAULT: NONE \ingroup Config*/
@@ -6175,6 +6180,10 @@ void CConfig::SetOutput(SU2_COMPONENT val_software, unsigned short val_izone) {
                 break;
             }
             cout << "." << endl;
+
+            if (sstParsedOptions.prodLim) cout << "Changing the value of the TKE production limiter constant to " << prodLimConst << endl;
+            if (sstParsedOptions.llt) cout << "Changing the value of the lower limits of TKE and Omega " << endl;
+
             break;
         }
         switch (Kind_Trans_Model) {

SU2_CFD/include/numerics/CNumerics.hpp

@@ -190,6 +190,8 @@ class CNumerics {
 
   bool bounded_scalar = false;    /*!< \brief Flag for bounded scalar problem */
 
+  su2double ProdDistr[5];
+
 public:
   /*!
    * \brief Return type used in some "ComputeResidual" overloads to give a
@@ -1605,6 +1607,8 @@ class CNumerics {
    * \return is_bounded_scalar : scalar solver uses bounded scalar convective transport
    */
   inline bool GetBoundedScalar() const { return bounded_scalar;}
+
+  inline su2double GetProdDest(int index) const {return ProdDistr[index];}
 };
 
 /*!

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -609,6 +609,8 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
   /*--- Ambient values for SST-SUST. ---*/
   const su2double kAmb, omegaAmb;
 
+  su2double Pk, Dk, Pw, Dw;
+
   su2double F1_i, F2_i, CDkw_i;
   su2double Residual[2];
   su2double* Jacobian_i[2];
@@ -805,6 +807,11 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       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);
+      if (!sstParsedOptions.modified) P -= Eddy_Viscosity_i * diverg*diverg * 2.0/3.0;
+      if (sstParsedOptions.fullProd) P -= Density_i * ScalarVar_i[0] * diverg * 2.0/3.0;
+
+      su2double PLim = 0.0;
+      if ( P > prod_limit ) PLim = 1.0;
       su2double pk = max(0.0, min(P, prod_limit));
 
       const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i;
@@ -843,18 +850,28 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       }
 
       /*--- Add the production terms to the residuals. ---*/
+      Pk = pk;
+      Pw = pw;
 
       Residual[0] += pk * Volume;
       Residual[1] += pw * Volume;
 
       /*--- Add the dissipation  terms to the residuals.---*/
+      Dk = dk;
+      Dw = dw;
 
       Residual[0] -= dk * Volume;
       Residual[1] -= dw * Volume;
 
-      /*--- Cross diffusion ---*/
+      ProdDistr[0] = Pk;
+      ProdDistr[1] = Dk;
+      ProdDistr[2] = Pw;
+      ProdDistr[3] = Dw;
+      ProdDistr[4] = PLim;
 
-      Residual[1] += (1.0 - F1_i) * CDkw_i * Volume;
+      /*--- Cross diffusion ---*/
+      const su2double CrossDiffusionTerm =  CDkw_i < 0.0 ? -min(Pw, abs(CrossDiffusionTerm)) : CrossDiffusionTerm;
+      Residual[1] += (1.0 - F1_i) * CrossDiffusionTerm * Volume;
 
       /*--- Contribution due to 2D axisymmetric formulation ---*/
 
@@ -863,6 +880,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Implicit part ---*/
 
       Jacobian_i[0][0] = -beta_star * ScalarVar_i[1] * Volume;
+      if (sstParsedOptions.fullProd) Jacobian_i[0][0] -= diverg * Volume*2.0/3.0;
       Jacobian_i[0][1] = -beta_star * ScalarVar_i[0] * Volume;
       Jacobian_i[1][0] = 0.0;
       Jacobian_i[1][1] = -2.0 * beta_blended * ScalarVar_i[1] * Volume;

SU2_CFD/include/variables/CTurbVariable.hpp

@@ -43,6 +43,11 @@ class CTurbVariable : public CScalarVariable {
   static constexpr size_t MAXNVAR = 2;
   VectorType turb_index;
   VectorType intermittency;         /*!< \brief Value of the intermittency for the trans. model. */
+  VectorType Pk;
+  VectorType Pw;
+  VectorType Dk;
+  VectorType Dw;
+  VectorType PkLim;
 
   /*!
    * \brief Constructor of the class.
@@ -100,5 +105,38 @@ class CTurbVariable : public CScalarVariable {
    */
   inline void SetIntermittency(unsigned long iPoint, su2double val_intermittency) final { intermittency(iPoint) = val_intermittency; }
 
+  /*!
+   * \brief Set the intermittency of the transition model.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_intermittency - New value of the intermittency.
+   */
+  inline void SetProdDestr(unsigned long iPoint, su2double* val_ProdDestr) final { 
+    Pk(iPoint) = val_ProdDestr[0];
+    Dk(iPoint) = val_ProdDestr[1];
+    Pw(iPoint) = val_ProdDestr[2];
+    Dw(iPoint) = val_ProdDestr[3]; 
+    PkLim(iPoint) = val_ProdDestr[4]; 
+  }
+
+  /*!
+   * \brief Set the intermittency of the transition model.
+   * \param[in] iPoint - Point index.
+   * \param[in] val_intermittency - New value of the intermittency.
+   */
+  inline su2double GetProdTKE(unsigned long iPoint) const final { 
+    return Pk(iPoint); 
+  }
+  inline su2double GetDestrTKE(unsigned long iPoint) const final { 
+    return Dk(iPoint); 
+  }
+  inline su2double GetProdW(unsigned long iPoint) const final { 
+    return Pw(iPoint); 
+  }
+  inline su2double GetDestrW(unsigned long iPoint) const final { 
+    return Dw(iPoint); 
+  }
+  inline su2double GetPkLim(unsigned long iPoint) const final { 
+    return PkLim(iPoint); 
+  }
 };
 

SU2_CFD/include/variables/CVariable.hpp

@@ -2340,4 +2340,11 @@ class CVariable {
 
   inline virtual const su2double *GetScalarSources(unsigned long iPoint) const { return nullptr; }
   inline virtual const su2double *GetScalarLookups(unsigned long iPoint) const { return nullptr; }
+
+  inline virtual void SetProdDestr(unsigned long iPoint, su2double* val_ProdDestr) { }
+  inline virtual su2double GetProdTKE(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double GetDestrTKE(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double GetProdW(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double GetDestrW(unsigned long iPoint) const { return 0.0; }
+  inline virtual su2double GetPkLim(unsigned long iPoint) const { return 0.0; }
 };

SU2_CFD/src/output/CFlowOutput.cpp

@@ -1252,6 +1252,13 @@ void CFlowOutput::SetVolumeOutputFieldsScalarSolution(const CConfig* config){
     case TURB_FAMILY::KW:
       AddVolumeOutput("TKE", "Turb_Kin_Energy", "SOLUTION", "Turbulent kinetic energy");
       AddVolumeOutput("DISSIPATION", "Omega", "SOLUTION", "Rate of dissipation");
+      AddVolumeOutput("PROD_TKE", "Prod_TKE", "SOLUTION", "Production of turbulent kinetic energy");
+      AddVolumeOutput("DESTR_TKE", "Destr_TKE", "SOLUTION", "Destruction of turbulent kinetic energy");
+      AddVolumeOutput("PROD_TKE_LIM", "Prod_TKE_Lim", "SOLUTION", "Check if production limiter has been used for TKE");
+      AddVolumeOutput("PROD_W", "Prod_W", "SOLUTION", "Production of rate of dissipation");
+      AddVolumeOutput("DESTR_W", "Destr_W", "SOLUTION", "Destruction of rate of dissipation");
+      AddVolumeOutput("CDkw", "CDkw", "SOLUTION", "Cross-Diffusion term");
+      AddVolumeOutput("F1", "F1", "SOLUTION", "F1 blending function");
       break;
 
     case TURB_FAMILY::NONE:
@@ -1479,6 +1486,8 @@ void CFlowOutput::SetVolumeOutputFieldsScalarMisc(const CConfig* config) {
       AddVolumeOutput("VORTICITY", "Vorticity", "VORTEX_IDENTIFICATION", "Value of the vorticity");
     }
     AddVolumeOutput("Q_CRITERION", "Q_Criterion", "VORTEX_IDENTIFICATION", "Value of the Q-Criterion");
+    AddVolumeOutput("WALL_DISTANCE", "Wall_Distance", "MISC", "Wall distance value");
+    AddVolumeOutput("STRAIN_MAG", "Strain_Magnitude", "MISC", "Strain magnitude value");
   }
 
   // Timestep info
@@ -1512,6 +1521,8 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
       SetVolumeOutputValue("VORTICITY", iPoint, Node_Flow->GetVorticity(iPoint)[2]);
     }
     SetVolumeOutputValue("Q_CRITERION", iPoint, GetQCriterion(Node_Flow->GetVelocityGradient(iPoint)));
+    SetVolumeOutputValue("WALL_DISTANCE", iPoint, Node_Geo->GetWall_Distance(iPoint));
+    SetVolumeOutputValue("STRAIN_MAG", iPoint, Node_Turb->GetStrainMag(iPoint));
   }
 
   const bool limiter = (config->GetKind_SlopeLimit_Turb() != LIMITER::NONE);
@@ -1528,6 +1539,13 @@ void CFlowOutput::LoadVolumeDataScalar(const CConfig* config, const CSolver* con
     case TURB_FAMILY::KW:
       SetVolumeOutputValue("TKE", iPoint, Node_Turb->GetSolution(iPoint, 0));
       SetVolumeOutputValue("DISSIPATION", iPoint, Node_Turb->GetSolution(iPoint, 1));
+      SetVolumeOutputValue("PROD_TKE", iPoint, Node_Turb->GetProdTKE(iPoint));
+      SetVolumeOutputValue("DESTR_TKE", iPoint, Node_Turb->GetDestrTKE(iPoint));
+      SetVolumeOutputValue("PROD_TKE_LIM", iPoint, Node_Turb->GetPkLim(iPoint));
+      SetVolumeOutputValue("PROD_W", iPoint, Node_Turb->GetProdW(iPoint));
+      SetVolumeOutputValue("DESTR_W", iPoint, Node_Turb->GetDestrW(iPoint));
+      SetVolumeOutputValue("CDkw", iPoint, Node_Turb->GetCrossDiff(iPoint));
+      SetVolumeOutputValue("F1", iPoint, Node_Turb->GetF1blending(iPoint));
       SetVolumeOutputValue("RES_TKE", iPoint, turb_solver->LinSysRes(iPoint, 0));
       SetVolumeOutputValue("RES_DISSIPATION", iPoint, turb_solver->LinSysRes(iPoint, 1));
       if (limiter) {

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -4768,6 +4768,8 @@ 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;
+  CVariable* turbNodes = nullptr;
+  if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   auto *Normal = new su2double[nDim];
 
@@ -4917,7 +4919,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();
+      // if (tkeNeeded) Energy += GetTke_Inf();
+      if (tkeNeeded) Energy += turbNodes->GetSolution(iPoint,0);
 
       /*--- Store new primitive state for computing the flux. ---*/
 
@@ -5011,6 +5014,8 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container,
              implicit     = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT),
              gravity      = config->GetGravityForce(),
              tkeNeeded    = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
+  CVariable* turbNodes = nullptr;
+  if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   su2double **P_Tensor = new su2double*[nVar],
             **invP_Tensor = new su2double*[nVar];
@@ -5110,7 +5115,8 @@ 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();
+        // if (tkeNeeded) Energy_e += GetTke_Inf();
+        if (tkeNeeded) Energy_e += turbNodes->GetSolution(iPoint,0);
         break;
 
       case STATIC_SUPERSONIC_INFLOW_PT:
@@ -7176,6 +7182,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 ---*/
@@ -7258,6 +7267,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
         }
         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);
@@ -7341,6 +7351,8 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con
   const bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   const auto Marker_Tag = config->GetMarker_All_TagBound(val_marker);
   const bool tkeNeeded = (config->GetKind_Turb_Model() == TURB_MODEL::SST);
+  CVariable* turbNodes = nullptr;
+  if (tkeNeeded) turbNodes = solver_container[TURB_SOL]->GetNodes();
 
   /*--- Supersonic inlet flow: there are no outgoing characteristics,
    so all flow variables can be imposed at the inlet.
@@ -7361,14 +7373,15 @@ 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) Energy += GetTke_Inf();
+  const su2double Energy_woTKE = Pressure / (Density * Gamma_Minus_One) + 0.5 * Velocity2;
+  su2double Energy = Energy_woTKE;
 
   /*--- Loop over all the vertices on this boundary marker ---*/
 
   SU2_OMP_FOR_DYN(OMP_MIN_SIZE)
   for (auto iVertex = 0ul; iVertex < geometry->nVertex[val_marker]; iVertex++) {
     const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
+    if (tkeNeeded) Energy = Energy_woTKE + turbNodes->GetSolution(iPoint,0);
 
     if (!geometry->nodes->GetDomain(iPoint)) continue;