diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index 50add943c1f..1c79535c5aa 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -48,7 +48,11 @@ class CLookUpTable { version_lut, /*!< \brief LUT version as specified in LUT file.*/ version_reader, /*!< \brief Reader version (should be equal or above LUT version).*/ name_CV1, /*!< \brief Name of controlling variable 1.*/ - name_CV2; /*!< \brief Name of xontrolling variable 2.*/ + name_CV2; /*!< \brief Name of controlling variable 2.*/ + + unsigned long idx_CV1, /*!< \brief Table variable index of controlling variable 1.*/ + idx_CV2, /*!< \brief Table variable index of controlling variable 2.*/ + idx_null; /*!< \brief Variable index corresponding to NULL variable. */ su2vector n_points, /*!< \brief Number of data poins per table level.*/ n_triangles, /*!< \brief Number of triangles per table level.*/ @@ -56,7 +60,7 @@ class CLookUpTable { unsigned long n_variables, n_table_levels = 1; - su2vector z_values_levels; /*!< \brief Constant z-values of each table level.*/ + std::vector z_values_levels; /*!< \brief Constant z-values of each table level.*/ unsigned short table_dim = 2; /*!< \brief Table dimension.*/ /*! @@ -105,21 +109,6 @@ class CLookUpTable { */ su2vector> interp_mat_inv_x_y; - /*! \brief - * Returns the index to the variable in the lookup table. - */ - inline unsigned int GetIndexOfVar(const std::string& nameVar) const { - int index = find(names_var.begin(), names_var.end(), nameVar) - names_var.begin(); - if (index == int(names_var.size())) { - index = -1; - std::string error_msg = "Variable '"; - error_msg.append(nameVar); - error_msg.append("' is not in the lookup table."); - SU2_MPI::Error(error_msg, CURRENT_FUNCTION); - } - return index; - } - /*! \brief * Returns true if the string is null or zero (ignores case). */ @@ -182,8 +171,8 @@ class CLookUpTable { * \param[in] interp_mat_inv - Inverse matrix for interpolation. * \param[out] interp_coeffs - Interpolation coefficients. */ - void GetInterpCoeffs(su2double val_CV1, su2double val_CV2, su2activematrix& interp_mat_inv, - std::array& interp_coeffs); + void GetInterpCoeffs(su2double val_CV1, su2double val_CV2, const su2activematrix& interp_mat_inv, + std::array& interp_coeffs) const; /*! * \brief Compute interpolated value of a point P in the triangle. @@ -239,6 +228,28 @@ class CLookUpTable { const std::vector& names_var, std::vector& var_vals, const unsigned long i_level = 0); + /*! + * \brief Interpolate data based on distance-weighted averaging on the nearest two table nodes. + * \param[in] val_CV1 - First coordinate of point P(val_CV1,val_CV2) to check. + * \param[in] val_CV2 - Second coordinate of point P(val_CV1,val_CV2) to check. + * \param[in] idx_var - Vector of table variable indices to be looked up. + * \param[out] val_vars - Pointer to the vector of stored values of the variables to look up. + */ + void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, + const std::vector& idx_var, std::vector& var_vals, + const unsigned long i_level = 0); + + /*! + * \brief Interpolate data based on distance-weighted averaging on the nearest two table nodes. + * \param[in] val_CV1 - First coordinate of point P(val_CV1,val_CV2) to check. + * \param[in] val_CV2 - Second coordinate of point P(val_CV1,val_CV2) to check. + * \param[in] idx_var - Vector of table variable indices to be looked up. + * \param[out] val_vars - Pointer to the vector of stored values of the variables to look up. + */ + void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, + const std::vector& idx_var, std::vector& var_vals, + const unsigned long i_level = 0); + /*! * \brief Interpolate data based on distance-weighted averaging on the nearest two table nodes for a single variable. * \param[in] val_CV1 - First coordinate of point P(val_CV1,val_CV2) to check. @@ -284,6 +295,27 @@ class CLookUpTable { const su2double val_CV1, const su2double val_CV2, const su2double val_CV3); + /*! + * \brief Find the triangle within which the query point (val_CV1, val_CV2) is located. + * \param[in] val_CV1 - First controlling variable value. + * \param[in] val_CV2 - Second controlling variable value. + * \param[in] id_triangle - Reference to inclusion triangle index. + * \param[in] iLevel - Table level index. + * \returns if query point is within data set. + */ + bool FindInclusionTriangle(const su2double val_CV1, const su2double val_CV2, unsigned long& id_triangle, + const unsigned long iLevel = 0); + + /*! + * \brief Identify the nearest second nearest hull nodes w.r.t. the query point (val_CV1, val_CV2). + * \param[in] val_CV1 - First controlling variable value. + * \param[in] val_CV2 - Second controlling variable value. + * \param[in] iLevel - Table level index. + * \returns pair with nearest node index(first) and second nearest node(second). + */ + std::pair FindNearestNeighbors(const su2double val_CV1, const su2double val_CV2, + const unsigned long iLevel = 0); + public: CLookUpTable(const std::string& file_name_lut, std::string name_CV1_in, std::string name_CV2_in); @@ -293,62 +325,66 @@ class CLookUpTable { void PrintTableInfo(); /*! - * \brief Lookup 1 value of the single variable "val_name_var" using controlling variable values(val_CV1,val_CV2). - * \param[in] val_name_var - String name of the variable to look up. + * \brief Lookup value of variable stored under idx_var using controlling variable values(val_CV1,val_CV2). + * \param[in] idx_var - Column index corresponding to look-up data. * \param[out] val_var - The stored value of the variable to look up. * \param[in] val_CV1 - Value of controlling variable 1. * \param[in] val_CV2 - Value of controlling variable 2. - * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. - */ - unsigned long LookUp_XY(const std::string& val_name_var, su2double* val_var, su2double val_CV1, su2double val_CV2, - unsigned long i_level = 0); - - /*! - * \brief Lookup 1 value for each of the variables in "val_name_var" using controlling variable - * values(val_CV1,val_CV2). - * \param[in] val_names_var - vector of string names of the variables to look up. - * \param[out] val_vars - pointer to the vector of stored values of the variables to look up. - * \param[in] val_CV1 - value of controlling variable 1. - * \param[in] val_CV2 - value of controlling variable 2. - * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + * \returns whether query is inside (true) or outside (false) data set. */ - unsigned long LookUp_XY(const std::vector& val_names_var, std::vector& val_vars, - su2double val_CV1, su2double val_CV2, unsigned long i_level = 0); + bool LookUp_XY(const unsigned long idx_var, su2double* val_var, const su2double val_CV1, const su2double val_CV2, + const unsigned long i_level = 0); /*! - * \brief Lookup the value of the variable "val_name_var" using controlling variable values(val_CV1,val_CV2). - * \param[in] val_name_var - String name of the variable to look up. + * \brief Lookup the value of the variable stored under idx_var using controlling variable values(val_CV1,val_CV2). + * \param[in] idx_var - Table data column indices corresponding to look-up variables. * \param[out] val_var - The stored value of the variable to look up. * \param[in] val_CV1 - Value of controlling variable 1. * \param[in] val_CV2 - Value of controlling variable 2. - * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + * \returns whether query is inside (true) or outside (false) data set. */ - unsigned long LookUp_XY(const std::vector& val_names_var, std::vector& val_vars, - su2double val_CV1, su2double val_CV2, unsigned long i_level = 0); + bool LookUp_XY(const std::vector& idx_var, std::vector& val_vars, const su2double val_CV1, + const su2double val_CV2, const unsigned long i_level = 0); /*! - * \brief Lookup the value of the variable "val_name_var" using controlling variable values(val_CV1,val_CV2,val_z). - * \param[in] val_name_var - String name of the variable to look up. + * \brief Lookup the values of the variables stored under idx_var using controlling variable values(val_CV1,val_CV2). + * \param[in] idx_var - Table data column indices corresponding to look-up variables. * \param[out] val_var - The stored value of the variable to look up. * \param[in] val_CV1 - Value of controlling variable 1. * \param[in] val_CV2 - Value of controlling variable 2. - * \param[in] val_CV3 - Value of controlling variable 3. - * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + * \returns whether query is inside (true) or outside (false) data set. + */ + bool LookUp_XY(const std::vector& idx_var, std::vector& val_vars, const su2double val_CV1, + su2double val_CV2, const unsigned long i_level = 0); + /*! + * \brief Lookup the value of the variable stored under idx_var using controlling variable values(val_CV1,val_CV2, + * val_CV3). \param[in] val_name_var - String name of the variable to look up. \param[out] val_var - The stored value + * of the variable to look up. \param[in] val_CV1 - Value of controlling variable 1. \param[in] val_CV2 - Value of + * controlling variable 2. \param[in] val_CV3 - Value of controlling variable 3. \returns whether query is inside + * (true) or outside (false) data set. */ - unsigned long LookUp_XYZ(const std::string& val_name_var, su2double* val_var, su2double val_CV1, su2double val_CV2, - su2double val_CV3); + bool LookUp_XYZ(const unsigned long idx_var, su2double* val_var, const su2double val_CV1, const su2double val_CV2, + const su2double val_CV3); /*! - * \brief Lookup the value of the variable "val_name_var" using controlling variable values(val_CV1,val_CV2,val_z). - * \param[in] val_name_var - String name of the variable to look up. - * \param[out] val_var - The stored value of the variable to look up. - * \param[in] val_CV1 - Value of controlling variable 1. - * \param[in] val_CV2 - Value of controlling variable 2. - * \param[in] val_CV3 - Value of controlling variable 3. - * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + * \brief Lookup the values of the variables stored under idx_var using controlling variable values + * (val_CV1,val_CV2,val_z). \param[in] idx_var - Table variable index to look up. \param[out] val_var - The stored + * value of the variable to look up. \param[in] val_CV1 - Value of controlling variable 1. \param[in] val_CV2 - Value + * of controlling variable 2. \param[in] val_CV3 - Value of controlling variable 3. \returns whether query is inside + * (true) or outside (false) data set. + */ + bool LookUp_XYZ(const std::vector& idx_var, std::vector& val_vars, const su2double val_CV1, + const su2double val_CV2, const su2double val_CV3 = 0); + + /*! + * \brief Lookup the values of the variables stored under idx_var using controlling variable values + * (val_CV1,val_CV2,val_z). \param[in] idx_var - Table variable index to look up. \param[out] val_var - The stored + * value of the variable to look up. \param[in] val_CV1 - Value of controlling variable 1. \param[in] val_CV2 - Value + * of controlling variable 2. \param[in] val_CV3 - Value of controlling variable 3. \returns whether query is inside + * (true) or outside (false) data set. */ - unsigned long LookUp_XYZ(const std::vector& val_names_var, std::vector& val_vars, - su2double val_CV1, su2double val_CV2, su2double val_CV3 = 0); + bool LookUp_XYZ(const std::vector& idx_var, std::vector& val_vars, const su2double val_CV1, + const su2double val_CV2, const su2double val_CV3 = 0); /*! * \brief Find the table levels with constant z-values directly above and below query val_z. @@ -362,7 +398,7 @@ class CLookUpTable { * \brief Determine the minimum and maximum value of the second controlling variable. * \returns Pair of minimum and maximum value of controlling variable 2. */ - inline std::pair GetTableLimitsY(unsigned long i_level = 0) const { + inline std::pair GetTableLimitsY(const unsigned long i_level = 0) const { return limits_table_y[i_level]; } @@ -370,7 +406,29 @@ class CLookUpTable { * \brief Determine the minimum and maximum value of the first controlling variable. * \returns Pair of minimum and maximum value of controlling variable 1. */ - inline std::pair GetTableLimitsX(unsigned long i_level = 0) const { + inline std::pair GetTableLimitsX(const unsigned long i_level = 0) const { return limits_table_x[i_level]; } + + /*! + * \brief Returns the index to the variable in the lookup table. + * \param[in] nameVar - Variable name for which to retrieve the column index. + * \returns Table data column index corresponding to variable. + */ + inline unsigned int GetIndexOfVar(const std::string& nameVar) const { + int index = find(names_var.begin(), names_var.end(), nameVar) - names_var.begin(); + if (index == int(names_var.size())) { + index = -1; + std::string error_msg = "Variable '"; + error_msg.append(nameVar); + error_msg.append("' is not in the lookup table."); + SU2_MPI::Error(error_msg, CURRENT_FUNCTION); + } + return index; + } + + /*! \brief + * Returns the table variable index which will always return zero when looked up. + */ + unsigned long GetNullIndex() const { return idx_null; } }; diff --git a/Common/include/containers/CTrapezoidalMap.hpp b/Common/include/containers/CTrapezoidalMap.hpp index c82d0281a2c..2225a9b718c 100644 --- a/Common/include/containers/CTrapezoidalMap.hpp +++ b/Common/include/containers/CTrapezoidalMap.hpp @@ -69,7 +69,7 @@ class CTrapezoidalMap { * \param[in] val_y - y-coordinate or second independent variable * \param[out] val_index - index to the triangle */ - unsigned long GetTriangle(su2double val_x, su2double val_y); + unsigned long GetTriangle(const su2double val_x, const su2double val_y); /*! * \brief get the indices of the vertical coordinate band (xmin,xmax) in the 2D search space @@ -78,7 +78,7 @@ class CTrapezoidalMap { * \param[out] val_band - a pair(i_low,i_up) , the lower index and upper index between which the value val_x * can be found */ - std::pair GetBand(su2double val_x); + std::pair GetBand(const su2double val_x); /*! * \brief for a given coordinate (val_x,value), known to be in the band (xmin,xmax) with band index (i_low,i_up), diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 206ea2ff648..6efa79bcef9 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -40,6 +40,10 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, LoadTableRaw(var_file_name_lut); + /* Store indices of controlling variables. */ + idx_CV1 = GetIndexOfVar(name_CV1); + idx_CV2 = GetIndexOfVar(name_CV2); + FindTableLimits(name_CV1, name_CV2); if (rank == MASTER_NODE) @@ -53,6 +57,9 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, PrintTableInfo(); + /* Add additional variable index which will always result in zero when looked up. */ + idx_null = names_var.size(); + if (rank == MASTER_NODE) switch (table_dim) { case 2: cout << "Building a trapezoidal map for the (" + name_CV1 + ", " + name_CV2 + @@ -158,8 +165,6 @@ void CLookUpTable::LoadTableRaw(const string& var_file_name_lut) { } void CLookUpTable::FindTableLimits(const string& name_cv1, const string& name_cv2) { - int idx_CV1 = GetIndexOfVar(name_cv1); - int idx_CV2 = GetIndexOfVar(name_cv2); limits_table_x.resize(n_table_levels); limits_table_y.resize(n_table_levels); @@ -214,7 +219,7 @@ void CLookUpTable::PrintTableInfo() { cout << "| Variable names: |" << endl; cout << "| |" << endl; - for (unsigned long i_var = 0; i_var < names_var.size(); i_var++) + for (auto i_var = 0u; i_var < names_var.size(); i_var++) cout << "| " << right << setw(3) << i_var << ": " << left << setw(59) << names_var[i_var] << " |" << endl; cout << "+------------------------------------------------------------------+" << endl; @@ -238,7 +243,7 @@ void CLookUpTable::PrintTableInfo() { cout << "| Variable names: |" << endl; cout << "| |" << endl; - for (unsigned long i_var = 0; i_var < names_var.size(); i_var++) + for (auto i_var = 0u; i_var < names_var.size(); i_var++) cout << "| " << right << setw(3) << i_var << ": " << left << setw(59) << names_var[i_var] << " |" << endl; cout << "+------------------------------------------------------------------+" << endl; @@ -257,11 +262,11 @@ void CLookUpTable::IdentifyUniqueEdges() { vector > neighborElemsOfPoint; neighborElemsOfPoint.resize(n_points[i_level]); - for (unsigned long iElem = 0; iElem < n_triangles[i_level]; iElem++) { + for (auto iElem = 0u; iElem < n_triangles[i_level]; iElem++) { /* loop over 3 points per triangle */ - for (unsigned long iPoint = 0; iPoint < N_POINTS_TRIANGLE; iPoint++) { + for (auto iPoint = 0u; iPoint < N_POINTS_TRIANGLE; iPoint++) { /* get the global ID of the current point */ - const unsigned long GlobalIndex = triangles[i_level][iElem][iPoint]; + const auto GlobalIndex = triangles[i_level][iElem][iPoint]; /* add the current element ID to the neighbor list for this point */ neighborElemsOfPoint[GlobalIndex].push_back(iElem); @@ -270,7 +275,7 @@ void CLookUpTable::IdentifyUniqueEdges() { /* remove duplicates from the neighboring element lists*/ vector::iterator vecIt; - for (unsigned long iPoint = 0; iPoint < n_points[i_level]; iPoint++) { + for (auto iPoint = 0u; iPoint < n_points[i_level]; iPoint++) { /* sort neighboring elements for each point */ sort(neighborElemsOfPoint[iPoint].begin(), neighborElemsOfPoint[iPoint].end()); @@ -285,12 +290,12 @@ void CLookUpTable::IdentifyUniqueEdges() { the point IDs that are neighboring points */ vector > neighborPointsOfPoint; neighborPointsOfPoint.resize(n_points[i_level]); - for (unsigned long iPoint = 0; iPoint < n_points[i_level]; iPoint++) { - for (unsigned long iElem = 0; iElem < neighborElemsOfPoint[iPoint].size(); iElem++) { + for (auto iPoint = 0u; iPoint < n_points[i_level]; iPoint++) { + for (auto iElem = 0u; iElem < neighborElemsOfPoint[iPoint].size(); iElem++) { /* loop over element points */ - for (unsigned long jPoint = 0; jPoint < N_POINTS_TRIANGLE; jPoint++) { + for (auto jPoint = 0u; jPoint < N_POINTS_TRIANGLE; jPoint++) { /* get the global ID of the current point */ - const unsigned long GlobalIndex = triangles[i_level][neighborElemsOfPoint[iPoint][iElem]][jPoint]; + const auto GlobalIndex = triangles[i_level][neighborElemsOfPoint[iPoint][iElem]][jPoint]; /* add the current element ID to the neighbor list for this point */ if (GlobalIndex != iPoint) neighborPointsOfPoint[iPoint].push_back(GlobalIndex); @@ -299,7 +304,7 @@ void CLookUpTable::IdentifyUniqueEdges() { } /* remove duplicates from the neighboring points lists */ - for (unsigned long iPoint = 0; iPoint < n_points[i_level]; iPoint++) { + for (auto iPoint = 0u; iPoint < n_points[i_level]; iPoint++) { /* sort neighboring points for each point */ sort(neighborPointsOfPoint[iPoint].begin(), neighborPointsOfPoint[iPoint].end()); @@ -315,10 +320,10 @@ void CLookUpTable::IdentifyUniqueEdges() { that the smaller global index is in the first position and the larger is in the second to make for a unique set, so there's no need to remove duplicates. */ - for (unsigned long iPoint = 0; iPoint < n_points[i_level]; iPoint++) { - for (unsigned long jPoint = 0; jPoint < neighborPointsOfPoint[iPoint].size(); jPoint++) { + for (auto iPoint = 0u; iPoint < n_points[i_level]; iPoint++) { + for (auto jPoint = 0u; jPoint < neighborPointsOfPoint[iPoint].size(); jPoint++) { /* Store the neighbor index more clearly. */ - const unsigned long GlobalIndex = neighborPointsOfPoint[iPoint][jPoint]; + const auto GlobalIndex = neighborPointsOfPoint[iPoint][jPoint]; /* Store the edge so that the lower index of the pair is always first. */ if (iPoint < GlobalIndex) { @@ -334,17 +339,17 @@ void CLookUpTable::IdentifyUniqueEdges() { pair, loop through the neighboring elements and store the two elements that contain the second point in the edge ('left' and 'right' of the edge). */ edge_to_triangle[i_level].resize(edges[i_level].size()); - for (unsigned long iEdge = 0; iEdge < edges[i_level].size(); iEdge++) { + for (auto iEdge = 0u; iEdge < edges[i_level].size(); iEdge++) { /* Store the two points of the edge more clearly. */ - const unsigned long iPoint = edges[i_level][iEdge][0]; - const unsigned long jPoint = edges[i_level][iEdge][1]; + const auto iPoint = edges[i_level][iEdge][0]; + const auto jPoint = edges[i_level][iEdge][1]; /* Loop over all neighboring elements to iPoint. */ - for (unsigned long iElem = 0; iElem < neighborElemsOfPoint[iPoint].size(); iElem++) { + for (auto iElem = 0u; iElem < neighborElemsOfPoint[iPoint].size(); iElem++) { /* loop over 3 points per triangle */ - for (unsigned long kPoint = 0; kPoint < N_POINTS_TRIANGLE; kPoint++) { + for (auto kPoint = 0u; kPoint < N_POINTS_TRIANGLE; kPoint++) { /* Get the global ID of the current point. */ - const unsigned long GlobalIndex = triangles[i_level][neighborElemsOfPoint[iPoint][iElem]][kPoint]; + const auto GlobalIndex = triangles[i_level][neighborElemsOfPoint[iPoint][iElem]][kPoint]; /* Add the current element ID to the neighbor list for this point. */ if (GlobalIndex == jPoint) edge_to_triangle[i_level][iEdge].push_back(neighborElemsOfPoint[iPoint][iElem]); @@ -360,18 +365,18 @@ void CLookUpTable::ComputeInterpCoeffs() { std::array next_triangle; - const su2double* val_CV1 = GetDataP(name_CV1, i_level); - const su2double* val_CV2 = GetDataP(name_CV2, i_level); + const su2double* val_CV1 = table_data[i_level][idx_CV1]; + const su2double* val_CV2 = table_data[i_level][idx_CV2]; /* calculate weights for each triangle (basically a distance function) and * build inverse interpolation matrices */ interp_mat_inv_x_y[i_level].resize(n_triangles[i_level]); - for (unsigned long i_triangle = 0; i_triangle < n_triangles[i_level]; i_triangle++) { - for (int p = 0; p < 3; p++) { + for (auto i_triangle = 0u; i_triangle < n_triangles[i_level]; i_triangle++) { + for (auto p = 0u; p < N_POINTS_TRIANGLE; p++) { next_triangle[p] = triangles[i_level][i_triangle][p]; } - su2activematrix x_interp_mat_inv(3, 3); + su2activematrix x_interp_mat_inv(N_POINTS_TRIANGLE, N_POINTS_TRIANGLE); GetInterpMatInv(val_CV1, val_CV2, next_triangle, x_interp_mat_inv); interp_mat_inv_x_y[i_level][i_triangle] = x_interp_mat_inv; } @@ -380,24 +385,23 @@ void CLookUpTable::ComputeInterpCoeffs() { void CLookUpTable::GetInterpMatInv(const su2double* vec_x, const su2double* vec_y, std::array& point_ids, su2activematrix& interp_mat_inv) { - const unsigned int M = 3; - CSquareMatrixCM global_M(3); + CSquareMatrixCM global_M(N_POINTS_TRIANGLE); /* setup LHM matrix for the interpolation */ - for (unsigned int i_point = 0; i_point < M; i_point++) { - su2double x = vec_x[point_ids[i_point]]; - su2double y = vec_y[point_ids[i_point]]; + for (auto i_vertex = 0u; i_vertex < N_POINTS_TRIANGLE; i_vertex++) { + su2double x = vec_x[point_ids[i_vertex]]; + su2double y = vec_y[point_ids[i_vertex]]; - global_M(i_point, 0) = SU2_TYPE::GetValue(1.0); - global_M(i_point, 1) = SU2_TYPE::GetValue(x); - global_M(i_point, 2) = SU2_TYPE::GetValue(y); + global_M(i_vertex, 0) = SU2_TYPE::GetValue(1.0); + global_M(i_vertex, 1) = SU2_TYPE::GetValue(x); + global_M(i_vertex, 2) = SU2_TYPE::GetValue(y); } global_M.Invert(); global_M.Transpose(); - for (unsigned int i = 0; i < M; i++) { - for (unsigned int j = 0; j < M; j++) { + for (auto i = 0u; i < N_POINTS_TRIANGLE; i++) { + for (auto j = 0u; j < N_POINTS_TRIANGLE; j++) { interp_mat_inv[i][j] = global_M(i, j); } } @@ -406,61 +410,47 @@ void CLookUpTable::GetInterpMatInv(const su2double* vec_x, const su2double* vec_ std::pair CLookUpTable::FindInclusionLevels(const su2double val_CV3) { /*--- Find the table levels with constant z-values directly below and above the query value val_CV3 ---*/ + su2double val_z = val_CV3; + unsigned long i_up{0}, i_low{0}; /* Check if val_CV3 lies outside table bounds */ - if (val_CV3 >= *limits_table_z.second) { - return std::make_pair(n_table_levels - 1, n_table_levels - 1); - } else if (val_CV3 <= *limits_table_z.first) { - return std::make_pair(0, 0); - } else { - std::pair inclusion_levels; - /* Loop over table levels to find the levels directly above and below the query value */ - for (auto i_level = 0ul; i_level < n_table_levels - 1; i_level++) { - if ((val_CV3 >= z_values_levels[i_level]) && (val_CV3 < z_values_levels[i_level + 1])) { - inclusion_levels = std::make_pair(i_level, i_level + 1); - } - } - return inclusion_levels; + if (val_z < z_values_levels.front()) { + val_z = z_values_levels.front(); + return make_pair(i_low, i_up); } -} - -unsigned long CLookUpTable::LookUp_XYZ(const std::string& val_name_var, su2double* val_var, su2double val_CV1, - su2double val_CV2, su2double val_CV3) { - /*--- Perform quasi-3D interpolation for a single variable named val_name_var on a query point - with coordinates val_CV1, val_CV2, and val_CV3 ---*/ - - /* 1: Find table levels directly above and below the query point (the levels that sandwhich val_CV3) */ - std::pair inclusion_levels = FindInclusionLevels(val_CV3); - bool within_z_limits = (inclusion_levels.first != inclusion_levels.second); - - if (within_z_limits) { - /* 2: Determine val_CV1 and val_CV2 for the inclusion table levels. */ - std::array, 2> lower_upper_CV1_2 = - ComputeNormalizedXY(inclusion_levels, val_CV1, val_CV2, val_CV3); - su2double val_CV1_lower = lower_upper_CV1_2[0][0], val_CV2_lower = lower_upper_CV1_2[0][1], - val_CV1_upper = lower_upper_CV1_2[1][0], val_CV2_upper = lower_upper_CV1_2[1][1]; + if (val_z > z_values_levels.back()) { + val_z = z_values_levels.back(); + i_low = n_table_levels - 1; + i_up = n_table_levels - 1; + return make_pair(i_low, i_up); + } + std::pair::iterator, std::vector::iterator> bounds; + bounds = std::equal_range(z_values_levels.begin(), z_values_levels.end(), val_z); - /* 3: Perform 2D interpolations on upper and lower inclusion levels */ - unsigned long lower_level = inclusion_levels.first; - unsigned long upper_level = inclusion_levels.second; + /*--- if upper bound = 0, then use the range [0,1] ---*/ + i_up = max(1, bounds.first - z_values_levels.begin()); + i_low = i_up - 1; + return make_pair(i_low, i_up); +} - su2double val_var_lower, val_var_upper; - unsigned long exit_code_lower = LookUp_XY(val_name_var, &val_var_lower, val_CV1_lower, val_CV2_lower, lower_level); - unsigned long exit_code_upper = LookUp_XY(val_name_var, &val_var_upper, val_CV1_upper, val_CV2_upper, upper_level); +bool CLookUpTable::LookUp_XYZ(const unsigned long idx_var, su2double* val_var, const su2double val_CV1, + const su2double val_CV2, const su2double val_CV3) { + vector vec_idx_var = {idx_var}; + vector vec_val_var = {val_var}; + return LookUp_XYZ(vec_idx_var, vec_val_var, val_CV1, val_CV2, val_CV3); +} - /* 4: Perform linear interpolation along the z-direction using the x-y interpolation results - from upper and lower trapezoidal maps */ - Linear_Interpolation(val_CV3, lower_level, upper_level, val_var_lower, val_var_upper, val_var); +bool CLookUpTable::LookUp_XYZ(const std::vector& idx_var, std::vector& val_vars, + const su2double val_CV1, const su2double val_CV2, const su2double val_CV3) { + vector var_vals_output; + var_vals_output.resize(val_vars.size()); + bool inside = LookUp_XYZ(idx_var, var_vals_output, val_CV1, val_CV2, val_CV3); + for (auto iVar = 0u; iVar < val_vars.size(); iVar++) *val_vars[iVar] = var_vals_output[iVar]; - return max(exit_code_lower, exit_code_upper); - } else { - /* Perform single, 2D interpolation when val_CV3 lies outside table bounds */ - unsigned long bound_level = inclusion_levels.first; - LookUp_XY(val_name_var, val_var, val_CV1, val_CV2, bound_level); - return 1; - } + return inside; } -unsigned long CLookUpTable::LookUp_XYZ(const std::vector& val_names_var, std::vector& val_vars, - su2double val_CV1, su2double val_CV2, su2double val_CV3) { + +bool CLookUpTable::LookUp_XYZ(const std::vector& idx_var, std::vector& val_vars, + const su2double val_CV1, const su2double val_CV2, const su2double val_CV3) { /*--- Perform quasi-3D interpolation for a vector of variables with names val_names_var on a query point with coordinates val_CV1, val_CV2, and val_CV3 ---*/ @@ -482,19 +472,19 @@ unsigned long CLookUpTable::LookUp_XYZ(const std::vector& val_names std::vector val_vars_lower, val_vars_upper; val_vars_lower.resize(val_vars.size()); val_vars_upper.resize(val_vars.size()); - unsigned long exit_code_lower = LookUp_XY(val_names_var, val_vars_lower, val_CV1_lower, val_CV2_lower, lower_level); - unsigned long exit_code_upper = LookUp_XY(val_names_var, val_vars_upper, val_CV1_upper, val_CV2_upper, upper_level); + auto inside_lower = LookUp_XY(idx_var, val_vars_lower, val_CV1_lower, val_CV2_lower, lower_level); + auto inside_upper = LookUp_XY(idx_var, val_vars_upper, val_CV1_upper, val_CV2_upper, upper_level); /* 4: Perform linear interpolation along the z-direction using the x-y interpolation results from upper and lower trapezoidal maps */ Linear_Interpolation(val_CV3, lower_level, upper_level, val_vars_lower, val_vars_upper, val_vars); - return max(exit_code_lower, exit_code_upper); + return (inside_upper && inside_lower); } else { /* Perform single, 2D interpolation when val_CV3 lies outside table bounds */ unsigned long bound_level = inclusion_levels.first; - LookUp_XY(val_names_var, val_vars, val_CV1, val_CV2, bound_level); - return 1; + LookUp_XY(idx_var, val_vars, val_CV1, val_CV2, bound_level); + return false; } } @@ -579,105 +569,76 @@ std::array, 2> CLookUpTable::ComputeNormalizedXY( return lower_upper_CVs; } -unsigned long CLookUpTable::LookUp_XY(const string& val_name_var, su2double* val_var, su2double val_CV1, - su2double val_CV2, unsigned long i_level) { - unsigned long exit_code = 1; - - if (noSource(val_name_var)) { - *val_var = 0.0; - exit_code = 0; - return exit_code; - } - - /* check if x and y value is in table range */ - if ((val_CV1 >= *limits_table_x[i_level].first && val_CV1 <= *limits_table_x[i_level].second) && - (val_CV2 >= *limits_table_y[i_level].first && val_CV2 <= *limits_table_y[i_level].second)) { - /* find the triangle that holds the (x, y) point */ - unsigned long id_triangle = trap_map_x_y[i_level].GetTriangle(val_CV1, val_CV2); +bool CLookUpTable::LookUp_XY(const vector& idx_var, vector& val_vars, const su2double val_CV1, + const su2double val_CV2, unsigned long i_level) { + unsigned long id_triangle; + bool inside = FindInclusionTriangle(val_CV1, val_CV2, id_triangle, i_level); - if (IsInTriangle(val_CV1, val_CV2, id_triangle, i_level)) { - /* get interpolation coefficients for point on triangle */ - std::array interp_coeffs{0}; - GetInterpCoeffs(val_CV1, val_CV2, interp_mat_inv_x_y[i_level][id_triangle], interp_coeffs); + /* loop over variable names and interpolate / get values */ + if (inside) { + std::array interp_coeffs{0}; + std::array triangle{0}; + for (auto iVertex = 0u; iVertex < N_POINTS_TRIANGLE; iVertex++) + triangle[iVertex] = triangles[i_level][id_triangle][iVertex]; - /* first, copy the single triangle from the large triangle list*/ - std::array triangle{0}; - for (int p = 0; p < 3; p++) triangle[p] = triangles[i_level][id_triangle][p]; + GetInterpCoeffs(val_CV1, val_CV2, interp_mat_inv_x_y[i_level][id_triangle], interp_coeffs); - *val_var = Interpolate(GetDataP(val_name_var, i_level), triangle, interp_coeffs); - exit_code = 0; + for (auto iVar = 0u; iVar < idx_var.size(); iVar++) { + if (idx_var[iVar] == idx_null) { + val_vars[iVar] = 0; + } else { + su2double interp_var = Interpolate(table_data[i_level][idx_var[iVar]], triangle, interp_coeffs); + val_vars[iVar] = interp_var; + } } - } - if (exit_code == 1) InterpolateToNearestNeighbors(val_CV1, val_CV2, val_name_var, val_var, i_level); + } else + InterpolateToNearestNeighbors(val_CV1, val_CV2, idx_var, val_vars, i_level); - return exit_code; + return inside; } -unsigned long CLookUpTable::LookUp_XY(const vector& val_names_var, vector& val_vars, - su2double val_CV1, su2double val_CV2, unsigned long i_level) { - vector look_up_data(val_names_var.size()); +bool CLookUpTable::LookUp_XY(const unsigned long idx_var, su2double* val_var, const su2double val_CV1, + const su2double val_CV2, const unsigned long i_level) { + vector vec_idx_var = {idx_var}; + vector val_vars = {val_var}; + return LookUp_XY(vec_idx_var, val_vars, val_CV1, val_CV2, i_level); +} - for (long unsigned int i_var = 0; i_var < val_vars.size(); ++i_var) { - look_up_data[i_var] = &val_vars[i_var]; - } +bool CLookUpTable::LookUp_XY(const vector& idx_var, vector& val_vars, + const su2double val_CV1, const su2double val_CV2, const unsigned long i_level) { + vector output_var_vals; + output_var_vals.resize(val_vars.size()); + bool inside = LookUp_XY(idx_var, output_var_vals, val_CV1, val_CV2, i_level); - unsigned long exit_code = LookUp_XY(val_names_var, look_up_data, val_CV1, val_CV2, i_level); + for (auto iVar = 0u; iVar < val_vars.size(); iVar++) *val_vars[iVar] = output_var_vals[iVar]; - return exit_code; + return inside; } -unsigned long CLookUpTable::LookUp_XY(const vector& val_names_var, vector& val_vars, - su2double val_CV1, su2double val_CV2, unsigned long i_level) { - unsigned long exit_code = 1; - unsigned long id_triangle = 0; - std::array interp_coeffs{0}; - std::array triangle{0}; - +bool CLookUpTable::FindInclusionTriangle(const su2double val_CV1, const su2double val_CV2, unsigned long& id_triangle, + const unsigned long iLevel) { /* check if x value is in table x-dimension range * and if y is in table y-dimension table range */ - if ((val_CV1 >= *limits_table_x[i_level].first && val_CV1 <= *limits_table_x[i_level].second) && - (val_CV2 >= *limits_table_y[i_level].first && val_CV2 <= *limits_table_y[i_level].second)) { + if ((val_CV1 >= *limits_table_x[iLevel].first && val_CV1 <= *limits_table_x[iLevel].second) && + (val_CV2 >= *limits_table_y[iLevel].first && val_CV2 <= *limits_table_y[iLevel].second)) { /* if so, try to find the triangle that holds the (prog, enth) point */ - id_triangle = trap_map_x_y[i_level].GetTriangle(val_CV1, val_CV2); + id_triangle = trap_map_x_y[iLevel].GetTriangle(val_CV1, val_CV2); /* check if point is inside a triangle (if table domain is non-rectangular, * the previous range check might be true but the point could still be outside of the domain) */ - if (IsInTriangle(val_CV1, val_CV2, id_triangle, i_level)) { - /* if so, get interpolation coefficients for point in the triangle */ - GetInterpCoeffs(val_CV1, val_CV2, interp_mat_inv_x_y[i_level][id_triangle], interp_coeffs); - - /* exit_code 0 means point was in triangle */ - exit_code = 0; - } + return IsInTriangle(val_CV1, val_CV2, id_triangle, iLevel); } - - /* loop over variable names and interpolate / get values */ - for (long unsigned int i_var = 0; i_var < val_names_var.size(); ++i_var) { - if (noSource(val_names_var[i_var])) { - *val_vars[i_var] = 0.0; - exit_code = 0; - } else { - if (exit_code == 0) { - /* first, copy the single triangle from the large triangle list*/ - for (int p = 0; p < 3; p++) triangle[p] = triangles[i_level][id_triangle][p]; - *val_vars[i_var] = Interpolate(GetDataP(val_names_var[i_var], i_level), triangle, interp_coeffs); - } else { - InterpolateToNearestNeighbors(val_CV1, val_CV2, val_names_var[i_var], val_vars[i_var], i_level); - } - } - } - - return exit_code; + return false; } -void CLookUpTable::GetInterpCoeffs(su2double val_CV1, su2double val_CV2, su2activematrix& interp_mat_inv, - std::array& interp_coeffs) { - std::array query_vector = {1, val_CV1, val_CV2}; +void CLookUpTable::GetInterpCoeffs(su2double val_CV1, su2double val_CV2, const su2activematrix& interp_mat_inv, + std::array& interp_coeffs) const { + std::array query_vector = {1, val_CV1, val_CV2}; su2double d; - for (int i = 0; i < 3; i++) { + for (auto i = 0u; i < N_POINTS_TRIANGLE; i++) { d = 0; - for (int j = 0; j < 3; j++) { + for (auto j = 0u; j < N_POINTS_TRIANGLE; j++) { d = d + interp_mat_inv[i][j] * query_vector[j]; } interp_coeffs[i] = d; @@ -689,7 +650,7 @@ su2double CLookUpTable::Interpolate(const su2double* val_samples, std::array CLookUpTable::FindNearestNeighbors(const su2double val_CV1, + const su2double val_CV2, + const unsigned long i_level) { + su2double min_distance = 1e99, second_distance = 1e99; + su2double norm_coeff_x = 1. / (*limits_table_x[i_level].second - *limits_table_x[i_level].first); + su2double norm_coeff_y = 1. / (*limits_table_y[i_level].second - *limits_table_y[i_level].first); + su2double val_CV1_norm = val_CV1 / (*limits_table_x[i_level].second - *limits_table_x[i_level].first); + su2double val_CV2_norm = val_CV2 / (*limits_table_y[i_level].second - *limits_table_y[i_level].first); + + const su2double* x_table = table_data[i_level][idx_CV1]; + const su2double* y_table = table_data[i_level][idx_CV2]; + unsigned long i_nearest = 0, i_second_nearest = 0; + + for (unsigned long i_point = 0; i_point < n_hull_points[i_level]; ++i_point) { + su2double next_x_norm = x_table[hull[i_level][i_point]] * norm_coeff_x; + su2double next_y_norm = y_table[hull[i_level][i_point]] * norm_coeff_y; + + su2double next_distance = pow(val_CV1_norm - next_x_norm, 2) + pow(val_CV2_norm - next_y_norm, 2); + + /* Check if distance to node is lower than current nearest or second nearest node. */ + if (next_distance < min_distance) { + second_distance = min_distance; + min_distance = next_distance; + + i_second_nearest = i_nearest; + i_nearest = hull[i_level][i_point]; + } else if ((next_distance > min_distance) && (next_distance < second_distance)) { + i_second_nearest = hull[i_level][i_point]; + + second_distance = next_distance; + } + } + return make_pair(i_nearest, i_second_nearest); +} + unsigned long CLookUpTable::FindNearestNeighborOnHull(su2double val_CV1, su2double val_CV2, unsigned long i_level) { su2double min_distance = 1.e99; su2double next_distance = 1.e99; - su2double next_x_norm; - su2double next_y_norm; unsigned long neighbor_id = 0; - const su2double* x_table = GetDataP(name_CV1, i_level); - const su2double* y_table = GetDataP(name_CV2, i_level); + const su2double* x_table = table_data[i_level][idx_CV1]; + const su2double* y_table = table_data[i_level][idx_CV2]; su2double norm_coeff_x = 1. / (limits_table_x[i_level].second - limits_table_x[i_level].first); su2double norm_coeff_y = 1. / (limits_table_y[i_level].second - limits_table_y[i_level].first); su2double val_CV1_norm = val_CV1 / (limits_table_x[i_level].second - limits_table_x[i_level].first); su2double val_CV2_norm = val_CV2 / (limits_table_y[i_level].second - limits_table_y[i_level].first); - for (unsigned long i_point = 0; i_point < n_hull_points[i_level]; ++i_point) { - next_x_norm = x_table[hull[i_level][i_point]] * norm_coeff_x; - next_y_norm = y_table[hull[i_level][i_point]] * norm_coeff_y; + for (auto i_point = 0u; i_point < n_hull_points[i_level]; ++i_point) { + su2double next_x_norm = x_table[hull[i_level][i_point]] * norm_coeff_x; + su2double next_y_norm = y_table[hull[i_level][i_point]] * norm_coeff_y; next_distance = sqrt(pow(val_CV1_norm - next_x_norm, 2) + pow(val_CV2_norm - next_y_norm, 2)); @@ -727,39 +721,25 @@ unsigned long CLookUpTable::FindNearestNeighborOnHull(su2double val_CV1, su2doub } void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, - const std::vector& names_var, - std::vector& var_vals, const unsigned long i_level) { + const std::vector& idx_var, + std::vector& var_vals, const unsigned long i_level) { /* Interpolate data using distance-weighted averaging on the two nearest table nodes. */ - su2double min_distance = 1e99, second_distance = 1e99; su2double norm_coeff_x = 1. / (*limits_table_x[i_level].second - *limits_table_x[i_level].first); su2double norm_coeff_y = 1. / (*limits_table_y[i_level].second - *limits_table_y[i_level].first); su2double val_CV1_norm = val_CV1 / (*limits_table_x[i_level].second - *limits_table_x[i_level].first); su2double val_CV2_norm = val_CV2 / (*limits_table_y[i_level].second - *limits_table_y[i_level].first); - const su2double* x_table = GetDataP(name_CV1, i_level); - const su2double* y_table = GetDataP(name_CV2, i_level); - unsigned long i_nearest = 0, i_second_nearest = 0; - - for (unsigned long i_point = 0; i_point < n_hull_points[i_level]; ++i_point) { - su2double next_x_norm = x_table[hull[i_level][i_point]] * norm_coeff_x; - su2double next_y_norm = y_table[hull[i_level][i_point]] * norm_coeff_y; - - su2double next_distance = pow(val_CV1_norm - next_x_norm, 2) + pow(val_CV2_norm - next_y_norm, 2); - - /* Check if distance to node is lower than current nearest or second nearest node. */ - if (next_distance < min_distance) { - second_distance = min_distance; - min_distance = next_distance; - - i_second_nearest = i_nearest; - i_nearest = hull[i_level][i_point]; - } else if ((next_distance > min_distance) && (next_distance < second_distance)) { - i_second_nearest = hull[i_level][i_point]; - - second_distance = next_distance; - } - } + std::pair nearest_IDs = FindNearestNeighbors(val_CV1, val_CV2, i_level); + unsigned long i_nearest = nearest_IDs.first, i_second_nearest = nearest_IDs.second; + su2double x_nearest = table_data[i_level][idx_CV1][nearest_IDs.first], + x_second_nearest = table_data[i_level][idx_CV1][nearest_IDs.second], + y_nearest = table_data[i_level][idx_CV2][nearest_IDs.first], + y_second_nearest = table_data[i_level][idx_CV2][nearest_IDs.second]; + su2double min_distance = + pow(val_CV1_norm - x_nearest * norm_coeff_x, 2) + pow(val_CV2_norm - y_nearest * norm_coeff_y, 2); + su2double second_distance = + pow(val_CV1_norm - x_second_nearest * norm_coeff_x, 2) + pow(val_CV2_norm - y_second_nearest * norm_coeff_y, 2); min_distance = sqrt(min_distance); second_distance = sqrt(second_distance); @@ -767,12 +747,34 @@ void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const /* Interpolate data using distance-weighted averaging */ su2double delimiter = (1.0 / min_distance) + (1.0 / second_distance); for (auto iVar = 0u; iVar < var_vals.size(); iVar++) { - su2double data_nearest = GetDataP(names_var[iVar], i_level)[i_nearest], - data_second_nearest = GetDataP(names_var[iVar], i_level)[i_second_nearest]; - *var_vals[iVar] = (data_nearest * (1.0 / min_distance) + data_second_nearest * (1.0 / second_distance)) / delimiter; + if (idx_var[iVar] == idx_null) { + var_vals[iVar] = 0; + } else { + su2double data_nearest = table_data[i_level][idx_var[iVar]][i_nearest], + data_second_nearest = table_data[i_level][idx_var[iVar]][i_second_nearest]; + var_vals[iVar] = + (data_nearest * (1.0 / min_distance) + data_second_nearest * (1.0 / second_distance)) / delimiter; + } } } +void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, + const std::vector& names_var, + std::vector& var_vals, const unsigned long i_level) { + vector idx_var; + idx_var.resize(names_var.size()); + for (auto iVar = 0u; iVar < names_var.size(); iVar++) idx_var[iVar] = GetIndexOfVar(names_var[iVar]); +} + +void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, + const std::vector& idx_var, + std::vector& var_vals, const unsigned long i_level) { + vector output_var_vals; + output_var_vals.resize(var_vals.size()); + InterpolateToNearestNeighbors(val_CV1, val_CV2, idx_var, output_var_vals, i_level); + for (auto iVar = 0u; iVar < var_vals.size(); iVar++) *var_vals[iVar] = output_var_vals[iVar]; +} + void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, const std::string& name_var, su2double* var_vals, const unsigned long i_level) { @@ -783,14 +785,14 @@ void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const bool CLookUpTable::IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, unsigned long i_level) { - su2double tri_x_0 = GetDataP(name_CV1, i_level)[triangles[i_level][val_id_triangle][0]]; - su2double tri_y_0 = GetDataP(name_CV2, i_level)[triangles[i_level][val_id_triangle][0]]; + su2double tri_x_0 = table_data[i_level][idx_CV1][triangles[i_level][val_id_triangle][0]]; + su2double tri_y_0 = table_data[i_level][idx_CV2][triangles[i_level][val_id_triangle][0]]; - su2double tri_x_1 = GetDataP(name_CV1, i_level)[triangles[i_level][val_id_triangle][1]]; - su2double tri_y_1 = GetDataP(name_CV2, i_level)[triangles[i_level][val_id_triangle][1]]; + su2double tri_x_1 = table_data[i_level][idx_CV1][triangles[i_level][val_id_triangle][1]]; + su2double tri_y_1 = table_data[i_level][idx_CV2][triangles[i_level][val_id_triangle][1]]; - su2double tri_x_2 = GetDataP(name_CV1, i_level)[triangles[i_level][val_id_triangle][2]]; - su2double tri_y_2 = GetDataP(name_CV2, i_level)[triangles[i_level][val_id_triangle][2]]; + su2double tri_x_2 = table_data[i_level][idx_CV1][triangles[i_level][val_id_triangle][2]]; + su2double tri_y_2 = table_data[i_level][idx_CV2][triangles[i_level][val_id_triangle][2]]; su2double area_tri = TriArea(tri_x_0, tri_y_0, tri_x_1, tri_y_1, tri_x_2, tri_y_2); diff --git a/Common/src/containers/CTrapezoidalMap.cpp b/Common/src/containers/CTrapezoidalMap.cpp index b0c1d972b7a..e5494437362 100644 --- a/Common/src/containers/CTrapezoidalMap.cpp +++ b/Common/src/containers/CTrapezoidalMap.cpp @@ -178,7 +178,7 @@ CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* sa } } -unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) { +unsigned long CTrapezoidalMap::GetTriangle(const su2double val_x, const su2double val_y) { /* find x band in which val_x sits */ pair band = GetBand(val_x); @@ -210,16 +210,16 @@ unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) { return triangle[0]; } -pair CTrapezoidalMap::GetBand(su2double val_x) { +pair CTrapezoidalMap::GetBand(const su2double val_x) { unsigned long i_low = 0; unsigned long i_up = 0; - + su2double val_x_sample = val_x; /* check if val_x is in x-bounds of the table, if not then project val_x to either x-min or x-max */ - if (val_x < unique_bands_x.front()) val_x = unique_bands_x.front(); - if (val_x > unique_bands_x.back()) val_x = unique_bands_x.back(); + if (val_x_sample < unique_bands_x.front()) val_x_sample = unique_bands_x.front(); + if (val_x_sample > unique_bands_x.back()) val_x_sample = unique_bands_x.back(); std::pair::iterator, std::vector::iterator> bounds; - bounds = std::equal_range(unique_bands_x.begin(), unique_bands_x.end(), val_x); + bounds = std::equal_range(unique_bands_x.begin(), unique_bands_x.end(), val_x_sample); /*--- if upper bound = 0, then use the range [0,1] ---*/ i_up = max(1, bounds.first - unique_bands_x.begin()); diff --git a/SU2_CFD/include/fluid/CDataDrivenFluid.hpp b/SU2_CFD/include/fluid/CDataDrivenFluid.hpp index 57ffca7cc88..1a524326c8a 100644 --- a/SU2_CFD/include/fluid/CDataDrivenFluid.hpp +++ b/SU2_CFD/include/fluid/CDataDrivenFluid.hpp @@ -95,7 +95,14 @@ class CDataDrivenFluid final : public CFluidModel { vector MLP_inputs; /*!< \brief Inputs for the multi-layer perceptron look-up operation. */ CLookUpTable* lookup_table; /*!< \brief Look-up table regression object. */ - + unsigned long LUT_idx_s, + LUT_idx_dsde_rho, + LUT_idx_dsdrho_e, + LUT_idx_d2sde2, + LUT_idx_d2sdedrho, + LUT_idx_d2sdrho2; + vector LUT_lookup_indices; + unsigned long outside_dataset, /*!< \brief Density-energy combination lies outside data set. */ nIter_Newton; /*!< \brief Number of Newton solver iterations. */ diff --git a/SU2_CFD/include/fluid/CFluidFlamelet.hpp b/SU2_CFD/include/fluid/CFluidFlamelet.hpp index b066b2d0353..12f370ad5eb 100644 --- a/SU2_CFD/include/fluid/CFluidFlamelet.hpp +++ b/SU2_CFD/include/fluid/CFluidFlamelet.hpp @@ -60,6 +60,10 @@ class CFluidFlamelet final : public CFluidModel { CLookUpTable* look_up_table; + vector LUT_idx_TD, + LUT_idx_Sources, + LUT_idx_LookUp; + /*--- Class variables for the multi-layer perceptron method ---*/ #ifdef USE_MLPCPP MLPToolbox::CLookUp_ANN* lookup_mlp; /*!< \brief Multi-layer perceptron collection. */ @@ -79,6 +83,17 @@ class CFluidFlamelet final : public CFluidModel { void PreprocessLookUp(CConfig* config); + /*! \brief + * Returns true if the string is null or zero (ignores case). + */ + inline bool noSource(const std::string& name_var) const { + if (name_var.compare("NULL") == 0 || name_var.compare("Null") == 0 || name_var.compare("null") == 0 || + name_var.compare("ZERO") == 0 || name_var.compare("Zero") == 0 || name_var.compare("zero") == 0) { + return true; + } else { + return false; + } + } public: CFluidFlamelet(CConfig* config, su2double value_pressure_operating); diff --git a/SU2_CFD/src/fluid/CDataDrivenFluid.cpp b/SU2_CFD/src/fluid/CDataDrivenFluid.cpp index fd11176eeba..09f187f2d1d 100644 --- a/SU2_CFD/src/fluid/CDataDrivenFluid.cpp +++ b/SU2_CFD/src/fluid/CDataDrivenFluid.cpp @@ -118,6 +118,21 @@ void CDataDrivenFluid::MapInputs_to_Outputs() { lookup_mlp->PairVariableswithMLPs(*iomap_rhoe); MLP_inputs.resize(2); #endif + } else { + /*--- Retrieve column indices of LUT output variables ---*/ + LUT_idx_s = lookup_table->GetIndexOfVar(output_names_rhoe[idx_s]); + LUT_idx_dsdrho_e = lookup_table->GetIndexOfVar(output_names_rhoe[idx_dsdrho_e]); + LUT_idx_dsde_rho = lookup_table->GetIndexOfVar(output_names_rhoe[idx_dsde_rho]); + LUT_idx_d2sde2 = lookup_table->GetIndexOfVar(output_names_rhoe[idx_d2sde2]); + LUT_idx_d2sdedrho= lookup_table->GetIndexOfVar(output_names_rhoe[idx_d2sdedrho]); + LUT_idx_d2sdrho2 = lookup_table->GetIndexOfVar(output_names_rhoe[idx_d2sdrho2]); + + LUT_lookup_indices.push_back(LUT_idx_s); + LUT_lookup_indices.push_back(LUT_idx_dsde_rho); + LUT_lookup_indices.push_back(LUT_idx_dsdrho_e); + LUT_lookup_indices.push_back(LUT_idx_d2sde2); + LUT_lookup_indices.push_back(LUT_idx_d2sdedrho); + LUT_lookup_indices.push_back(LUT_idx_d2sdrho2); } } @@ -233,21 +248,10 @@ unsigned long CDataDrivenFluid::Predict_MLP(su2double rho, su2double e) { } unsigned long CDataDrivenFluid::Predict_LUT(su2double rho, su2double e) { - unsigned long exit_code; - std::vector output_names_rhoe_LUT; - std::vector outputs_LUT; - output_names_rhoe_LUT.resize(output_names_rhoe.size()); - for (auto iOutput = 0u; iOutput < output_names_rhoe.size(); iOutput++) { - output_names_rhoe_LUT[iOutput] = output_names_rhoe[iOutput]; - } - - outputs_LUT.resize(outputs_rhoe.size()); - for (auto iOutput = 0u; iOutput < outputs_rhoe.size(); iOutput++) { - outputs_LUT[iOutput] = outputs_rhoe[iOutput]; - } - - exit_code = lookup_table->LookUp_XY(output_names_rhoe_LUT, outputs_LUT, rho, e); - return exit_code; + bool inside = lookup_table->LookUp_XY(LUT_lookup_indices, outputs_rhoe, rho, e); + if (inside) + return 0; + return 1; } void CDataDrivenFluid::Evaluate_Dataset(su2double rho, su2double e) { diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 569c4ae6a2a..b4991e96c21 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -196,6 +196,22 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { lookup_mlp->PairVariableswithMLPs(*iomap_Sources); lookup_mlp->PairVariableswithMLPs(*iomap_LookUp); #endif + } else { + for (auto iVar=0u; iVar < varnames_TD.size(); iVar++) { + LUT_idx_TD.push_back(look_up_table->GetIndexOfVar(varnames_TD[iVar])); + } + for (auto iVar=0u; iVar < varnames_Sources.size(); iVar++) { + unsigned long LUT_idx; + if (noSource(varnames_Sources[iVar])) { + LUT_idx = look_up_table->GetNullIndex(); + } else { + LUT_idx = look_up_table->GetIndexOfVar(varnames_Sources[iVar]); + } + LUT_idx_Sources.push_back(LUT_idx); + } + for (auto iVar=0u; iVar < varnames_LookUp.size(); iVar++) { + LUT_idx_LookUp.push_back(look_up_table->GetIndexOfVar(varnames_LookUp[iVar])); + } } } @@ -207,21 +223,25 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca vector varnames; vector val_vars; vector refs_vars; + vector LUT_idx; switch (lookup_type) { case FLAMELET_LOOKUP_OPS::TD: varnames = varnames_TD; + LUT_idx = LUT_idx_TD; #ifdef USE_MLPCPP iomap_Current = iomap_TD; #endif break; case FLAMELET_LOOKUP_OPS::SOURCES: varnames = varnames_Sources; + LUT_idx = LUT_idx_Sources; #ifdef USE_MLPCPP iomap_Current = iomap_Sources; #endif break; case FLAMELET_LOOKUP_OPS::LOOKUP: varnames = varnames_LookUp; + LUT_idx = LUT_idx_LookUp; #ifdef USE_MLPCPP iomap_Current = iomap_LookUp; #endif @@ -233,13 +253,16 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); /*--- Add all quantities and their names to the look up vectors. ---*/ + bool inside; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: if (include_mixture_fraction) { - extrapolation = look_up_table->LookUp_XYZ(varnames, output_refs, val_prog, val_enth, val_mixfrac); + inside = look_up_table->LookUp_XYZ(LUT_idx, output_refs, val_prog, val_enth, val_mixfrac); } else { - extrapolation = look_up_table->LookUp_XY(varnames, output_refs, val_prog, val_enth); + inside = look_up_table->LookUp_XY(LUT_idx, output_refs, val_prog, val_enth); } + if (inside) extrapolation = 0; + else extrapolation = 1; break; case ENUM_DATADRIVEN_METHOD::MLP: refs_vars.resize(output_refs.size()); diff --git a/UnitTests/Common/containers/CLookupTable_tests.cpp b/UnitTests/Common/containers/CLookupTable_tests.cpp index 487ecdc481b..7192d529398 100644 --- a/UnitTests/Common/containers/CLookupTable_tests.cpp +++ b/UnitTests/Common/containers/CLookupTable_tests.cpp @@ -50,8 +50,9 @@ TEST_CASE("LUTreader", "[tabulated chemistry]") { su2double prog = 0.55; su2double enth = -0.5; string look_up_tag = "Density"; + unsigned long idx_tag = look_up_table.GetIndexOfVar(look_up_tag); su2double look_up_dat; - look_up_table.LookUp_XY(look_up_tag, &look_up_dat, prog, enth); + look_up_table.LookUp_XY(idx_tag, &look_up_dat, prog, enth); CHECK(look_up_dat == Approx(1.02)); /*--- look up a single value for viscosity ---*/ @@ -59,7 +60,8 @@ TEST_CASE("LUTreader", "[tabulated chemistry]") { prog = 0.6; enth = 0.9; look_up_tag = "Viscosity"; - look_up_table.LookUp_XY(look_up_tag, &look_up_dat, prog, enth); + idx_tag = look_up_table.GetIndexOfVar(look_up_tag); + look_up_table.LookUp_XY(idx_tag, &look_up_dat, prog, enth); CHECK(look_up_dat == Approx(0.0000674286)); /* find the table limits */ @@ -77,7 +79,8 @@ TEST_CASE("LUTreader", "[tabulated chemistry]") { prog = 1.10; enth = 1.1; look_up_tag = "Density"; - look_up_table.LookUp_XY(look_up_tag, &look_up_dat, prog, enth); + idx_tag = look_up_table.GetIndexOfVar(look_up_tag); + look_up_table.LookUp_XY(idx_tag, &look_up_dat, prog, enth); CHECK(look_up_dat == Approx(1.1738796125)); } @@ -98,8 +101,9 @@ TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { su2double enth = -0.5; su2double mfrac = 0.5; string look_up_tag = "Density"; + unsigned long idx_tag = look_up_table.GetIndexOfVar(look_up_tag); su2double look_up_dat; - look_up_table.LookUp_XYZ(look_up_tag, &look_up_dat, prog, enth, mfrac); + look_up_table.LookUp_XYZ(idx_tag, &look_up_dat, prog, enth, mfrac); CHECK(look_up_dat == Approx(1.02)); /*--- look up a single value for viscosity ---*/ @@ -108,7 +112,8 @@ TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { enth = 0.9; mfrac = 0.8; look_up_tag = "Viscosity"; - look_up_table.LookUp_XYZ(look_up_tag, &look_up_dat, prog, enth, mfrac); + idx_tag = look_up_table.GetIndexOfVar(look_up_tag); + look_up_table.LookUp_XYZ(idx_tag, &look_up_dat, prog, enth, mfrac); CHECK(look_up_dat == Approx(0.0000674286)); /* find the table limits */ @@ -127,6 +132,7 @@ TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { enth = 1.1; mfrac = 2.0; look_up_tag = "Density"; - look_up_table.LookUp_XYZ(look_up_tag, &look_up_dat, prog, enth, mfrac); + idx_tag = look_up_table.GetIndexOfVar(look_up_tag); + look_up_table.LookUp_XYZ(idx_tag, &look_up_dat, prog, enth, mfrac); CHECK(look_up_dat == Approx(1.1738796125)); }