Fix treatment of reference temperature (#2688)

bigfooted · web-flow · commit 02b07fb7f162 · 2026-01-09T23:33:03.000+01:00
* update treatment of reference temperature
* update pywrapper testcase
diff --git a/SU2_CFD/include/fluid/CIncIdealGas.hpp b/SU2_CFD/include/fluid/CIncIdealGas.hpp
@@ -39,7 +39,7 @@ class CIncIdealGas final : public CFluidModel {
   /*!
    * \brief Constructor of the class.
    */
-  CIncIdealGas(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure, su2double val_Temperature_Ref) {
+  CIncIdealGas(su2double val_Cp, su2double val_gas_constant, su2double val_operating_pressure, su2double val_Std_Ref_Temp_ND) {
     /*--- 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. ---*/
@@ -48,7 +48,7 @@ class CIncIdealGas final : public CFluidModel {
     Gamma = 1.0;
     Cp = val_Cp;
     Cv = Cp;
-    Std_Ref_Temp_ND = STD_REF_TEMP / val_Temperature_Ref;
+    Std_Ref_Temp_ND = val_Std_Ref_Temp_ND;
   }
 
   /*!
diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp
@@ -297,7 +297,7 @@ 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 CIncIdealGas(config->GetSpecific_Heat_Cp(), config->GetGas_Constant(), Pressure_Thermodynamic, Temperature_FreeStream);
+      auxFluidModel = new CIncIdealGas(config->GetSpecific_Heat_Cp(), config->GetGas_Constant(), Pressure_Thermodynamic, STD_REF_TEMP);
       auxFluidModel->SetTDState_T(Temperature_FreeStream);
       Pressure_Thermodynamic = auxFluidModel->GetPressure();
       config->SetPressure_Thermodynamic(Pressure_Thermodynamic);
@@ -307,7 +307,7 @@ 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, Temperature_FreeStream);
+      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++)
@@ -383,10 +383,10 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
     Pressure_Ref    = 1.0;
   }
   else if (config->GetRef_Inc_NonDim() == INITIAL_VALUES) {
-Density_Ref     = Density_FreeStream;
-Velocity_Ref    = ModVel_FreeStream;
-Temperature_Ref = Temperature_FreeStream;
-Pressure_Ref    = Density_Ref*Velocity_Ref*Velocity_Ref;
+    Density_Ref     = Density_FreeStream;
+    Velocity_Ref    = ModVel_FreeStream;
+    Temperature_Ref = Temperature_FreeStream;
+    Pressure_Ref    = Density_Ref*Velocity_Ref*Velocity_Ref;
   }
   else if (config->GetRef_Inc_NonDim() == REFERENCE_VALUES) {
     Density_Ref     = config->GetInc_Density_Ref();
@@ -484,7 +484,7 @@ Pressure_Ref    = Density_Ref*Velocity_Ref*Velocity_Ref;
         break;
 
       case INC_IDEAL_GAS:
-        fluidModel = new CIncIdealGas(Specific_Heat_CpND, Gas_ConstantND, Pressure_ThermodynamicND, Temperature_FreeStreamND);
+        fluidModel = new CIncIdealGas(Specific_Heat_CpND, Gas_ConstantND, Pressure_ThermodynamicND, STD_REF_TEMP / config->GetTemperature_Ref());
         fluidModel->SetTDState_T(Temperature_FreeStreamND);
         break;
 
diff --git a/TestCases/hybrid_regression.py b/TestCases/hybrid_regression.py
@@ -375,7 +375,7 @@ def main():
     inc_buoyancy.cfg_dir   = "incomp_navierstokes/buoyancy_cavity"
     inc_buoyancy.cfg_file  = "lam_buoyancy_cavity.cfg"
     inc_buoyancy.test_iter = 20
-    inc_buoyancy.test_vals = [-3.860443, -3.122399, 4.074529, -23.413000]
+    inc_buoyancy.test_vals = [-3.860465, -3.122402, 1.322722, -23.413000]
     test_list.append(inc_buoyancy)
 
     # Laminar heated cylinder with polynomial fluid model
@@ -665,8 +665,8 @@ def main():
     slinc_steady.cfg_dir   = "sliding_interface/incompressible_steady"
     slinc_steady.cfg_file  = "config.cfg"
     slinc_steady.test_iter = 19
-    slinc_steady.test_vals         = [19.000000, -1.048972, -1.324971]
-    slinc_steady.test_vals_aarch64 = [19.000000, -1.048972, -1.324971]
+    slinc_steady.test_vals         = [19.000000, -1.154874, -1.378120]
+    slinc_steady.test_vals_aarch64 = [19.000000, -1.154874, -1.378120]
     slinc_steady.multizone = True
     test_list.append(slinc_steady)
 
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
@@ -619,7 +619,7 @@ def main():
     inc_buoyancy.cfg_dir   = "incomp_navierstokes/buoyancy_cavity"
     inc_buoyancy.cfg_file  = "lam_buoyancy_cavity.cfg"
     inc_buoyancy.test_iter = 20
-    inc_buoyancy.test_vals = [-3.860462, -3.122401, 4.074510, -23.414000]
+    inc_buoyancy.test_vals = [-3.860462, -3.122401, 1.322722, -23.414000]
     test_list.append(inc_buoyancy)
 
     # Laminar heated cylinder with polynomial fluid model
@@ -1225,7 +1225,7 @@ def main():
     slinc_steady.cfg_dir   = "sliding_interface/incompressible_steady"
     slinc_steady.cfg_file  = "config.cfg"
     slinc_steady.test_iter = 19
-    slinc_steady.test_vals = [19.000000, -1.013427, -1.222965]
+    slinc_steady.test_vals = [19.000000, -1.131557, -1.370471]
     slinc_steady.timeout   = 100
     slinc_steady.tol       = 0.00002
     slinc_steady.multizone = True
@@ -1352,7 +1352,7 @@ def main():
     p1rad.cfg_dir   = "radiation/p1model"
     p1rad.cfg_file  = "configp1.cfg"
     p1rad.test_iter = 50
-    p1rad.test_vals = [-7.700974, -7.734765, -2.061661, 0.382811, -57.861000]
+    p1rad.test_vals = [-8.197112, -7.856394, -2.380743, 0.381122, -57.186000]
     test_list.append(p1rad)
 
 
@@ -1520,8 +1520,8 @@ def main():
     pywrapper_buoyancy.cfg_dir = "py_wrapper/custom_source_buoyancy"
     pywrapper_buoyancy.cfg_file = "lam_buoyancy_cavity.cfg"
     pywrapper_buoyancy.test_iter = 0
-    pywrapper_buoyancy.test_vals = [-13.314107, -13.764877, -13.136439, -5.379183]
-    pywrapper_buoyancy.test_vals_aarch64 = [-17.746018, -17.460693, -17.430708, -12.260624]
+    pywrapper_buoyancy.test_vals = [-13.227985, -13.678478, -13.050758, -7.777151]
+    pywrapper_buoyancy.test_vals_aarch64 = [-13.227985, -13.678478, -13.050758, -7.777151]
     pywrapper_buoyancy.command = TestCase.Command("mpirun -np 2", "python", "run.py")
     test_list.append(pywrapper_buoyancy)
 
@@ -1530,7 +1530,7 @@ def main():
     pywrapper_zimont.cfg_dir = "py_wrapper/turbulent_premixed_psi"
     pywrapper_zimont.cfg_file = "psi.cfg"
     pywrapper_zimont.test_iter = 0
-    pywrapper_zimont.test_vals = [-9.317943, -0.377235, -1.091632, -0.502950, 2.968250, -1.242417]
+    pywrapper_zimont.test_vals = [-3.415653, -2.221781, -3.107330, -2.569248, 0.531838, -3.796447]
     pywrapper_zimont.command = TestCase.Command("mpirun -np 2", "python", "run.py")
     test_list.append(pywrapper_zimont)
 
diff --git a/TestCases/parallel_regression_AD.py b/TestCases/parallel_regression_AD.py
@@ -290,7 +290,7 @@ def main():
     discadj_cht.cfg_dir   = "coupled_cht/disc_adj_incomp_2d"
     discadj_cht.cfg_file  = "cht_2d_3cylinders.cfg"
     discadj_cht.test_iter = 10
-    discadj_cht.test_vals = [-6.111499, -2.668835, -2.669557, -2.669594]
+    discadj_cht.test_vals = [-6.308199, -3.086021, -3.086022, -3.086022]
     test_list.append(discadj_cht)
 
     # 2D DA cht streamwise periodic case, 2 zones, avg temp objective
diff --git a/TestCases/py_wrapper/turbulent_premixed_psi/run.py b/TestCases/py_wrapper/turbulent_premixed_psi/run.py
@@ -50,10 +50,9 @@
 k_u = 0.0523
 # unburnt heat capacity of methane-air (phi=0.5, P=5)
 cp_u = 1350.0
-
-# P = rho*R*T
-# 5 = 2.55 * R * 673
-# R = 0.0029
+# laminar burning velocity of methane-air at phi=0.5, P=5
+Slu = 0.232
+Le = 1.0
 
 
 # ################################################################## #
@@ -94,18 +93,6 @@ def update_temperature(SU2Driver, iPoint):
     C = SU2Driver.Solution(iSPECIESSOLVER)(iPoint,0)
     T = Tu*(1-C) + Tf*C
 
-    #iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW']
-    #solindex = getsolvar(SU2Driver)
-    #iTEMP = solindex.get("TEMPERATURE")
-    #T = SU2Driver.Solution(iFLOWSOLVER)(iPoint,iTEMP)
-    #SU2Driver.Solution(iFLOWSOLVER).Set(iPoint,iTEMP,cp_u*T)
-
-    prim_indices = SU2Driver.GetPrimitiveIndices()
-    iTemp = prim_indices['TEMPERATURE']
-    SU2Driver.Primitives().Set(iPoint,iTemp, T)
-    #ih = prim_indices['ENTHALPY']
-    #SU2Driver.Primitives().Set(iPoint,ih, cp_u*(T-Tref))
-        
     iFLOWSOLVER = SU2Driver.GetSolverIndices()['INC.FLOW']
     iENTH = 3
     #h = 
@@ -115,7 +102,7 @@ def update_temperature(SU2Driver, iPoint):
 # ################################################################## #
 # Source term according to Zimont
 # ################################################################## #
-def zimont(SU2Driver, iPoint):
+def zimont(SU2Driver, iPoint, nDim):
 
     iSSTSOLVER = SU2Driver.GetSolverIndices()['SST']
     tke, dissipation = SU2Driver.Solution(iSSTSOLVER)(iPoint)
@@ -127,19 +114,20 @@ def zimont(SU2Driver, iPoint):
     iDENSITY = primindex.get("DENSITY")
     iMU = primindex.get("LAMINAR_VISCOSITY")
 
-    # laminar burning velocity of methane-air at phi=0.5, P=5
-    Slu = 0.232
-
     rho = SU2Driver.Primitives()(iPoint,iDENSITY)
     mu = SU2Driver.Primitives()(iPoint,iMU)
     nu=mu/rho
     # Turbulent Flamespeed Closure with Dinkelacker correction
     up = np.sqrt((2.0/3.0) * tke )
     lt = (0.09**0.75) * (tke**1.5) / dissipation
     Re = up*lt/nu
-    Le = 1.0
     Ut = Slu * (1.0 + (0.46/Le) * np.power(Re,0.25) * np.power(up/Slu,0.3) * np.power(Pu,0.2) )
-    norm_gradc = np.sqrt(gradc[0]*gradc[0] + gradc[1]*gradc[1])
+
+    norm_gradc = 0.0
+    for idim in range(nDim):
+      norm_gradc += gradc[idim]*gradc[idim]
+    norm_gradc = np.sqrt(norm_gradc)
+
     Sc = rho_u * Ut * norm_gradc
 
     return Sc
@@ -233,7 +221,7 @@ def main():
     for i_node in range(driver.GetNumberNodes() - driver.GetNumberHaloNodes()):
       # add source term:
       # default TFC of Zimont: rho*Sc = rho_u * U_t * grad(c)
-      S = zimont(driver,i_node)
+      S = zimont(driver,i_node, nDim)
       Source.Set(i_node,0,S)
 
     # for the update of temperature, we need to update also the halo nodes
diff --git a/TestCases/radiation/p1adjoint/of_grad_cd.csv.ref b/TestCases/radiation/p1adjoint/of_grad_cd.csv.ref
@@ -1,6 +1,6 @@
 "VARIABLE"      , "GRADIENT"      , "FINDIFF_STEP"
-0               , -0.000958205    , 0.001
-1               , -0.00432733     , 0.001
+0               , -0.000958193    , 0.001
+1               , -0.00432732     , 0.001
 2               , -0.0111443      , 0.001
 3               , -0.0214515      , 0.001
 4               , -0.0343997      , 0.001
@@ -20,16 +20,16 @@
 18              , -0.042129       , 0.001
 19              , -0.0303137      , 0.001
 20              , -0.0203813      , 0.001
-21              , -0.0126226      , 0.001
-22              , -0.00700441     , 0.001
-23              , -0.00322726     , 0.001
-24              , -0.000951904    , 0.001
-25              , 0.00128803      , 0.001
-26              , 0.0051774       , 0.001
+21              , -0.0126225      , 0.001
+22              , -0.0070044      , 0.001
+23              , -0.00322725     , 0.001
+24              , -0.000951896    , 0.001
+25              , 0.00128802      , 0.001
+26              , 0.00517739      , 0.001
 27              , 0.0125547       , 0.001
 28              , 0.0234674       , 0.001
 29              , 0.0370465       , 0.001
-30              , 0.0521956       , 0.001
+30              , 0.0521957       , 0.001
 31              , 0.0678379       , 0.001
 32              , 0.0827257       , 0.001
 33              , 0.0955787       , 0.001
@@ -47,5 +47,5 @@
 45              , 0.0192316       , 0.001
 46              , 0.0124655       , 0.001
 47              , 0.00752232      , 0.001
-48              , 0.0039072       , 0.001
+48              , 0.00390719      , 0.001
 49              , 0.00131195      , 0.001
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
@@ -415,7 +415,7 @@ def main():
     inc_buoyancy.cfg_dir   = "incomp_navierstokes/buoyancy_cavity"
     inc_buoyancy.cfg_file  = "lam_buoyancy_cavity.cfg"
     inc_buoyancy.test_iter = 20
-    inc_buoyancy.test_vals = [-3.860463, -3.122402, 4.074509, -23.414000]
+    inc_buoyancy.test_vals = [-3.860466, -3.122403, 1.322721, -23.413000]
     test_list.append(inc_buoyancy)
 
     # Laminar heated cylinder with polynomial fluid model
@@ -994,7 +994,7 @@ def main():
     slinc_steady.cfg_dir   = "sliding_interface/incompressible_steady"
     slinc_steady.cfg_file  = "config.cfg"
     slinc_steady.test_iter = 19
-    slinc_steady.test_vals = [19.000000, -1.050011, -1.303591]
+    slinc_steady.test_vals = [19.000000, -1.148249, -1.398402]
     slinc_steady.timeout   = 100
     slinc_steady.multizone = True
     test_list.append(slinc_steady)
@@ -1106,7 +1106,7 @@ def main():
     p1rad.cfg_dir   = "radiation/p1model"
     p1rad.cfg_file  = "configp1.cfg"
     p1rad.test_iter = 50
-    p1rad.test_vals = [-7.939639, -7.898740, -2.238329, 0.388591, -56.886000]
+    p1rad.test_vals = [-8.284674, -8.008659, -2.424204, 0.389030, -56.560000]
     test_list.append(p1rad)
 
     # ###############################
diff --git a/TestCases/serial_regression_AD.py b/TestCases/serial_regression_AD.py
@@ -219,7 +219,7 @@ def main():
     discadj_cht.cfg_dir   = "coupled_cht/disc_adj_incomp_2d"
     discadj_cht.cfg_file  = "cht_2d_3cylinders.cfg"
     discadj_cht.test_iter = 10
-    discadj_cht.test_vals = [-6.111521, -2.668836, -2.669558, -2.669595]
+    discadj_cht.test_vals = [-6.308199, -3.086023, -3.086024, -3.086024]
     test_list.append(discadj_cht)
 
     ######################################
diff --git a/TestCases/tutorials.py b/TestCases/tutorials.py
@@ -49,7 +49,7 @@ def main():
     cht_incompressible_unsteady.cfg_dir   = "../Tutorials/multiphysics/unsteady_cht/"
     cht_incompressible_unsteady.cfg_file  = "cht_2d_3cylinders.cfg"
     cht_incompressible_unsteady.test_iter = 2
-    cht_incompressible_unsteady.test_vals = [-2.536012, -0.080399, -0.080399, -0.080399, -12.421963, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 238.240000]
+    cht_incompressible_unsteady.test_vals = [-3.075372, -0.080399, -0.080399, -0.080399, -11.163219, 0.000000, 0.000000, 0.000000, 0.000000, 0.000000, 238.240000]
     cht_incompressible_unsteady.multizone = True
     cht_incompressible_unsteady.unsteady  = True
     test_list.append(cht_incompressible_unsteady)
