Merge pull request #2062 from su2code/sst_consistency

Make SST-1994m consistent with literature and use 2003m by default

Author: pcarruscag

Common/include/option_structure.hpp

@@ -1020,8 +1020,9 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
 
   const bool default_version = !found_1994 && !found_1994m && !found_2003 && !found_2003m;
 
-  const bool sst_1994 = found_1994 || found_1994m || default_version;
-  const bool sst_2003 = found_2003 || found_2003m;
+  const bool sst_1994 = found_1994 || found_1994m;
+  /*--- Default version since v8. ---*/
+  const bool sst_2003 = found_2003 || found_2003m || default_version;
 
   /*--- When V2003m or V1994m is selected, we automatically select sst_m. ---*/
   const bool sst_m = found_1994m || found_2003m || default_version;
@@ -1037,12 +1038,6 @@ inline SST_ParsedOptions ParseSSTOptions(const SST_OPTIONS *SST_Options, unsigne
     SSTParsedOptions.version = SST_OPTIONS::V2003;
   } else {
     SSTParsedOptions.version = SST_OPTIONS::V1994;
-
-    if (rank==MASTER_NODE) {
-      std::cout <<
-        "WARNING: The current SST-1994m model is inconsistent with literature. We recommend using the SST-2003m model.\n"
-        "In SU2 v8 the 2003m model will become default, and the inconsistency will be fixed." << std::endl;
-    }
   }
 
   // Parse production modifications

SU2_CFD/include/numerics/turbulent/turb_sources.hpp

@@ -49,7 +49,7 @@ struct CSAVariables {
   const su2double cw1 = cb1 / k2 + (1 + cb2) / sigma;
   const su2double cr1 = 0.5;
   const su2double CRot = 1.0;
-  const su2double c2 = 0.7, c3 = 0.9; 
+  const su2double c2 = 0.7, c3 = 0.9;
 
   /*--- List of auxiliary functions ---*/
   su2double ft2, d_ft2, r, d_r, g, d_g, glim, fw, d_fw, Ji, d_Ji, S, Shat, d_Shat, fv1, d_fv1, fv2, d_fv2;
@@ -805,11 +805,6 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       const su2double prod_limit = prod_lim_const * beta_star * Density_i * ScalarVar_i[1] * ScalarVar_i[0];
 
       su2double P = Eddy_Viscosity_i * pow(P_Base, 2);
-
-      if (sstParsedOptions.version == SST_OPTIONS::V1994) {
-        /*--- INTRODUCE THE SST-V1994m BUG WHERE DIVERGENCE TERM WILL BE REMOVED ---*/
-        P -= 2.0 / 3.0 * Density_i * ScalarVar_i[0] * diverg;
-      }
       su2double pk = max(0.0, min(P, prod_limit));
 
       const auto& eddy_visc_var = sstParsedOptions.version == SST_OPTIONS::V1994 ? VorticityMag : StrainMag_i;
@@ -818,8 +813,7 @@ class CSourcePieceWise_TurbSST final : public CNumerics {
       /*--- Production limiter only for V2003, recompute for V1994. ---*/
       su2double pw;
       if (sstParsedOptions.version == SST_OPTIONS::V1994) {
-        /*--- INTRODUCE THE SST-V1994m BUG WHERE DIVERGENCE TERM WILL BE REMOVED ---*/
-        pw = alfa_blended * Density_i * max(pow(P_Base, 2) - 2.0 / 3.0 * zeta * diverg, 0.0);
+        pw = alfa_blended * Density_i * pow(P_Base, 2);
       } else {
         pw = (alfa_blended * Density_i / Eddy_Viscosity_i) * pk;
       }

SU2_CFD/src/solvers/CAdjNSSolver.cpp

@@ -766,12 +766,8 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
             Enthalpy = solver_container[FLOW_SOL]->GetNodes()->GetEnthalpy(iPoint);
 
             /*--- Turbulent kinetic energy ---*/
-            /// TODO: This does not seem to be consistent with the primal treatment.
-            if (config->GetKind_Turb_Model() == TURB_MODEL::SST)
-              val_turb_ke = solver_container[TURB_SOL]->GetNodes()->GetSolution(iPoint,0);
-            else
-              val_turb_ke = 0.0;
-
+            // turb_ke is not considered in the stress tensor, see #797
+            val_turb_ke = 0.0;
             CNumerics::ComputeStressTensor(nDim, tau, PrimVar_Grad+1, Laminar_Viscosity, Density, val_turb_ke);
 
             /*--- Form normal_grad_gridvel = \partial_n (u_omega) ---*/

TestCases/axisymmetric_rans/air_nozzle/air_nozzle_restart.cfg

@@ -93,4 +93,4 @@ SOLUTION_FILENAME= solution_flow.dat
 RESTART_FILENAME= restart_flow.dat
 OUTPUT_WRT_FREQ= 1000
 SCREEN_OUTPUT= (INNER_ITER, RMS_DENSITY, RMS_ENERGY, RMS_TKE, RMS_DISSIPATION, TOTAL_HEATFLUX, \
-                RMS_ADJ_DENSITY, RMS_ADJ_ENERGY, RMS_ADJ_TKE, RMS_ADJ_DISSIPATION, SENS_AOA, SENS_MACH)
+                RMS_ADJ_DENSITY, RMS_ADJ_ENERGY, RMS_ADJ_TKE, RMS_ADJ_DISSIPATION)

TestCases/disc_adj_incomp_rans/naca0012/turb_naca0012_sst.cfg

@@ -13,7 +13,6 @@
 %
 SOLVER= INC_RANS
 KIND_TURB_MODEL= SST
-MATH_PROBLEM= DISCRETE_ADJOINT
 RESTART_SOL= NO
 
 % -------------------- INCOMPRESSIBLE FREE-STREAM DEFINITION ------------------%
@@ -45,11 +44,9 @@ MARKER_MONITORING= ( airfoil )
 % ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
 %
 NUM_METHOD_GRAD= GREEN_GAUSS
-CFL_NUMBER= 25.0
+CFL_NUMBER= 1500.0
 CFL_ADAPT= NO
-CFL_ADAPT_PARAM= ( 1.5, 0.5, 25.0, 10000.0 )
-MAX_DELTA_TIME= 1E6
-RK_ALPHA_COEFF= ( 0.66667, 0.66667, 1.000000 )
+QUASI_NEWTON_NUM_SAMPLES= 30
 ITER= 99999
 OBJECTIVE_FUNCTION= DRAG
 OBJECTIVE_WEIGHT = 1.0
@@ -67,7 +64,7 @@ LINEAR_SOLVER= FGMRES
 LINEAR_SOLVER_PREC= ILU
 LINEAR_SOLVER_ILU_FILL_IN= 0
 LINEAR_SOLVER_ERROR= 1E-10
-LINEAR_SOLVER_ITER= 10
+LINEAR_SOLVER_ITER= 15
 
 % -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
 %
@@ -115,7 +112,7 @@ GRAD_OBJFUNC_FILENAME= of_grad.dat
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
 OUTPUT_WRT_FREQ= 250
-SCREEN_OUTPUT= (INNER_ITER, RMS_ADJ_PRESSURE, RMS_ADJ_TKE, RMS_ADJ_DISSIPATION, SENS_PRESS, SENS_AOA)
+SCREEN_OUTPUT= (INNER_ITER, RMS_ADJ_PRESSURE, RMS_ADJ_TKE, RMS_ADJ_DISSIPATION, SENS_PRESS, SENS_AOA, LINSOL_RESIDUAL)
 
 % --------------------- OPTIMAL SHAPE DESIGN DEFINITION -----------------------%
 %

TestCases/disc_adj_turbomachinery/transonic_stator_2D/transonic_stator.cfg

@@ -14,7 +14,7 @@
 %
 SOLVER= RANS
 KIND_TURB_MODEL= SST
-MATH_PROBLEM= DISCRETE_ADJOINT
+SST_OPTIONS= V1994m
 RESTART_SOL= NO
 NZONES= 1
 
@@ -185,4 +185,4 @@ GRAD_OBJFUNC_FILENAME= of_grad.dat
 SURFACE_FILENAME= surface_flow
 SURFACE_ADJ_FILENAME= surface_adjoint
 OUTPUT_WRT_FREQ= 500
-SCREEN_OUTPUT=(INNER_ITER, RMS_ADJ_DENSITY, RMS_ADJ_TKE)
+SCREEN_OUTPUT=(INNER_ITER, RMS_RES)