Fix NASA Polynomial Tests and CConfig Integration

Common/include/CConfig.hpp

@@ -1120,6 +1120,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 {
@@ -2591,6 +2594,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.
@@ -4254,6 +4272,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,172 @@
+/*!
+ * \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 7-coefficient polynomial format for thermodynamic properties:
+ *   Cp/R = a1 + a2*T + a3*T^2 + a4*T^3 + a5*T^4
+ *   H/(R*T) = a1 + a2*T/2 + a3*T^2/3 + a4*T^3/4 + a5*T^4/5 + a6/T
+ *   S/R = a1*ln(T) + a2*T + a3*T^2/2 + a4*T^3/3 + a5*T^4/4 + a7
+ * 
+ * Uses a single temperature range provided via CP_POLYCOEFFS (indices 0-6).
+ */
+template <int N_COEFFS = 7>
+class CIncIdealGasNASA final : public CFluidModel {
+ public:
+  /*!
+   * \brief Constructor of the class.
+   */
+  CIncIdealGasNASA(su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref) {
+    /*--- 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;
+  }
+
+  /*!
+   * \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).
+          Indices 0-4: Cp coefficients (a1-a5)
+          Index 5: Enthalpy constant (a6)
+          Index 6: Entropy constant (a7) ---*/
+    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);
+    Temperature_Max = config->GetTemperatureLimits(1);
+  }
+
+  /*!
+   * \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 {
+    Temperature = t;
+    Density = Pressure / (Temperature * Gas_Constant);
+
+    su2double Cp_over_R = coeffs_[0] + coeffs_[1]*t + coeffs_[2]*t*t + 
+                          coeffs_[3]*t*t*t + coeffs_[4]*t*t*t*t;
+
+    Cp = Cp_over_R * Gas_Constant;
+
+    su2double H_over_RT = coeffs_[0] + coeffs_[1]*t/2.0 + coeffs_[2]*t*t/3.0 + 
+                          coeffs_[3]*t*t*t/4.0 + coeffs_[4]*t*t*t*t/5.0 + coeffs_[5]/t;
+
+    Enthalpy = H_over_RT * Gas_Constant * t;
+    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 {
+    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 */
+
+    su2double Cp_iter = 0.0;
+    su2double delta_temp_iter = 1e10;
+    su2double delta_enthalpy_iter;
+    const int counter_limit = 50;
+    int counter = 0;
+
+    while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) {
+
+      su2double Cp_over_R = coeffs_[0] + coeffs_[1]*temp_iter + coeffs_[2]*temp_iter*temp_iter + 
+                            coeffs_[3]*temp_iter*temp_iter*temp_iter + 
+                            coeffs_[4]*temp_iter*temp_iter*temp_iter*temp_iter;
+
+      Cp_iter = Cp_over_R * Gas_Constant;
+
+      su2double H_over_RT = coeffs_[0] + coeffs_[1]*temp_iter/2.0 + coeffs_[2]*temp_iter*temp_iter/3.0 + 
+                            coeffs_[3]*temp_iter*temp_iter*temp_iter/4.0 + 
+                            coeffs_[4]*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + coeffs_[5]/temp_iter;
+
+      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;
+
+      if (temp_iter < Temperature_Min) {
+        // Clamp to min but don't break immediately to allow recovery if simple overshoot
+         if (abs(delta_temp_iter) < toll) break; 
+         temp_iter = Temperature_Min;
+      }
+      if (temp_iter > Temperature_Max) {
+         if (abs(delta_temp_iter) < toll) break;
+         temp_iter = Temperature_Max;
+      }
+    }
+
+    Temperature = temp_iter;
+    Cp = Cp_iter;
+
+    if (counter == counter_limit) {
+      if (SU2_MPI::GetRank() == MASTER_NODE)
+          std::cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation (NASA Model)." << std::endl;
+    }
+
+    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. */
+
+  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 (requires 7 coefficients) ---*/
+        auxFluidModel = new CIncIdealGasNASA<7>(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<7>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / 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);