diff --git a/Common/include/toolboxes/C1DInterpolation.hpp b/Common/include/toolboxes/C1DInterpolation.hpp index 348a9dddefa..703f208ea86 100644 --- a/Common/include/toolboxes/C1DInterpolation.hpp +++ b/Common/include/toolboxes/C1DInterpolation.hpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2019, SU2 Contributors (cf. AUTHORS.md) @@ -38,130 +38,131 @@ using namespace std; class C1DInterpolation{ protected: - bool Point_Match = false; /*!< \brief to make sure all points from the inlet data match the vertex coordinates */ - vector X; /*!< \brief x for inlet data */ - vector Data; /*!< \brief f(x) for inlet data */ + bool Point_Match = false; /*!< \brief to make sure all points from the inlet data match the vertex coordinates */ + vector X; /*!< \brief x for inlet data */ + vector Data; /*!< \brief f(x) for inlet data */ public: - /*! - * \brief Virtual destructor of the C1DInterpolation class. - */ - virtual ~C1DInterpolation() = default; - - /*! - * \brief virtual method for setting the cofficients for the respective spline spline. - * \param[in] X - the x values. - * \param[in] Data - the f(x) values. - */ - virtual void SetSpline(const vector &X, const vector &Data){} - - /*! - * \brief virtual method for evaluating the value of the respective Spline. - * \param[in] Point_Interp - the point where interpolation value is required. - * \returns the interpolated value. - */ - virtual su2double EvaluateSpline(su2double Point_Interp){return 0;} - - /*! - * \brief bool variable to make sure all vertex points fell in range of the inlet data. - * \returns the bool variable. - */ - bool GetPointMatch() const {return Point_Match;} + /*! + * \brief Virtual destructor of the C1DInterpolation class. + */ + virtual ~C1DInterpolation() = default; + + /*! + * \brief virtual method for setting the cofficients for the respective spline spline. + * \param[in] X - the x values. + * \param[in] Data - the f(x) values. + */ + virtual void SetSpline(const vector &X, const vector &Data){} + + /*! + * \brief virtual method for evaluating the value of the respective Spline. + * \param[in] Point_Interp - the point where interpolation value is required. + * \returns the interpolated value. + */ + virtual su2double EvaluateSpline(su2double Point_Interp){return 0;} + + /*! + * \brief bool variable to make sure all vertex points fell in range of the inlet data. + * \returns the bool variable. + */ + bool GetPointMatch() const {return Point_Match;} }; -class CAkimaInterpolation final: public C1DInterpolation{ +class CAkimaInterpolation final: public C1DInterpolation{ private: - vector x,y,b,c,d; /*!< \brief local variables for Akima spline cooefficients */ - int n; /*!< \brief local variable for holding the size of the vector */ + vector x,y,b,c,d; /*!< \brief local variables for Akima spline cooefficients */ + int n; /*!< \brief local variable for holding the size of the vector */ + public: - - /*! - * \brief Constructor of the CAkimaInterpolation class. - * \param[in] X - the x values. - * \param[in] Data - the f(x) values. - */ - CAkimaInterpolation(vector &X, vector &Data){ - SetSpline(X,Data); - } - - /*! - * \brief Destructor of the CAkimaInterpolation class. - */ - ~CAkimaInterpolation(){} - - /*! - * \brief for setting the cofficients for the Akima spline. - * \param[in] X - the x values. - * \param[in] Data - the f(x) values. - */ - void SetSpline(const vector &X, const vector &Data) override; - - /*! - * \brief For evaluating the value of Akima Spline. - * \param[in] Point_Interp - the point where interpolation value is required. - * \returns the interpolated value. - */ - su2double EvaluateSpline(su2double Point_Interp) override; + /*! + * \brief Constructor of the CAkimaInterpolation class. + * \param[in] X - the x values. + * \param[in] Data - the f(x) values. + */ + CAkimaInterpolation(vector &X, vector &Data){ + SetSpline(X,Data); + } + + /*! + * \brief Destructor of the CAkimaInterpolation class. + */ + ~CAkimaInterpolation(){} + + /*! + * \brief for setting the cofficients for the Akima spline. + * \param[in] X - the x values. + * \param[in] Data - the f(x) values. + */ + void SetSpline(const vector &X, const vector &Data) override; + + /*! + * \brief For evaluating the value of Akima Spline. + * \param[in] Point_Interp - the point where interpolation value is required. + * \returns the interpolated value. + */ + su2double EvaluateSpline(su2double Point_Interp) override; }; class CLinearInterpolation final: public C1DInterpolation{ - private: - vector x,y,dydx; /*!< \brief local variables for linear 'spline' cooefficients */ - int n; /*!< \brief local variable for holding the size of the vector */ - public: - - /*! - * \brief Constructor of the CLinearInterpolation class. - * \param[in] X - the x values. - * \param[in] Data - the f(x) values. - */ - CLinearInterpolation(vector &X, vector &Data){ - SetSpline(X,Data); - } - - /*! - * \brief Destructor of the CInletInterpolation class. - */ - ~CLinearInterpolation(){} - - /*! - * \brief for setting the cofficients for Linear 'spline'. - * \param[in] X - the x values. - * \param[in] Data - the f(x) values. - */ - void SetSpline(const vector &X, const vector &Data) override; - - /*! - * \brief For evaluating the value for Linear 'spline'. - * \param[in] Point_Interp - the point where interpolation value is required. - * \returns the interpolated value. - */ - su2double EvaluateSpline(su2double Point_Interp) override; +private: + vector x,y,dydx; /*!< \brief local variables for linear 'spline' cooefficients */ + int n; /*!< \brief local variable for holding the size of the vector */ + +public: + /*! + * \brief Constructor of the CLinearInterpolation class. + * \param[in] X - the x values. + * \param[in] Data - the f(x) values. + */ + CLinearInterpolation(vector &X, vector &Data){ + SetSpline(X,Data); + } + + /*! + * \brief Destructor of the CInletInterpolation class. + */ + ~CLinearInterpolation(){} + + /*! + * \brief for setting the cofficients for Linear 'spline'. + * \param[in] X - the x values. + * \param[in] Data - the f(x) values. + */ + void SetSpline(const vector &X, const vector &Data) override; + + /*! + * \brief For evaluating the value for Linear 'spline'. + * \param[in] Point_Interp - the point where interpolation value is required. + * \returns the interpolated value. + */ + su2double EvaluateSpline(su2double Point_Interp) override; }; /*! -* \brief to correct for interpolation type. -* \param[in] Inlet_Interpolated - the interpolated data after spline evaluation. -* \param[in] Theta - the angle of the vertex (in xy plane). -* \param[in] nDim - the dimensions of the case. -* \param[in] Coord - the coordinates of the vertex. -* \param[in] nVar_Turb - the number of turbulence variables as defined by turbulence model -* \param[in] ENUM_INLET_INTERPOLATIONTYPE - enum of the interpolation type to be done -* \returns the corrected Inlet Interpolated Data. -*/ -vector CorrectedInletValues(const vector &Inlet_Interpolated, - su2double Theta , - unsigned short nDim, - su2double *Coord, - unsigned short nVar_Turb, - ENUM_INLET_INTERPOLATIONTYPE Interpolation_Type); + * \brief to correct for interpolation type. + * \param[in] Inlet_Interpolated - the interpolated data after spline evaluation. + * \param[in] Theta - the angle of the vertex (in xy plane). + * \param[in] nDim - the dimensions of the case. + * \param[in] Coord - the coordinates of the vertex. + * \param[in] nVar_Turb - the number of turbulence variables as defined by turbulence model + * \param[in] ENUM_INLET_INTERPOLATIONTYPE - enum of the interpolation type to be done + * \returns the corrected Inlet Interpolated Data. + */ +vector CorrectedInletValues(const vector &Inlet_Interpolated, + su2double Theta , + unsigned short nDim, + const su2double *Coord, + unsigned short nVar_Turb, + ENUM_INLET_INTERPOLATIONTYPE Interpolation_Type); /*! -* \brief to print the Inlet Interpolated Data -* \param[in] Inlet_Interpolated_Interpolated - the final vector for the interpolated data -* \param[in] Marker - name of the inlet marker -* \param[in] nVertex - total number of vertexes. -* \param[in] nDim - the dimensions of the problem. -* \param[in] nColumns - the number of columns in the final interpolated data -*/ -void PrintInletInterpolatedData(const vector& Inlet_Data_Interpolated, string Marker, unsigned long nVertex, unsigned short nDim, unsigned short nColumns); + * \brief to print the Inlet Interpolated Data + * \param[in] Inlet_Interpolated_Interpolated - the final vector for the interpolated data + * \param[in] Marker - name of the inlet marker + * \param[in] nVertex - total number of vertexes. + * \param[in] nDim - the dimensions of the problem. + * \param[in] nColumns - the number of columns in the final interpolated data + */ +void PrintInletInterpolatedData(const vector& Inlet_Data_Interpolated, string Marker, + unsigned long nVertex, unsigned short nDim, unsigned short nColumns); diff --git a/Common/src/toolboxes/C1DInterpolation.cpp b/Common/src/toolboxes/C1DInterpolation.cpp index 52b63408ac0..c7d2c31bbd1 100644 --- a/Common/src/toolboxes/C1DInterpolation.cpp +++ b/Common/src/toolboxes/C1DInterpolation.cpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2019, SU2 Contributors (cf. AUTHORS.md) @@ -24,167 +24,177 @@ * You should have received a copy of the GNU Lesser General Public * License along with SU2. If not, see . */ - + #include "../../include/toolboxes/C1DInterpolation.hpp" su2double CAkimaInterpolation::EvaluateSpline(su2double Point_Interp){ - Point_Match = false; - for (int i=0; i=x[i] && Point_Interp<=x[i+1]) - Point_Match = true; + Point_Match = false; - if (Point_Match == false){ - cout<<"WARNING: Extrapolating data for radius: "<=x[i] && Point_Interp<=x[i+1]) + Point_Match = true; - int i=0, j=n-1; + if (Point_Match == false){ + cout<<"WARNING: Extrapolating data for radius: "<1){ - int m=(i+j) / 2 ; - if (Point_Interp>x[m]) {i=m;} - else {j=m;} - } + int i=0, j=n-1; + + while (j-i>1){ + int m = (i+j) / 2; + if (Point_Interp>x[m]) {i = m;} + else {j = m;} + } - su2double h=Point_Interp-x[i] ; - return y[i]+h*(b[i]+h*(c[i]+h*d[i])) ; + su2double h=Point_Interp-x[i]; + + return y[i]+h*(b[i]+h*(c[i]+h*d[i])); } su2double CLinearInterpolation::EvaluateSpline(su2double Point_Interp){ - Point_Match = false; - for (int i=0; i=x[i] && Point_Interp<=x[i+1]){ - Point_Match = true; - return (Point_Interp-x[i])*dydx[i]+y[i]; - } + Point_Match = false; + + for (int i=0; i=x[i] && Point_Interp<=x[i+1]){ + Point_Match = true; + return (Point_Interp-x[i])*dydx[i]+y[i]; } - return 0; + } + return 0; } void CLinearInterpolation::SetSpline(const vector &X, const vector &Data){ - n = X.size(); - su2double h; - x.resize(n); - y.resize(n); - dydx.resize(n-1); - - x=X; - y=Data; - - for (int i=0; i &X,const vector &Data){ - - n = X.size(); - vector h (n-1); - vector p (n-1); - - /*---calculating finite differences (h) and gradients (p) ---*/ - for (int i=0; i &X,const vector &Data){ + + n = X.size(); + vector h (n-1); + vector p (n-1); + + /*---calculating finite differences (h) and gradients (p) ---*/ + for (int i=0; i CorrectedInletValues(const vector &Inlet_Interpolated , - su2double Theta , - unsigned short nDim, - su2double *Coord, - unsigned short nVar_Turb, - ENUM_INLET_INTERPOLATIONTYPE Interpolation_Type){ +vector CorrectedInletValues(const vector &Inlet_Interpolated , + su2double Theta , + unsigned short nDim, + const su2double *Coord, + unsigned short nVar_Turb, + ENUM_INLET_INTERPOLATIONTYPE Interpolation_Type){ - unsigned short size_columns=Inlet_Interpolated.size()+nDim; - vector Inlet_Values(size_columns); - su2double unit_r, unit_Theta, unit_m, Alpha, Phi; + unsigned short size_columns=Inlet_Interpolated.size()+nDim; + vector Inlet_Values(size_columns); + su2double unit_r, unit_Theta, unit_m, Alpha, Phi; - /*---For x,y,z,T,P columns---*/ - for (int i=0;i& Inlet_Data_Interpolated, string Marker, + unsigned long nVertex, unsigned short nDim, unsigned short nColumns){ - return Inlet_Values; -} + ofstream myfile; + myfile.precision(16); + myfile.open("Interpolated_Data_"+Marker+".dat",ios_base::out); -void PrintInletInterpolatedData(const vector& Inlet_Data_Interpolated, string Marker, unsigned long nVertex, unsigned short nDim, unsigned short nColumns){ - ofstream myfile; - myfile.precision(16); - myfile.open("Interpolated_Data_"+Marker+".dat",ios_base::out); - - if (myfile.is_open()){ - for (unsigned long iVertex = 0; iVertex < nVertex; iVertex++) { - for (unsigned short iVar=0; iVar < nColumns; iVar++){ - myfile<GetTime_Marching() == DT_STEPPING_1ST) || (config->GetTime_Marching() == DT_STEPPING_2ND)); bool time_stepping = config->GetTime_Marching() == TIME_STEPPING; @@ -4177,9 +4178,9 @@ void CSolver::LoadInletProfile(CGeometry **geometry, string profile_filename = config->GetInlet_FileName(); ifstream inlet_file; string Interpolation_Function, Interpolation_Type; - bool Interpolate; + bool Interpolate = false; - su2double *Normal = new su2double[nDim]; + su2double *Normal = new su2double[nDim]; unsigned long Marker_Counter = 0; @@ -4236,207 +4237,218 @@ void CSolver::LoadInletProfile(CGeometry **geometry, const su2double tolerance = config->GetInlet_Profile_Matching_Tolerance(); for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) { - if (config->GetMarker_All_KindBC(iMarker) == KIND_MARKER) { - /*--- Get tag in order to identify the correct inlet data. ---*/ + /*--- Skip if this is the wrong type of marker. ---*/ - Marker_Tag = config->GetMarker_All_TagBound(iMarker); + if (config->GetMarker_All_KindBC(iMarker) != KIND_MARKER) continue; - for (jMarker = 0; jMarker < profileReader.GetNumberOfProfiles(); jMarker++) { + /*--- Get tag in order to identify the correct inlet data. ---*/ - /*--- If we have found the matching marker string, continue. ---*/ + Marker_Tag = config->GetMarker_All_TagBound(iMarker); - if (profileReader.GetTagForProfile(jMarker) == Marker_Tag) { + for (jMarker = 0; jMarker < profileReader.GetNumberOfProfiles(); jMarker++) { - /*--- Increment our counter for marker matches. ---*/ + /*--- If we have not found the matching marker string, continue to next marker. ---*/ - Marker_Counter++; + if (profileReader.GetTagForProfile(jMarker) != Marker_Tag) continue; - /*--- Get data for this profile. ---*/ + /*--- Increment our counter for marker matches. ---*/ - vector Inlet_Data = profileReader.GetDataForProfile(jMarker); - unsigned short nColumns = profileReader.GetNumberOfColumnsInProfile(jMarker); - vector Inlet_Data_Interpolated ((nCol_InletFile+nDim)*geometry[MESH_0]->nVertex[iMarker]); + Marker_Counter++; - /*--- Define Inlet Values vectors before and after interpolation (if needed) ---*/ - vector Inlet_Values(nCol_InletFile+nDim); - vector Inlet_Interpolated(nColumns); + /*--- Get data for this profile. ---*/ - unsigned long nRows = profileReader.GetNumberOfRowsInProfile(jMarker); + vector Inlet_Data = profileReader.GetDataForProfile(jMarker); + unsigned short nColumns = profileReader.GetNumberOfColumnsInProfile(jMarker); + vector Inlet_Data_Interpolated ((nCol_InletFile+nDim)*geometry[MESH_0]->nVertex[iMarker]); - /*--- Pointer to call Set and Evaluate functions. ---*/ - vector interpolator (nColumns); - string interpolation_function, interpolation_type; + /*--- Define Inlet Values vectors before and after interpolation (if needed) ---*/ + vector Inlet_Values(nCol_InletFile+nDim); + vector Inlet_Interpolated(nColumns); - /*--- Define the reference for interpolation. ---*/ - unsigned short radius_index=0; - vector InletRadii = profileReader.GetColumnForProfile(jMarker, radius_index); - vector Interpolation_Column (nRows); + unsigned long nRows = profileReader.GetNumberOfRowsInProfile(jMarker); - switch(config->GetKindInletInterpolationFunction()){ + /*--- Pointer to call Set and Evaluate functions. ---*/ + vector interpolator (nColumns); + string interpolation_function, interpolation_type; - case (NONE): - Interpolate = false; - break; + /*--- Define the reference for interpolation. ---*/ + unsigned short radius_index=0; + vector InletRadii = profileReader.GetColumnForProfile(jMarker, radius_index); + vector Interpolation_Column (nRows); - case (AKIMA_1D): - for (unsigned short iCol=0; iCol < nColumns; iCol++){ - Interpolation_Column = profileReader.GetColumnForProfile(jMarker, iCol); - interpolator[iCol] = new CAkimaInterpolation(InletRadii,Interpolation_Column); - interpolation_function = "AKIMA"; - Interpolate = true; - } - break; - - case (LINEAR_1D): - for (unsigned short iCol=0; iCol < nColumns; iCol++){ - Interpolation_Column = profileReader.GetColumnForProfile(jMarker, iCol); - interpolator[iCol] = new CLinearInterpolation(InletRadii,Interpolation_Column); - interpolation_function = "LINEAR"; - Interpolate = true; - } - break; - - default: - SU2_MPI::Error("Error in the Kind_InletInterpolation Marker\n",CURRENT_FUNCTION); - break; - } + switch(config->GetKindInletInterpolationFunction()){ - if (Interpolate == true){ - switch(config->GetKindInletInterpolationType()){ - case(VR_VTHETA): - interpolation_type="VR_VTHETA"; - break; - case(ALPHA_PHI): - interpolation_type="ALPHA_PHI"; - break; - } - cout<<"Inlet Interpolation being done using "<nVertex[iMarker]; iVertex++) { - - iPoint = geometry[MESH_0]->vertex[iMarker][iVertex]->GetNode(); - Coord = geometry[MESH_0]->node[iPoint]->GetCoord(); + default: + SU2_MPI::Error("Error in the Kind_InletInterpolation Marker\n",CURRENT_FUNCTION); + break; + } - if(Interpolate == false){ + if (Interpolate == true){ + switch(config->GetKindInletInterpolationType()){ + case(VR_VTHETA): + interpolation_type="VR_VTHETA"; + break; + case(ALPHA_PHI): + interpolation_type="ALPHA_PHI"; + break; + } + cout<<"Inlet Interpolation being done using "<nVertex[iMarker]; iVertex++) { - /*--- Get the coords for this data point. ---*/ + iPoint = geometry[MESH_0]->vertex[iMarker][iVertex]->GetNode(); + Coord = geometry[MESH_0]->node[iPoint]->GetCoord(); - index = iRow*nColumns; + if(Interpolate == false) { - dist = 0.0; - for (unsigned short iDim = 0; iDim < nDim; iDim++) - dist += pow(Inlet_Data[index+iDim] - Coord[iDim], 2); - dist = sqrt(dist); + min_dist = 1e16; - /*--- Check is this is the closest point and store data if so. ---*/ + /*--- Find the distance to the closest point in our inlet profile data. ---*/ - if (dist < min_dist) { - min_dist = dist; - for (iVar = 0; iVar < nColumns; iVar++) - Inlet_Values[iVar] = Inlet_Data[index+iVar]; - } + for (iRow = 0; iRow < nRows; iRow++) { - } + /*--- Get the coords for this data point. ---*/ - /*--- If the diff is less than the tolerance, match the two. - We could modify this to simply use the nearest neighbor, or - eventually add something more elaborate here for interpolation. ---*/ + index = iRow*nColumns; - if (min_dist < tolerance) { + dist = 0.0; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + dist += pow(Inlet_Data[index+iDim] - Coord[iDim], 2); + dist = sqrt(dist); - solver[MESH_0][KIND_SOLVER]->SetInletAtVertex(Inlet_Values.data(), iMarker, iVertex); + /*--- Check is this is the closest point and store data if so. ---*/ - } else { + if (dist < min_dist) { + min_dist = dist; + for (iVar = 0; iVar < nColumns; iVar++) + Inlet_Values[iVar] = Inlet_Data[index+iVar]; + } - unsigned long GlobalIndex = geometry[MESH_0]->node[iPoint]->GetGlobalIndex(); - cout << "WARNING: Did not find a match between the points in the inlet file" << endl; - cout << "and point " << GlobalIndex; - cout << std::scientific; - cout << " at location: [" << Coord[0] << ", " << Coord[1]; - if (nDim ==3) error_msg << ", " << Coord[2]; - cout << "]" << endl; - cout << "Distance to closest point: " << min_dist << endl; - cout << "Current tolerance: " << tolerance << endl; - cout << endl; - cout << "You can widen the tolerance for point matching by changing the value" << endl; - cout << "of the option INLET_MATCHING_TOLERANCE in your *.cfg file." << endl; - local_failure++; - break; - } + } - } + /*--- If the diff is less than the tolerance, match the two. + We could modify this to simply use the nearest neighbor, or + eventually add something more elaborate here for interpolation. ---*/ + + if (min_dist < tolerance) { + + solver[MESH_0][KIND_SOLVER]->SetInletAtVertex(Inlet_Values.data(), iMarker, iVertex); + + } else { + + unsigned long GlobalIndex = geometry[MESH_0]->node[iPoint]->GetGlobalIndex(); + cout << "WARNING: Did not find a match between the points in the inlet file" << endl; + cout << "and point " << GlobalIndex; + cout << std::scientific; + cout << " at location: [" << Coord[0] << ", " << Coord[1]; + if (nDim ==3) error_msg << ", " << Coord[2]; + cout << "]" << endl; + cout << "Distance to closest point: " << min_dist << endl; + cout << "Current tolerance: " << tolerance << endl; + cout << endl; + cout << "You can widen the tolerance for point matching by changing the value" << endl; + cout << "of the option INLET_MATCHING_TOLERANCE in your *.cfg file." << endl; + local_failure++; + break; + } - else if(Interpolate == true){ - - /* --- Calculating the radius and angle of the vertex ---*/ - /* --- Flow should be in z direction for 3D cases ---*/ - /* --- Or in x direction for 2D cases ---*/ - Interp_Radius = sqrt(pow(Coord[0],2)+ pow(Coord[1],2)); - Theta = atan2(Coord[1],Coord[0]); - - /* --- Evaluating and saving the final spline data ---*/ - for (unsigned short iVar=0; iVar < nColumns; iVar++){ - - /*---Evaluate spline will get the respective value of the Data set (column) specified - for that interpolator[iVar], cycling through all columns to get all the - data for that vertex ---*/ - Inlet_Interpolated[iVar]=interpolator[iVar]->EvaluateSpline(Interp_Radius); - if (interpolator[iVar]->GetPointMatch() == false){ - cout << "WARNING: Did not find a match between the radius in the inlet file " ; - cout << std::scientific; - cout << "at location: [" << Coord[0] << ", " << Coord[1]; - if (nDim == 3) {cout << ", " << Coord[2];} - cout << "]"; - cout << " with Radius: "<< Interp_Radius << endl; - cout << "You can add a row for Radius: " << Interp_Radius <<" in the inlet file "; - cout << "to eliminate this issue or give proper data" << endl; - local_failure++; - break; - } - } + } - /* --- Correcting for Interpolation Type ---*/ - switch(config->GetKindInletInterpolationType()){ - case(VR_VTHETA): - Inlet_Values = CorrectedInletValues(Inlet_Interpolated, Theta, nDim, Coord, nVar_Turb, VR_VTHETA); - break; - case(ALPHA_PHI): - Inlet_Values = CorrectedInletValues(Inlet_Interpolated, Theta, nDim, Coord, nVar_Turb, ALPHA_PHI); + else if(Interpolate == true) { + + /* --- Calculating the radius and angle of the vertex ---*/ + /* --- Flow should be in z direction for 3D cases ---*/ + /* --- Or in x direction for 2D cases ---*/ + Interp_Radius = sqrt(pow(Coord[0],2)+ pow(Coord[1],2)); + Theta = atan2(Coord[1],Coord[0]); + + /* --- Evaluating and saving the final spline data ---*/ + for (unsigned short iVar=0; iVar < nColumns; iVar++){ + + /*---Evaluate spline will get the respective value of the Data set (column) specified + for that interpolator[iVar], cycling through all columns to get all the + data for that vertex ---*/ + Inlet_Interpolated[iVar]=interpolator[iVar]->EvaluateSpline(Interp_Radius); + if (interpolator[iVar]->GetPointMatch() == false){ + cout << "WARNING: Did not find a match between the radius in the inlet file " ; + cout << std::scientific; + cout << "at location: [" << Coord[0] << ", " << Coord[1]; + if (nDim == 3) {cout << ", " << Coord[2];} + cout << "]"; + cout << " with Radius: "<< Interp_Radius << endl; + cout << "You can add a row for Radius: " << Interp_Radius <<" in the inlet file "; + cout << "to eliminate this issue or give proper data" << endl; + local_failure++; break; - } - - solver[MESH_0][KIND_SOLVER]->SetInletAtVertex(Inlet_Values.data(), iMarker, iVertex); - - for (unsigned short iVar=0; iVar < (nCol_InletFile+nDim); iVar++) - Inlet_Data_Interpolated[iVertex*(nCol_InletFile+nDim)+iVar] = Inlet_Values[iVar]; } } - if(config->GetPrintInlet_InterpolatedData() == true) - PrintInletInterpolatedData(Inlet_Data_Interpolated,profileReader.GetTagForProfile(jMarker),geometry[MESH_0]->nVertex[iMarker],nDim, nCol_InletFile+nDim); - - for (int i=0; iGetKindInletInterpolationType()){ + case(VR_VTHETA): + Inlet_Values = CorrectedInletValues(Inlet_Interpolated, Theta, nDim, Coord, nVar_Turb, VR_VTHETA); + break; + case(ALPHA_PHI): + Inlet_Values = CorrectedInletValues(Inlet_Interpolated, Theta, nDim, Coord, nVar_Turb, ALPHA_PHI); + break; + } + + solver[MESH_0][KIND_SOLVER]->SetInletAtVertex(Inlet_Values.data(), iMarker, iVertex); + + for (unsigned short iVar=0; iVar < (nCol_InletFile+nDim); iVar++) + Inlet_Data_Interpolated[iVertex*(nCol_InletFile+nDim)+iVar] = Inlet_Values[iVar]; + } + + } // end iVertex loop + + if(config->GetPrintInlet_InterpolatedData() == true) { + PrintInletInterpolatedData(Inlet_Data_Interpolated, profileReader.GetTagForProfile(jMarker), + geometry[MESH_0]->nVertex[iMarker], nDim, nCol_InletFile+nDim); } - if (local_failure > 0) break; - } - } + + for (int i=0; i 0) break; + + } // end iMarker loop SU2_MPI::Allreduce(&local_failure, &global_failure, 1, MPI_UNSIGNED_SHORT, MPI_SUM, MPI_COMM_WORLD); @@ -4459,8 +4471,6 @@ void CSolver::LoadInletProfile(CGeometry **geometry, for (jMarker = 0; jMarker < profileReader.GetNumberOfProfiles(); jMarker++) { if (profileReader.GetTagForProfile(jMarker) == Marker_Tag) { nColumns = profileReader.GetNumberOfColumnsInProfile(jMarker); - if(Interpolate == true) - nColumns+=nDim; } } vector Inlet_Values(nColumns);