Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix dissipation in transition model and update inlet profile (initial profile from config) #1268

Merged
merged 13 commits into from
Apr 21, 2021
Merged
3 changes: 3 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -94,3 +94,6 @@ build/

# ninja binary (build system)
ninja

# Ignore vscode folder
.vscode/
7 changes: 6 additions & 1 deletion SU2_CFD/include/CMarkerProfileReaderFVM.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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. */
Expand Down Expand Up @@ -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.
Expand Down
57 changes: 44 additions & 13 deletions SU2_CFD/src/CMarkerProfileReaderFVM.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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. ---*/

Expand All @@ -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;
Comment on lines 47 to +51
Copy link
Contributor

@TobiKattmann TobiKattmann Apr 28, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess this could be done with an initializer list for the ctor .. but I am unsure if one is to be preferred over the other. The core guidelines dont seem to have a strong opinion either. CE Playground


/* Attempt to open the specified file. */
ifstream profile_file;
Expand All @@ -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++) {
Expand All @@ -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;
}
}
Comment on lines 130 to 134
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

in that case I would remove the whole else statement ... unless there are further plans


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);
Expand All @@ -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. ---*/

Expand Down Expand Up @@ -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;

Expand All @@ -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;

Expand All @@ -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;

}

}
Expand All @@ -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);
Expand Down
7 changes: 5 additions & 2 deletions SU2_CFD/src/numerics/transition.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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) ---*/
Expand Down Expand Up @@ -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) --*/
Expand Down
3 changes: 2 additions & 1 deletion SU2_CFD/src/numerics/turbulent/turb_sources.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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();
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

👍 every now and then CFD-online helps us find some bugs :)


su2double nu_t = (TurbVar_i[0]*fv1); //S-A variable

Expand Down
9 changes: 6 additions & 3 deletions SU2_CFD/src/output/CFlowIncOutput.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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));
Comment on lines 596 to +601
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I guess one could move this config->GetEnergy_Equation() out of the if nDim and make the acces solver[FLOW_SOL]->LinSysRes(iPoint, nDim+1)

}

switch(config->GetKind_Turb_Model()){
Expand Down
90 changes: 87 additions & 3 deletions SU2_CFD/src/solvers/CSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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);
Comment on lines +3359 to +3360
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is going to sound like neat picking in this particular context but bear with me (you never know the next piece of code you will work on).
From an efficiency standpoint you should always move as little data as possible (talking about the spaces in setw and passing numbers as strings).
This is particularly relevant for strings because of something called "small string optimization", where up to a certain number of characters, strings do not allocate memory on the heap.
In this example, this also means that you are forcing the file format from the outside, this should be the responsibility of the class (profile reader in this case) consider for example that we may want to use it for CSV formats.


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. ---*/

Expand Down
9 changes: 6 additions & 3 deletions SU2_CFD/src/solvers/CTransLMSolver.cpp
Original file line number Diff line number Diff line change
@@ -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"
*
Expand Down Expand Up @@ -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() {
Expand Down Expand Up @@ -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) {
Expand Down
2 changes: 1 addition & 1 deletion TestCases/parallel_regression.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
Loading