From 87b3fdc103100b75a380d5e76e5d41d490148af8 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Wed, 14 Jan 2026 12:35:55 +0530 Subject: [PATCH 01/11] Fix NASA Polynomial Tests and CConfig Integration --- Common/include/CConfig.hpp | 49 +++++ Common/src/CConfig.cpp | 11 + SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 178 ++++++++++++++++ .../SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp | 196 ++++++++++++++++++ UnitTests/meson.build | 1 + 5 files changed, 435 insertions(+) create mode 100644 SU2_CFD/include/fluid/CIncIdealGasNASA.hpp create mode 100644 UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index f198df434761..0c553e9cae30 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1120,6 +1120,13 @@ class CConfig { array cp_polycoeffs{{0.0}}; /*!< \brief Array for specific heat polynomial coefficients. */ array mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ array kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ + + array NASA_Coeffs_Low{{0.0}}; /*!< \brief NASA polynomial coefficients (low temp range). */ + array NASA_Coeffs_High{{0.0}}; /*!< \brief NASA polynomial coefficients (high temp range). */ + su2double NASA_Temp_Low; /*!< \brief NASA polynomial low temperature bound. */ + su2double NASA_Temp_Mid; /*!< \brief NASA polynomial mid temperature break. */ + su2double NASA_Temp_High; /*!< \brief NASA polynomial high temperature bound. */ + bool Body_Force; /*!< \brief Flag to know if a body force is included in the formulation. */ struct CMUSCLRampParam { @@ -2591,6 +2598,48 @@ class CConfig { */ void SetEnergy_Ref(su2double val_energy_ref) { Energy_Ref = val_energy_ref; } + /*! + * \brief Get the value of the NASA polynomial coefficient (low temp range). + * \param[in] val_index - Index of the coefficient. + * \return Value of the coefficient. + */ + su2double GetNASA_CoeffLowND(unsigned short val_index) const { return NASA_Coeffs_Low[val_index]; } + + /*! + * \brief Get the value of the NASA polynomial coefficient (high temp range). + * \param[in] val_index - Index of the coefficient. + * \return Value of the coefficient. + */ + su2double GetNASA_CoeffHighND(unsigned short val_index) const { return NASA_Coeffs_High[val_index]; } + + /*! + * \brief Get the value of the NASA polynomial low temperature bound. + * \return Value of the temperature. + */ + su2double GetNASA_TempLow() const { return NASA_Temp_Low; } + + /*! + * \brief Get the value of the NASA polynomial mid temperature break. + * \return Value of the temperature. + */ + su2double GetNASA_TempMid() const { return NASA_Temp_Mid; } + + /*! + * \brief Get the value of the NASA polynomial high temperature bound. + * \return Value of the temperature. + */ + su2double GetNASA_TempHigh() const { return NASA_Temp_High; } + + // Setters for NASA polynomials (useful for testing) + void SetNASA_CoeffLowND(unsigned short val_index, su2double val) { NASA_Coeffs_Low[val_index] = val; } + void SetNASA_CoeffHighND(unsigned short val_index, su2double val) { NASA_Coeffs_High[val_index] = val; } + void SetNASA_TempLow(su2double val) { NASA_Temp_Low = val; } + void SetNASA_TempMid(su2double val) { NASA_Temp_Mid = val; } + void SetNASA_TempHigh(su2double val) { NASA_Temp_High = 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. diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index a9ecb3a94b24..fca77dd55bc9 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1339,6 +1339,17 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data()); + /* DESCRIPTION: NASA polynomial coefficients (low temp range) */ + addDoubleArrayOption("NASA_POLYCOEFFS_LOW", 7, false, NASA_Coeffs_Low.data()); + /* DESCRIPTION: NASA polynomial coefficients (high temp range) */ + addDoubleArrayOption("NASA_POLYCOEFFS_HIGH", 7, false, NASA_Coeffs_High.data()); + /* DESCRIPTION: NASA polynomial low temperature bound */ + addDoubleOption("NASA_TEMP_LOW", NASA_Temp_Low, 0.0); + /* DESCRIPTION: NASA polynomial mid temperature break */ + addDoubleOption("NASA_TEMP_MID", NASA_Temp_Mid, 0.0); + /* DESCRIPTION: NASA polynomial high temperature bound */ + addDoubleOption("NASA_TEMP_HIGH", NASA_Temp_High, 0.0); + /*!\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); /*!\brief REYNOLDS_LENGTH \n DESCRIPTION: Reynolds length (1 m by default). Used for compressible solver: incompressible solver will use 1.0. \ingroup Config */ diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp new file mode 100644 index 000000000000..11243cf7b522 --- /dev/null +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -0,0 +1,178 @@ +/*! + * \file CIncIdealGasNASA.hpp + * \brief Defines the incompressible Ideal Gas model with NASA polynomials for Cp. + * \author Pratyksh Gupta (based on CIncIdealGasPolynomial by T. Economon) + * \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 . + */ + +#pragma once + +#include +#include + +#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 + * + * Supports dual temperature ranges (low/high) with separate coefficients. + */ +template +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 { + T_mid = config->GetNASA_TempMid(); + T_low = config->GetNASA_TempLow(); + T_high = config->GetNASA_TempHigh(); + + for (int i = 0; i < N_COEFFS; ++i) { + coeffs_low_[i] = config->GetNASA_CoeffLowND(i); + coeffs_high_[i] = config->GetNASA_CoeffHighND(i); + } + + Temperature_Min = config->GetTemperatureLimits(0); + Temperature_Max = config->GetTemperatureLimits(1); + + if (T_mid <= T_low || T_high <= T_mid) { + SU2_MPI::Error("Invalid NASA polynomial temperature ranges. Must have T_low < T_mid < T_high.", + CURRENT_FUNCTION); + } + } + + /*! + * \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); + + const su2double* coeffs = (t < T_mid) ? coeffs_low_.data() : coeffs_high_.data(); + + 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 = (T_low + T_high) / 2.0; + 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)) { + const su2double* coeffs = (temp_iter < T_mid) ? coeffs_low_.data() : coeffs_high_.data(); + + 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) { + cout << "Warning: Negative temperature found during Newton-Raphson." << endl; + temp_iter = Temperature_Min; + break; + } + if (temp_iter > Temperature_Max) { + cout << "Warning: Temperature exceeds maximum during Newton-Raphson." << endl; + temp_iter = Temperature_Max; + break; + } + } + + Temperature = temp_iter; + Cp = Cp_iter; + + if (counter == counter_limit) { + cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation." << 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 coeffs_low_; /*!< \brief NASA polynomial coefficients for low temperature range. */ + std::array coeffs_high_; /*!< \brief NASA polynomial coefficients for high temperature range. */ + + su2double T_low{200.0}; /*!< \brief Lower temperature bound. */ + su2double T_mid{1000.0}; /*!< \brief Mid-point temperature (split between low/high ranges). */ + su2double T_high{6000.0}; /*!< \brief Upper temperature bound. */ + + su2double Temperature_Min{0.0}; /*!< \brief Minimum temperature allowed in Newton-Raphson iterations. */ + su2double Temperature_Max{1e10}; /*!< \brief Maximum temperature allowed in Newton-Raphson iterations. */ +}; diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp new file mode 100644 index 000000000000..f1a43f2b43a0 --- /dev/null +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -0,0 +1,196 @@ +/*! + * \file CIncIdealGasNASA_tests.cpp + * \brief Unit tests for the NASA polynomial fluid model class. + * \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 . + */ + +#include "catch.hpp" +#include +#include +#include "../../../SU2_CFD/include/fluid/CIncIdealGasNASA.hpp" + +// Helper to config an already constructed CConfig object +void ConfigureNASA(CConfig& config) { + su2double nasa_low[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372}; + su2double nasa_high[7] = {2.92664, 1.4879768e-3, -5.68476e-7, 1.0097038e-10, -6.753351e-15, -922.7977, 5.980528}; + + for(int i=0; i<7; ++i) { + config.SetNASA_CoeffLowND(i, nasa_low[i]); + config.SetNASA_CoeffHighND(i, nasa_high[i]); + } + config.SetNASA_TempLow(200.0); + config.SetNASA_TempMid(1000.0); + config.SetNASA_TempHigh(6000.0); + // Ensure we set the temperature limits that SetCpModel reads + config.SetTemperatureLimits(0, 200.0); + config.SetTemperatureLimits(1, 6000.0); +} + +// Helper to populate the stringstream with minimal required options +void SetupConfigStream(std::stringstream& ss) { + ss << "SOLVER= EULER" << std::endl; + ss << "MATH_PROBLEM= DIRECT" << std::endl; + ss << "KIND_TURB_MODEL= NONE" << std::endl; + ss << "KIND_TRANS_MODEL= NONE" << std::endl; +} + +TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { + const su2double R_N2 = 296.8; + const su2double P = 101325.0; + const su2double T_ref = 298.15; + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); + nasa_model.SetCpModel(&config, T_ref); + + SECTION("Low temperature range (300K)") { + nasa_model.SetTDState_T(300.0); + const su2double cp = nasa_model.GetCp(); + CHECK(cp == Approx(1040.0).epsilon(0.01)); + } + + SECTION("High temperature range (1500K)") { + nasa_model.SetTDState_T(1500.0); + const su2double cp = nasa_model.GetCp(); + CHECK(cp == Approx(1240.0).epsilon(0.01)); + } + + SECTION("Temperature range transition") { + nasa_model.SetTDState_T(999.0); + const su2double cp_below = nasa_model.GetCp(); + + nasa_model.SetTDState_T(1001.0); + const su2double cp_above = nasa_model.GetCp(); + + const su2double discontinuity = std::abs(cp_above - cp_below); + CHECK(discontinuity < 50.0); + } +} + +TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]") { + const su2double R_N2 = 296.8; + const su2double P = 101325.0; + const su2double T_ref = 298.15; + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); + nasa_model.SetCpModel(&config, T_ref); + + SECTION("Round-trip T->H->T at 500K") { + const su2double T_target = 500.0; + + nasa_model.SetTDState_T(T_target); + const su2double H = nasa_model.GetEnthalpy(); + + nasa_model.SetTDState_h(H); + const su2double T_recovered = nasa_model.GetTemperature(); + + CHECK(T_recovered == Approx(T_target).epsilon(1e-6)); + } + + SECTION("Multiple temperature points") { + const su2double test_temps[] = {250.0, 500.0, 1000.0, 2000.0, 4000.0}; + + for (const auto T_in : test_temps) { + nasa_model.SetTDState_T(T_in); + const su2double H = nasa_model.GetEnthalpy(); + + nasa_model.SetTDState_h(H); + const su2double T_out = nasa_model.GetTemperature(); + + CHECK(T_out == Approx(T_in).epsilon(1e-6)); + } + } +} + +TEST_CASE("NASA polynomial boundary conditions", "[CIncIdealGasNASA]") { + const su2double R_N2 = 296.8; + const su2double P = 101325.0; + const su2double T_ref = 298.15; + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); + nasa_model.SetCpModel(&config, T_ref); + + SECTION("Lower temperature bound (200K)") { + nasa_model.SetTDState_T(200.0); + const su2double cp = nasa_model.GetCp(); + + CHECK(cp > 0.0); + CHECK(cp < 2000.0); + } + + SECTION("Upper temperature bound (6000K)") { + nasa_model.SetTDState_T(6000.0); + const su2double cp = nasa_model.GetCp(); + + CHECK(cp > 0.0); + CHECK(cp < 3000.0); + } +} + +TEST_CASE("NASA polynomial thermodynamic consistency", "[CIncIdealGasNASA]") { + const su2double R_N2 = 296.8; + const su2double P = 101325.0; + const su2double T_ref = 298.15; + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); + nasa_model.SetCpModel(&config, T_ref); + + SECTION("Density from ideal gas law") { + const su2double T = 300.0; + nasa_model.SetTDState_T(T); + + const su2double rho = nasa_model.GetDensity(); + const su2double rho_expected = P / (R_N2 * T); + + CHECK(rho == Approx(rho_expected).epsilon(1e-10)); + } + + SECTION("Cv from Cp (incompressible: gamma = 1)") { + nasa_model.SetTDState_T(500.0); + + const su2double cp = nasa_model.GetCp(); + const su2double cv = nasa_model.GetCv(); + + CHECK(cv == Approx(cp).epsilon(1e-10)); + } +} diff --git a/UnitTests/meson.build b/UnitTests/meson.build index 3e64cfb9aa3d..6b92257f0b0f 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -15,6 +15,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'Common/toolboxes/multilayer_perceptron/CLookUp_ANN_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', 'SU2_CFD/fluid/CFluidModel_tests.cpp', + 'SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/windowing.cpp']) From 238142de967070dba254a6b273c18784dcb32011 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Wed, 14 Jan 2026 14:03:46 +0530 Subject: [PATCH 02/11] Apply pre-commit style fixes --- Common/include/CConfig.hpp | 6 +- .../SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp | 88 +++++++++---------- 2 files changed, 47 insertions(+), 47 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index 0c553e9cae30..ff6953989cc6 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1120,7 +1120,7 @@ class CConfig { array cp_polycoeffs{{0.0}}; /*!< \brief Array for specific heat polynomial coefficients. */ array mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ array kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ - + array NASA_Coeffs_Low{{0.0}}; /*!< \brief NASA polynomial coefficients (low temp range). */ array NASA_Coeffs_High{{0.0}}; /*!< \brief NASA polynomial coefficients (high temp range). */ su2double NASA_Temp_Low; /*!< \brief NASA polynomial low temperature bound. */ @@ -2629,14 +2629,14 @@ class CConfig { * \return Value of the temperature. */ su2double GetNASA_TempHigh() const { return NASA_Temp_High; } - + // Setters for NASA polynomials (useful for testing) void SetNASA_CoeffLowND(unsigned short val_index, su2double val) { NASA_Coeffs_Low[val_index] = val; } void SetNASA_CoeffHighND(unsigned short val_index, su2double val) { NASA_Coeffs_High[val_index] = val; } void SetNASA_TempLow(su2double val) { NASA_Temp_Low = val; } void SetNASA_TempMid(su2double val) { NASA_Temp_Mid = val; } void SetNASA_TempHigh(su2double val) { NASA_Temp_High = val; } - + // Setter for TemperatureLimits (for testing) void SetTemperatureLimits(unsigned short val_index, su2double val) { TemperatureLimits[val_index] = val; } diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp index f1a43f2b43a0..83b474787fc5 100644 --- a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -32,34 +32,34 @@ // Helper to config an already constructed CConfig object void ConfigureNASA(CConfig& config) { - su2double nasa_low[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372}; - su2double nasa_high[7] = {2.92664, 1.4879768e-3, -5.68476e-7, 1.0097038e-10, -6.753351e-15, -922.7977, 5.980528}; - - for(int i=0; i<7; ++i) { - config.SetNASA_CoeffLowND(i, nasa_low[i]); - config.SetNASA_CoeffHighND(i, nasa_high[i]); - } - config.SetNASA_TempLow(200.0); - config.SetNASA_TempMid(1000.0); - config.SetNASA_TempHigh(6000.0); - // Ensure we set the temperature limits that SetCpModel reads - config.SetTemperatureLimits(0, 200.0); - config.SetTemperatureLimits(1, 6000.0); + su2double nasa_low[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372}; + su2double nasa_high[7] = {2.92664, 1.4879768e-3, -5.68476e-7, 1.0097038e-10, -6.753351e-15, -922.7977, 5.980528}; + + for (int i = 0; i < 7; ++i) { + config.SetNASA_CoeffLowND(i, nasa_low[i]); + config.SetNASA_CoeffHighND(i, nasa_high[i]); + } + config.SetNASA_TempLow(200.0); + config.SetNASA_TempMid(1000.0); + config.SetNASA_TempHigh(6000.0); + // Ensure we set the temperature limits that SetCpModel reads + config.SetTemperatureLimits(0, 200.0); + config.SetTemperatureLimits(1, 6000.0); } // Helper to populate the stringstream with minimal required options void SetupConfigStream(std::stringstream& ss) { - ss << "SOLVER= EULER" << std::endl; - ss << "MATH_PROBLEM= DIRECT" << std::endl; - ss << "KIND_TURB_MODEL= NONE" << std::endl; - ss << "KIND_TRANS_MODEL= NONE" << std::endl; + ss << "SOLVER= EULER" << std::endl; + ss << "MATH_PROBLEM= DIRECT" << std::endl; + ss << "KIND_TURB_MODEL= NONE" << std::endl; + ss << "KIND_TRANS_MODEL= NONE" << std::endl; } TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { const su2double R_N2 = 296.8; const su2double P = 101325.0; const su2double T_ref = 298.15; - + std::stringstream ss; SetupConfigStream(ss); CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); @@ -67,26 +67,26 @@ TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); - + SECTION("Low temperature range (300K)") { nasa_model.SetTDState_T(300.0); const su2double cp = nasa_model.GetCp(); CHECK(cp == Approx(1040.0).epsilon(0.01)); } - + SECTION("High temperature range (1500K)") { nasa_model.SetTDState_T(1500.0); const su2double cp = nasa_model.GetCp(); CHECK(cp == Approx(1240.0).epsilon(0.01)); } - + SECTION("Temperature range transition") { nasa_model.SetTDState_T(999.0); const su2double cp_below = nasa_model.GetCp(); - + nasa_model.SetTDState_T(1001.0); const su2double cp_above = nasa_model.GetCp(); - + const su2double discontinuity = std::abs(cp_above - cp_below); CHECK(discontinuity < 50.0); } @@ -96,7 +96,7 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]" const su2double R_N2 = 296.8; const su2double P = 101325.0; const su2double T_ref = 298.15; - + std::stringstream ss; SetupConfigStream(ss); CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); @@ -104,29 +104,29 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]" CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); - + SECTION("Round-trip T->H->T at 500K") { const su2double T_target = 500.0; - + nasa_model.SetTDState_T(T_target); const su2double H = nasa_model.GetEnthalpy(); - + nasa_model.SetTDState_h(H); const su2double T_recovered = nasa_model.GetTemperature(); - + CHECK(T_recovered == Approx(T_target).epsilon(1e-6)); } - + SECTION("Multiple temperature points") { const su2double test_temps[] = {250.0, 500.0, 1000.0, 2000.0, 4000.0}; - + for (const auto T_in : test_temps) { nasa_model.SetTDState_T(T_in); const su2double H = nasa_model.GetEnthalpy(); - + nasa_model.SetTDState_h(H); const su2double T_out = nasa_model.GetTemperature(); - + CHECK(T_out == Approx(T_in).epsilon(1e-6)); } } @@ -136,7 +136,7 @@ TEST_CASE("NASA polynomial boundary conditions", "[CIncIdealGasNASA]") { const su2double R_N2 = 296.8; const su2double P = 101325.0; const su2double T_ref = 298.15; - + std::stringstream ss; SetupConfigStream(ss); CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); @@ -144,19 +144,19 @@ TEST_CASE("NASA polynomial boundary conditions", "[CIncIdealGasNASA]") { CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); - + SECTION("Lower temperature bound (200K)") { nasa_model.SetTDState_T(200.0); const su2double cp = nasa_model.GetCp(); - + CHECK(cp > 0.0); CHECK(cp < 2000.0); } - + SECTION("Upper temperature bound (6000K)") { nasa_model.SetTDState_T(6000.0); const su2double cp = nasa_model.GetCp(); - + CHECK(cp > 0.0); CHECK(cp < 3000.0); } @@ -166,7 +166,7 @@ TEST_CASE("NASA polynomial thermodynamic consistency", "[CIncIdealGasNASA]") { const su2double R_N2 = 296.8; const su2double P = 101325.0; const su2double T_ref = 298.15; - + std::stringstream ss; SetupConfigStream(ss); CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); @@ -174,23 +174,23 @@ TEST_CASE("NASA polynomial thermodynamic consistency", "[CIncIdealGasNASA]") { CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); - + SECTION("Density from ideal gas law") { const su2double T = 300.0; nasa_model.SetTDState_T(T); - + const su2double rho = nasa_model.GetDensity(); const su2double rho_expected = P / (R_N2 * T); - + CHECK(rho == Approx(rho_expected).epsilon(1e-10)); } - + SECTION("Cv from Cp (incompressible: gamma = 1)") { nasa_model.SetTDState_T(500.0); - + const su2double cp = nasa_model.GetCp(); const su2double cv = nasa_model.GetCv(); - + CHECK(cv == Approx(cp).epsilon(1e-10)); } } From 56c0f4d6ff7c370a98fef8a936f57b69d76111cb Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Wed, 14 Jan 2026 16:33:07 +0530 Subject: [PATCH 03/11] Refactor NASA polynomials to use CP_POLYCOEFFS and single CP_NASA flag --- Common/include/CConfig.hpp | 46 +----- Common/include/option_structure.hpp | 2 +- Common/src/CConfig.cpp | 12 +- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 76 +++++----- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 32 +++- .../SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp | 141 ++++-------------- 6 files changed, 96 insertions(+), 213 deletions(-) diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp index ff6953989cc6..a2c4b4ec4504 100644 --- a/Common/include/CConfig.hpp +++ b/Common/include/CConfig.hpp @@ -1121,11 +1121,7 @@ class CConfig { array mu_polycoeffs{{0.0}}; /*!< \brief Array for viscosity polynomial coefficients. */ array kt_polycoeffs{{0.0}}; /*!< \brief Array for thermal conductivity polynomial coefficients. */ - array NASA_Coeffs_Low{{0.0}}; /*!< \brief NASA polynomial coefficients (low temp range). */ - array NASA_Coeffs_High{{0.0}}; /*!< \brief NASA polynomial coefficients (high temp range). */ - su2double NASA_Temp_Low; /*!< \brief NASA polynomial low temperature bound. */ - su2double NASA_Temp_Mid; /*!< \brief NASA polynomial mid temperature break. */ - su2double NASA_Temp_High; /*!< \brief NASA polynomial high temperature bound. */ + 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. */ @@ -2599,43 +2595,16 @@ class CConfig { void SetEnergy_Ref(su2double val_energy_ref) { Energy_Ref = val_energy_ref; } /*! - * \brief Get the value of the NASA polynomial coefficient (low temp range). - * \param[in] val_index - Index of the coefficient. - * \return Value of the coefficient. + * \brief Get the flag for using NASA polynomial format for Cp. + * \return True if NASA format is used. */ - su2double GetNASA_CoeffLowND(unsigned short val_index) const { return NASA_Coeffs_Low[val_index]; } + bool GetCp_NASA_Format(void) const { return Cp_NASA_Format; } /*! - * \brief Get the value of the NASA polynomial coefficient (high temp range). - * \param[in] val_index - Index of the coefficient. - * \return Value of the coefficient. + * \brief Set the flag for using NASA polynomial format for Cp. + * \param[in] val - Value of the flag. */ - su2double GetNASA_CoeffHighND(unsigned short val_index) const { return NASA_Coeffs_High[val_index]; } - - /*! - * \brief Get the value of the NASA polynomial low temperature bound. - * \return Value of the temperature. - */ - su2double GetNASA_TempLow() const { return NASA_Temp_Low; } - - /*! - * \brief Get the value of the NASA polynomial mid temperature break. - * \return Value of the temperature. - */ - su2double GetNASA_TempMid() const { return NASA_Temp_Mid; } - - /*! - * \brief Get the value of the NASA polynomial high temperature bound. - * \return Value of the temperature. - */ - su2double GetNASA_TempHigh() const { return NASA_Temp_High; } - - // Setters for NASA polynomials (useful for testing) - void SetNASA_CoeffLowND(unsigned short val_index, su2double val) { NASA_Coeffs_Low[val_index] = val; } - void SetNASA_CoeffHighND(unsigned short val_index, su2double val) { NASA_Coeffs_High[val_index] = val; } - void SetNASA_TempLow(su2double val) { NASA_Temp_Low = val; } - void SetNASA_TempMid(su2double val) { NASA_Temp_Mid = val; } - void SetNASA_TempHigh(su2double val) { NASA_Temp_High = val; } + 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; } @@ -4303,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. diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 690350fae050..7bfff7fb5c9c 100644 --- a/Common/include/option_structure.hpp +++ b/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 diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index fca77dd55bc9..516d1c0eb966 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1339,16 +1339,8 @@ void CConfig::SetConfig_Options() { /* DESCRIPTION: Definition of the temperature polynomial coefficients for specific heat Cp. */ addDoubleArrayOption("KT_POLYCOEFFS", N_POLY_COEFFS, false, kt_polycoeffs.data()); - /* DESCRIPTION: NASA polynomial coefficients (low temp range) */ - addDoubleArrayOption("NASA_POLYCOEFFS_LOW", 7, false, NASA_Coeffs_Low.data()); - /* DESCRIPTION: NASA polynomial coefficients (high temp range) */ - addDoubleArrayOption("NASA_POLYCOEFFS_HIGH", 7, false, NASA_Coeffs_High.data()); - /* DESCRIPTION: NASA polynomial low temperature bound */ - addDoubleOption("NASA_TEMP_LOW", NASA_Temp_Low, 0.0); - /* DESCRIPTION: NASA polynomial mid temperature break */ - addDoubleOption("NASA_TEMP_MID", NASA_Temp_Mid, 0.0); - /* DESCRIPTION: NASA polynomial high temperature bound */ - addDoubleOption("NASA_TEMP_HIGH", NASA_Temp_High, 0.0); + /* 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); diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index 11243cf7b522..7398e3abde6a 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -1,7 +1,7 @@ /*! * \file CIncIdealGasNASA.hpp * \brief Defines the incompressible Ideal Gas model with NASA polynomials for Cp. - * \author Pratyksh Gupta (based on CIncIdealGasPolynomial by T. Economon) + * \author Pratyksh Gupta * \version 8.4.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -40,8 +40,9 @@ * 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 * - * Supports dual temperature ranges (low/high) with separate coefficients. + * Uses a single temperature range provided via CP_POLYCOEFFS (indices 0-6). */ template class CIncIdealGasNASA final : public CFluidModel { @@ -64,22 +65,21 @@ class CIncIdealGasNASA final : public CFluidModel { * \param[in] config - configuration container for the problem. */ void SetCpModel(const CConfig* config, su2double val_Temperature_Ref) override { - T_mid = config->GetNASA_TempMid(); - T_low = config->GetNASA_TempLow(); - T_high = config->GetNASA_TempHigh(); + /*--- 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) { - coeffs_low_[i] = config->GetNASA_CoeffLowND(i); - coeffs_high_[i] = config->GetNASA_CoeffHighND(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); - - if (T_mid <= T_low || T_high <= T_mid) { - SU2_MPI::Error("Invalid NASA polynomial temperature ranges. Must have T_low < T_mid < T_high.", - CURRENT_FUNCTION); - } } /*! @@ -90,15 +90,13 @@ class CIncIdealGasNASA final : public CFluidModel { Temperature = t; Density = Pressure / (Temperature * Gas_Constant); - const su2double* coeffs = (t < T_mid) ? coeffs_low_.data() : coeffs_high_.data(); - - su2double Cp_over_R = coeffs[0] + coeffs[1]*t + coeffs[2]*t*t + - coeffs[3]*t*t*t + coeffs[4]*t*t*t*t; + 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; + 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; @@ -112,7 +110,9 @@ class CIncIdealGasNASA final : public CFluidModel { Enthalpy = val_enthalpy; const su2double toll = 1e-5; - su2double temp_iter = (T_low + T_high) / 2.0; + 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; @@ -120,17 +120,16 @@ class CIncIdealGasNASA final : public CFluidModel { int counter = 0; while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { - const su2double* coeffs = (temp_iter < T_mid) ? coeffs_low_.data() : coeffs_high_.data(); - 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; + 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 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; @@ -139,14 +138,13 @@ class CIncIdealGasNASA final : public CFluidModel { temp_iter += delta_temp_iter; if (temp_iter < Temperature_Min) { - cout << "Warning: Negative temperature found during Newton-Raphson." << endl; - temp_iter = Temperature_Min; - break; + // 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) { - cout << "Warning: Temperature exceeds maximum during Newton-Raphson." << endl; - temp_iter = Temperature_Max; - break; + if (abs(delta_temp_iter) < toll) break; + temp_iter = Temperature_Max; } } @@ -154,7 +152,8 @@ class CIncIdealGasNASA final : public CFluidModel { Cp = Cp_iter; if (counter == counter_limit) { - cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation." << endl; + 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); @@ -166,13 +165,8 @@ class CIncIdealGasNASA final : public CFluidModel { 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 coeffs_low_; /*!< \brief NASA polynomial coefficients for low temperature range. */ - std::array coeffs_high_; /*!< \brief NASA polynomial coefficients for high temperature range. */ - - su2double T_low{200.0}; /*!< \brief Lower temperature bound. */ - su2double T_mid{1000.0}; /*!< \brief Mid-point temperature (split between low/high ranges). */ - su2double T_high{6000.0}; /*!< \brief Upper temperature bound. */ + std::array coeffs_; /*!< \brief NASA polynomial coefficients. */ - su2double Temperature_Min{0.0}; /*!< \brief Minimum temperature allowed in Newton-Raphson iterations. */ - su2double Temperature_Max{1e10}; /*!< \brief Maximum temperature allowed in Newton-Raphson iterations. */ + su2double Temperature_Min{0.0}; /*!< \brief Minimum temperature allowed. */ + su2double Temperature_Max{1e10}; /*!< \brief Maximum temperature allowed. */ }; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index c99f87896c2d..5f8790052a4a 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/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(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(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(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(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); diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp index 83b474787fc5..7bd7583d5cee 100644 --- a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -1,6 +1,6 @@ /*! * \file CIncIdealGasNASA_tests.cpp - * \brief Unit tests for the NASA polynomial fluid model class. + * \brief Unit tests for the NASA polynomial fluid model class (Single Range). * \author Pratyksh Gupta * \version 8.4.0 "Harrier" * @@ -32,27 +32,26 @@ // Helper to config an already constructed CConfig object void ConfigureNASA(CConfig& config) { - su2double nasa_low[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372}; - su2double nasa_high[7] = {2.92664, 1.4879768e-3, -5.68476e-7, 1.0097038e-10, -6.753351e-15, -922.7977, 5.980528}; + // Nitrogen (N2) Low Range Coefficients (200K - 1000K) + su2double nasa_coeffs[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372}; + config.SetCp_NASA_Format(true); for (int i = 0; i < 7; ++i) { - config.SetNASA_CoeffLowND(i, nasa_low[i]); - config.SetNASA_CoeffHighND(i, nasa_high[i]); + config.SetCp_PolyCoeff(i, nasa_coeffs[i]); } - config.SetNASA_TempLow(200.0); - config.SetNASA_TempMid(1000.0); - config.SetNASA_TempHigh(6000.0); + // Ensure we set the temperature limits that SetCpModel reads + // Note: These limits are usually set via TEMPERATURE_LIMITS option, here manually. config.SetTemperatureLimits(0, 200.0); - config.SetTemperatureLimits(1, 6000.0); + config.SetTemperatureLimits(1, 1000.0); } // Helper to populate the stringstream with minimal required options void SetupConfigStream(std::stringstream& ss) { - ss << "SOLVER= EULER" << std::endl; - ss << "MATH_PROBLEM= DIRECT" << std::endl; - ss << "KIND_TURB_MODEL= NONE" << std::endl; - ss << "KIND_TRANS_MODEL= NONE" << std::endl; + ss << "SOLVER= EULER\n"; + ss << "MATH_PROBLEM= DIRECT\n"; + ss << "KIND_TURB_MODEL= NONE\n"; + ss << "KIND_TRANS_MODEL= NONE\n"; } TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { @@ -68,28 +67,12 @@ TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); - SECTION("Low temperature range (300K)") { + SECTION("Valid temperature range (300K)") { nasa_model.SetTDState_T(300.0); const su2double cp = nasa_model.GetCp(); + // Cp of N2 at 300K is approx 1040 J/kgK. CHECK(cp == Approx(1040.0).epsilon(0.01)); } - - SECTION("High temperature range (1500K)") { - nasa_model.SetTDState_T(1500.0); - const su2double cp = nasa_model.GetCp(); - CHECK(cp == Approx(1240.0).epsilon(0.01)); - } - - SECTION("Temperature range transition") { - nasa_model.SetTDState_T(999.0); - const su2double cp_below = nasa_model.GetCp(); - - nasa_model.SetTDState_T(1001.0); - const su2double cp_above = nasa_model.GetCp(); - - const su2double discontinuity = std::abs(cp_above - cp_below); - CHECK(discontinuity < 50.0); - } } TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]") { @@ -105,92 +88,18 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]" CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); - SECTION("Round-trip T->H->T at 500K") { - const su2double T_target = 500.0; - + SECTION("Inversion Consistency") { + su2double T_target = 800.0; nasa_model.SetTDState_T(T_target); - const su2double H = nasa_model.GetEnthalpy(); - - nasa_model.SetTDState_h(H); - const su2double T_recovered = nasa_model.GetTemperature(); - - CHECK(T_recovered == Approx(T_target).epsilon(1e-6)); - } - - SECTION("Multiple temperature points") { - const su2double test_temps[] = {250.0, 500.0, 1000.0, 2000.0, 4000.0}; - - for (const auto T_in : test_temps) { - nasa_model.SetTDState_T(T_in); - const su2double H = nasa_model.GetEnthalpy(); - - nasa_model.SetTDState_h(H); - const su2double T_out = nasa_model.GetTemperature(); - - CHECK(T_out == Approx(T_in).epsilon(1e-6)); - } - } -} - -TEST_CASE("NASA polynomial boundary conditions", "[CIncIdealGasNASA]") { - const su2double R_N2 = 296.8; - const su2double P = 101325.0; - const su2double T_ref = 298.15; - - std::stringstream ss; - SetupConfigStream(ss); - CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); - ConfigureNASA(config); - - CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); - nasa_model.SetCpModel(&config, T_ref); + su2double H_target = nasa_model.GetEnthalpy(); - SECTION("Lower temperature bound (200K)") { - nasa_model.SetTDState_T(200.0); - const su2double cp = nasa_model.GetCp(); - - CHECK(cp > 0.0); - CHECK(cp < 2000.0); - } - - SECTION("Upper temperature bound (6000K)") { - nasa_model.SetTDState_T(6000.0); - const su2double cp = nasa_model.GetCp(); - - CHECK(cp > 0.0); - CHECK(cp < 3000.0); - } -} - -TEST_CASE("NASA polynomial thermodynamic consistency", "[CIncIdealGasNASA]") { - const su2double R_N2 = 296.8; - const su2double P = 101325.0; - const su2double T_ref = 298.15; - - std::stringstream ss; - SetupConfigStream(ss); - CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); - ConfigureNASA(config); - - CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); - nasa_model.SetCpModel(&config, T_ref); - - SECTION("Density from ideal gas law") { - const su2double T = 300.0; - nasa_model.SetTDState_T(T); - - const su2double rho = nasa_model.GetDensity(); - const su2double rho_expected = P / (R_N2 * T); - - CHECK(rho == Approx(rho_expected).epsilon(1e-10)); - } - - SECTION("Cv from Cp (incompressible: gamma = 1)") { - nasa_model.SetTDState_T(500.0); - - const su2double cp = nasa_model.GetCp(); - const su2double cv = nasa_model.GetCv(); - - CHECK(cv == Approx(cp).epsilon(1e-10)); + // Reset state to something else + nasa_model.SetTDState_T(300.0); + + // Invert from H_target + nasa_model.SetTDState_h(H_target); + su2double T_recovered = nasa_model.GetTemperature(); + + CHECK(T_recovered == Approx(T_target).margin(0.001)); } } From 448854f64d6e81f3387d7d5b0a2a19b4d61e7cb8 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Wed, 14 Jan 2026 16:37:09 +0530 Subject: [PATCH 04/11] Apply clang-format fixes --- UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp index 7bd7583d5cee..9f22c8987305 100644 --- a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -39,7 +39,7 @@ void ConfigureNASA(CConfig& config) { for (int i = 0; i < 7; ++i) { config.SetCp_PolyCoeff(i, nasa_coeffs[i]); } - + // Ensure we set the temperature limits that SetCpModel reads // Note: These limits are usually set via TEMPERATURE_LIMITS option, here manually. config.SetTemperatureLimits(0, 200.0); @@ -95,11 +95,11 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]" // Reset state to something else nasa_model.SetTDState_T(300.0); - + // Invert from H_target nasa_model.SetTDState_h(H_target); su2double T_recovered = nasa_model.GetTemperature(); - + CHECK(T_recovered == Approx(T_target).margin(0.001)); } } From 8ce89c99128b229b367a106c3d84ad258d512764 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Wed, 14 Jan 2026 17:54:59 +0530 Subject: [PATCH 05/11] Allow variable number of coeffs for polynomial CP, MU, KT --- Common/src/CConfig.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp index 516d1c0eb966..a113f2dcd09f 100644 --- a/Common/src/CConfig.cpp +++ b/Common/src/CConfig.cpp @@ -1333,11 +1333,11 @@ 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); From 3701aa27cf71b9a34017e7756a861d090040a232 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Thu, 15 Jan 2026 14:19:10 +0530 Subject: [PATCH 06/11] Add NASA polynomial reference citation and improve code clarity with explicit coefficient variables --- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 34 +++++++++++++++------- 1 file changed, 23 insertions(+), 11 deletions(-) diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index 7398e3abde6a..9c1157a46c90 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -37,7 +37,8 @@ * \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: + * Implements NASA 7-coefficient polynomial format (NASA SP-273) for thermodynamic properties: + * 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. * 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 @@ -90,13 +91,18 @@ class CIncIdealGasNASA final : public CFluidModel { 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; + 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]; + + su2double Cp_over_R = a1 + a2*t + a3*t*t + a4*t*t*t + a5*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; + su2double H_over_RT = a1 + a2*t/2.0 + a3*t*t/3.0 + a4*t*t*t/4.0 + a5*t*t*t*t/5.0 + a6/t; Enthalpy = H_over_RT * Gas_Constant * t; Cv = Cp / Gamma; @@ -113,6 +119,13 @@ class CIncIdealGasNASA final : public CFluidModel { 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]; + su2double Cp_iter = 0.0; su2double delta_temp_iter = 1e10; su2double delta_enthalpy_iter; @@ -121,15 +134,14 @@ class CIncIdealGasNASA final : public CFluidModel { 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; + su2double Cp_over_R = a1 + a2*temp_iter + a3*temp_iter*temp_iter + + a4*temp_iter*temp_iter*temp_iter + a5*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 H_over_RT = a1 + a2*temp_iter/2.0 + a3*temp_iter*temp_iter/3.0 + + a4*temp_iter*temp_iter*temp_iter/4.0 + + a5*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + a6/temp_iter; su2double Enthalpy_iter = H_over_RT * Gas_Constant * temp_iter; From eca925a242ba2c4d6ab9a0f27ae82ccc6d6668cf Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Thu, 15 Jan 2026 21:59:28 +0530 Subject: [PATCH 07/11] Implement full NASA-9 polynomial format with NASA-7 backward compatibility --- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 55 +++++++++++++------ SU2_CFD/src/solvers/CIncEulerSolver.cpp | 6 +- .../SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp | 21 +++++-- 3 files changed, 58 insertions(+), 24 deletions(-) diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index 9c1157a46c90..a1b6264042e9 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -37,15 +37,19 @@ * \brief Child class for defining an incompressible ideal gas model with NASA polynomials. * \author Pratyksh Gupta * - * Implements NASA 7-coefficient polynomial format (NASA SP-273) for thermodynamic properties: + * 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. - * 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). + * 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 +template class CIncIdealGasNASA final : public CFluidModel { public: /*! @@ -68,9 +72,13 @@ class CIncIdealGasNASA final : public CFluidModel { 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) ---*/ + 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); @@ -97,12 +105,20 @@ class CIncIdealGasNASA final : public CFluidModel { 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]; - su2double Cp_over_R = a1 + a2*t + a3*t*t + a4*t*t*t + a5*t*t*t*t; + const su2double t_inv = 1.0 / t; + const su2double t_inv2 = t_inv * t_inv; + + // 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*t_inv2 + a2*t_inv + a3 + a4*t + a5*t*t + a6*t*t*t + a7*t*t*t*t; Cp = Cp_over_R * Gas_Constant; - su2double H_over_RT = a1 + a2*t/2.0 + a3*t*t/3.0 + a4*t*t*t/4.0 + a5*t*t*t*t/5.0 + a6/t; + // 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*t_inv2 + a2*std::log(t)*t_inv + a3 + a4*t/2.0 + a5*t*t/3.0 + + a6*t*t*t/4.0 + a7*t*t*t*t/5.0 + a8*t_inv; Enthalpy = H_over_RT * Gas_Constant * t; Cv = Cp / Gamma; @@ -125,6 +141,8 @@ class CIncIdealGasNASA final : public CFluidModel { 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]; su2double Cp_iter = 0.0; su2double delta_temp_iter = 1e10; @@ -134,14 +152,19 @@ class CIncIdealGasNASA final : public CFluidModel { while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { - su2double Cp_over_R = a1 + a2*temp_iter + a3*temp_iter*temp_iter + - a4*temp_iter*temp_iter*temp_iter + a5*temp_iter*temp_iter*temp_iter*temp_iter; + const su2double t_inv = 1.0 / temp_iter; + const su2double t_inv2 = t_inv * t_inv; + + // 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*t_inv2 + a2*t_inv + a3 + a4*temp_iter + a5*temp_iter*temp_iter + + a6*temp_iter*temp_iter*temp_iter + a7*temp_iter*temp_iter*temp_iter*temp_iter; Cp_iter = Cp_over_R * Gas_Constant; - su2double H_over_RT = a1 + a2*temp_iter/2.0 + a3*temp_iter*temp_iter/3.0 + - a4*temp_iter*temp_iter*temp_iter/4.0 + - a5*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + a6/temp_iter; + // 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*t_inv2 + a2*std::log(temp_iter)*t_inv + a3 + a4*temp_iter/2.0 + + a5*temp_iter*temp_iter/3.0 + a6*temp_iter*temp_iter*temp_iter/4.0 + + a7*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + a8*t_inv; su2double Enthalpy_iter = H_over_RT * Gas_Constant * temp_iter; diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 5f8790052a4a..fa3cc5e900c4 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -309,8 +309,8 @@ 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(); 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); + /*--- 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(config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP); @@ -511,7 +511,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i case INC_IDEAL_GAS_POLY: if (config->GetCp_NASA_Format()) { - fluidModel = new CIncIdealGasNASA<7>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); + fluidModel = new CIncIdealGasNASA<9>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); } else { fluidModel = new CIncIdealGasPolynomial(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); } diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp index 9f22c8987305..304ac41511ef 100644 --- a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -32,11 +32,22 @@ // Helper to config an already constructed CConfig object void ConfigureNASA(CConfig& config) { - // Nitrogen (N2) Low Range Coefficients (200K - 1000K) - su2double nasa_coeffs[7] = {3.298677, 1.4082404e-3, -3.963222e-6, 5.641515e-9, -2.444854e-12, -1020.8999, 3.950372}; + // Nitrogen (N2) Low Range Coefficients (200K - 1000K) in NASA-7 format + // NASA-9 indexing: a1=0, a2=0 (no inverse T terms), a3-a7 (polynomial), a8-a9 (integration constants) + su2double nasa_coeffs[9] = { + 0.0, // a1 (T^-2 term, not used in NASA-7) + 0.0, // a2 (T^-1 term, not used in NASA-7) + 3.298677, // a3 (constant term) + 1.4082404e-3, // a4 (T term) + -3.963222e-6, // a5 (T^2 term) + 5.641515e-9, // a6 (T^3 term) + -2.444854e-12, // a7 (T^4 term) + -1020.8999, // a8 (H integration constant) + 3.950372 // a9 (S integration constant) + }; config.SetCp_NASA_Format(true); - for (int i = 0; i < 7; ++i) { + for (int i = 0; i < 9; ++i) { config.SetCp_PolyCoeff(i, nasa_coeffs[i]); } @@ -64,7 +75,7 @@ TEST_CASE("NASA polynomial Cp calculation", "[CIncIdealGasNASA]") { CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); ConfigureNASA(config); - CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); + CIncIdealGasNASA<9> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); SECTION("Valid temperature range (300K)") { @@ -85,7 +96,7 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]" CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); ConfigureNASA(config); - CIncIdealGasNASA<7> nasa_model(R_N2, P, T_ref); + CIncIdealGasNASA<9> nasa_model(R_N2, P, T_ref); nasa_model.SetCpModel(&config, T_ref); SECTION("Inversion Consistency") { From 7fbbf18c352bb6bfc871c02f52e960ff6efcc1f9 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Fri, 16 Jan 2026 10:45:30 +0530 Subject: [PATCH 08/11] Fix NASA polynomial temperature scaling for non-dimensional simulations --- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 31 +++++++++------- SU2_CFD/src/solvers/CIncEulerSolver.cpp | 2 +- .../SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp | 37 +++++++++++++++++++ config_template.cfg | 5 +++ 4 files changed, 61 insertions(+), 14 deletions(-) diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index a1b6264042e9..7f8028d47829 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -55,7 +55,7 @@ class CIncIdealGasNASA final : public CFluidModel { /*! * \brief Constructor of the class. */ - CIncIdealGasNASA(su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref) { + 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). ---*/ @@ -63,6 +63,7 @@ class CIncIdealGasNASA final : public CFluidModel { Pressure = val_operating_pressure; Gamma = 1.0; Std_Ref_Temp_ND = val_Temperature_Ref; + Ref_Temp_Dim = val_Ref_Temp_Dim; } /*! @@ -87,8 +88,8 @@ class CIncIdealGasNASA final : public CFluidModel { } } - Temperature_Min = config->GetTemperatureLimits(0); - Temperature_Max = config->GetTemperatureLimits(1); + Temperature_Min = config->GetTemperatureLimits(0) / Ref_Temp_Dim; + Temperature_Max = config->GetTemperatureLimits(1) / Ref_Temp_Dim; } /*! @@ -108,17 +109,19 @@ class CIncIdealGasNASA final : public CFluidModel { const su2double a7 = coeffs_[6]; const su2double a8 = coeffs_[7]; - const su2double t_inv = 1.0 / t; + // Convert to dimensional temperature for polynomial evaluation (NASA coeffs expect Kelvin) + const su2double T_dim = t * Ref_Temp_Dim; + const su2double t_inv = 1.0 / T_dim; const su2double t_inv2 = t_inv * t_inv; // 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*t_inv2 + a2*t_inv + a3 + a4*t + a5*t*t + a6*t*t*t + a7*t*t*t*t; + su2double Cp_over_R = a1*t_inv2 + a2*t_inv + a3 + a4*T_dim + a5*T_dim*T_dim + a6*T_dim*T_dim*T_dim + a7*T_dim*T_dim*T_dim*T_dim; 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*t_inv2 + a2*std::log(t)*t_inv + a3 + a4*t/2.0 + a5*t*t/3.0 + - a6*t*t*t/4.0 + a7*t*t*t*t/5.0 + a8*t_inv; + su2double H_over_RT = -a1*t_inv2 + a2*std::log(T_dim)*t_inv + a3 + a4*T_dim/2.0 + a5*T_dim*T_dim/3.0 + + a6*T_dim*T_dim*T_dim/4.0 + a7*T_dim*T_dim*T_dim*T_dim/5.0 + a8*t_inv; Enthalpy = H_over_RT * Gas_Constant * t; Cv = Cp / Gamma; @@ -152,19 +155,20 @@ class CIncIdealGasNASA final : public CFluidModel { while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { - const su2double t_inv = 1.0 / temp_iter; + const su2double T_dim = temp_iter * Ref_Temp_Dim; + const su2double t_inv = 1.0 / T_dim; const su2double t_inv2 = t_inv * t_inv; // 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*t_inv2 + a2*t_inv + a3 + a4*temp_iter + a5*temp_iter*temp_iter + - a6*temp_iter*temp_iter*temp_iter + a7*temp_iter*temp_iter*temp_iter*temp_iter; + su2double Cp_over_R = a1*t_inv2 + a2*t_inv + a3 + a4*T_dim + a5*T_dim*T_dim + + a6*T_dim*T_dim*T_dim + a7*T_dim*T_dim*T_dim*T_dim; 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*t_inv2 + a2*std::log(temp_iter)*t_inv + a3 + a4*temp_iter/2.0 + - a5*temp_iter*temp_iter/3.0 + a6*temp_iter*temp_iter*temp_iter/4.0 + - a7*temp_iter*temp_iter*temp_iter*temp_iter/5.0 + a8*t_inv; + su2double H_over_RT = -a1*t_inv2 + a2*std::log(T_dim)*t_inv + a3 + a4*T_dim/2.0 + + a5*T_dim*T_dim/3.0 + a6*T_dim*T_dim*T_dim/4.0 + + a7*T_dim*T_dim*T_dim*T_dim/5.0 + a8*t_inv; su2double Enthalpy_iter = H_over_RT * Gas_Constant * temp_iter; @@ -199,6 +203,7 @@ class CIncIdealGasNASA final : public CFluidModel { 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 coeffs_; /*!< \brief NASA polynomial coefficients. */ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index fa3cc5e900c4..03daa5915bd2 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -511,7 +511,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i case INC_IDEAL_GAS_POLY: if (config->GetCp_NASA_Format()) { - fluidModel = new CIncIdealGasNASA<9>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); + fluidModel = new CIncIdealGasNASA<9>(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref(), config->GetTemperature_Ref()); } else { fluidModel = new CIncIdealGasPolynomial(Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref()); } diff --git a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp index 304ac41511ef..1c4ea8d7c578 100644 --- a/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp +++ b/UnitTests/SU2_CFD/fluid/CIncIdealGasNASA_tests.cpp @@ -114,3 +114,40 @@ TEST_CASE("NASA polynomial enthalpy-temperature inversion", "[CIncIdealGasNASA]" CHECK(T_recovered == Approx(T_target).margin(0.001)); } } + +TEST_CASE("NASA polynomial non-dimensional scaling", "[CIncIdealGasNASA]") { + const su2double R_N2_dim = 296.8; + const su2double P_dim = 101325.0; + const su2double T_ref_dim = 288.15; // Reference Temperature + + // Non-dimensional Gas Constant set to 1.0 for simplicity (Cp output will be Cp/R) + const su2double R_nd = 1.0; + const su2double P_nd = 1.0; // Arbitrary for Cp check + const su2double T_ref_nd_for_class = 298.15 / T_ref_dim; // Std_Ref_Temp_ND + + std::stringstream ss; + SetupConfigStream(ss); + CConfig config(ss, SU2_COMPONENT::SU2_CFD, false); + ConfigureNASA(config); + + // Initialize with dimensional T_ref passed as 4th argument + CIncIdealGasNASA<9> nasa_model(R_nd, P_nd, T_ref_nd_for_class, T_ref_dim); + nasa_model.SetCpModel(&config, T_ref_nd_for_class); + + SECTION("Correct scaling for T = 300K") { + const su2double T_target_dim = 300.0; + const su2double T_target_nd = T_target_dim / T_ref_dim; + + nasa_model.SetTDState_T(T_target_nd); + + const su2double cp_nd = nasa_model.GetCp(); + // Expected Cp/R at 300K for N2 from coeffs: + // a1=a2=0, a3=3.298677, a4=1.4082404e-3, a5=-3.963222e-6, a6=5.641515e-9, a7=-2.444854e-12 + // Cp/R = 3.298677 + 1.4082404e-3*300 - 3.963222e-6*300^2 + 5.641515e-9*300^3 - 2.444854e-12*300^4 + // = 3.298677 + 0.422472 - 0.356690 + 0.152321 - 0.019803 + // = ~3.497 + + // Check if Cp output (which is Cp/R * 1.0) matches approx 3.5 + CHECK(cp_nd == Approx(3.497).epsilon(0.01)); + } +} diff --git a/config_template.cfg b/config_template.cfg index e519f5971aaf..12d01629282b 100644 --- a/config_template.cfg +++ b/config_template.cfg @@ -426,8 +426,13 @@ MOLECULAR_WEIGHT= 28.96, 16.043 % % Temperature polynomial coefficients (up to quartic) for specific heat Cp. % Format -> Cp(T) : b0 + b1*T + b2*T^2 + b3*T^3 + b4*T^4 +% For NASA 9-coefficient format (see CP_NASA option), use 9 coefficients: +% a1 (T^-2), a2 (T^-1), a3 (T^0), a4 (T), a5 (T^2), a6 (T^3), a7 (T^4), a8 (H const), a9 (S const) CP_POLYCOEFFS= (0.0, 0.0, 0.0, 0.0, 0.0) % +% Flag to use NASA polynomial format for Cp in incompressible ideal gas model. +CP_NASA= NO +% % --- Nonequilibrium fluid options % % Gas model - mixture From 0cc637ab6eb4a11812a7646d66b78fdc19ca3dd0 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Sat, 17 Jan 2026 11:34:27 +0530 Subject: [PATCH 09/11] Address code review for NASA polynomials logic: use pow and simplify clamping --- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 32 +++++++--------------- 1 file changed, 10 insertions(+), 22 deletions(-) diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index 7f8028d47829..6889fee7dc68 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -111,17 +111,15 @@ class CIncIdealGasNASA final : public CFluidModel { // Convert to dimensional temperature for polynomial evaluation (NASA coeffs expect Kelvin) const su2double T_dim = t * Ref_Temp_Dim; - const su2double t_inv = 1.0 / T_dim; - const su2double t_inv2 = t_inv * t_inv; - + // 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*t_inv2 + a2*t_inv + a3 + a4*T_dim + a5*T_dim*T_dim + a6*T_dim*T_dim*T_dim + a7*T_dim*T_dim*T_dim*T_dim; + 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*t_inv2 + a2*std::log(T_dim)*t_inv + a3 + a4*T_dim/2.0 + a5*T_dim*T_dim/3.0 + - a6*T_dim*T_dim*T_dim/4.0 + a7*T_dim*T_dim*T_dim*T_dim/5.0 + a8*t_inv; + su2double H_over_RT = -a1 * pow(T_dim, -2.0) + a2 * std::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; Cv = Cp / Gamma; @@ -156,19 +154,17 @@ class CIncIdealGasNASA final : public CFluidModel { while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { const su2double T_dim = temp_iter * Ref_Temp_Dim; - const su2double t_inv = 1.0 / T_dim; - const su2double t_inv2 = t_inv * t_inv; // 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*t_inv2 + a2*t_inv + a3 + a4*T_dim + a5*T_dim*T_dim + - a6*T_dim*T_dim*T_dim + a7*T_dim*T_dim*T_dim*T_dim; + 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*t_inv2 + a2*std::log(T_dim)*t_inv + a3 + a4*T_dim/2.0 + - a5*T_dim*T_dim/3.0 + a6*T_dim*T_dim*T_dim/4.0 + - a7*T_dim*T_dim*T_dim*T_dim/5.0 + a8*t_inv; + su2double H_over_RT = -a1 * pow(T_dim, -2.0) + a2 * std::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; @@ -176,15 +172,7 @@ class CIncIdealGasNASA final : public CFluidModel { 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; - } + temp_iter = std::fmin(std::fmax(Temperature_Min, temp_iter), Temperature_Max); } Temperature = temp_iter; From db284ac4757ac728dae8a2820dbaa706c16510b7 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Sat, 17 Jan 2026 11:36:34 +0530 Subject: [PATCH 10/11] Code cleanup: check linter suggestions --- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index 6889fee7dc68..744fe745ef2c 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -180,7 +180,7 @@ class CIncIdealGasNASA final : public CFluidModel { 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; + std::cout << "Warning: Newton-Raphson exceeds max. iterations in temperature computation (NASA Model).\n"; } Density = Pressure / (Temperature * Gas_Constant); From 9e203c8ad9c80f3d358d815094ccd57e11b59913 Mon Sep 17 00:00:00 2001 From: guptapratykshh Date: Sun, 18 Jan 2026 20:35:11 +0530 Subject: [PATCH 11/11] Complete NASA-9 model: add entropy, improve AD compatibility and robustness --- SU2_CFD/include/fluid/CIncIdealGasNASA.hpp | 31 +++++++++++++++++++--- 1 file changed, 27 insertions(+), 4 deletions(-) diff --git a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp index 744fe745ef2c..93693f910dec 100644 --- a/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp +++ b/SU2_CFD/include/fluid/CIncIdealGasNASA.hpp @@ -97,6 +97,8 @@ class CIncIdealGasNASA final : public CFluidModel { * \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); @@ -108,6 +110,7 @@ class CIncIdealGasNASA final : public CFluidModel { 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; @@ -118,10 +121,17 @@ class CIncIdealGasNASA final : public CFluidModel { 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 * std::log(T_dim) / T_dim + a3 + a4 * T_dim / 2.0 + a5 * pow(T_dim, 2.0) / 3.0 + + 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; } @@ -130,6 +140,11 @@ class CIncIdealGasNASA final : public CFluidModel { * \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; @@ -144,6 +159,7 @@ class CIncIdealGasNASA final : public CFluidModel { 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; @@ -151,7 +167,7 @@ class CIncIdealGasNASA final : public CFluidModel { const int counter_limit = 50; int counter = 0; - while ((abs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { + while ((fabs(delta_temp_iter) > toll) && (counter++ < counter_limit)) { const su2double T_dim = temp_iter * Ref_Temp_Dim; @@ -162,7 +178,7 @@ class CIncIdealGasNASA final : public CFluidModel { 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 * std::log(T_dim) / T_dim + a3 + a4 * T_dim / 2.0 + + 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; @@ -172,11 +188,18 @@ class CIncIdealGasNASA final : public CFluidModel { delta_temp_iter = delta_enthalpy_iter / Cp_iter; temp_iter += delta_temp_iter; - temp_iter = std::fmin(std::fmax(Temperature_Min, temp_iter), Temperature_Max); + 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)