Fix dissipation in transition model and update inlet profile (initial…

… profile from config) (#1268)

* fixed turbulence dissipation bug factor 100 in transition turbulence model
* do not write residual temperature to file if energy eq. is off
* inlet profile is now initialized from config options and added column names
* changed residuals of schubauer-klebanoff case
* change output profile file to example_filename

.gitignore

@@ -94,3 +94,6 @@ build/
 
 # ninja binary (build system)
 ninja
+
+# Ignore vscode folder
+.vscode/

SU2_CFD/include/CMarkerProfileReaderFVM.hpp

@@ -74,6 +74,9 @@ class CMarkerProfileReaderFVM {
 
   string filename;  /*!< \brief File name of the marker profile file. */
 
+  vector<string> columnNames; /*!< \brief string containing all the names of the columns, one for each marker */
+  vector<string> columnValues; /*!< \brief initial values for the profile, in string format */
+
   vector<string> profileTags;  /*!< \brief Auxiliary structure for holding the string names of the markers in a profile file. */
 
   vector<unsigned long> numberOfRowsInProfile;  /*!< \brief Auxiliary structure for holding the number of rows for a particular marker in a profile file. */
@@ -113,7 +116,9 @@ class CMarkerProfileReaderFVM {
                           CConfig        *val_config,
                           string         val_filename,
                           unsigned short val_kind_marker,
-                          unsigned short val_number_vars);
+                          unsigned short val_number_vars,
+                          vector<string> val_columnNames,
+                          vector<string> val_columnValues);
 
   /*!
    * \brief Destructor of the CMeshReaderFVM class.

SU2_CFD/src/CMarkerProfileReaderFVM.cpp

@@ -31,7 +31,9 @@ CMarkerProfileReaderFVM::CMarkerProfileReaderFVM(CGeometry      *val_geometry,
                                                  CConfig        *val_config,
                                                  string         val_filename,
                                                  unsigned short val_kind_marker,
-                                                 unsigned short val_number_vars) {
+                                                 unsigned short val_number_vars,
+                                                 vector<string> val_columnNames,
+                                                 vector<string> val_columnValues) {
 
   /*--- Store input values and pointers to class data. ---*/
 
@@ -45,6 +47,8 @@ CMarkerProfileReaderFVM::CMarkerProfileReaderFVM(CGeometry      *val_geometry,
   filename     = val_filename;
   markerType   = val_kind_marker;
   numberOfVars = val_number_vars;
+  columnNames  = val_columnNames;
+  columnValues = val_columnValues;
 
   /* Attempt to open the specified file. */
   ifstream profile_file;
@@ -71,21 +75,27 @@ void CMarkerProfileReaderFVM::ReadMarkerProfile() {
 
   ifstream profile_file;
   profile_file.open(filename.data(), ios::in);
+  bool nmarkFound = false; 
+  unsigned long skip = 0;
 
   /*--- Identify the markers and data set in the profile file ---*/
 
   string text_line;
+  /*--- We search the file until we find the keyword NMARK=. Data before NMARK= will be ignored. 
+        This allows for some information in a header that will be ignored by the profile reader. ---*/
   while (getline (profile_file, text_line)) {
-
+    /*--- read NMARK ---*/
     string::size_type position = text_line.find ("NMARK=",0);
     if (position != string::npos) {
+      nmarkFound = true;
       text_line.erase (0,6); numberOfProfiles = atoi(text_line.c_str());
 
       numberOfRowsInProfile.resize(numberOfProfiles);
       numberOfColumnsInProfile.resize(numberOfProfiles);
 
       for (unsigned short iMarker = 0 ; iMarker < numberOfProfiles; iMarker++) {
 
+        /*--- read MARKER_TAG ---*/
         getline (profile_file, text_line);
         text_line.erase (0,11);
         for (unsigned short iChar = 0; iChar < 20; iChar++) {
@@ -95,24 +105,40 @@ void CMarkerProfileReaderFVM::ReadMarkerProfile() {
         }
         profileTags.push_back(text_line.c_str());
 
+        /*--- read NROW ---*/
         getline (profile_file, text_line);
         text_line.erase (0,5); numberOfRowsInProfile[iMarker] = atoi(text_line.c_str());
 
+        /*--- read NCOL ---*/
         getline (profile_file, text_line);
         text_line.erase (0,5); numberOfColumnsInProfile[iMarker] = atoi(text_line.c_str());
 
-        /*--- Skip the data. This is read in the next loop. ---*/
+        /*--- read the column format description. This line is not required, so if we cannot find it, we just continue  ---*/
+        getline (profile_file, text_line);
+        string::size_type dataheader = text_line.find ("# COORD",0);
+        if (dataheader == 0) {
+          skip = 0;
+        } else {
+          /*--- no header, but we have read a line, so we have to read one line less data ---*/
+          skip = 1;
+        }
 
-        for (unsigned long iRow = 0; iRow < numberOfRowsInProfile[iMarker]; iRow++) getline (profile_file, text_line);
+        /*--- Skip the data. This is read in the next loop. ---*/
+
+        for (unsigned long iRow = 0; iRow < (numberOfRowsInProfile[iMarker]-skip); iRow++) getline (profile_file, text_line);
 
       }
     } else {
-      SU2_MPI::Error("While opening profile file, no \"NMARK=\" specification was found", CURRENT_FUNCTION);
+      //cout << "inlet profile reader is ignoring line: " << text_line << endl;
     }
   }
 
   profile_file.close();
 
+  if (nmarkFound==false) {
+    SU2_MPI::Error("While opening profile file, no \"NMARK=\" specification was found", CURRENT_FUNCTION);
+  }
+
   /*--- Compute array bounds and offsets. Allocate data structure. ---*/
 
   profileData.resize(numberOfProfiles);
@@ -132,11 +158,14 @@ void CMarkerProfileReaderFVM::ReadMarkerProfile() {
 
       for (unsigned short iMarker = 0; iMarker < numberOfProfiles; iMarker++) {
 
-        /*--- Skip the tag, nRow, and nCol lines. ---*/
+        /*--- Skip the tag, nRow and nCol lines. ---*/
 
         getline (profile_file, text_line);
         getline (profile_file, text_line);
         getline (profile_file, text_line);
+
+        /*--- if skip=0 then we can expect column format description ---*/
+        if (skip == 0) getline (profile_file, text_line);
 
         /*--- Now read the data for each row and store. ---*/
 
@@ -419,7 +448,7 @@ void CMarkerProfileReaderFVM::WriteMarkerProfileTemplate() {
 
   if (rank == MASTER_NODE) {
 
-    ofstream node_file("profile_example.dat");
+    ofstream node_file("example_"+filename);
 
     node_file << "NMARK= " << numberOfProfiles << endl;
 
@@ -435,6 +464,9 @@ void CMarkerProfileReaderFVM::WriteMarkerProfileTemplate() {
       node_file << "NROW="        << numberOfRowsInProfile[iMarker] << endl;
       node_file << "NCOL="        << nColumns          << endl;
 
+      /*--- header line (names of the columns) --- */
+      node_file << columnNames[iMarker] << endl;
+
       node_file << setprecision(15);
       node_file << std::scientific;
 
@@ -445,10 +477,9 @@ void CMarkerProfileReaderFVM::WriteMarkerProfileTemplate() {
         for (unsigned short iDim = 0; iDim < dimension; iDim++) {
           node_file << profileCoords[iMarker][iDim][iPoint] << "\t";
         }
-        for (unsigned short iDim = 0; iDim < numberOfVars; iDim++) {
-          node_file << 0.0 << "\t";
-        }
-        node_file << endl;
+
+        node_file << columnValues[iMarker] << endl;
+
       }
 
     }
@@ -460,8 +491,8 @@ void CMarkerProfileReaderFVM::WriteMarkerProfileTemplate() {
     err << endl;
     err << "  Could not find the input file for the marker profile." << endl;
     err << "  Looked for: " << filename << "." << endl;
-    err << "  Created a template profile file with node coordinates" << endl;
-    err << "  and correct number of columns at `profile_example.dat`." << endl;
+    err << "  Created a template profile file with default values" << endl;
+    err << "  named example_" << filename << endl;
     err << "  You can use this file as a guide for making your own profile" << endl;
     err << "  specification." << endl << endl;
     SU2_MPI::Error(err.str(), CURRENT_FUNCTION);

SU2_CFD/src/numerics/transition.cpp

@@ -399,7 +399,9 @@ void CSourcePieceWise_TransLM::ComputeResidual_TransLM(su2double *val_residual,
   /* -- These lines must be manually reinserted into the differentiated routine! --*/
 //  rey  = config->GetReynolds();
 //  mach = config->GetMach();
-  tu   = config->GetTurbulenceIntensity_FreeStream();
+
+  /*--- turbulence intensity in the transition model is a percentage, so multiply by 100 ---*/
+  tu   = 100.0 * config->GetTurbulenceIntensity_FreeStream();
   //SU2_CPP2C COMMENT END
 
   /*--- Compute vorticity and strain (TODO: Update for 3D) ---*/
@@ -586,7 +588,8 @@ void CSourcePieceWise_TransLM::CSourcePieceWise_TransLM__ComputeResidual_TransLM
 //  tu = 0.0;
 //  rey  = config->GetReynolds();
 //  mach = config->GetMach();
-  tu   = config->GetTurbulenceIntensity_FreeStream();
+  /* --- turbulence intensity for transition models is in percentages, so multiply by 100 ---*/ 
+  tu   = 100.0 * config->GetTurbulenceIntensity_FreeStream();
   /*--- Compute vorticity and strain (TODO: Update for 3D) ---*/
   Vorticity = fabs(PrimVar_Grad_i[1][1] - PrimVar_Grad_i[2][0]);
   /*-- Strain = sqrt(2*Sij*Sij) --*/

SU2_CFD/src/numerics/turbulent/turb_sources.cpp

@@ -139,7 +139,8 @@ CNumerics::ResidualType<> CSourcePieceWise_TurbSA::ComputeResidual(const CConfig
       const su2double chi_1 = 0.002;
       const su2double chi_2 = 50.0;
 
-      su2double tu = config->GetTurbulenceIntensity_FreeStream();
+      /*--- turbulence intensity is u'/U so we multiply by 100 to get percentage ---*/
+      su2double tu = 100.0 * config->GetTurbulenceIntensity_FreeStream();
 
       su2double nu_t = (TurbVar_i[0]*fv1); //S-A variable
 

SU2_CFD/src/output/CFlowIncOutput.cpp

@@ -435,7 +435,8 @@ void CFlowIncOutput::SetVolumeOutputFields(CConfig *config){
   AddVolumeOutput("RES_VELOCITY-Y", "Residual_Velocity_y", "RESIDUAL", "Residual of the y-velocity component");
   if (nDim == 3)
     AddVolumeOutput("RES_VELOCITY-Z", "Residual_Velocity_z", "RESIDUAL", "Residual of the z-velocity component");
-  AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature");
+  if (config->GetEnergy_Equation())  
+    AddVolumeOutput("RES_TEMPERATURE", "Residual_Temperature", "RESIDUAL", "Residual of the temperature");
 
   switch(config->GetKind_Turb_Model()){
   case SST: case SST_SUST:
@@ -593,9 +594,11 @@ void CFlowIncOutput::LoadVolumeData(CConfig *config, CGeometry *geometry, CSolve
   SetVolumeOutputValue("RES_VELOCITY-Y", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 2));
   if (nDim == 3){
     SetVolumeOutputValue("RES_VELOCITY-Z", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 3));
-    SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 4));
+    if (config->GetEnergy_Equation())  
+      SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 4));
   } else {
-    SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 3));
+    if (config->GetEnergy_Equation())  
+      SetVolumeOutputValue("RES_TEMPERATURE", iPoint, solver[FLOW_SOL]->LinSysRes(iPoint, 3));
   }
 
   switch(config->GetKind_Turb_Model()){

SU2_CFD/src/solvers/CSolver.cpp

@@ -3294,8 +3294,6 @@ void CSolver::LoadInletProfile(CGeometry **geometry,
   const auto KIND_SOLVER = val_kind_solver;
   const auto KIND_MARKER = val_kind_marker;
 
-  /*--- Local variables ---*/
-
   const bool time_stepping = (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_1ST) ||
                              (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND) ||
                              (config->GetTime_Marching() == TIME_MARCHING::TIME_STEPPING);
@@ -3308,6 +3306,10 @@ void CSolver::LoadInletProfile(CGeometry **geometry,
   unsigned short nVar_Turb = 0;
   if (config->GetKind_Turb_Model() != NONE) nVar_Turb = solver[MESH_0][TURB_SOL]->GetnVar();
 
+  /*--- names of the columns in the profile ---*/
+  vector<string> columnNames;
+  vector<string> columnValues;
+
   /*--- Count the number of columns that we have for this flow case,
    excluding the coordinates. Here, we have 2 entries for the total
    conditions or mass flow, another nDim for the direction vector, and
@@ -3317,6 +3319,12 @@ void CSolver::LoadInletProfile(CGeometry **geometry,
 
   const unsigned short nCol_InletFile = 2 + nDim + nVar_Turb;
 
+  /*--- for incompressible flow, we can switch the energy equation off ---*/
+  /*--- for now, we write the temperature even if we are not using it ---*/
+  /*--- because a number of routines depend on the presence of the temperature field ---*/
+  //if (config->GetEnergy_Equation() ==false)
+  //nCol_InletFile = nCol_InletFile -1;
+
   /*--- Multizone problems require the number of the zone to be appended. ---*/
 
   if (nZone > 1)
@@ -3327,9 +3335,85 @@ void CSolver::LoadInletProfile(CGeometry **geometry,
   if (time_stepping)
     profile_filename = config->GetUnsteady_FileName(profile_filename, val_iter, ".dat");
 
+
+  // create vector of column names
+  for (unsigned short iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {
+
+    /*--- Skip if this is the wrong type of marker. ---*/
+    if (config->GetMarker_All_KindBC(iMarker) != KIND_MARKER) continue;
+
+    string Marker_Tag = config->GetMarker_All_TagBound(iMarker);
+    su2double p_total   = config->GetInlet_Ptotal(Marker_Tag);
+    su2double t_total   = config->GetInlet_Ttotal(Marker_Tag);
+    auto flow_dir = config->GetInlet_FlowDir(Marker_Tag);
+    std::stringstream columnName,columnValue;    
+
+    columnValue << setprecision(15);
+    columnValue << std::scientific;
+
+    columnValue << t_total << "\t" << p_total <<"\t";
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      columnValue << flow_dir[iDim] <<"\t";
+    }
+
+    columnName << "# COORD-X  " << setw(24) << "COORD-Y    " << setw(24);
+    if(nDim==3) columnName << "COORD-Z    " << setw(24); 
+
+    if (config->GetKind_Regime()==ENUM_REGIME::COMPRESSIBLE){
+      switch (config->GetKind_Inlet()) {
+        /*--- compressible conditions ---*/
+        case INLET_TYPE::TOTAL_CONDITIONS:
+          columnName << "TEMPERATURE" << setw(24) << "PRESSURE   " << setw(24); 
+          break;
+        case INLET_TYPE::MASS_FLOW:
+          columnName << "DENSITY    " << setw(24) << "VELOCITY   " << setw(24); 
+          break;
+        default:
+          SU2_MPI::Error("Unsupported INLET_TYPE.", CURRENT_FUNCTION);
+          break;        }
+    } else {
+      switch (config->GetKind_Inc_Inlet(Marker_Tag)) {
+        /*--- incompressible conditions ---*/  
+        case INLET_TYPE::VELOCITY_INLET:
+          columnName << "TEMPERATURE" << setw(24) << "VELOCITY   " << setw(24); 
+          break;
+        case INLET_TYPE::PRESSURE_INLET:
+          columnName << "TEMPERATURE" << setw(24) << "PRESSURE   " << setw(24); 
+          break;
+        default:
+          SU2_MPI::Error("Unsupported INC_INLET_TYPE.", CURRENT_FUNCTION);
+          break;
+      }
+    }
+
+    columnName << "NORMAL-X   " << setw(24) << "NORMAL-Y   " << setw(24);
+    if(nDim==3)  columnName << "NORMAL-Z   " << setw(24);
+
+    switch (config->GetKind_Turb_Model()) {
+      case NO_TURB_MODEL:
+        /*--- no turbulence model---*/
+        break;
+      case SA: case SA_NEG: case SA_E: case SA_COMP: case SA_E_COMP:
+        /*--- 1-equation turbulence model: SA ---*/
+        columnName << "NU_TILDE   " << setw(24);
+        columnValue << config->GetNuFactor_FreeStream()*config->GetViscosity_FreeStreamND()/config->GetDensity_FreeStreamND() <<"\t";
+        break;
+      case SST: case SST_SUST:
+        /*--- 2-equation turbulence model (SST) ---*/
+        columnName << "TKE        " << setw(24) << "DISSIPATION";
+        columnValue << config->GetTke_FreeStream() << "\t" << config->GetOmega_FreeStream() <<"\t";
+        break;
+    }
+
+    columnNames.push_back(columnName.str());   
+    columnValues.push_back(columnValue.str());
+
+  }
+
+
   /*--- Read the profile data from an ASCII file. ---*/
 
-  CMarkerProfileReaderFVM profileReader(geometry[MESH_0], config, profile_filename, KIND_MARKER, nCol_InletFile);
+  CMarkerProfileReaderFVM profileReader(geometry[MESH_0], config, profile_filename, KIND_MARKER, nCol_InletFile, columnNames,columnValues);
 
   /*--- Load data from the restart into correct containers. ---*/
 

SU2_CFD/src/solvers/CTransLMSolver.cpp

@@ -1,6 +1,6 @@
 /*!
  * \file CTransLMSolver.cpp
- * \brief Main subrotuines for Transition model solver.
+ * \brief Main subroutines for Langtry-Menter Transition model solver.
  * \author A. Aranake
  * \version 7.1.1 "Blackbird"
  *
@@ -30,7 +30,10 @@
 #include "../../include/variables/CTransLMVariable.hpp"
 #include "../../include/variables/CTurbSAVariable.hpp"
 
-
+/*---  This is the implementation of the Langtry-Menter transition model.
+       The main reference for this model is:Langtry, Menter, AIAA J. 47(12) 2009 
+       DOI: https://doi.org/10.2514/1.42362 ---*/
+
 CTransLMSolver::CTransLMSolver(void) : CTurbSolver() {}
 
 CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned short iMesh) : CTurbSolver() {
@@ -104,7 +107,7 @@ CTransLMSolver::CTransLMSolver(CGeometry *geometry, CConfig *config, unsigned sh
 
     /*--- Read farfield conditions from config ---*/
     Intermittency_Inf = config->GetIntermittency_FreeStream();
-    tu_Inf            = config->GetTurbulenceIntensity_FreeStream();
+    tu_Inf            = 100.0 * config->GetTurbulenceIntensity_FreeStream();
 
     /*-- Initialize REth from correlation --*/
     if (tu_Inf <= 1.3) {

TestCases/parallel_regression.py

@@ -607,7 +607,7 @@ def main():
     schubauer_klebanoff_transition.cfg_dir      = "transition/Schubauer_Klebanoff"
     schubauer_klebanoff_transition.cfg_file     = "transitional_BC_model_ConfigFile.cfg"
     schubauer_klebanoff_transition.test_iter    = 10
-    schubauer_klebanoff_transition.test_vals    = [-7.994740, -14.268433, 0.000046, 0.007987]
+    schubauer_klebanoff_transition.test_vals    = [-7.994740, -13.240225, 0.000046, 0.007987]
     schubauer_klebanoff_transition.su2_exec     = "parallel_computation.py -f"
     schubauer_klebanoff_transition.timeout      = 1600
     schubauer_klebanoff_transition.tol          = 0.00001