From b174af6ab12513e20da199764000a1099aa4a8bc Mon Sep 17 00:00:00 2001
From: vdweide <e.t.a.vanderweide@utwente.nl>
Date: Tue, 20 Aug 2019 16:47:36 +0200
Subject: [PATCH 1/6] Added the SST turbulence model with sustaining terms

---
 Common/include/option_structure.hpp       | 20 ++++----
 Common/src/config_structure.cpp           | 17 ++++---
 Common/src/geometry_structure.cpp         |  6 ++-
 SU2_CFD/include/numerics_structure.hpp    |  8 ++-
 SU2_CFD/src/drivers/CDriver.cpp           | 62 ++++++++++++++---------
 SU2_CFD/src/numerics_direct_turbulent.cpp | 31 ++++++++++--
 SU2_CFD/src/output_structure.cpp          | 40 +++++++++------
 SU2_CFD/src/solver_adjoint_mean.cpp       |  4 +-
 SU2_CFD/src/solver_adjoint_turbulent.cpp  | 10 ++--
 SU2_CFD/src/solver_direct_mean.cpp        | 59 +++++++++++----------
 SU2_CFD/src/solver_direct_mean_fem.cpp    |  2 +-
 SU2_CFD/src/solver_direct_mean_inc.cpp    | 23 ++++-----
 SU2_CFD/src/solver_direct_turbulent.cpp   | 15 +++---
 config_template.cfg                       |  2 +-
 14 files changed, 179 insertions(+), 120 deletions(-)

diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp
index b0a54b8fb3f..94ef18344dc 100644
--- a/Common/include/option_structure.hpp
+++ b/Common/include/option_structure.hpp
@@ -819,28 +819,30 @@ static const map<string, ENUM_LIMITER> Limiter_Map = CCreateMap<string, ENUM_LIM
  */
 enum ENUM_TURB_MODEL {
   NO_TURB_MODEL = 0, /*!< \brief No turbulence model. */
-  SA      = 1, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
-  SA_NEG  = 2, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
-  SST     = 3, /*!< \brief Kind of Turbulence model (Menter SST). */
-  SA_E    = 4, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards). */
-  SA_COMP = 5, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Compressibility Correction). */
-  SA_E_COMP = 6, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards with Compressibility Correction). */
+  SA        = 1, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
+  SA_NEG    = 2, /*!< \brief Kind of Turbulent model (Spalart-Allmaras). */
+  SA_E      = 3, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards). */
+  SA_COMP   = 4, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Compressibility Correction). */
+  SA_E_COMP = 5, /*!< \brief Kind of Turbulent model (Spalart-Allmaras Edwards with Compressibility Correction). */
+  SST       = 6, /*!< \brief Kind of Turbulence model (Menter SST). */
+  SST_SUST  = 7  /*!< \brief Kind of Turbulence model (Menter SST with sustaining terms for free-stream preservation). */
 };
 static const map<string, ENUM_TURB_MODEL> Turb_Model_Map = CCreateMap<string, ENUM_TURB_MODEL>
 ("NONE", NO_TURB_MODEL)
 ("SA", SA)
 ("SA_NEG", SA_NEG)
-("SST", SST)
 ("SA_E", SA_E)
 ("SA_COMP", SA_COMP)
-("SA_E_COMP", SA_E_COMP);
+("SA_E_COMP", SA_E_COMP)
+("SST", SST)
+("SST_SUST", SST_SUST);
 
 /*!
  * \brief types of transition models
  */
 enum ENUM_TRANS_MODEL {
   NO_TRANS_MODEL = 0,            /*!< \brief No transition model. */
-  LM = 1,	/*!< \brief Kind of transition model (LM for Spalart-Allmaras). */
+  LM = 1,	/*!< \brief Kind of transition model (Langtry-Menter (LM) for SST and Spalart-Allmaras). */
   BC = 2	/*!< \brief Kind of transition model (BAS-CAKMAKCIOGLU (BC) for Spalart-Allmaras). */
 };
 static const map<string, ENUM_TRANS_MODEL> Trans_Model_Map = CCreateMap<string, ENUM_TRANS_MODEL>
diff --git a/Common/src/config_structure.cpp b/Common/src/config_structure.cpp
index a30427de3f3..31abef80b5c 100644
--- a/Common/src/config_structure.cpp
+++ b/Common/src/config_structure.cpp
@@ -4121,7 +4121,7 @@ void CConfig::SetPostprocessing(unsigned short val_software, unsigned short val_
 
   /* --- Throw error if UQ used for any turbulence model other that SST --- */
 
-  if (Kind_Solver == RANS && Kind_Turb_Model != SST && using_uq){
+  if (Kind_Solver == RANS && Kind_Turb_Model != SST && Kind_Turb_Model != SST_SUST && using_uq){
     SU2_MPI::Error("UQ capabilities only implemented for SST turbulence model", CURRENT_FUNCTION);
   }
 
@@ -5010,21 +5010,22 @@ void CConfig::SetOutput(unsigned short val_software, unsigned short val_izone) {
         if (Kind_Regime == INCOMPRESSIBLE) cout << "Incompressible RANS equations." << endl;
         cout << "Turbulence model: ";
         switch (Kind_Turb_Model) {
-          case SA:     cout << "Spalart Allmaras" << endl; break;
-          case SA_NEG: cout << "Negative Spalart Allmaras" << endl; break;
-          case SST:    cout << "Menter's SST"     << endl; break;
-          case SA_E:   cout << "Edwards Spalart Allmaras" << endl; break;
+          case SA:        cout << "Spalart Allmaras" << endl; break;
+          case SA_NEG:    cout << "Negative Spalart Allmaras" << endl; break;
+          case SA_E:      cout << "Edwards Spalart Allmaras" << endl; break;
           case SA_COMP:   cout << "Compressibility Correction Spalart Allmaras" << endl; break;
-          case SA_E_COMP:   cout << "Compressibility Correction Edwards Spalart Allmaras" << endl; break;
+          case SA_E_COMP: cout << "Compressibility Correction Edwards Spalart Allmaras" << endl; break;
+          case SST:       cout << "Menter's SST"     << endl; break;
+          case SST_SUST:  cout << "Menter's SST with sustaining terms" << endl; break;
         }
         if (QCR) cout << "Using Quadratic Constitutive Relation, 2000 version (QCR2000)" << endl;
         cout << "Hybrid RANS/LES: ";
         switch (Kind_HybridRANSLES){
           case NO_HYBRIDRANSLES: cout <<  "No Hybrid RANS/LES" << endl; break;
-          case SA_DES:  cout << "Detached Eddy Simulation (DES97) " << endl; break;
+          case SA_DES:   cout << "Detached Eddy Simulation (DES97) " << endl; break;
           case SA_DDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Standard SGS" << endl; break;
           case SA_ZDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Vorticity-based SGS" << endl; break;
-          case SA_EDDES:  cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
+          case SA_EDDES: cout << "Delayed Detached Eddy Simulation (DDES) with Shear-layer Adapted SGS" << endl; break;
         }
         if (using_uq){
           cout << "Perturbing Reynold's Stress Matrix towards "<< eig_val_comp << " component turbulence"<< endl;
diff --git a/Common/src/geometry_structure.cpp b/Common/src/geometry_structure.cpp
index 60b251f9393..8acb723b17d 100644
--- a/Common/src/geometry_structure.cpp
+++ b/Common/src/geometry_structure.cpp
@@ -16024,8 +16024,10 @@ void CPhysicalGeometry::SetSensitivity(CConfig *config) {
   string filename = config->GetSolution_AdjFileName();
   bool compressible = (config->GetKind_Regime() == COMPRESSIBLE);
   bool incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
-  bool sst = config->GetKind_Turb_Model() == SST;
-  bool sa = (config->GetKind_Turb_Model() == SA) || (config->GetKind_Turb_Model() == SA_NEG);
+  bool sst = (config->GetKind_Turb_Model() == SST)  || (config->GetKind_Turb_Model() == SST_SUST);
+  bool sa  = (config->GetKind_Turb_Model() == SA)   || (config->GetKind_Turb_Model() == SA_NEG)  ||
+             (config->GetKind_Turb_Model() == SA_E) || (config->GetKind_Turb_Model() == SA_COMP) ||
+             (config->GetKind_Turb_Model() == SA_E_COMP);
   bool grid_movement = config->GetGrid_Movement();
   bool frozen_visc = config->GetFrozen_Visc_Disc();
   unsigned short Kind_Solver = config->GetKind_Solver();
diff --git a/SU2_CFD/include/numerics_structure.hpp b/SU2_CFD/include/numerics_structure.hpp
index 591192ac9dd..1fe60acddaf 100644
--- a/SU2_CFD/include/numerics_structure.hpp
+++ b/SU2_CFD/include/numerics_structure.hpp
@@ -5074,8 +5074,11 @@ class CSourcePieceWise_TurbSST : public CNumerics {
   a1;
   
   su2double CDkw_i, CDkw_j;
-  
+
+  su2double kAmb, omegaAmb;
+ 
   bool incompressible;
+  bool sustaining_terms;
   
 public:
   
@@ -5085,7 +5088,8 @@ class CSourcePieceWise_TurbSST : public CNumerics {
    * \param[in] val_nVar - Number of variables of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, su2double* constants, CConfig *config);
+  CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, su2double* constants,
+                           su2double val_kine_Inf, su2double val_omega_Inf, CConfig *config);
   
   /*!
    * \brief Destructor of the class.
diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index 1bec79bf50b..af13ab35b72 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -1091,12 +1091,13 @@ void CDriver::Solver_Preprocessing(CConfig* config, CGeometry** geometry, CSolve
   
   if (turbulent || fem_turbulent)
     switch (config->GetKind_Turb_Model()) {
-      case SA:     spalart_allmaras = true;     break;
-      case SA_NEG: neg_spalart_allmaras = true; break;
-      case SST:    menter_sst = true;           break;
-      case SA_E:   e_spalart_allmaras = true;   break;
-      case SA_COMP: comp_spalart_allmaras = true; break;
+      case SA:        spalart_allmaras = true;        break;
+      case SA_NEG:    neg_spalart_allmaras = true;    break;
+      case SA_E:      e_spalart_allmaras = true;      break;
+      case SA_COMP:   comp_spalart_allmaras = true;   break;
       case SA_E_COMP: e_comp_spalart_allmaras = true; break;
+      case SST:       menter_sst = true;              break;
+      case SST_SUST:  menter_sst = true;              break;
       default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
     }
   
@@ -1566,12 +1567,14 @@ void CDriver::Solver_Postprocessing(CSolver ****solver, CGeometry **geometry,
   
   if (turbulent)
     switch (config->GetKind_Turb_Model()) {
-    case SA:     spalart_allmaras = true;     break;
-    case SA_NEG: neg_spalart_allmaras = true; break;
-    case SST:    menter_sst = true;           break;
-    case SA_E: e_spalart_allmaras = true; break;
-    case SA_COMP: comp_spalart_allmaras = true; break;
+    case SA:        spalart_allmaras = true;        break;
+    case SA_NEG:    neg_spalart_allmaras = true;    break;
+    case SA_E:      e_spalart_allmaras = true;      break;
+    case SA_COMP:   comp_spalart_allmaras = true;   break;
     case SA_E_COMP: e_comp_spalart_allmaras = true; break;
+    case SST:       menter_sst = true;              break;
+    case SST_SUST:  menter_sst = true;              break;
+    default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
     }
   
   /*--- Definition of the Class for the solution: solver_container[DOMAIN][MESH_LEVEL][EQUATION]. Note that euler, ns
@@ -1804,6 +1807,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CSolver ***solver, CNumeri
   numerics = new CNumerics***[config->GetnMGLevels()+1];
   
   su2double *constants = NULL;
+  su2double kine_Inf = 0.0, omega_Inf = 0.0;
   
   bool
   euler, adj_euler,
@@ -1854,14 +1858,24 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CSolver ***solver, CNumeri
 
   if (turbulent || fem_turbulent)
     switch (config->GetKind_Turb_Model()) {
-      case SA:     spalart_allmaras = true;     break;
-      case SA_NEG: neg_spalart_allmaras = true; break;
-      case SA_E:   e_spalart_allmaras = true; break;
-      case SA_COMP:   comp_spalart_allmaras = true; break;
-      case SA_E_COMP:   e_comp_spalart_allmaras = true; break;
-      case SST:    menter_sst = true; constants = solver[MESH_0][TURB_SOL]->GetConstants(); break;
+      case SA:        spalart_allmaras = true;        break;
+      case SA_NEG:    neg_spalart_allmaras = true;    break;
+      case SA_E:      e_spalart_allmaras = true;      break;
+      case SA_COMP:   comp_spalart_allmaras = true;   break;
+      case SA_E_COMP: e_comp_spalart_allmaras = true; break;
+      case SST:       menter_sst = true;              break;
+      case SST_SUST:  menter_sst = true;              break;
       default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
     }
+
+  /*--- If the Menter SST model is used, store the constants of the model and determine the
+        free stream values of the turbulent kinetic energy and dissipation rate. ---*/
+
+  if (menter_sst) {
+    constants = solver[MESH_0][TURB_SOL]->GetConstants();
+    kine_Inf  = solver[MESH_0][TURB_SOL]->GetTke_Inf();
+    omega_Inf = solver[MESH_0][TURB_SOL]->GetOmega_Inf();
+  }
   
   /*--- Number of variables for the template ---*/
   
@@ -2258,7 +2272,7 @@ void CDriver::Numerics_Preprocessing(CConfig *config, CSolver ***solver, CNumeri
       else if (comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSA_COMP(nDim, nVar_Turb, config);
       else if (e_comp_spalart_allmaras) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSA_E_COMP(nDim, nVar_Turb, config);
       else if (neg_spalart_allmaras) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSA_Neg(nDim, nVar_Turb, config);
-      else if (menter_sst) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, config);
+      else if (menter_sst) numerics[iMGlevel][TURB_SOL][SOURCE_FIRST_TERM] = new CSourcePieceWise_TurbSST(nDim, nVar_Turb, constants, kine_Inf, omega_Inf, config);
       numerics[iMGlevel][TURB_SOL][SOURCE_SECOND_TERM] = new CSourceNothing(nDim, nVar_Turb, config);
     }
     
@@ -2666,13 +2680,14 @@ void CDriver::Numerics_Postprocessing(CNumerics *****numerics,
 
   if (turbulent || fem_turbulent)
     switch (config->GetKind_Turb_Model()) {
-      case SA:     spalart_allmaras = true;     break;
-      case SA_NEG: neg_spalart_allmaras = true; break;
-      case SST:    menter_sst = true;  break;
-      case SA_COMP: comp_spalart_allmaras = true; break;
-      case SA_E: e_spalart_allmaras = true; break;
+      case SA:        spalart_allmaras = true;        break;
+      case SA_NEG:    neg_spalart_allmaras = true;    break;
+      case SA_COMP:   comp_spalart_allmaras = true;   break;
+      case SA_E:      e_spalart_allmaras = true;      break;
       case SA_E_COMP: e_comp_spalart_allmaras = true; break;
-
+      case SST:       menter_sst = true;              break;
+      case SST_SUST:  menter_sst = true;              break;
+      default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
     }
   
   /*--- Solver definition for the template problem ---*/
@@ -2715,6 +2730,7 @@ void CDriver::Numerics_Postprocessing(CNumerics *****numerics,
           
           /*--- Definition of the boundary condition method ---*/
           for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++)
+      default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
             delete numerics[val_iInst][iMGlevel][FLOW_SOL][CONV_BOUND_TERM];
           
         }
diff --git a/SU2_CFD/src/numerics_direct_turbulent.cpp b/SU2_CFD/src/numerics_direct_turbulent.cpp
index 533704fb331..0d8d069855b 100644
--- a/SU2_CFD/src/numerics_direct_turbulent.cpp
+++ b/SU2_CFD/src/numerics_direct_turbulent.cpp
@@ -1166,9 +1166,11 @@ void CAvgGrad_TurbSST::FinishResidualCalc(su2double *val_residual, su2double **J
 }
 
 CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, unsigned short val_nVar, su2double *constants,
-                                                   CConfig *config) : CNumerics(val_nDim, val_nVar, config) {
+                                                   su2double val_kine_Inf, su2double val_omega_Inf, CConfig *config)
+  : CNumerics(val_nDim, val_nVar, config) {
   
-  incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
+  incompressible   = (config->GetKind_Regime()     == INCOMPRESSIBLE);
+  sustaining_terms = (config->GetKind_Turb_Model() == SST_SUST);
   
   /*--- Closure constants ---*/
   beta_star     = constants[6];
@@ -1179,6 +1181,10 @@ CSourcePieceWise_TurbSST::CSourcePieceWise_TurbSST(unsigned short val_nDim, unsi
   alfa_1        = constants[8];
   alfa_2        = constants[9];
   a1            = constants[7];
+
+  /*--- Set the ambient values of k and omega to the free stream values. ---*/
+  kAmb     = val_kine_Inf;
+  omegaAmb = val_omega_Inf;
 }
 
 CSourcePieceWise_TurbSST::~CSourcePieceWise_TurbSST(void) { }
@@ -1256,10 +1262,27 @@ void CSourcePieceWise_TurbSST::ComputeResidual(su2double *val_residual, su2doubl
    else {
      pw = StrainMag_i*StrainMag_i - 2.0/3.0*zeta*diverg;
    }
-   pw = max(pw,0.0);
+   pw = alfa_blended*Density_i*max(pw,0.0);
+
+   /*--- Sustaining terms, if desired. Note that if the production terms are
+         larger equal than the sustaining terms, the original formulation is
+         obtained again. This is in contrast to the version in literature
+         where the sustaining terms are simply added. This latter approach could
+         lead to problems for very big values of the free-stream turbulence
+         intensity. ---*/
+
+   if ( sustaining_terms ) {
+     const su2double sust_k = beta_star*Density_i*kAmb*omegaAmb;
+     const su2double sust_w = beta_blended*Density_i*omegaAmb*omegaAmb;
+
+     pk = max(pk, sust_k);
+     pw = max(pw, sust_w);
+   }
+
+   /*--- Add the production terms to the residuals. ---*/
 
    val_residual[0] += pk*Volume;
-   val_residual[1] += alfa_blended*Density_i*pw*Volume;
+   val_residual[1] += pw*Volume;
 
    /*--- Dissipation ---*/
 
diff --git a/SU2_CFD/src/output_structure.cpp b/SU2_CFD/src/output_structure.cpp
index ba2623996f7..59a4bd3a43b 100644
--- a/SU2_CFD/src/output_structure.cpp
+++ b/SU2_CFD/src/output_structure.cpp
@@ -4460,13 +4460,17 @@ void COutput::SetConvHistory_Header(ofstream *ConvHist_file, CConfig *config, un
     case SA:case SA_NEG:case SA_E: case SA_COMP: case SA_E_COMP: 
       SPRINTF (turb_resid, ",\"Res_Turb[0]\"");
       break;
-    case SST:   	SPRINTF (turb_resid, ",\"Res_Turb[0]\",\"Res_Turb[1]\""); break;
+    case SST:case SST_SUST:
+      SPRINTF (turb_resid, ",\"Res_Turb[0]\",\"Res_Turb[1]\"");
+      break;
   }
   switch (config->GetKind_Turb_Model()) {
     case SA:case SA_NEG:case SA_E: case SA_COMP: case SA_E_COMP:
       SPRINTF (adj_turb_resid, ",\"Res_AdjTurb[0]\"");
       break;
-    case SST:   	SPRINTF (adj_turb_resid, ",\"Res_AdjTurb[0]\",\"Res_AdjTurb[1]\""); break;
+    case SST:case SST_SUST:
+      SPRINTF (adj_turb_resid, ",\"Res_AdjTurb[0]\",\"Res_AdjTurb[1]\"");
+      break;
   }
   char fem_resid[]= ",\"Res_FEM[0]\",\"Res_FEM[1]\",\"Res_FEM[2]\"";
   char heat_resid[]= ",\"Res_Heat\"";
@@ -4827,7 +4831,7 @@ void COutput::SetConvHistory_Body(ofstream *ConvHist_file,
     if (turbulent) {
       switch (config[val_iZone]->GetKind_Turb_Model()) {
         case SA: case SA_NEG: case SA_E: case SA_E_COMP: case SA_COMP: nVar_Turb = 1; break;
-        case SST:    nVar_Turb = 2; break;
+        case SST: case SST_SUST: nVar_Turb = 2; break;
       }
     }
     if (transition) nVar_Trans = 2;
@@ -4846,7 +4850,7 @@ void COutput::SetConvHistory_Body(ofstream *ConvHist_file,
     if (turbulent) {
       switch (config[val_iZone]->GetKind_Turb_Model()) {
         case SA: case SA_NEG: case SA_E: case SA_E_COMP: case SA_COMP: nVar_AdjTurb = 1; break;
-        case SST:    nVar_AdjTurb = 2; break;
+        case SST: case SST_SUST: nVar_AdjTurb = 2; break;
       }
     }
     if (weakly_coupled_heat) nVar_AdjHeat = 1;
@@ -5709,7 +5713,7 @@ void COutput::SetConvHistory_Body(ofstream *ConvHist_file,
 
             switch (config[val_iZone]->GetKind_Turb_Model()) {
               case SA: case SA_NEG: case SA_E: case SA_E_COMP: case SA_COMP:        cout << "       Res[nu]"; break;
-              case SST:	      cout << "     Res[kine]" << "     Res[omega]"; break;
+              case SST:	case SST_SUST: cout << "     Res[kine]" << "     Res[omega]"; break;
             }
 
             if (weakly_coupled_heat) {
@@ -6839,7 +6843,11 @@ void COutput::SpecialOutput_ForcesBreakdown(CSolver *****solver, CGeometry ****g
         switch (Kind_Turb_Model) {
           case SA:        Breakdown_file << "Spalart Allmaras" << "\n"; break;
           case SA_NEG:    Breakdown_file << "Negative Spalart Allmaras" << "\n"; break;
+          case SA_E:      Breakdown_file << "Edwards Spalart Allmaras" << "\n"; break;
+          case SA_COMP:   Breakdown_file << "Compressibility Correction Spalart Allmaras" << "\n"; break;
+          case SA_E_COMP: Breakdown_file << "Compressibility Correction Edwards Spalart Allmaras" << "\n"; break;
           case SST:       Breakdown_file << "Menter's SST"     << "\n"; break;
+          case SST_SUST:  Breakdown_file << "Menter's SST with sustaining terms" << "\n"; break;
         }
         break;
     }
@@ -11670,7 +11678,7 @@ void COutput::WriteTurboPerfConvHistory(CConfig *config){
   string inMarker_Tag, outMarker_Tag, inMarkerTag_Mix;
   unsigned short nZone       = config->GetnZone();
   bool turbulent = ((config->GetKind_Solver() == RANS) || (config->GetKind_Solver() == DISC_ADJ_RANS));
-  bool menter_sst       = (config->GetKind_Turb_Model() == SST);
+  bool menter_sst = (config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST);
 
   unsigned short nBladesRow, nStages;
   unsigned short iStage;
@@ -12714,7 +12722,7 @@ void COutput::LoadLocalData_Flow(CConfig *config, CGeometry *geometry, CSolver *
   Variable_Names.push_back("Energy");
   
   if (SecondIndex != NONE) {
-    if (config->GetKind_Turb_Model() == SST) {
+    if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
       Variable_Names.push_back("TKE");
       Variable_Names.push_back("Omega");
     } else {
@@ -12740,7 +12748,7 @@ void COutput::LoadLocalData_Flow(CConfig *config, CGeometry *geometry, CSolver *
       Variable_Names.push_back("Limiter_Energy");
       
       if (SecondIndex != NONE) {
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
           Variable_Names.push_back("Limiter_TKE");
           Variable_Names.push_back("Limiter_Omega");
         } else {
@@ -12762,7 +12770,7 @@ void COutput::LoadLocalData_Flow(CConfig *config, CGeometry *geometry, CSolver *
       Variable_Names.push_back("Residual_Energy");
       
       if (SecondIndex != NONE) {
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
           Variable_Names.push_back("Residual_TKE");
           Variable_Names.push_back("Residual_Omega");
         } else {
@@ -13357,7 +13365,7 @@ void COutput::LoadLocalData_IncFlow(CConfig *config, CGeometry *geometry, CSolve
   if (energy || weakly_coupled_heat) Variable_Names.push_back("Temperature");
 
   if (SecondIndex != NONE) {
-    if (config->GetKind_Turb_Model() == SST) {
+    if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
       Variable_Names.push_back("TKE");
       Variable_Names.push_back("Omega");
     } else {
@@ -13384,7 +13392,7 @@ void COutput::LoadLocalData_IncFlow(CConfig *config, CGeometry *geometry, CSolve
         Variable_Names.push_back("Limiter_Temperature");
 
       if (SecondIndex != NONE) {
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
           Variable_Names.push_back("Limiter_TKE");
           Variable_Names.push_back("Limiter_Omega");
         } else {
@@ -13407,7 +13415,7 @@ void COutput::LoadLocalData_IncFlow(CConfig *config, CGeometry *geometry, CSolve
         Variable_Names.push_back("Residual_Temperature");
 
       if (SecondIndex != NONE) {
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
           Variable_Names.push_back("Residual_TKE");
           Variable_Names.push_back("Residual_Omega");
         } else {
@@ -13950,7 +13958,7 @@ void COutput::LoadLocalData_AdjFlow(CConfig *config, CGeometry *geometry, CSolve
     Variable_Names.push_back("Adjoint_Energy");
   }
   if (SecondIndex != NONE) {
-    if (config->GetKind_Turb_Model() == SST) {
+    if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
       Variable_Names.push_back("Adjoint_TKE");
       Variable_Names.push_back("Adjoint_Omega");
     } else {
@@ -13996,7 +14004,7 @@ void COutput::LoadLocalData_AdjFlow(CConfig *config, CGeometry *geometry, CSolve
         Variable_Names.push_back("Limiter_Adjoint_Energy");
       }
       if (SecondIndex != NONE) {
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
           Variable_Names.push_back("Limiter_Adjoint_TKE");
           Variable_Names.push_back("Limiter_Adjoint_Omega");
         } else {
@@ -14025,7 +14033,7 @@ void COutput::LoadLocalData_AdjFlow(CConfig *config, CGeometry *geometry, CSolve
         Variable_Names.push_back("Residual_Adjoint_Energy");
       }
       if (SecondIndex != NONE) {
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)) {
           Variable_Names.push_back("Residual_Adjoint_TKE");
           Variable_Names.push_back("Residual_Adjoint_Omega");
         } else {
@@ -18812,7 +18820,7 @@ void COutput::Write_InletFile_Flow(CConfig *config, CGeometry *geometry, CSolver
         nVar_Turb = 1;
         turb_val[0] = solver[TURB_SOL]->GetNuTilde_Inf();
         break;
-      case SST:
+      case SST: case SST_SUST:
         nVar_Turb = 2;
         turb_val[0] = solver[TURB_SOL]->GetTke_Inf();
         turb_val[1] = solver[TURB_SOL]->GetOmega_Inf();
diff --git a/SU2_CFD/src/solver_adjoint_mean.cpp b/SU2_CFD/src/solver_adjoint_mean.cpp
index e3c398e6846..8d33d44db44 100644
--- a/SU2_CFD/src/solver_adjoint_mean.cpp
+++ b/SU2_CFD/src/solver_adjoint_mean.cpp
@@ -4216,7 +4216,7 @@ void CAdjEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
         visc_numerics->SetAdjointVar(Psi_domain, Psi_outlet);
 
         /*--- Turbulent kinetic energy ---*/
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 
 
@@ -5613,7 +5613,7 @@ void CAdjNSSolver::Viscous_Sensitivity(CGeometry *geometry, CSolver **solver_con
             
             /*--- Turbulent kinetic energy ---*/
             
-            if (config->GetKind_Turb_Model() == SST)
+            if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
               val_turb_ke = solver_container[TURB_SOL]->node[iPoint]->GetSolution(0);
             else
               val_turb_ke = 0.0;
diff --git a/SU2_CFD/src/solver_adjoint_turbulent.cpp b/SU2_CFD/src/solver_adjoint_turbulent.cpp
index 82fc48c8587..0e2cf8e5eaf 100644
--- a/SU2_CFD/src/solver_adjoint_turbulent.cpp
+++ b/SU2_CFD/src/solver_adjoint_turbulent.cpp
@@ -50,9 +50,13 @@ CAdjTurbSolver::CAdjTurbSolver(CGeometry *geometry, CConfig *config, unsigned sh
   
   /*--- Dimension of the problem  ---*/
   switch (config->GetKind_Turb_Model()) {
-    case SA :     nVar = 1; break;
-    case SA_NEG : nVar = 1; break;
-    case SST :    nVar = 2; break;
+    case SA :        nVar = 1; break;
+    case SA_NEG :    nVar = 1; break;
+    case SA_E :      nVar = 1; break;
+    case SA_COMP :   nVar = 1; break;
+    case SA_E_COMP : nVar = 1; break;
+    case SST :       nVar = 2; break;
+    case SST_SUST :  nVar = 2; break;
   }
   
   /*--- Initialize nVarGrad for deallocation ---*/
diff --git a/SU2_CFD/src/solver_direct_mean.cpp b/SU2_CFD/src/solver_direct_mean.cpp
index 0bf1abc05f1..17506a1b735 100644
--- a/SU2_CFD/src/solver_direct_mean.cpp
+++ b/SU2_CFD/src/solver_direct_mean.cpp
@@ -2230,7 +2230,7 @@ void CEulerSolver::SetNondimensionalization(CConfig *config, unsigned short iMes
   bool grid_movement      = config->GetGrid_Movement();
   bool gravity            = config->GetGravityForce();
   bool turbulent          = (config->GetKind_Solver() == RANS) || (config->GetKind_Solver() == DISC_ADJ_RANS);
-  bool tkeNeeded          = ((turbulent) && (config->GetKind_Turb_Model() == SST));
+  bool tkeNeeded          = (turbulent && (config->GetKind_Turb_Model() == SST || config->GetKind_Turb_Model() == SST_SUST));
   bool free_stream_temp   = (config->GetKind_FreeStreamOption() == TEMPERATURE_FS);
   bool reynolds_init      = (config->GetKind_InitOption() == REYNOLDS);
   bool aeroelastic        = config->GetAeroelastic_Simulation();
@@ -7834,7 +7834,7 @@ void CEulerSolver::BC_Euler_Wall(CGeometry *geometry, CSolver **solver_container
   bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
   bool grid_movement = config->GetGrid_Movement();
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   
   Normal = new su2double[nDim];
   NormalArea = new su2double[nDim];
@@ -8021,9 +8021,8 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container,
   bool implicit       = config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT;
   bool grid_movement  = config->GetGrid_Movement();
   bool viscous        = config->GetViscous();
-  bool tkeNeeded = (((config->GetKind_Solver() == RANS ) ||
-                     (config->GetKind_Solver() == DISC_ADJ_RANS))
-                    && (config->GetKind_Turb_Model() == SST));
+  bool tkeNeeded = (((config->GetKind_Solver() == RANS ) || (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
     
   su2double *Normal = new su2double[nDim];
   
@@ -8238,7 +8237,7 @@ void CEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_container,
         
         /*--- Turbulent kinetic energy ---*/
         
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                               solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
@@ -8285,7 +8284,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container,
   bool viscous              = config->GetViscous();
   bool gravity = (config->GetGravityForce());
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   
   su2double *Normal, *FlowDirMix, TangVelocity, NormalVelocity;
   Normal = new su2double[nDim];
@@ -8726,7 +8725,7 @@ void CEulerSolver::BC_Riemann(CGeometry *geometry, CSolver **solver_container,
         
         /*--- Turbulent kinetic energy ---*/
 
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
         /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -8796,7 +8795,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain
   bool viscous              = config->GetViscous();
   bool gravity = (config->GetGravityForce());
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-      (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
 
   su2double *Normal, *turboNormal, *UnitNormal, *FlowDirMix, FlowDirMixMag, *turboVelocity;
   Normal = new su2double[nDim];
@@ -9234,7 +9233,7 @@ void CEulerSolver::BC_TurboRiemann(CGeometry *geometry, CSolver **solver_contain
 
           /*--- Turbulent kinetic energy ---*/
           
-          if (config->GetKind_Turb_Model() == SST)
+          if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
             visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 
           /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -10129,7 +10128,7 @@ void CEulerSolver::BC_Giles(CGeometry *geometry, CSolver **solver_container,
         
         /*--- Turbulent kinetic energy ---*/
         
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
         /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -10199,7 +10198,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
   unsigned short Kind_Inlet = config->GetKind_Inlet();
   string Marker_Tag         = config->GetMarker_All_TagBound(val_marker);
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   su2double *Normal = new su2double[nDim];
     
   /*--- Loop over all the vertices on this boundary marker ---*/
@@ -10457,7 +10456,7 @@ void CEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
 //
 //        /*--- Turbulent kinetic energy ---*/
 //
-//        if (config->GetKind_Turb_Model() == SST)
+//        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
 //          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -10500,7 +10499,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
   string Marker_Tag       = config->GetMarker_All_TagBound(val_marker);
   bool gravity = (config->GetGravityForce());
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   su2double *Normal = new su2double[nDim];
   
   /*--- Loop over all the vertices on this boundary marker ---*/
@@ -10637,7 +10636,7 @@ void CEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
 //
 //        /*--- Turbulent kinetic energy ---*/
 //
-//        if (config->GetKind_Turb_Model() == SST)
+//        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
 //          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -10676,7 +10675,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con
   bool grid_movement  = config->GetGrid_Movement();
   string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   su2double *Normal = new su2double[nDim];
   su2double *Velocity = new su2double[nDim];
   
@@ -10779,7 +10778,7 @@ void CEulerSolver::BC_Supersonic_Inlet(CGeometry *geometry, CSolver **solver_con
 //
 //        /*--- Turbulent kinetic energy ---*/
 //
-//        if (config->GetKind_Turb_Model() == SST)
+//        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
 //          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -10898,7 +10897,7 @@ void CEulerSolver::BC_Supersonic_Outlet(CGeometry *geometry, CSolver **solver_co
  //
  //       /*--- Turbulent kinetic energy ---*/
  //
- //       if (config->GetKind_Turb_Model() == SST)
+ //       if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
  //         visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -10941,7 +10940,7 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai
   su2double Gas_Constant = config->GetGas_ConstantND();
   string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   su2double Baseline_Press = 0.75 * config->GetPressure_FreeStreamND();
   bool Engine_HalfModel = config->GetEngine_HalfModel();
 
@@ -11122,7 +11121,7 @@ void CEulerSolver::BC_Engine_Inflow(CGeometry *geometry, CSolver **solver_contai
 //
 //        /*--- Turbulent kinetic energy ---*/
 //
-//        if (config->GetKind_Turb_Model() == SST)
+//        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
 //          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -11161,7 +11160,7 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta
   bool grid_movement        = config->GetGrid_Movement();
   string Marker_Tag = config->GetMarker_All_TagBound(val_marker);
   bool tkeNeeded = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   su2double DampingFactor = config->GetDamp_Engine_Exhaust();
   su2double Baseline_Press = 0.75 * config->GetPressure_FreeStreamND();
   
@@ -11377,7 +11376,7 @@ void CEulerSolver::BC_Engine_Exhaust(CGeometry *geometry, CSolver **solver_conta
 //
 //        /*--- Turbulent kinetic energy ---*/
 //
-//        if (config->GetKind_Turb_Model() == SST)
+//        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
 //          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -11626,7 +11625,7 @@ void CEulerSolver::BC_Sym_Plane(CGeometry      *geometry,
         visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected);
         
         /*--- Turbulent kinetic energy. ---*/
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                               solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
@@ -11778,7 +11777,7 @@ void CEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_cont
 
               /*--- Turbulent kinetic energy ---*/
 
-              if (config->GetKind_Turb_Model() == SST)
+              if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
                 visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 
               /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -11975,7 +11974,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C
   su2double Gas_Constant  = config->GetGas_ConstantND();
   bool grid_movement      = config->GetGrid_Movement();
   bool tkeNeeded          = (((config->GetKind_Solver() == RANS )|| (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                            (config->GetKind_Turb_Model() == SST));
+                             ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   bool ratio              = (config->GetActDisk_Jump() == RATIO);
   su2double SecondaryFlow = config->GetSecondaryFlow_ActDisk();
   
@@ -12337,7 +12336,7 @@ void CEulerSolver::BC_ActDisk(CGeometry *geometry, CSolver **solver_container, C
 //
 //        /*--- Turbulent kinetic energy ---*/
 //
-//        if (config->GetKind_Turb_Model() == SST)
+//        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
 //          visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 //
 //        /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -12835,7 +12834,7 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig
    (that could appear in the restart file before the grid velocities). ---*/
   unsigned short turbVars = 0;
   if (turbulent){
-    if (turb_model == SST) turbVars = 2;  
+    if ((turb_model == SST) || (turb_model == SST_SUST)) turbVars = 2;
     else turbVars = 1;
   }
 
@@ -13298,7 +13297,7 @@ void CEulerSolver::TurboAverageProcess(CSolver **solver, CGeometry *geometry, CC
   su2double MachTest, soundSpeed;
   bool turbulent = ((config->GetKind_Solver() == RANS) || (config->GetKind_Solver() == DISC_ADJ_RANS));
   bool spalart_allmaras = (config->GetKind_Turb_Model() == SA);
-  bool menter_sst       = (config->GetKind_Turb_Model() == SST);
+  bool menter_sst       = ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST));
 
   /*-- Variables declaration and allocation ---*/
   Velocity            = new su2double[nDim];
@@ -15198,7 +15197,7 @@ unsigned long CNSSolver::SetPrimitive_Variables(CSolver **solver_container, CCon
   unsigned short turb_model = config->GetKind_Turb_Model();
   bool RightSol = true;
   
-  bool tkeNeeded            = (turb_model == SST);
+  bool tkeNeeded = ((turb_model == SST) || (turb_model == SST_SUST)) ;
 
   for (iPoint = 0; iPoint < nPoint; iPoint ++) {
     
@@ -15469,7 +15468,7 @@ void CNSSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_container
     
     /*--- Turbulent kinetic energy ---*/
     
-    if (config->GetKind_Turb_Model() == SST)
+    if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
       numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                      solver_container[TURB_SOL]->node[jPoint]->GetSolution(0));
     
diff --git a/SU2_CFD/src/solver_direct_mean_fem.cpp b/SU2_CFD/src/solver_direct_mean_fem.cpp
index 4ba280de9d5..4d9ef83e13a 100644
--- a/SU2_CFD/src/solver_direct_mean_fem.cpp
+++ b/SU2_CFD/src/solver_direct_mean_fem.cpp
@@ -803,7 +803,7 @@ void CFEM_DG_EulerSolver::SetNondimensionalization(CConfig        *config,
   bool viscous            = config->GetViscous();
   bool grid_movement      = config->GetGrid_Movement();
   bool turbulent          = (config->GetKind_Solver() == FEM_RANS) || (config->GetKind_Solver() == FEM_LES);
-  bool tkeNeeded          = ((turbulent) && (config->GetKind_Turb_Model() == SST));
+  bool tkeNeeded          = ((turbulent) && ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   bool free_stream_temp   = (config->GetKind_FreeStreamOption() == TEMPERATURE_FS);
   bool reynolds_init      = (config->GetKind_InitOption() == REYNOLDS);
 
diff --git a/SU2_CFD/src/solver_direct_mean_inc.cpp b/SU2_CFD/src/solver_direct_mean_inc.cpp
index 9670cd473f5..b0029b0ee70 100644
--- a/SU2_CFD/src/solver_direct_mean_inc.cpp
+++ b/SU2_CFD/src/solver_direct_mean_inc.cpp
@@ -824,7 +824,7 @@ void CIncEulerSolver::SetNondimensionalization(CConfig *config, unsigned short i
   bool grid_movement = config->GetGrid_Movement();
   bool turbulent     = ((config->GetKind_Solver() == RANS) ||
                         (config->GetKind_Solver() == DISC_ADJ_RANS));
-  bool tkeNeeded     = ((turbulent) && (config->GetKind_Turb_Model() == SST));
+  bool tkeNeeded     = ((turbulent) && ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   bool energy        = config->GetEnergy_Equation();
   bool boussinesq    = (config->GetKind_DensityModel() == BOUSSINESQ);
 
@@ -4608,9 +4608,8 @@ void CIncEulerSolver::BC_Euler_Wall(CGeometry *geometry, CSolver **solver_contai
   su2double Density = 0.0, Pressure = 0.0, *Normal = NULL, Area, *NormalArea, turb_ke;
   
   bool implicit = (config->GetKind_TimeIntScheme_Flow() == EULER_IMPLICIT);
-  bool tkeNeeded = (((config->GetKind_Solver() == RANS ) ||
-                     (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
-                    (config->GetKind_Turb_Model() == SST));
+  bool tkeNeeded = (((config->GetKind_Solver() == RANS ) || (config->GetKind_Solver() == DISC_ADJ_RANS)) &&
+                    ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST)));
   
   Normal     = new su2double[nDim];
   NormalArea = new su2double[nDim];
@@ -4793,7 +4792,7 @@ void CIncEulerSolver::BC_Far_Field(CGeometry *geometry, CSolver **solver_contain
         
         /*--- Turbulent kinetic energy ---*/
         
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                               solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
@@ -5046,7 +5045,7 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container,
         
         /*--- Turbulent kinetic energy ---*/
         
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                               solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
@@ -5256,7 +5255,7 @@ void CIncEulerSolver::BC_Outlet(CGeometry *geometry, CSolver **solver_container,
         
         /*--- Turbulent kinetic energy ---*/
         
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                               solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
@@ -5503,7 +5502,7 @@ void CIncEulerSolver::BC_Sym_Plane(CGeometry      *geometry,
         visc_numerics->SetPrimVarGradient(node[iPoint]->GetGradient_Primitive(), Grad_Reflected);
         
         /*--- Turbulent kinetic energy. ---*/
-        if (config->GetKind_Turb_Model() == SST)
+        if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
           visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                               solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
         
@@ -5641,7 +5640,7 @@ void CIncEulerSolver::BC_Fluid_Interface(CGeometry *geometry, CSolver **solver_c
 
               /*--- Turbulent kinetic energy ---*/
 
-              if (config->GetKind_Turb_Model() == SST)
+              if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
                 visc_numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0), solver_container[TURB_SOL]->node[iPoint]->GetSolution(0));
 
               /*--- Set the wall shear stress values (wall functions) to -1 (no evaluation using wall functions) ---*/
@@ -6544,7 +6543,7 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf
    (that could appear in the restart file before the grid velocities). ---*/
   unsigned short turbVars = 0;
   if (turbulent){
-    if (turb_model == SST) turbVars = 2;
+    if ((turb_model == SST) || (turb_model == SST_SUST)) turbVars = 2;
     else turbVars = 1;
   }
   
@@ -7553,7 +7552,7 @@ unsigned long CIncNSSolver::SetPrimitive_Variables(CSolver **solver_container, C
   unsigned short turb_model = config->GetKind_Turb_Model();
   bool physical = true;
   
-  bool tkeNeeded = (turb_model == SST);
+  bool tkeNeeded = ((turb_model == SST) || (turb_model == SST_SUST));
   
   for (iPoint = 0; iPoint < nPoint; iPoint++) {
     
@@ -7836,7 +7835,7 @@ void CIncNSSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contai
     
     /*--- Turbulent kinetic energy ---*/
     
-    if (config->GetKind_Turb_Model() == SST)
+    if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
       numerics->SetTurbKineticEnergy(solver_container[TURB_SOL]->node[iPoint]->GetSolution(0),
                                      solver_container[TURB_SOL]->node[jPoint]->GetSolution(0));
     
diff --git a/SU2_CFD/src/solver_direct_turbulent.cpp b/SU2_CFD/src/solver_direct_turbulent.cpp
index a92e8b8f03c..0d7dccd1e99 100644
--- a/SU2_CFD/src/solver_direct_turbulent.cpp
+++ b/SU2_CFD/src/solver_direct_turbulent.cpp
@@ -235,7 +235,7 @@ void CTurbSolver::Viscous_Residual(CGeometry *geometry, CSolver **solver_contain
     numerics->SetTurbVarGradient(node[iPoint]->GetGradient(), node[jPoint]->GetGradient());
     
     /*--- Menter's first blending function (only SST)---*/
-    if (config->GetKind_Turb_Model() == SST)
+    if ((config->GetKind_Turb_Model() == SST) || (config->GetKind_Turb_Model() == SST_SUST))
       numerics->SetF1blending(node[iPoint]->GetF1blending(), node[jPoint]->GetF1blending());
     
     /*--- Compute residual, and Jacobians ---*/
@@ -420,7 +420,7 @@ void CTurbSolver::ImplicitEuler_Iteration(CGeometry *geometry, CSolver **solver_
         
         break;
 
-      case SST:
+      case SST: case SST_SUST:
         
         for (iPoint = 0; iPoint < nPointDomain; iPoint++) {
           
@@ -477,6 +477,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
   bool grid_movement = config->GetGrid_Movement();
   
   bool incompressible = (config->GetKind_Regime() == INCOMPRESSIBLE);
+  unsigned short turbModel = config->GetKind_Turb_Model();
 
   /*--- Store the physical time step ---*/
   
@@ -506,7 +507,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
       /*--- Compute the dual time-stepping source term based on the chosen
        time discretization scheme (1st- or 2nd-order).---*/
       
-      if (config->GetKind_Turb_Model() == SST) {
+      if ((turbModel == SST) || (turbModel == SST_SUST)) {
         
         /*--- If this is the SST model, we need to multiply by the density
          in order to get the conservative variables ---*/
@@ -596,7 +597,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
       
       /*--- Multiply by density at node i for the SST model ---*/
       
-      if (config->GetKind_Turb_Model() == SST) {
+      if ((turbModel == SST) || (turbModel == SST_SUST)) {
         if (incompressible) Density_n = solver_container[FLOW_SOL]->node[iPoint]->GetDensity(); // Temporary fix
         else Density_n = solver_container[FLOW_SOL]->node[iPoint]->GetSolution_time_n()[0];
         for (iVar = 0; iVar < nVar; iVar++)
@@ -613,7 +614,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
       
       /*--- Multiply by density at node j for the SST model ---*/
       
-      if (config->GetKind_Turb_Model() == SST) {
+      if ((turbModel == SST) || (turbModel == SST_SUST)) {
         if (incompressible) Density_n = solver_container[FLOW_SOL]->node[jPoint]->GetDensity(); // Temporary fix
         else Density_n = solver_container[FLOW_SOL]->node[jPoint]->GetSolution_time_n()[0];
         for (iVar = 0; iVar < nVar; iVar++)
@@ -655,7 +656,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
         
         /*--- Multiply by density at node i for the SST model ---*/
         
-        if (config->GetKind_Turb_Model() == SST) {
+        if ((turbModel == SST) || (turbModel == SST_SUST)) {
           if (incompressible) Density_n = solver_container[FLOW_SOL]->node[iPoint]->GetDensity(); // Temporary fix
           else Density_n = solver_container[FLOW_SOL]->node[iPoint]->GetSolution_time_n()[0];
           for (iVar = 0; iVar < nVar; iVar++)
@@ -694,7 +695,7 @@ void CTurbSolver::SetResidual_DualTime(CGeometry *geometry, CSolver **solver_con
        introduction of the GCL term above, the remainder of the source residual
        due to the time discretization has a new form.---*/
       
-      if (config->GetKind_Turb_Model() == SST) {
+      if ((turbModel == SST) || (turbModel == SST_SUST)) {
         
         /*--- If this is the SST model, we need to multiply by the density
          in order to get the conservative variables ---*/
diff --git a/config_template.cfg b/config_template.cfg
index 14578550bc4..2160c5b8a9d 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -17,7 +17,7 @@
 %                               POISSON_EQUATION)
 PHYSICAL_PROBLEM= EULER
 %
-% Specify turbulence model (NONE, SA, SA_NEG, SST, SA_E, SA_COMP, SA_E_COMP)
+% Specify turbulence model (NONE, SA, SA_NEG, SST, SA_E, SA_COMP, SA_E_COMP, SST_SUST)
 KIND_TURB_MODEL= NONE
 %
 % Specify subgrid scale model(NONE, IMPLICIT_LES, SMAGORINSKY, WALE, VREMAN)

From b19fadc4146f19f39c97d91775a59ff8127b5100 Mon Sep 17 00:00:00 2001
From: vdweide <e.t.a.vanderweide@utwente.nl>
Date: Tue, 20 Aug 2019 22:56:10 +0200
Subject: [PATCH 2/6] Fixed a bug and added a couple of regression test cases
 for SST_SUST

---
 SU2_CFD/src/drivers/CDriver.cpp               |   1 -
 .../naca0012/naca0012_SST_SUST.cfg            | 257 +++++++++++++++++
 TestCases/parallel_regression.py              |  33 +++
 .../rans/naca0012/turb_NACA0012_sst_sust.cfg  | 264 ++++++++++++++++++
 .../rans/rae2822/turb_SST_SUST_RAE2822.cfg    | 242 ++++++++++++++++
 TestCases/serial_regression.py                |  38 ++-
 6 files changed, 832 insertions(+), 3 deletions(-)
 create mode 100644 TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg
 create mode 100644 TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg
 create mode 100644 TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg

diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp
index af13ab35b72..1a59bf76839 100644
--- a/SU2_CFD/src/drivers/CDriver.cpp
+++ b/SU2_CFD/src/drivers/CDriver.cpp
@@ -2730,7 +2730,6 @@ void CDriver::Numerics_Postprocessing(CNumerics *****numerics,
           
           /*--- Definition of the boundary condition method ---*/
           for (iMGlevel = 0; iMGlevel <= config->GetnMGLevels(); iMGlevel++)
-      default: SU2_MPI::Error("Specified turbulence model unavailable or none selected", CURRENT_FUNCTION); break;
             delete numerics[val_iInst][iMGlevel][FLOW_SOL][CONV_BOUND_TERM];
           
         }
diff --git a/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg b/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg
new file mode 100644
index 00000000000..015b1173aba
--- /dev/null
+++ b/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg
@@ -0,0 +1,257 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: 2D NACA 0012 Airfoil Validation Case (incompressible)      %
+%                   http://turbmodels.larc.nasa.gov/naca0012_val_sa.html       %
+% Author: Francisco Palacios	                                               %
+% Institution: Stanford University                                             %
+% Date: Feb 18th, 2013                                                         %
+% File Version 6.2.0 "Falcon"                                                %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+% Physical governing equations (EULER, NAVIER_STOKES,
+%                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
+%                               POISSON_EQUATION)
+PHYSICAL_PROBLEM= NAVIER_STOKES
+%
+% Specify turbulent model (NONE, SA, SA_NEG, SST)
+KIND_TURB_MODEL= SST_SUST
+%
+% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
+MATH_PROBLEM= DIRECT
+%
+% Regime type (COMPRESSIBLE, INCOMPRESSIBLE, FREESURFACE)
+REGIME_TYPE= INCOMPRESSIBLE
+%
+% Restart solution (NO, YES)
+RESTART_SOL= NO
+
+% ------------------------- UNSTEADY SIMULATION -------------------------------%
+%
+% Unsteady simulation (NO, TIME_STEPPING, DUAL_TIME_STEPPING-1ST_ORDER, 
+%                      DUAL_TIME_STEPPING-2ND_ORDER)
+UNSTEADY_SIMULATION= NO
+
+% ---------------- INCOMPRESSIBLE FLOW CONDITION DEFINITION -------------------%
+%
+% Initial density for incompressible flows (1.2886 kg/m^3 by default)
+INC_DENSITY_INIT= 2.13163
+%
+% Initial velocity for incompressible flows (1.0,0,0 m/s by default)
+INC_VELOCITY_INIT= ( 52.1572, 0.0, 0.0 )
+
+% --------------------------- VISCOSITY MODEL ---------------------------------%
+%
+% Viscosity model (SUTHERLAND, CONSTANT_VISCOSITY).
+VISCOSITY_MODEL= CONSTANT_VISCOSITY
+%
+% Molecular Viscosity that would be constant (1.716E-5 by default)
+MU_CONSTANT= 1.853e-05
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+% Reference origin for moment computation
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+%
+% Reference length for pitching, rolling, and yawing non-dimensional moment
+REF_LENGTH= 1.0
+%
+% Reference area for force coefficients (0 implies automatic calculation)
+REF_AREA= 1.0
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+% Navier-Stokes wall boundary marker(s) (NONE = no marker)
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+%
+% Farfield boundary marker(s) (NONE = no marker)
+MARKER_FAR= ( farfield )
+%
+% Marker(s) of the surface to be plotted or designed
+MARKER_PLOTTING= ( airfoil )
+%
+% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+%
+% Courant-Friedrichs-Lewy condition of the finest grid
+CFL_NUMBER= 10.0
+%
+% Adaptive CFL number (NO, YES)
+CFL_ADAPT= NO
+%
+% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
+%                                        CFL max value )
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+%
+% Number of total iterations
+EXT_ITER= 2500
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+%
+% Coefficient for the limiter
+VENKAT_LIMITER_COEFF= 0.01
+%
+% Coefficient for the sharp edges limiter
+ADJ_SHARP_LIMITER_COEFF= 3.0
+%
+% Reference coefficient (sensitivity) for detecting sharp edges.
+REF_SHARP_EDGES= 3.0
+%
+% Remove sharp edges from the sensitivity evaluation (NO, YES)
+SENS_REMOVE_SHARP= NO
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+% Linear solver for implicit formulations (BCGSTAB, FGMRES)
+LINEAR_SOLVER= FGMRES
+%
+% Preconditioner of the Krylov linear solver (JACOBI, LINELET, LU_SGS)
+LINEAR_SOLVER_PREC= LU_SGS
+%
+% Minimum error of the linear solver for implicit formulations
+LINEAR_SOLVER_ERROR= 1E-4
+%
+% Max number of iterations of the linear solver for the implicit formulation
+LINEAR_SOLVER_ITER= 5
+
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+% Multi-Grid Levels (0 = no multi-grid)
+MGLEVEL= 0
+%
+% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE)
+MGCYCLE= V_CYCLE
+%
+% Multi-grid pre-smoothing level
+MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
+%
+% Multi-grid post-smoothing level
+MG_POST_SMOOTH= ( 0, 0, 0, 0 )
+%
+% Jacobi implicit smoothing of the correction
+MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
+%
+% Damping factor for the residual restriction
+MG_DAMP_RESTRICTION= 0.75
+%
+% Damping factor for the correction prolongation
+MG_DAMP_PROLONGATION= 0.75
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC,
+%                              TURKEL_PREC, MSW)
+CONV_NUM_METHOD_FLOW= FDS
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_FLOW= YES
+%
+% Slope limiter (VENKATAKRISHNAN, MINMOD)
+SLOPE_LIMITER_FLOW= VENKATAKRISHNAN
+%
+% 2nd and 4th order artificial dissipation coefficients
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+%
+% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+% Convective numerical method (SCALAR_UPWIND)
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_TURB= NO
+%
+% Slope limiter (VENKATAKRISHNAN, MINMOD)
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+%
+% Time discretization (EULER_IMPLICIT)
+TIME_DISCRE_TURB= EULER_IMPLICIT
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence criteria (CAUCHY, RESIDUAL)
+%
+CONV_CRITERIA= RESIDUAL
+%
+% Residual reduction (order of magnitude with respect to the initial value)
+RESIDUAL_REDUCTION= 6
+%
+% Min value of the residual (log10 of the residual)
+RESIDUAL_MINVAL= -10
+%
+% Start convergence criteria at iteration number
+STARTCONV_ITER= 10
+%
+% Number of elements to apply the criteria
+CAUCHY_ELEMS= 100
+%
+% Epsilon to control the series convergence
+CAUCHY_EPS= 1E-6
+%
+% Function to apply the criteria (LIFT, DRAG, NEARFIELD_PRESS, SENS_GEOMETRY, 
+% 	      	    		 SENS_MACH, DELTA_LIFT, DELTA_DRAG)
+CAUCHY_FUNC_FLOW= DRAG
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+% Mesh input file
+MESH_FILENAME= mesh_NACA0012_turb_897x257.su2
+%
+% Mesh input file format (SU2, CGNS, NETCDF_ASCII)
+MESH_FORMAT= SU2
+%
+% Mesh output file
+MESH_OUT_FILENAME= mesh_out.su2
+%
+% Restart flow input file
+SOLUTION_FLOW_FILENAME= solution_flow.dat
+%
+% Restart adjoint input file
+SOLUTION_ADJ_FILENAME= solution_adj.dat
+%
+% Output file format (PARAVIEW, TECPLOT, STL)
+OUTPUT_FORMAT= TECPLOT
+%
+% Output file convergence history (w/o extension) 
+CONV_FILENAME= history
+%
+% Output file restart flow
+RESTART_FLOW_FILENAME= restart_flow.dat
+%
+% Output file restart adjoint
+RESTART_ADJ_FILENAME= restart_adj.dat
+%
+% Output file flow (w/o extension) variables
+VOLUME_FLOW_FILENAME= flow
+%
+% Output file adjoint (w/o extension) variables
+VOLUME_ADJ_FILENAME= adjoint
+%
+% Output objective function gradient (using continuous adjoint)
+GRAD_OBJFUNC_FILENAME= of_grad.dat
+%
+% Output file surface flow coefficient (w/o extension)
+SURFACE_FLOW_FILENAME= surface_flow
+%
+% Output file surface adjoint coefficient (w/o extension)
+SURFACE_ADJ_FILENAME= surface_adjoint
+%
+% Writing solution file frequency
+WRT_SOL_FREQ= 100
+%
+% Writing convergence history frequency
+WRT_CON_FREQ= 1
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
index e7494c9b6fb..18318dd89d2 100644
--- a/TestCases/parallel_regression.py
+++ b/TestCases/parallel_regression.py
@@ -215,6 +215,17 @@ def main():
     rae2822_sst.tol       = 0.00001
     test_list.append(rae2822_sst)
 
+    # RAE2822 SST_SUST
+    rae2822_sst_sust           = TestCase('rae2822_sst_sust')
+    rae2822_sst_sust.cfg_dir   = "rans/rae2822"
+    rae2822_sst_sust.cfg_file  = "turb_SST_SUST_RAE2822.cfg"
+    rae2822_sst_sust.test_iter = 20
+    rae2822_sst_sust.test_vals = [-2.407053, 4.916635, 0.827705, 0.053331] #last 4 columns
+    rae2822_sst_sust.su2_exec  = "parallel_computation.py -f"
+    rae2822_sst_sust.timeout   = 1600
+    rae2822_sst_sust.tol       = 0.00001
+    test_list.append(rae2822_sst_sust)
+
     # Flat plate
     turb_flatplate           = TestCase('turb_flatplate')
     turb_flatplate.cfg_dir   = "rans/flatplate"
@@ -270,6 +281,17 @@ def main():
     turb_naca0012_sst.tol       = 0.00001
     test_list.append(turb_naca0012_sst)
 
+    # NACA0012 (SST_SUST, FUN3D finest grid results: CL=1.0840, CD=0.01253)
+    turb_naca0012_sst_sust           = TestCase('turb_naca0012_sst_sust')
+    turb_naca0012_sst_sust.cfg_dir   = "rans/naca0012"
+    turb_naca0012_sst_sust.cfg_file  = "turb_NACA0012_sst_sust.cfg"
+    turb_naca0012_sst_sust.test_iter = 10
+    turb_naca0012_sst_sust.test_vals = [-13.280836, -5.645158, 1.022304, 0.019539] #last 4 columns
+    turb_naca0012_sst_sust.su2_exec  = "parallel_computation.py -f"
+    turb_naca0012_sst_sust.timeout   = 3200
+    turb_naca0012_sst_sust.tol       = 0.00001
+    test_list.append(turb_naca0012_sst_sust)
+
     # PROPELLER
     propeller           = TestCase('propeller')
     propeller.cfg_dir   = "rans/propeller"
@@ -374,6 +396,17 @@ def main():
     inc_turb_naca0012.timeout   = 1600
     inc_turb_naca0012.tol       = 0.00001
     test_list.append(inc_turb_naca0012)
+
+    # NACA0012, SST_SUST
+    inc_turb_naca0012_sst_sust           = TestCase('inc_turb_naca0012_sst_sust')
+    inc_turb_naca0012_sst_sust.cfg_dir   = "incomp_rans/naca0012"
+    inc_turb_naca0012_sst_sust.cfg_file  = "naca0012_SST_SUST.cfg"
+    inc_turb_naca0012_sst_sust.test_iter = 20
+    inc_turb_naca0012_sst_sust.test_vals = [-7.277551, 0.147212, -0.000000, 0.311977] #last 4 columns
+    inc_turb_naca0012_sst_sust.su2_exec  = "parallel_computation.py -f"
+    inc_turb_naca0012_sst_sust.timeout   = 1600
+    inc_turb_naca0012_sst_sust.tol       = 0.00001
+    test_list.append(inc_turb_naca0012_sst_sust)
     
     ####################
     ### DG-FEM Euler ###
diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg
new file mode 100644
index 00000000000..f6c0013f09f
--- /dev/null
+++ b/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg
@@ -0,0 +1,264 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: 2D NACA 0012 Airfoil Validation Case (compressible)        %
+%                   http://turbmodels.larc.nasa.gov/naca0012_val_sa.html       %
+% Author: Thomas D. Economon                                                   %
+% Institution: Stanford University                                             %
+% Date: Feb 18th, 2013                                                         %
+% File Version 3.2.1 "eagle"                                                   %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+% Physical governing equations (EULER, NAVIER_STOKES,
+%                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
+%                               POISSON_EQUATION)                             
+PHYSICAL_PROBLEM= NAVIER_STOKES
+%
+% Specify turbulent model (NONE, SA, SA_NEG, SST)
+KIND_TURB_MODEL= SST_SUST
+%
+% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
+MATH_PROBLEM= DIRECT
+%
+% Restart solution (NO, YES)
+RESTART_SOL= YES
+%
+% Write binary restart files (YES, NO)
+WRT_BINARY_RESTART= NO
+%
+% Read binary restart files (YES, NO)
+READ_BINARY_RESTART= NO
+%
+% Minimize the required output memory
+LOW_MEMORY_OUTPUT= NO
+
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+% Mach number (non-dimensional, based on the free-stream values)
+MACH_NUMBER= 0.15
+%
+% Angle of attack (degrees, only for compressible flows)
+AOA= 10.0
+%
+% Free-stream temperature (288.15 K by default)
+FREESTREAM_TEMPERATURE= 300.0
+%
+% Reynolds number (non-dimensional, based on the free-stream values)
+REYNOLDS_NUMBER= 6.0E6
+%
+% Reynolds length (1 m by default)
+REYNOLDS_LENGTH= 1.0
+%
+% Free-stream Turbulence Intensity, sqrt(2 k_inf/3)/U_inf (0.0005)
+FREESTREAM_TURBULENCEINTENSITY = 0.001
+%
+% Free-stream Turbulent to Laminar viscosity ratio
+FREESTREAM_TURB2LAMVISCRATIO = 0.1
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+% Reference origin for moment computation
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+%
+% Reference length for pitching, rolling, and yawing non-dimensional moment
+REF_LENGTH= 1.0
+%
+% Reference area for force coefficients (0 implies automatic calculation)
+REF_AREA= 1.0
+%
+% Flow non-dimensionalization (DIMENSIONAL, FREESTREAM_PRESS_EQ_ONE,
+%                              FREESTREAM_VEL_EQ_MACH, FREESTREAM_VEL_EQ_ONE)
+REF_DIMENSIONALIZATION= FREESTREAM_PRESS_EQ_ONE
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+% Navier-Stokes wall boundary marker(s) (NONE = no marker)
+MARKER_HEATFLUX= ( airfoil, 0.0 )
+%
+% Farfield boundary marker(s) (NONE = no marker)
+MARKER_FAR= ( farfield )
+%
+% Marker(s) of the surface to be plotted or designed
+MARKER_PLOTTING= ( airfoil )
+%
+% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated
+MARKER_MONITORING= ( airfoil )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+%
+% Courant-Friedrichs-Lewy condition of the finest grid
+CFL_NUMBER= 10.0
+%
+% Max Delta time
+MAX_DELTA_TIME= 1E10
+%
+% Adaptive CFL number (NO, YES)
+CFL_ADAPT= NO
+%
+% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
+%                                        CFL max value )
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+%
+% Number of total iterations
+EXT_ITER= 99999
+
+% ----------------------- SLOPE LIMITER DEFINITION ----------------------------%
+%
+% Coefficient for the limiter
+VENKAT_LIMITER_COEFF= 0.03
+%
+% Freeze the value of the limiter after a number of iterations
+LIMITER_ITER= 99999
+
+% ------------------------ LINEAR SOLVER DEFINITION ---------------------------%
+%
+% Linear solver or smoother for implicit formulations (BCGSTAB, FGMRES, SMOOTHER)
+LINEAR_SOLVER= FGMRES
+%
+% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI)
+LINEAR_SOLVER_PREC= LU_SGS
+%
+% Minimum error of the linear solver for implicit formulations
+LINEAR_SOLVER_ERROR= 1E-10
+%
+% Max number of iterations of the linear solver for the implicit formulation
+LINEAR_SOLVER_ITER= 2
+
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+% Multi-Grid Levels (0 = no multi-grid)
+MGLEVEL= 0
+%
+% Multigrid pre-smoothing level
+MG_PRE_SMOOTH= ( 1, 1, 1, 1, 1, 1 )
+%
+% Multigrid post-smoothing level
+MG_POST_SMOOTH= ( 0, 0, 0, 0, 0, 0 )
+%
+% Jacobi implicit smoothing of the correction
+MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0, 0, 0 )
+%
+% Damping factor for the residual restriction
+MG_DAMP_RESTRICTION= 0.75
+%
+% Damping factor for the correction prolongation
+MG_DAMP_PROLONGATION= 0.75
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC,
+%                              TURKEL_PREC, MSW)
+CONV_NUM_METHOD_FLOW= ROE
+%
+% Spatial numerical order integration (1ST_ORDER, 2ND_ORDER, 2ND_ORDER_LIMITER)
+MUSCL_FLOW= YES
+%
+% Slope limiter (VENKATAKRISHNAN, MINMOD)
+SLOPE_LIMITER_FLOW= VENKATAKRISHNAN
+%
+% 2nd and 4th order artificial dissipation coefficients
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+%
+% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+% Convective numerical method (SCALAR_UPWIND)
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_TURB= NO
+%
+% Slope limiter (VENKATAKRISHNAN, MINMOD)
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+%
+% Time discretization (EULER_IMPLICIT)
+TIME_DISCRE_TURB= EULER_IMPLICIT
+%
+% Reduction factor of the CFL coefficient in the turbulence problem
+CFL_REDUCTION_TURB= 1.0
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence criteria (CAUCHY, RESIDUAL)
+%
+CONV_CRITERIA= RESIDUAL
+%
+% Residual reduction (order of magnitude with respect to the initial value)
+RESIDUAL_REDUCTION= 10
+%
+% Min value of the residual (log10 of the residual)
+RESIDUAL_MINVAL= -12
+%
+% Start convergence criteria at iteration number
+STARTCONV_ITER= 10
+%
+% Number of elements to apply the criteria
+CAUCHY_ELEMS= 100
+%
+% Epsilon to control the series convergence
+CAUCHY_EPS= 1E-6
+%
+% Function to apply the criteria (LIFT, DRAG, NEARFIELD_PRESS, SENS_GEOMETRY, 
+% 	      	    		 SENS_MACH, DELTA_LIFT, DELTA_DRAG)
+CAUCHY_FUNC_FLOW= DRAG
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+% Mesh input file
+MESH_FILENAME= n0012_225-65.su2
+%
+% Mesh input file format (SU2, CGNS, NETCDF_ASCII)
+MESH_FORMAT= SU2
+%
+% Mesh output file
+MESH_OUT_FILENAME= mesh_out.su2
+%
+% Restart flow input file
+SOLUTION_FLOW_FILENAME= solution_flow_sst_sust.dat
+%
+% Restart adjoint input file
+SOLUTION_ADJ_FILENAME= solution_adj.dat
+%
+% Output file format (PARAVIEW, TECPLOT, STL)
+OUTPUT_FORMAT= TECPLOT
+%
+% Output file convergence history (w/o extension) 
+CONV_FILENAME= history
+%
+% Output file restart flow
+RESTART_FLOW_FILENAME= restart_flow.dat
+%
+% Output file restart adjoint
+RESTART_ADJ_FILENAME= restart_adj.dat
+%
+% Output file flow (w/o extension) variables
+VOLUME_FLOW_FILENAME= flow
+%
+% Output file adjoint (w/o extension) variables
+VOLUME_ADJ_FILENAME= adjoint
+%
+% Output objective function gradient (using continuous adjoint)
+GRAD_OBJFUNC_FILENAME= of_grad.dat
+%
+% Output file surface flow coefficient (w/o extension)
+SURFACE_FLOW_FILENAME= surface_flow
+%
+% Output file surface adjoint coefficient (w/o extension)
+SURFACE_ADJ_FILENAME= surface_adjoint
+%
+% Writing solution file frequency
+WRT_SOL_FREQ= 10000
+%
+% Writing convergence history frequency
+WRT_CON_FREQ= 1
diff --git a/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg b/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg
new file mode 100644
index 00000000000..ba9da316c35
--- /dev/null
+++ b/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg
@@ -0,0 +1,242 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                                                                              %
+% SU2 configuration file                                                       %
+% Case description: Transonic simulation RAE2822 (RANS)                        %
+% Author: Francisco Palacios                                                   %
+% Institution: Stanford University                                             %
+% Date: 5/15/2013                                                              %
+% File Version 6.2.0 "Falcon"                                                %
+%                                                                              %
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ------------- DIRECT, ADJOINT, AND LINEARIZED PROBLEM DEFINITION ------------%
+%
+% Physical governing equations (EULER, NAVIER_STOKES,
+%                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
+%                               POISSON_EQUATION)                          
+PHYSICAL_PROBLEM= NAVIER_STOKES
+%
+% Specify turbulent model (NONE, SA, SA_NEG, SST, SST_SUST)
+KIND_TURB_MODEL= SST_SUST
+%
+% Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
+MATH_PROBLEM= DIRECT
+%
+% Restart solution (NO, YES)
+RESTART_SOL= NO
+
+% -------------------- COMPRESSIBLE FREE-STREAM DEFINITION --------------------%
+%
+% Mach number (non-dimensional, based on the free-stream values)
+MACH_NUMBER= 0.729
+%
+% Angle of attack (degrees, only for compressible flows)
+AOA= 2.31
+%
+% Free-stream temperature (288.15 K by default)
+FREESTREAM_TEMPERATURE= 288.15
+%
+% Reynolds number (non-dimensional, based on the free-stream values)
+REYNOLDS_NUMBER= 6.5E6
+%
+% Reynolds length (1 m by default)
+REYNOLDS_LENGTH= 1.0
+
+% Free-stream Turbulence Intensity, sqrt(2 k_inf/3)/U_inf (0.0005)
+FREESTREAM_TURBULENCEINTENSITY = 0.001
+
+% Free-stream Turbulent to Laminar viscosity ratio
+FREESTREAM_TURB2LAMVISCRATIO = 100.0
+
+% ---------------------- REFERENCE VALUE DEFINITION ---------------------------%
+%
+% Reference origin for moment computation
+REF_ORIGIN_MOMENT_X = 0.25
+REF_ORIGIN_MOMENT_Y = 0.00
+REF_ORIGIN_MOMENT_Z = 0.00
+%
+% Reference length for pitching, rolling, and yawing non-dimensional moment
+REF_LENGTH= 1.0
+%
+% Reference area for force coefficients (0 implies automatic calculation)
+REF_AREA= 1.0
+
+% -------------------- BOUNDARY CONDITION DEFINITION --------------------------%
+%
+% Navier-Stokes wall boundary marker(s) (NONE = no marker)
+MARKER_HEATFLUX= ( AIRFOIL, 0.0 )
+%
+% Farfield boundary marker(s) (NONE = no marker)
+MARKER_FAR= ( FARFIELD )
+%
+% Marker(s) of the surface to be plotted or designed
+MARKER_PLOTTING= ( AIRFOIL )
+%
+% Marker(s) of the surface where the functional (Cd, Cl, etc.) will be evaluated
+MARKER_MONITORING= ( AIRFOIL )
+
+% ------------- COMMON PARAMETERS DEFINING THE NUMERICAL METHOD ---------------%
+%
+% Numerical method for spatial gradients (GREEN_GAUSS, WEIGHTED_LEAST_SQUARES)
+NUM_METHOD_GRAD= WEIGHTED_LEAST_SQUARES
+%
+% Courant-Friedrichs-Lewy condition of the finest grid
+CFL_NUMBER= 2.5
+%
+% Adaptive CFL number (NO, YES)
+CFL_ADAPT= NO
+%
+% Parameters of the adaptive CFL number (factor down, factor up, CFL min value,
+%                                        CFL max value )
+CFL_ADAPT_PARAM= ( 1.5, 0.5, 1.0, 100.0 )
+%
+% Number of total iterations
+EXT_ITER= 99999
+%
+% Linear solver for the implicit formulation (BCGSTAB, FGMRES)
+LINEAR_SOLVER= FGMRES
+%
+% Preconditioner of the Krylov linear solver (ILU, LU_SGS, LINELET, JACOBI)
+LINEAR_SOLVER_PREC= LU_SGS
+%
+% Min error of the linear solver for the implicit formulation
+LINEAR_SOLVER_ERROR= 1E-6
+%
+% Max number of iterations of the linear solver for the implicit formulation
+LINEAR_SOLVER_ITER= 2
+
+% -------------------------- MULTIGRID PARAMETERS -----------------------------%
+%
+% Multi-Grid Levels (0 = no multi-grid)
+MGLEVEL= 3
+%
+% Multi-grid cycle (V_CYCLE, W_CYCLE, FULLMG_CYCLE)
+MGCYCLE= W_CYCLE
+%
+% Multi-grid pre-smoothing level
+MG_PRE_SMOOTH= ( 1, 2, 3, 3 )
+%
+% Multi-grid post-smoothing level
+MG_POST_SMOOTH= ( 0, 0, 0, 0 )
+%
+% Jacobi implicit smoothing of the correction
+MG_CORRECTION_SMOOTH= ( 0, 0, 0, 0 )
+%
+% Damping factor for the residual restriction
+MG_DAMP_RESTRICTION= 0.95
+%
+% Damping factor for the correction prolongation
+MG_DAMP_PROLONGATION= 0.95
+
+% -------------------- FLOW NUMERICAL METHOD DEFINITION -----------------------%
+%
+% Convective numerical method (JST, LAX-FRIEDRICH, CUSP, ROE, AUSM, HLLC,
+%                              TURKEL_PREC, MSW)
+CONV_NUM_METHOD_FLOW= JST
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the flow equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_FLOW= YES
+%
+% Slope limiter (NONE, VENKATAKRISHNAN, VENKATAKRISHNAN_WANG,
+%                BARTH_JESPERSEN, VAN_ALBADA_EDGE)
+SLOPE_LIMITER_FLOW= NONE
+%
+% Coefficient for the limiter (smooth regions)
+VENKAT_LIMITER_COEFF= 0.03
+%
+% 2nd and 4th order artificial dissipation coefficients
+JST_SENSOR_COEFF= ( 0.5, 0.02 )
+%
+% Time discretization (RUNGE-KUTTA_EXPLICIT, EULER_IMPLICIT, EULER_EXPLICIT)
+TIME_DISCRE_FLOW= EULER_IMPLICIT
+
+% -------------------- TURBULENT NUMERICAL METHOD DEFINITION ------------------%
+%
+% Convective numerical method (SCALAR_UPWIND)
+CONV_NUM_METHOD_TURB= SCALAR_UPWIND
+%
+% Monotonic Upwind Scheme for Conservation Laws (TVD) in the turbulence equations.
+%           Required for 2nd order upwind schemes (NO, YES)
+MUSCL_TURB= NO
+%
+% Slope limiter (VENKATAKRISHNAN, BARTH)
+SLOPE_LIMITER_TURB= VENKATAKRISHNAN
+%
+% Time discretization (EULER_IMPLICIT)
+TIME_DISCRE_TURB= EULER_IMPLICIT
+
+% --------------------------- CONVERGENCE PARAMETERS --------------------------%
+%
+% Convergence criteria (CAUCHY, RESIDUAL)
+%
+CONV_CRITERIA= RESIDUAL
+%
+% Residual reduction (order of magnitude with respect to the initial value)
+RESIDUAL_REDUCTION= 6
+%
+% Min value of the residual (log10 of the residual)
+RESIDUAL_MINVAL= -8
+%
+% Start convergence criteria at iteration number
+STARTCONV_ITER= 10
+%
+% Number of elements to apply the criteria
+CAUCHY_ELEMS= 100
+%
+% Epsilon to control the series convergence
+CAUCHY_EPS= 1E-6
+%
+% Function to apply the criteria (LIFT, DRAG, NEARFIELD_PRESS, SENS_GEOMETRY, 
+% 	      	    		 SENS_MACH, DELTA_LIFT, DELTA_DRAG)
+CAUCHY_FUNC_FLOW= DRAG
+
+% ------------------------- INPUT/OUTPUT INFORMATION --------------------------%
+%
+% Mesh input file
+MESH_FILENAME= mesh_RAE2822_turb.su2
+%
+% Mesh input file format (SU2, CGNS, NETCDF_ASCII)
+MESH_FORMAT= SU2
+%
+% Mesh output file
+MESH_OUT_FILENAME= mesh_out.su2
+%
+% Restart flow input file
+SOLUTION_FLOW_FILENAME= solution_flow.dat
+%
+% Restart adjoint input file
+SOLUTION_ADJ_FILENAME= solution_adj.dat
+%
+% Output file format (PARAVIEW, TECPLOT, STL)
+OUTPUT_FORMAT= TECPLOT
+%
+% Output file convergence history (w/o extension) 
+CONV_FILENAME= history
+%
+% Output file restart flow
+RESTART_FLOW_FILENAME= restart_flow.dat
+%
+% Output file restart adjoint
+RESTART_ADJ_FILENAME= restart_adj.dat
+%
+% Output file flow (w/o extension) variables
+VOLUME_FLOW_FILENAME= flow
+%
+% Output file adjoint (w/o extension) variables
+VOLUME_ADJ_FILENAME= adjoint
+%
+% Output objective function gradient (using continuous adjoint)
+GRAD_OBJFUNC_FILENAME= of_grad.dat
+%
+% Output file surface flow coefficient (w/o extension)
+SURFACE_FLOW_FILENAME= surface_flow
+%
+% Output file surface adjoint coefficient (w/o extension)
+SURFACE_ADJ_FILENAME= surface_adjoint
+%
+% Writing solution file frequency
+WRT_SOL_FREQ= 250
+%
+% Writing convergence history frequency
+WRT_CON_FREQ= 1
diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index ee6b0197975..a876c454cb3 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -213,6 +213,17 @@ def main():
     rae2822_sst.tol       = 0.00001
     test_list.append(rae2822_sst)
 
+    # RAE2822 SST_SUST
+    rae2822_sst_sust           = TestCase('rae2822_sst_sust')
+    rae2822_sst_sust.cfg_dir   = "rans/rae2822"
+    rae2822_sst_sust.cfg_file  = "turb_SST_SUST_RAE2822.cfg"
+    rae2822_sst_sust.test_iter = 20
+    rae2822_sst_sust.test_vals = [-2.401009, 4.909241, 0.825023, 0.052675] #last 4 columns
+    rae2822_sst_sust.su2_exec  = "SU2_CFD"
+    rae2822_sst_sust.timeout   = 1600
+    rae2822_sst_sust.tol       = 0.00001
+    test_list.append(rae2822_sst_sust)
+
     # Flat plate
     turb_flatplate           = TestCase('turb_flatplate')
     turb_flatplate.cfg_dir   = "rans/flatplate"
@@ -268,6 +279,18 @@ def main():
     turb_naca0012_sst.tol       = 0.00001
     test_list.append(turb_naca0012_sst)
 
+    # NACA0012 (SST, FUN3D results for finest grid: CL=1.0840, CD=0.01253)
+    turb_naca0012_sst_sust           = TestCase('turb_naca0012_sst_sust')
+    turb_naca0012_sst_sust.cfg_dir   = "rans/naca0012"
+    turb_naca0012_sst_sust.cfg_file  = "turb_NACA0012_sst_sust.cfg"
+    turb_naca0012_sst_sust.test_iter = 10
+    turb_naca0012_sst_sust.test_vals = [-13.280975, -5.646508, 1.022304, 0.019539] #last 4 columns
+    turb_naca0012_sst_sust.su2_exec  = "SU2_CFD"
+    turb_naca0012_sst_sust.timeout   = 3200
+    turb_naca0012_sst_sust.tol       = 0.00001
+    test_list.append(turb_naca0012_sst_sust)
+
+
     # PROPELLER 
     propeller           = TestCase('propeller')
     propeller.cfg_dir   = "rans/propeller"
@@ -362,7 +385,7 @@ def main():
     ### Incompressible RANS  ###
     ############################
 
-    # NACA0012
+    # NACA0012, SA
     inc_turb_naca0012           = TestCase('inc_turb_naca0012')
     inc_turb_naca0012.cfg_dir   = "incomp_rans/naca0012"
     inc_turb_naca0012.cfg_file  = "naca0012.cfg"
@@ -372,7 +395,18 @@ def main():
     inc_turb_naca0012.timeout   = 1600
     inc_turb_naca0012.tol       = 0.00001
     test_list.append(inc_turb_naca0012)
-    
+
+    # NACA0012, SST_SUST
+    inc_turb_naca0012_sst_sust           = TestCase('inc_turb_naca0012_sst_sust')
+    inc_turb_naca0012_sst_sust.cfg_dir   = "incomp_rans/naca0012"
+    inc_turb_naca0012_sst_sust.cfg_file  = "naca0012_SST_SUST.cfg"
+    inc_turb_naca0012_sst_sust.test_iter = 20
+    inc_turb_naca0012_sst_sust.test_vals = [-7.277591, 0.146922, 0.000021, 0.311963] #last 4 columns
+    inc_turb_naca0012_sst_sust.su2_exec  = "SU2_CFD"
+    inc_turb_naca0012_sst_sust.timeout   = 1600
+    inc_turb_naca0012_sst_sust.tol       = 0.00001
+    test_list.append(inc_turb_naca0012_sst_sust)
+
     ####################
     ### DG-FEM Euler ###
     ####################

From 1ca22aea0b2c7aa0d6ea5ea56a9cc8daa2a61793 Mon Sep 17 00:00:00 2001
From: vdweide <e.t.a.vanderweide@utwente.nl>
Date: Tue, 20 Aug 2019 23:01:24 +0200
Subject: [PATCH 3/6] A small change in the comments of the serial regression
 test

---
 TestCases/serial_regression.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py
index a876c454cb3..9a2bd71e296 100644
--- a/TestCases/serial_regression.py
+++ b/TestCases/serial_regression.py
@@ -279,7 +279,7 @@ def main():
     turb_naca0012_sst.tol       = 0.00001
     test_list.append(turb_naca0012_sst)
 
-    # NACA0012 (SST, FUN3D results for finest grid: CL=1.0840, CD=0.01253)
+    # NACA0012 (SST_SUST, FUN3D results for finest grid: CL=1.0840, CD=0.01253)
     turb_naca0012_sst_sust           = TestCase('turb_naca0012_sst_sust')
     turb_naca0012_sst_sust.cfg_dir   = "rans/naca0012"
     turb_naca0012_sst_sust.cfg_file  = "turb_NACA0012_sst_sust.cfg"

From efb7e999ff25e6913a1919014b7f07d3358ad75c Mon Sep 17 00:00:00 2001
From: vdweide <e.t.a.vanderweide@utwente.nl>
Date: Thu, 22 Aug 2019 07:38:25 +0200
Subject: [PATCH 4/6] Changed travis such that only I get an error message when
 something goes south

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 428f6d61f7f..23500b6cf8e 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,7 +18,7 @@ compiler:
 notifications:
     email:
         recipients:
-            - su2code-dev@lists.stanford.edu
+            - e.t.a.vanderweide@utwente.nl
   
 branches:
     only:

From 0fa9f534907ae9d937fd9535c18118ff4414f5da Mon Sep 17 00:00:00 2001
From: vdweide <e.t.a.vanderweide@utwente.nl>
Date: Thu, 22 Aug 2019 11:50:49 +0200
Subject: [PATCH 5/6] Reverted travis

---
 .travis.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.travis.yml b/.travis.yml
index 23500b6cf8e..428f6d61f7f 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -18,7 +18,7 @@ compiler:
 notifications:
     email:
         recipients:
-            - e.t.a.vanderweide@utwente.nl
+            - su2code-dev@lists.stanford.edu
   
 branches:
     only:

From 7bbb549da11e4f0584e9a3581f9114f1f584c813 Mon Sep 17 00:00:00 2001
From: vdweide <e.t.a.vanderweide@utwente.nl>
Date: Sun, 25 Aug 2019 15:43:16 +0200
Subject: [PATCH 6/6] Updated the config files for the SST_SUST test cases

---
 TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg | 5 +----
 TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg   | 2 +-
 TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg     | 2 +-
 3 files changed, 3 insertions(+), 6 deletions(-)

diff --git a/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg b/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg
index 015b1173aba..7f26864c322 100644
--- a/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg
+++ b/TestCases/incomp_rans/naca0012/naca0012_SST_SUST.cfg
@@ -15,7 +15,7 @@
 % Physical governing equations (EULER, NAVIER_STOKES,
 %                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
 %                               POISSON_EQUATION)
-PHYSICAL_PROBLEM= NAVIER_STOKES
+SOLVER= INC_NAVIER_STOKES
 %
 % Specify turbulent model (NONE, SA, SA_NEG, SST)
 KIND_TURB_MODEL= SST_SUST
@@ -23,9 +23,6 @@ KIND_TURB_MODEL= SST_SUST
 % Mathematical problem (DIRECT, CONTINUOUS_ADJOINT)
 MATH_PROBLEM= DIRECT
 %
-% Regime type (COMPRESSIBLE, INCOMPRESSIBLE, FREESURFACE)
-REGIME_TYPE= INCOMPRESSIBLE
-%
 % Restart solution (NO, YES)
 RESTART_SOL= NO
 
diff --git a/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg b/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg
index f6c0013f09f..6badd339536 100644
--- a/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg
+++ b/TestCases/rans/naca0012/turb_NACA0012_sst_sust.cfg
@@ -15,7 +15,7 @@
 % Physical governing equations (EULER, NAVIER_STOKES,
 %                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
 %                               POISSON_EQUATION)                             
-PHYSICAL_PROBLEM= NAVIER_STOKES
+SOLVER= NAVIER_STOKES
 %
 % Specify turbulent model (NONE, SA, SA_NEG, SST)
 KIND_TURB_MODEL= SST_SUST
diff --git a/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg b/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg
index ba9da316c35..00de695c341 100644
--- a/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg
+++ b/TestCases/rans/rae2822/turb_SST_SUST_RAE2822.cfg
@@ -14,7 +14,7 @@
 % Physical governing equations (EULER, NAVIER_STOKES,
 %                               WAVE_EQUATION, HEAT_EQUATION, FEM_ELASTICITY,
 %                               POISSON_EQUATION)                          
-PHYSICAL_PROBLEM= NAVIER_STOKES
+SOLVER= NAVIER_STOKES
 %
 % Specify turbulent model (NONE, SA, SA_NEG, SST, SST_SUST)
 KIND_TURB_MODEL= SST_SUST