Fix NASA Polynomial Tests and CConfig Integration

Common/include/CConfig.hpp

@@ -1121,6 +1121,9 @@ class CConfig {
   array<su2double, N_POLY_COEFFS> cp_polycoeffs{{0.0}};  /*!< \brief Array for specific heat polynomial coefficients. */
   array<su2double, N_POLY_COEFFS> mu_polycoeffs{{0.0}};  /*!< \brief Array for viscosity polynomial coefficients. */
   array<su2double, N_POLY_COEFFS> kt_polycoeffs{{0.0}};  /*!< \brief Array for thermal conductivity polynomial coefficients. */
+
+  bool Cp_NASA_Format{false}; /*!< \brief Flag to use NASA polynomial format for Cp. */
+
   bool Body_Force;                      /*!< \brief Flag to know if a body force is included in the formulation. */
 
   struct CMUSCLRampParam {
@@ -2592,6 +2595,21 @@ class CConfig {
    */
   void SetEnergy_Ref(su2double val_energy_ref) { Energy_Ref = val_energy_ref; }
 
+  /*!
+   * \brief Get the flag for using NASA polynomial format for Cp.
+   * \return True if NASA format is used.
+   */
+  bool GetCp_NASA_Format(void) const { return Cp_NASA_Format; }
+
+  /*!
+   * \brief Set the flag for using NASA polynomial format for Cp.
+   * \param[in] val - Value of the flag.
+   */
+  void SetCp_NASA_Format(bool val) { Cp_NASA_Format = val; }
+
+  // Setter for TemperatureLimits (for testing)
+  void SetTemperatureLimits(unsigned short val_index, su2double val) { TemperatureLimits[val_index] = val; }
+
   /*!
    * \brief Set the reference Omega for nondimensionalization.
    * \param[in] val_omega_ref - Value of the reference omega.
@@ -4255,6 +4273,7 @@ class CConfig {
    * \param[in] val_index - Index of the array with all polynomial coefficients.
    */
   void SetCp_PolyCoeffND(su2double val_coeff, unsigned short val_index) { CpPolyCoefficientsND[val_index] = val_coeff; }
+  void SetCp_PolyCoeff(unsigned short val_index, su2double val) { cp_polycoeffs[val_index] = val; }
 
   /*!
    * \brief Set the temperature polynomial coefficient for viscosity.

Common/include/option_structure.hpp

@@ -201,7 +201,7 @@ constexpr passivedouble COLORING_EFF_THRESH = 0.875;  /*!< \brief Below this val
    assume a quartic form (5 coefficients). For example,
    Cp(T) = b0 + b1*T + b2*T^2 + b3*T^3 + b4*T^4. By default, all coeffs
    are set to zero and will be properly non-dim. in the solver. ---*/
-constexpr int N_POLY_COEFFS = 5; /*!< \brief Number of coefficients in temperature polynomial fits for fluid models. */
+constexpr int N_POLY_COEFFS = 10; /*!< \brief Number of coefficients in temperature polynomial fits for fluid models. */
 
 /*!
  * \brief Boolean answers

Common/src/CConfig.cpp

@@ -1333,11 +1333,14 @@ void CConfig::SetConfig_Options() {
   /*--- Options related to temperature polynomial coefficients for fluid models. ---*/
 
   /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */
-  addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, false, cp_polycoeffs.data());
+  addDoubleArrayOption("CP_POLYCOEFFS", N_POLY_COEFFS, true, cp_polycoeffs.data());
   /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */
-  addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, false, mu_polycoeffs.data());
+  addDoubleArrayOption("MU_POLYCOEFFS", N_POLY_COEFFS, true, mu_polycoeffs.data());
   /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */
-  addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data());
+  addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, true, kt_polycoeffs.data());
+
+  /* DESCRIPTION: Flag to use NASA polynomial format for Cp in incompressible ideal gas model. */
+  addBoolOption("CP_NASA", Cp_NASA_Format, false);
 
   /*!\brief REYNOLDS_NUMBER \n DESCRIPTION: Reynolds number (non-dimensional, based on the free-stream values). Needed for viscous solvers. For incompressible solvers the Reynolds length will always be 1.0 \n DEFAULT: 0.0 \ingroup Config */
   addDoubleOption("REYNOLDS_NUMBER", Reynolds, 0.0);

SU2_CFD/include/fluid/CIncIdealGasNASA.hpp

@@ -0,0 +1,223 @@
+/*!
+ * \file CIncIdealGasNASA.hpp
+ * \brief Defines the incompressible Ideal Gas model with NASA polynomials for Cp.
+ * \author Pratyksh Gupta
+ * \version 8.4.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 <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <array>
+#include <iostream>
+
+#include "CFluidModel.hpp"
+
+/*!
+ * \class CIncIdealGasNASA
+ * \brief Child class for defining an incompressible ideal gas model with NASA polynomials.
+ * \author Pratyksh Gupta
+ * 
+ * Implements NASA 9-coefficient polynomial format for thermodynamic properties with backward compatibility for NASA-7.
+ * Ref: McBride, B.J., Zehe, M.J., and Gordon, S., "NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species", NASA/TP-2002-211556, 2002.
+ * 
+ * NASA-9 format (full):
+ *   Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
+ *   H/(RT) = -a1*T^-2 + a2*ln(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T
+ *   S/R = -a1*T^-2/2 - a2*T^-1 + a3*ln(T) + a4*T + a5*T^2/2 + a6*T^3/3 + a7*T^4/4 + a9
+ * 
+ * NASA-7 format (subset): Set a1=a2=0 to recover the traditional 7-coefficient format.
+ * 
+ * Uses a single temperature range provided via CP_POLYCOEFFS (indices 0-8).
+ */
+template <int N_COEFFS = 9>
+class CIncIdealGasNASA final : public CFluidModel {
+ public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CIncIdealGasNASA(su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref, su2double val_Ref_Temp_Dim = 1.0) {
+    /*--- In the incompressible ideal gas model, the thermodynamic pressure is decoupled
+    from the governing equations and held constant. The density is therefore only a
+    function of temperature variations. The gas is incompressible, so Cp = Cv (gamma = 1). ---*/
+    Gas_Constant = val_gas_constant;
+    Pressure = val_operating_pressure;
+    Gamma = 1.0;
+    Std_Ref_Temp_ND = val_Temperature_Ref;
+    Ref_Temp_Dim = val_Ref_Temp_Dim;
+  }
+
+  /*!
+   * \brief Set the NASA polynomial coefficients for variable Cp.
+   * \param[in] config - configuration container for the problem.
+   */
+  void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override {
+
+    /*--- Read NASA coefficients from the standard polynomial coefficient array (CP_POLYCOEFFS).
+          NASA-9 format uses indices 0-8:
+          Indices 0-1: Inverse temperature terms (a1*T^-2, a2*T^-1)
+          Indices 2-6: Polynomial terms (a3, a4*T, a5*T^2, a6*T^3, a7*T^4)
+          Index 7: Enthalpy constant (a8)
+          Index 8: Entropy constant (a9)
+
+          For NASA-7 compatibility: Set indices 0-1 to zero. ---*/
+    for (int i = 0; i < N_COEFFS; ++i) {
+      if (i < config->GetnPolyCoeffs()) {
+        coeffs_[i] = config->GetCp_PolyCoeff(i);
+      } else {
+        coeffs_[i] = 0.0; 
+      }
+    }
+
+    Temperature_Min = config->GetTemperatureLimits(0) / Ref_Temp_Dim;
+    Temperature_Max = config->GetTemperatureLimits(1) / Ref_Temp_Dim;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using Temperature.
+   * \param[in] t - Temperature value at the point.
+   */
+  void SetTDState_T(su2double t, const su2double *val_scalars = nullptr) override {
+    using std::pow;
+    using std::log;
+    Temperature = t;
+    Density = Pressure / (Temperature * Gas_Constant);
+
+    const su2double a1 = coeffs_[0];
+    const su2double a2 = coeffs_[1];
+    const su2double a3 = coeffs_[2];
+    const su2double a4 = coeffs_[3];
+    const su2double a5 = coeffs_[4];
+    const su2double a6 = coeffs_[5];
+    const su2double a7 = coeffs_[6];
+    const su2double a8 = coeffs_[7];
+    const su2double a9 = coeffs_[8];
+
+    // Convert to dimensional temperature for polynomial evaluation (NASA coeffs expect Kelvin)
+    const su2double T_dim = t * Ref_Temp_Dim;
+
+    // NASA-9: Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
+    su2double Cp_over_R = a1 * pow(T_dim, -2.0) + a2 * pow(T_dim, -1.0) + a3 + a4 * T_dim + a5 * pow(T_dim, 2.0) + a6 * pow(T_dim, 3.0) + a7 * pow(T_dim, 4.0);
+
+    Cp = Cp_over_R * Gas_Constant;
+
+    // NASA-9: H/(RT) = -a1*T^-2 + a2*ln(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T
+    su2double H_over_RT = -a1 * pow(T_dim, -2.0) + a2 * log(T_dim) / T_dim + a3 + a4 * T_dim / 2.0 + a5 * pow(T_dim, 2.0) / 3.0 + 
+                          a6 * pow(T_dim, 3.0) / 4.0 + a7 * pow(T_dim, 4.0) / 5.0 + a8 / T_dim;
+
+    Enthalpy = H_over_RT * Gas_Constant * t;
+
+    // NASA-9: S/R = -a1*T^-2/2 - a2*T^-1 + a3*ln(T) + a4*T + a5*T^2/2 + a6*T^3/3 + a7*T^4/4 + a9
+    su2double S_over_R = -a1 * pow(T_dim, -2.0) / 2.0 - a2 * pow(T_dim, -1.0) + a3 * log(T_dim) + a4 * T_dim + 
+                         a5 * pow(T_dim, 2.0) / 2.0 + a6 * pow(T_dim, 3.0) / 3.0 + a7 * pow(T_dim, 4.0) / 4.0 + a9;
+
+    Entropy = S_over_R * Gas_Constant;
+
+    Cv = Cp / Gamma;
+  }
+
+  /*!
+   * \brief Set the Dimensionless State using enthalpy.
+   * \param[in] val_enthalpy - Enthalpy value at the point.
+   */
+  void SetTDState_h(su2double val_enthalpy, const su2double* val_scalars = nullptr) override {
+    using std::pow;
+    using std::log;
+    using std::fmin;
+    using std::fmax;
+    using std::fabs;
+    Enthalpy = val_enthalpy;
+
+    const su2double toll = 1e-5;
+    su2double temp_iter = (Temperature_Min + Temperature_Max) / 2.0; /* Start in middle of allowed range */
+    if (temp_iter < 1.0) temp_iter = 300.0; /* Fallback if limits are not set or zero */
+
+    const su2double a1 = coeffs_[0];
+    const su2double a2 = coeffs_[1];
+    const su2double a3 = coeffs_[2];
+    const su2double a4 = coeffs_[3];
+    const su2double a5 = coeffs_[4];
+    const su2double a6 = coeffs_[5];
+    const su2double a7 = coeffs_[6];
+    const su2double a8 = coeffs_[7];
+    const su2double a9 = coeffs_[8];
+
+    su2double Cp_iter = 0.0;
+    su2double delta_temp_iter = 1e10;
+    su2double delta_enthalpy_iter;
+    const int counter_limit = 50;
+    int counter = 0;
+
+    while ((fabs(delta_temp_iter) > toll) && (counter++ < counter_limit)) {
+
+      const su2double T_dim = temp_iter * Ref_Temp_Dim;
+
+      // NASA-9: Cp/R = a1*T^-2 + a2*T^-1 + a3 + a4*T + a5*T^2 + a6*T^3 + a7*T^4
+      su2double Cp_over_R = a1 * pow(T_dim, -2.0) + a2 * pow(T_dim, -1.0) + a3 + a4 * T_dim + a5 * pow(T_dim, 2.0) + 
+                            a6 * pow(T_dim, 3.0) + a7 * pow(T_dim, 4.0);
+
+      Cp_iter = Cp_over_R * Gas_Constant;
+
+      // NASA-9: H/(RT) = -a1*T^-2 + a2*ln(T)/T + a3 + a4*T/2 + a5*T^2/3 + a6*T^3/4 + a7*T^4/5 + a8/T
+      su2double H_over_RT = -a1 * pow(T_dim, -2.0) + a2 * log(T_dim) / T_dim + a3 + a4 * T_dim / 2.0 + 
+                            a5 * pow(T_dim, 2.0) / 3.0 + a6 * pow(T_dim, 3.0) / 4.0 + 
+                            a7 * pow(T_dim, 4.0) / 5.0 + a8 / T_dim;
+
+      su2double Enthalpy_iter = H_over_RT * Gas_Constant * temp_iter;
+
+      delta_enthalpy_iter = Enthalpy - Enthalpy_iter;
+      delta_temp_iter = delta_enthalpy_iter / Cp_iter;
+      temp_iter += delta_temp_iter;
+
+      temp_iter = fmin(fmax(Temperature_Min, temp_iter), Temperature_Max);
+    }
+
+    Temperature = temp_iter;
+    Cp = Cp_iter;
+
+    // Evaluate Entropy at final state
+    const su2double T_dim_final = Temperature * Ref_Temp_Dim;
+    su2double S_over_R = -a1 * pow(T_dim_final, -2.0) / 2.0 - a2 * pow(T_dim_final, -1.0) + a3 * log(T_dim_final) + 
+                         a4 * T_dim_final + a5 * pow(T_dim_final, 2.0) / 2.0 + a6 * pow(T_dim_final, 3.0) / 3.0 + 
+                         a7 * pow(T_dim_final, 4.0) / 4.0 + a9;
+    Entropy = S_over_R * Gas_Constant;
+
+    if (counter == counter_limit) {
+      if (SU2_MPI::GetRank() == MASTER_NODE)
+          std::cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation (NASA Model).\n";
+    }
+
+    Density = Pressure / (Temperature * Gas_Constant);
+    Cv = Cp / Gamma;
+  }
+
+ private:
+  su2double Gas_Constant{0.0};     /*!< \brief Specific Gas Constant. */
+  su2double Gamma{0.0};            /*!< \brief Ratio of specific heats. */
+  su2double Std_Ref_Temp_ND{0.0};  /*!< \brief Nondimensional standard reference temperature for enthalpy. */
+  su2double Ref_Temp_Dim{1.0};     /*!< \brief Dimensional reference temperature for evaluating polynomials. */
+
+  std::array<su2double, N_COEFFS> coeffs_;  /*!< \brief NASA polynomial coefficients. */
+
+  su2double Temperature_Min{0.0};  /*!< \brief Minimum temperature allowed. */
+  su2double Temperature_Max{1e10}; /*!< \brief Maximum temperature allowed. */
+};

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -30,6 +30,7 @@
 #include "../../include/fluid/CConstantDensity.hpp"
 #include "../../include/fluid/CIncIdealGas.hpp"
 #include "../../include/fluid/CIncIdealGasPolynomial.hpp"
+#include "../../include/fluid/CIncIdealGasNASA.hpp"
 #include "../../include/variables/CIncNSVariable.hpp"
 #include "../../../Common/include/toolboxes/geometry_toolbox.hpp"
 #include "../../include/fluid/CFluidScalar.hpp"
@@ -307,11 +308,21 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
 
       config->SetGas_Constant(UNIVERSAL_GAS_CONSTANT/(config->GetMolecular_Weight()/1000.0));
       Pressure_Thermodynamic = Density_FreeStream*Temperature_FreeStream*config->GetGas_Constant();
-      auxFluidModel = new CIncIdealGasPolynomial<N_POLY_COEFFS>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
+      if (config->GetCp_NASA_Format()) {
+        /*--- Use the NASA polynomial model (supports NASA-9 format with 9 coefficients, backward compatible with NASA-7) ---*/
+        auxFluidModel = new CIncIdealGasNASA<9>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
+      } else {
+        /*--- Use the standard polynomial model ---*/
+        auxFluidModel = new CIncIdealGasPolynomial<N_POLY_COEFFS>(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
+      }
+
       if (viscous) {
         /*--- Variable Cp model via polynomial. ---*/
-        for (iVar = 0; iVar < config->GetnPolyCoeffs(); iVar++)
-          config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar), iVar);
+        if (!config->GetCp_NASA_Format()) {
+             for (iVar = 0; iVar < config->GetnPolyCoeffs(); iVar++)
+                config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar), iVar);
+        }
+        /*--- Both models refer to config for coeffs, so setup is handled inside SetCpModel ---*/
         auxFluidModel->SetCpModel(config, Temperature_FreeStream);
       }
       auxFluidModel->SetTDState_T(Temperature_FreeStream);
@@ -499,12 +510,19 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
         break;
 
       case INC_IDEAL_GAS_POLY:
-        fluidModel = new CIncIdealGasPolynomial<N_POLY_COEFFS>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref());
+        if (config->GetCp_NASA_Format()) {
+             fluidModel = new CIncIdealGasNASA<9>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref(), config->GetTemperature_Ref());
+        } else {
+             fluidModel = new CIncIdealGasPolynomial<N_POLY_COEFFS>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref());
+        }
+
         if (viscous) {
           /*--- Variable Cp model via polynomial. ---*/
-          config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(0)/Gas_Constant_Ref, 0);
-          for (iVar = 1; iVar < config->GetnPolyCoeffs(); iVar++)
-            config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar)*pow(Temperature_Ref,iVar)/Gas_Constant_Ref, iVar);
+          if (!config->GetCp_NASA_Format()) {
+              config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(0)/Gas_Constant_Ref, 0);
+              for (iVar = 1; iVar < config->GetnPolyCoeffs(); iVar++)
+                config->SetCp_PolyCoeffND(config->GetCp_PolyCoeff(iVar)*pow(Temperature_Ref,iVar)/Gas_Constant_Ref, iVar);
+          }
           fluidModel->SetCpModel(config, Temperature_FreeStreamND);
         }
         fluidModel->SetTDState_T(Temperature_FreeStreamND);