fix sensible enthalpy

bigfooted · bigfooted · commit 0f6df0ecc005 · 2025-12-23T20:25:33.000+01:00
diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
@@ -919,7 +919,6 @@ class CConfig {
   Initial_BCThrust,                /*!< \brief Ratio of turbulent to laminar viscosity at the actuator disk. */
   Pressure_FreeStream,             /*!< \brief Total pressure of the fluid. */
   Pressure_Thermodynamic,          /*!< \brief Thermodynamic pressure of the fluid. */
-  Standard_Ref_Temperature,        /*!< \brief Standard reference temperature for multicomponent flows. */
   Temperature_FreeStream,          /*!< \brief Total temperature of the fluid.  */
   Temperature_ve_FreeStream;       /*!< \brief Total vibrational-electronic temperature of the fluid.  */
   unsigned short wallModel_MaxIter; /*!< \brief maximum number of iterations for the Newton method for the wall model */
@@ -1959,12 +1958,6 @@ class CConfig {
    */
   su2double GetPressure_Thermodynamic(void) const { return Pressure_Thermodynamic; }
 
-  /*!
-   * \brief Get the value of the standard reference temperature for multicomponent flows.
-   * \return Standard reference temperature, Non-dimensionalized if it is needed for Non-Dimensional problems.
-   */
-  su2double GetStandard_RefTemperatureND(void) const { return Standard_Ref_Temperature / Temperature_Ref; }
-
   /*!
    * \brief Get the value of the non-dimensionalized thermodynamic pressure.
    * \return Non-dimensionalized thermodynamic pressure.
diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
@@ -96,7 +96,7 @@ const su2double UNIVERSAL_GAS_CONSTANT = 8.3144598;   /*!< \brief Universal gas
 const su2double BOLTZMANN_CONSTANT = 1.3806503E-23;   /*!< \brief Boltzmann's constant [J K^-1] */
 const su2double AVOGAD_CONSTANT = 6.0221415E26; /*!< \brief Avogadro's constant, number of particles in one kmole. */
 const su2double FUND_ELEC_CHARGE_CGS = 4.8032047E-10; /*!< \brief Fundamental electric charge in CGS units, cm^(3/2) g^(1/2) s^(-1). */
-
+const su2double STD_REF_TEMP = 298.15;  /*!< \brief Standard reference temperature for enthalpy in Kelvin. */
 const su2double EPS = 1.0E-16;        /*!< \brief Error scale. */
 const su2double TURB_EPS = 1.0E-16;   /*!< \brief Turbulent Error scale. */
 
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
@@ -1229,9 +1229,6 @@ void CConfig::SetConfig_Options() {
   addDoubleOption("GAMMA_VALUE", Gamma, 1.4);
   /*!\brief THERMODYNAMIC_PRESSURE  \n DESCRIPTION: Thermodynamics(operating) Pressure (101325 Pa), only for incompressible flows) \ingroup Config*/
   addDoubleOption("THERMODYNAMIC_PRESSURE", Pressure_Thermodynamic, 101325.0);
-  /*!\brief STANDARD_REFERENCE_TEMPERATURE  \n DESCRIPTION: Standard reference temperature (298.15K), only for
-   * multicomponent incompressible flows) \ingroup Config*/
-  addDoubleOption("STANDARD_REFERENCE_TEMPERATURE", Standard_Ref_Temperature, 298.15);
   /*!\brief CP_VALUE  \n DESCRIPTION: Specific heat at constant pressure, Cp (1004.703 J/kg*K (air), constant density incompressible fluids only) \ingroup Config*/
   addDoubleListOption("SPECIFIC_HEAT_CP", nSpecific_Heat_Cp, Specific_Heat_Cp);
   /*!\brief THERMAL_EXPANSION_COEFF  \n DESCRIPTION: Thermal expansion coefficient (0.00347 K^-1 (air), used for Boussinesq approximation for liquids/non-ideal gases) \ingroup Config*/
diff --git a/SU2_CFD/include/fluid/CConstantDensity.hpp b/SU2_CFD/include/fluid/CConstantDensity.hpp
@@ -56,16 +56,16 @@ class CConstantDensity final : public CFluidModel {
        decoupled equation. Hence, we update the value.
        Note Cp = Cv, (gamma = 1).*/
     Temperature = t;
-    Enthalpy = Cp * Temperature;
+    Enthalpy = Cp * (Temperature - STD_REF_TEMP);  // Sensible enthalpy relative to STD_REF_TEMP
   }
 
   /*!
    * \brief Set the Dimensionless State using Enthalpy.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - not used here. 
+   * \param[in] val_scalars - not used here.
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
     Enthalpy = val_enthalpy;
-    Temperature = Enthalpy / Cp;
+    Temperature = Enthalpy / Cp + STD_REF_TEMP;  // Temperature from sensible enthalpy
   }
 };
diff --git a/SU2_CFD/include/fluid/CFluidScalar.hpp b/SU2_CFD/include/fluid/CFluidScalar.hpp
@@ -42,7 +42,6 @@ class CFluidScalar final : public CFluidModel {
   const int n_species_mixture;            /*!< \brief Number of species in mixture. */
   su2double Gas_Constant;                 /*!< \brief Specific gas constant. */
   const su2double Pressure_Thermodynamic; /*!< \brief Constant pressure thermodynamic. */
-  const su2double Ref_Temperature;        /*!< \brief Standard Reference temperature, usually set to 298.15 K. */
   const su2double GasConstant_Ref;        /*!< \brief Gas constant reference needed for Nondimensional problems. */
   const su2double Prandtl_Turb_Number;    /*!< \brief Prandlt turbulent number.*/
   const su2double Schmidt_Turb_Number;    /*!< \brief Schmidt turbulent number.*/
@@ -177,7 +176,7 @@ class CFluidScalar final : public CFluidModel {
   /*!
    * \brief Virtual member.
    * \param[in] val_enthalpy - Enthalpy value at the point.
-   * \param[in] val_scalars - Scalar mass fractions. 
+   * \param[in] val_scalars - Scalar mass fractions.
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override;
 };
diff --git a/SU2_CFD/include/fluid/CIncIdealGas.hpp b/SU2_CFD/include/fluid/CIncIdealGas.hpp
@@ -58,7 +58,7 @@ class CIncIdealGas final : public CFluidModel {
     /*--- The EoS only depends upon temperature. ---*/
     Temperature = t;
     Density = Pressure / (Temperature * Gas_Constant);
-    Enthalpy = Cp * Temperature;
+    Enthalpy = Cp * (Temperature - STD_REF_TEMP);  // Sensible enthalpy relative to REF_TEMP
   }
 
   /*!
@@ -67,7 +67,7 @@ class CIncIdealGas final : public CFluidModel {
    */
   void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
     Enthalpy = val_enthalpy;
-    Temperature = Enthalpy / Cp;
+    Temperature = Enthalpy / Cp + STD_REF_TEMP;  // Temperature from sensible enthalpy
     Density = Pressure / (Temperature * Gas_Constant);
   }
 
diff --git a/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp b/SU2_CFD/include/fluid/CIncIdealGasPolynomial.hpp
@@ -58,11 +58,10 @@ class CIncIdealGasPolynomial final : public CFluidModel {
    * \param[in] config - configuration container for the problem.
    */
   void SetCpModel(const CConfig* config) override {
-    const su2double t_ref = config->GetStandard_RefTemperatureND();
     Enthalpy_Ref = 0.0;
     su2double t_i = 1.0;
     for (int i = 0; i < N; ++i) {
-      t_i *= t_ref;
+      t_i *= STD_REF_TEMP;
       coeffs_[i] = config->GetCp_PolyCoeffND(i);
       Enthalpy_Ref += coeffs_[i] * t_i / (i + 1);
     }
diff --git a/SU2_CFD/src/fluid/CFluidScalar.cpp b/SU2_CFD/src/fluid/CFluidScalar.cpp
@@ -45,7 +45,6 @@ CFluidScalar::CFluidScalar(su2double value_pressure_operating, const CConfig* co
     : CFluidModel(),
       n_species_mixture(config->GetnSpecies() + 1),
       Pressure_Thermodynamic(value_pressure_operating),
-      Ref_Temperature(config->GetStandard_RefTemperatureND()),
       GasConstant_Ref(config->GetGas_Constant_Ref()),
       Prandtl_Turb_Number(config->GetPrandtl_Turb()),
       Schmidt_Turb_Number(config->GetSchmidt_Number_Turbulent()),
@@ -219,14 +218,14 @@ su2double CFluidScalar::ComputeEnthalpyFromT(const su2double val_temperature, co
    * depend on temperature, but it does depend on mixture composition, enthalpy is directly computed from
    * the expression h_s = Cp(T - T_ref).
    */
-  su2double val_Enthalpy = Cp * (val_temperature - Ref_Temperature);
+  su2double val_Enthalpy = Cp * (val_temperature - STD_REF_TEMP);
   return val_Enthalpy;
 }
 
 void CFluidScalar::GetEnthalpyDiffusivity(su2double* enthalpy_diffusions) const {
-  const su2double enthalpy_species_N = specificHeat[n_species_mixture - 1] * (Temperature - Ref_Temperature);
+  const su2double enthalpy_species_N = specificHeat[n_species_mixture - 1] * (Temperature - STD_REF_TEMP);
   for (int iVar = 0; iVar < n_species_mixture - 1; iVar++) {
-    const su2double enthalpy_species_i = specificHeat[iVar] * (Temperature - Ref_Temperature);
+    const su2double enthalpy_species_i = specificHeat[iVar] * (Temperature - STD_REF_TEMP);
     enthalpy_diffusions[iVar] = Density * (enthalpy_species_i * massDiffusivity[iVar] -
                                            enthalpy_species_N * massDiffusivity[n_species_mixture - 1]);
     enthalpy_diffusions[iVar] += Mu_Turb * (enthalpy_species_i - enthalpy_species_N) / Schmidt_Turb_Number;
@@ -271,7 +270,7 @@ void CFluidScalar::SetTDState_h(const su2double val_enthalpy, const su2double* v
    * depend on temperature, but it does depend on mixture composition, temperature is directly solved from the
    * expression h_s = Cp(T - T_ref).
    */
-  Temperature = val_enthalpy / Cp + Ref_Temperature;
+  Temperature = val_enthalpy / Cp + STD_REF_TEMP;
   Density = Pressure_Thermodynamic / (Temperature * Gas_Constant);
   Cv = Cp - Gas_Constant;
 
@@ -282,4 +281,5 @@ void CFluidScalar::SetTDState_h(const su2double val_enthalpy, const su2double* v
   }
 
   Kt = WilkeConductivity(val_scalars);
+  ComputeMassDiffusivity();
 }
diff --git a/SU2_CFD/src/output/CFlowIncOutput.cpp b/SU2_CFD/src/output/CFlowIncOutput.cpp
@@ -330,7 +330,8 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
     AddVolumeOutput("LAMINAR_VISCOSITY", "Laminar_Viscosity", "PRIMITIVE", "Laminar viscosity");
     AddVolumeOutput("HEAT_CAPACITY", "Heat_Capacity", "PRIMITIVE", "Heat capacity");
     AddVolumeOutput("THERMAL_CONDUCTIVITY", "Thermal_Conductivity", "PRIMITIVE", "Thermal conductivity");
-    if (heat || flamelet) AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature");
+    if (!weakly_coupled_heat)
+      AddVolumeOutput("TEMPERATURE", "Temperature", "PRIMITIVE", "Temperature");
 
     AddVolumeOutput("SKIN_FRICTION-X", "Skin_Friction_Coefficient_x", "PRIMITIVE", "x-component of the skin friction vector");
     AddVolumeOutput("SKIN_FRICTION-Y", "Skin_Friction_Coefficient_y", "PRIMITIVE", "y-component of the skin friction vector");
@@ -354,7 +355,7 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("RES_VELOCITY-Y", "Residual_Velocity_y", "RESIDUAL", "Residual of the y-velocity component");
   if (nDim == 3)
     AddVolumeOutput("RES_VELOCITY-Z", "Residual_Velocity_z", "RESIDUAL", "Residual of the z-velocity component");
-  if (config->GetEnergy_Equation()){
+  if (heat){
     AddVolumeOutput("RES_ENTHALPY", "Residual_Enthalpy", "RESIDUAL", "Residual of the enthalpy");
   }
   SetVolumeOutputFieldsScalarResidual(config);
@@ -440,7 +441,8 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
     SetVolumeOutputValue("LAMINAR_VISCOSITY", iPoint, Node_Flow->GetLaminarViscosity(iPoint));
     SetVolumeOutputValue("HEAT_CAPACITY", iPoint, Node_Flow->GetSpecificHeatCp(iPoint));
     SetVolumeOutputValue("THERMAL_CONDUCTIVITY", iPoint, Node_Flow->GetThermalConductivity(iPoint));
-    if (heat || flamelet) SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint));
+    if (!weakly_coupled_heat)
+      SetVolumeOutputValue("TEMPERATURE", iPoint, Node_Flow->GetTemperature(iPoint));
   }
 
   SetVolumeOutputValue("RES_PRESSURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 0));
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -2145,9 +2145,9 @@ void CIncEulerSolver::SetPreconditioner(const CConfig *config, unsigned long iPo
      Therefore, we build inv(Precon) here and multiply by the residual
      later in the R-K and Euler Explicit time integration schemes. ---*/
 
-    
+
     Preconditioner[0][0] = Enthalpy * BetaInc2 * dRhodh / Density + BetaInc2;
-    
+
     for (iDim = 0; iDim < nDim; iDim++) Preconditioner[iDim + 1][0] = -1.0 * Velocity[iDim] / Density;
 
     if (energy) {
diff --git a/SU2_CFD/src/solvers/CIncNSSolver.cpp b/SU2_CFD/src/solvers/CIncNSSolver.cpp
@@ -892,14 +892,15 @@ void CIncNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_containe
       SU2_MPI::Allreduce(&globalCounter1, &notConvergedCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
       SU2_MPI::Allreduce(&globalCounter2, &smallYPlusCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
 
-      if (rank == MASTER_NODE) {
+      const bool write_iter = (config->GetTimeIter() % config->GetVolumeOutputFrequency(0) == 0);
+      if (rank == MASTER_NODE && write_iter) {
         if (notConvergedCounter)
           cout << "Warning: Computation of wall coefficients (y+) did not converge in "
                << notConvergedCounter << " points." << endl;
 
         if (smallYPlusCounter)
-          cout << "Warning: y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
-               << " points, for which the wall model is not active." << endl;
+          cout << "y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
+               << " points. No problem, but you can increase your near-wall mesh size." << endl;
       }
     }
     END_SU2_OMP_SAFE_GLOBAL_ACCESS
diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp
@@ -1028,14 +1028,16 @@ void CNSSolver::SetTau_Wall_WF(CGeometry *geometry, CSolver **solver_container,
       SU2_MPI::Allreduce(&globalCounter1, &notConvergedCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
       SU2_MPI::Allreduce(&globalCounter2, &smallYPlusCounter, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
 
-      if (rank == MASTER_NODE) {
+      const bool write_iter = (config->GetInnerIter() % config->GetVolumeOutputFrequency(0) == 0);
+
+      if (rank == MASTER_NODE && write_iter) {
         if (notConvergedCounter)
           cout << "Warning: Computation of wall coefficients (y+) did not converge in "
                << notConvergedCounter << " points." << endl;
 
         if (smallYPlusCounter)
-          cout << "Warning: y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
-               << " points, for which the wall model is not active." << endl;
+          cout << "y+ < " << config->GetwallModel_MinYPlus() << " in " << smallYPlusCounter
+               << " points. No problem, but you can increase your near-wallmesh size." << endl;
       }
     }
     END_SU2_OMP_SAFE_GLOBAL_ACCESS
diff --git a/SU2_CFD/src/solvers/CSpeciesSolver.cpp b/SU2_CFD/src/solvers/CSpeciesSolver.cpp
@@ -207,12 +207,11 @@ void CSpeciesSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfi
 
     const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE);
     const bool energy = config->GetEnergy_Equation();
-    const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET);
     const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat();
 
     /*--- for the flamelet model, the temperature is saved to file, but the energy equation is off ---*/
 
-   if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--;
+   if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--;
 
     /*--- Load data from the restart into correct containers. ---*/
 
@@ -303,6 +302,8 @@ void CSpeciesSolver::Preprocessing(CGeometry* geometry, CSolver** solver_contain
     const su2double* scalar = solver_container[SPECIES_SOL]->GetNodes()->GetSolution(iPoint);
     solver_container[FLOW_SOL]->GetFluidModel()->SetMassDiffusivityModel(config);
     solver_container[FLOW_SOL]->GetFluidModel()->SetTDState_T(temperature, scalar);
+    /*--- Recompute viscosity, important  to get diffusivity correct across MPI ranks. ---*/
+    const su2double viscosity =  solver_container[FLOW_SOL]->GetFluidModel()->GetLaminarViscosity();
     for (auto iVar = 0u; iVar <= nVar; iVar++) {
       const su2double mass_diffusivity = solver_container[FLOW_SOL]->GetFluidModel()->GetMassDiffusivity(iVar);
       nodes->SetDiffusivity(iPoint, mass_diffusivity, iVar);
@@ -579,4 +580,4 @@ void CSpeciesSolver::SetInitialCondition(CGeometry **geometry, CSolver ***solver
   const bool restart = config->GetRestart() || config->GetRestart_Flow();
 
   PushSolutionBackInTime(TimeIter, restart, solver_container, geometry, config);
-}
+}
diff --git a/SU2_CFD/src/solvers/CTurbSolver.cpp b/SU2_CFD/src/solvers/CTurbSolver.cpp
@@ -131,9 +131,8 @@ void CTurbSolver::LoadRestart(CGeometry** geometry, CSolver*** solver, CConfig*
     const bool incompressible = (config->GetKind_Regime() == ENUM_REGIME::INCOMPRESSIBLE);
     const bool energy = config->GetEnergy_Equation();
     const bool weakly_coupled_heat = config->GetWeakly_Coupled_Heat();
-    const bool flamelet = (config->GetKind_FluidModel() == FLUID_FLAMELET);
 
-    if (incompressible && ((!energy) && (!weakly_coupled_heat) && (!flamelet))) skipVars--;
+    if (incompressible && ((!energy) && (!weakly_coupled_heat))) skipVars--;
 
     /*--- Load data from the restart into correct containers. ---*/
 
@@ -252,7 +251,7 @@ void CTurbSolver::ComputeUnderRelaxationFactorHelper(su2double allowableRatio) {
 
   /* Loop over the solution update given by relaxing the linear
    system for this nonlinear iteration. */
-  
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++) {
     su2double localUnderRelaxation = 1.0;
@@ -265,16 +264,16 @@ void CTurbSolver::ComputeUnderRelaxationFactorHelper(su2double allowableRatio) {
       turbulence variables can change over a nonlinear iteration. */
       if (ratio > allowableRatio) {
         localUnderRelaxation = min(allowableRatio / ratio, localUnderRelaxation);
-      
+
       }
     }
 
     /* Threshold the relaxation factor in the event that there is
      a very small value. This helps avoid catastrophic crashes due
      to non-realizable states by canceling the update. */
-    
+
     if (localUnderRelaxation < 1e-10) localUnderRelaxation = 0.0;
-    
+
     /* Store the under-relaxation factor for this point. */
 
     nodes->SetUnderRelaxation(iPoint, localUnderRelaxation);
diff --git a/SU2_CFD/src/variables/CIncEulerVariable.cpp b/SU2_CFD/src/variables/CIncEulerVariable.cpp
@@ -75,11 +75,10 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
   SetPressure(iPoint);
 
-  /*--- Set the value of the enthalpy directly ---*/
-
   su2double Enthalpy = Solution(iPoint, nDim +1);
   FluidModel->SetTDState_h(Enthalpy);
   su2double Temperature = FluidModel->GetTemperature();
+
   const auto check_temp = SetTemperature(iPoint, Temperature, TemperatureLimits);
 
   /*--- Use the fluid model to compute the new value of density.
@@ -103,8 +102,7 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
     Enthalpy = Solution(iPoint, nDim+1);
     FluidModel->SetTDState_h(Enthalpy);
-    Temperature = FluidModel->GetTemperature();
-    SetTemperature(iPoint, Temperature, TemperatureLimits);
+    SetTemperature(iPoint, FluidModel->GetTemperature(), TemperatureLimits);
     SetDensity(iPoint, FluidModel->GetDensity());
 
     /*--- Flag this point as non-physical. ---*/
@@ -117,10 +115,13 @@ bool CIncEulerVariable::SetPrimVar(unsigned long iPoint, CFluidModel *FluidModel
 
   SetVelocity(iPoint);
 
-  /*--- Set enthalpy and specific heats (only necessary for consistency with preconditioning). ---*/
+  /*--- Set specific heats ---*/
 
   SetSpecificHeatCp(iPoint, FluidModel->GetCp());
   SetSpecificHeatCv(iPoint, FluidModel->GetCv());
+
+  /*--- Set enthalpy ---*/
+
   SetEnthalpy(iPoint, FluidModel->GetEnthalpy());
 
   return physical;
diff --git a/SU2_CFD/src/variables/CIncNSVariable.cpp b/SU2_CFD/src/variables/CIncNSVariable.cpp
@@ -60,21 +60,16 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do
 
   SetPressure(iPoint);
 
-  if (Energy) {
-    Enthalpy = Solution(iPoint, nDim + 1);
-    FluidModel->SetTDState_h(Enthalpy, scalar);
-    Temperature = FluidModel->GetTemperature();
-  } else {
-    /*--- When energy equation is switch off, a constant temperature is imposed, and enthalpy is recomputed based on
-     * this temperature. As in the fluid flamelet model, the temperature is retrieved from a look-up table, then the
-     * temperature is obtained directly from the fluidmodel. For the other fluid models, GetTemperature provides the
-     * same value as TemperatureInc ---*/
-    FluidModel->SetTDState_T(TemperatureInc, scalar);
-    Enthalpy = Solution(iPoint, nDim + 1) = FluidModel->GetEnthalpy();
-    Temperature = FluidModel->GetTemperature();
-  }
+  Enthalpy = Solution(iPoint, nDim + 1);
+  FluidModel->SetTDState_h(Enthalpy, scalar);
+  Temperature = FluidModel->GetTemperature();
+
   auto check_temp = SetTemperature(iPoint, Temperature, TemperatureLimits);
 
+  /*--- Use the fluid model to compute the new value of density.
+  Note that the thermodynamic pressure is constant and decoupled
+  from the dynamic pressure being iterated. ---*/
+
   /*--- Set the value of the density ---*/
 
   const auto check_dens = SetDensity(iPoint, FluidModel->GetDensity());
@@ -93,7 +88,6 @@ bool CIncNSVariable::SetPrimVar(unsigned long iPoint, su2double eddy_visc, su2do
     Enthalpy = Solution(iPoint, nDim + 1);
     FluidModel->SetTDState_h(Enthalpy, scalar);
     SetTemperature(iPoint, FluidModel->GetTemperature(), TemperatureLimits);
-
     SetDensity(iPoint, FluidModel->GetDensity());
 
     /*--- Flag this point as non-physical. ---*/
