diff --git a/.github/workflows/regression.yml b/.github/workflows/regression.yml index 8f7c35c7389..98cf7ed4bd8 100644 --- a/.github/workflows/regression.yml +++ b/.github/workflows/regression.yml @@ -129,7 +129,7 @@ jobs: with: # -t -c args: -b ${{github.ref}} -t develop -c develop -s ${{matrix.testscript}} - - name: Cleanup + - name: Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:220614-1237 with: entrypoint: /bin/rm @@ -192,12 +192,31 @@ jobs: done find $BIN_FOLDER -type f -exec chmod a+x '{}' \; ls -lahR $BIN_FOLDER + echo "cloning branch" + branch="${{github.ref}}" + name="SU2" + SRC_FOLDER="$PWD/src" + mkdir -p $SRC_FOLDER + cd $SRC_FOLDER + git clone --recursive --depth=1 --shallow-submodules https://github.com/su2code/SU2 $name + cd $name + git config --add remote.origin.fetch '+refs/pull/*/merge:refs/remotes/origin/refs/pull/*/merge' + git config --add remote.origin.fetch '+refs/heads/*:refs/remotes/origin/refs/heads/*' + git fetch origin --depth=1 $branch:$branch + git checkout $branch + git submodule update + echo $PWD + ls -lahR + cd .. + echo "done cloning" + echo $PWD + ls -lahR - name: Run Unit Tests uses: docker://ghcr.io/su2code/su2/test-su2:220614-1237 with: entrypoint: install/bin/${{matrix.testdriver}} - - name: Post Cleanup + - name: Post Cleanup uses: docker://ghcr.io/su2code/su2/test-su2:220614-1237 with: entrypoint: /bin/rm - args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} \ No newline at end of file + args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }} diff --git a/Common/include/containers/CFileReaderLUT.hpp b/Common/include/containers/CFileReaderLUT.hpp new file mode 100644 index 00000000000..a80f4261655 --- /dev/null +++ b/Common/include/containers/CFileReaderLUT.hpp @@ -0,0 +1,83 @@ +/*! + * \file CFileReaderLUT.hpp + * \brief reading lookup table for tabulated fluid properties + * \author D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ +#pragma once + +#include +#include +#include + +#include "../../Common/include/parallelization/mpi_structure.hpp" +#include "../../../Common/include/linear_algebra/blas_structure.hpp" +#include "../../../Common/include/toolboxes/CSquareMatrixCM.hpp" + +class CFileReaderLUT { + protected: + int rank; + + std::string version_lut; + std::string version_reader; + unsigned long n_points; + unsigned long n_triangles; + unsigned long n_hull_points; + unsigned long n_variables; + + /*! \brief Holds the variable names stored in the table file. + * Order is in sync with tableFlamelet. + */ + std::vector names_var; + + /*! \brief Holds all data stored in the table. + * First index addresses the variable while second index addresses the point. + */ + su2activematrix table_data; + + su2matrix triangles; + + std::vector hull; + + std::string SkipToFlag(std::ifstream* file_stream, const std::string& flag); + + public: + CFileReaderLUT(){}; + + inline const std::string& GetVersionLUT() const { return version_lut; } + inline const std::string& GetVersionReader() const { return version_reader; } + inline unsigned long GetNPoints() const { return n_points; } + inline unsigned long GetNTriangles() const { return n_triangles; } + inline unsigned long GetNHullPoints() const { return n_hull_points; } + inline unsigned long GetNVariables() const { return n_variables; } + + inline const std::vector& GetNamesVar() const { return names_var; } + + inline const su2activematrix& GetTableData() const { return table_data; } + + inline const su2matrix& GetTriangles() const { return triangles; }; + + inline const std::vector& GetHull() const { return hull; }; + + void ReadRawDRG(const std::string& file_name); +}; diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp new file mode 100644 index 00000000000..8bb7872108e --- /dev/null +++ b/Common/include/containers/CLookUpTable.hpp @@ -0,0 +1,267 @@ +/*! + * \file CLookupTable.hpp + * \brief tabulation of fluid properties + * \author D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include +#include + +#include "../../Common/include/option_structure.hpp" +#include "CFileReaderLUT.hpp" +#include "CTrapezoidalMap.hpp" + +class CLookUpTable { + protected: + int rank; /*!< \brief MPI Rank. */ + + std::string file_name_lut; + std::string version_lut; + std::string version_reader; + unsigned long n_points; + unsigned long n_triangles; + unsigned long n_variables; + unsigned long n_hull_points; + + /*! + * \brief the lower and upper limits of the enthalpy and progress variable. + */ + su2double limits_table_enth[2]; + su2double limits_table_prog[2]; + + /*! \brief Holds the variable names stored in the table file. + * Order is in sync with data + */ + std::vector names_var; + + /*! \brief + * Holds all data stored in the table. First index addresses the variable + * while second index addresses the point. + */ + su2activematrix table_data; + + su2matrix triangles; + + /* we do not know this size in advance until we go through the entire lookup table */ + std::vector > edges; + std::vector > edge_to_triangle; + + /*! \brief + * The hull contains the boundary of the lookup table. + */ + std::vector hull; + + CTrapezoidalMap trap_map_prog_enth; + + /*! \brief + * vector of all the weight factors for the interpolation. + */ + std::vector interp_mat_inv_prog_enth; + + /*! \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 Get the pointer to the column data of the table (density, temperature, source terms, ...). + * \returns pointer to the column data. + */ + inline const su2double* GetDataP(const std::string& name_var) const { + return table_data[GetIndexOfVar(name_var)]; + } + + /*! + * \brief find the table limits, i.e. the minimum and maximum values of the 2 independent. + * controlling variables (progress variable and enthalpy). We put the values in the variables. + * limits_table_prog[2] and limit_table_enth[2]. + * \param[in] name_prog - the string name for the first controlling variable. + * \param[in] name_enth - the string name of the second controlling variable. + */ + void FindTableLimits(const std::string& name_prog, const std::string& name_enth); + + /*! + * \brief construct a list of all the edges and a list of the pair of elements left and right of the edge. + */ + void IdentifyUniqueEdges(); + + /*! + * \brief read the lookup table from file and store the data. + * \param[in] file_name_lut - the filename of the lookup table. + */ + void LoadTableRaw(const std::string& file_name_lut); + + /*! + * \brief compute vector of all (inverse) interpolation coefficients "interp_mat_inv_prog_enth" of all triangles. + * \param[in] name_prog - the string name of the first controlling variable (progress variable). + * \param[in] name_enth - the string name of the second controlling variable (enthalpy). + */ + void ComputeInterpCoeffs(const std::string& name_prog, const std::string& name_enth); + + /*! + * \brief compute the inverse matrix for interpolation. + * \param[in] vec_x - pointer to first coordinate (progress variable). + * \param[in] vec_y - pointer to second coordinate (enthalpy). + * \param[in] point_ids - single triangle data. + * \param[out] interp_mat_inv - inverse matrix for interpolation. + */ + void GetInterpMatInv(const su2double* vec_x, const su2double* vec_y, std::array& point_ids, + su2activematrix& interp_mat_inv); + + /*! + * \brief compute the interpolation coefficients for the triangular interpolation. + * \param[in] val_x - value of first coordinate (progress variable). + * \param[in] val_y - value of second coordinate (enthalpy). + * \param[in] interp_mat_inv - inverse matrix for interpolation. + * \param[out] interp_coeffs - interpolation coefficients. + */ + void GetInterpCoeffs(su2double val_x, su2double val_y, su2activematrix& interp_mat_inv, + std::array& interp_coeffs); + + /*! + * \brief compute interpolated value of a point P in the triangle. + * \param[in] val_samples - pointer to the variable data. + * \param[in] val_triangle - ID to the triangle. + * \param[in] val_interp_coeffs - interpolation coefficients using the point data in P. + * \returns resulting value of the interpolation. + */ + su2double Interpolate(const su2double* val_samples, std::array& val_triangle, + std::array& val_interp_coeffs); + + /*! + * \brief find the point on the hull (boundary of the table) that is closest to the point P(val_prog,val_enth). + * \param[in] val_x - first coordinate of point P(val_x,val_y) to check. + * \param[in] val_y - second coordinate of point P(val_x,val_y) to check. + * \param[in] name_prog - string name of the first controlling variable. + * \param[in] name_enth - string name of the second controlling variable. + * \returns point id of the nearest neighbor on the hull. + */ + unsigned long FindNearestNeighborOnHull(su2double val_prog, su2double val_enth, const std::string& name_prog, const std::string& name_enth); + + /*! + * \brief determine if a point P(val_x,val_y) is inside the triangle val_id_triangle. + * \param[in] val_x - first coordinate of point P(val_x,val_y) to check. + * \param[in] val_y - second coordinate of point P(val_x,val_y) to check. + * \param[in] val_id_triangle - ID of the triangle to check. + * \param[in] name_prog - string name of the first controlling variable. + * \param[in] name_enth - string name of the second controlling variable. + * \returns true if the point is in the triangle, false if it is outside. + */ + bool IsInTriangle(su2double val_x, su2double val_y, unsigned long val_id_triangle, const std::string& name_prog, + const std::string& name_enth); + + /*! + * \brief compute the area of a triangle given the 3 points of the triangle. + * \param[in] x1 - the coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). + * \param[in] y1 - the coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). + * \param[in] x2 - the coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). + * \param[in] y2 - the coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). + * \param[in] x3 - the coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). + * \param[in] y3 - the coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). + * \returns the absolute value of the area of the triangle. + */ + inline su2double TriArea(su2double x1, su2double y1, su2double x2, su2double y2, su2double x3, su2double y3) { + return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) * 0.5); + } + + public: + CLookUpTable(const std::string& file_name_lut, const std::string& name_prog, const std::string& name_enth); + + + /*! + * \brief print information to screen. + */ + void PrintTableInfo(); + + /*! + * \brief lookup 1 value of the single variable "val_name_var" using controlling variable values(val_prog,val_enth) + * whose controlling variable names are "name_prog" and "name_enth". + * \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_prog - value of controlling variable 1 (progress variable). + * \param[in] val_enth - value of controlling variable 2 (enthalpy). + * \param[in] name_prog - string name of controlling variable 1 (progress variable). + * \param[in] name_enth - string name of controlling variable 2 (enthalpy). + * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + */ + unsigned long LookUp_ProgEnth(const std::string& val_name_var, su2double* val_var, su2double val_prog, su2double val_enth, + const std::string& name_prog, const std::string& name_enth); + + /*! + * \brief lookup 1 value for each of the variables in "val_name_var" using controlling variable values(val_prog,val_enth) + * whose controlling variable names are "name_prog" and "name_enth". + * \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_prog - value of controlling variable 1 (progress variable). + * \param[in] val_enth - value of controlling variable 2 (enthalpy). + * \param[in] name_prog - string name of controlling variable 1 (progress variable). + * \param[in] name_enth - string name of controlling variable 2 (enthalpy). + * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + */ + unsigned long LookUp_ProgEnth(const std::vector& val_names_var, std::vector& val_vars, su2double val_prog, + su2double val_enth, const std::string& name_prog, const std::string& name_enth); + + /*! + * \brief lookup the value of the variable "val_name_var" using controlling variable values(val_prog,val_enth) + * whose controlling variable names are "name_prog" and "name_enth". + * \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_prog - value of controlling variable 1 (progress variable). + * \param[in] val_enth - value of controlling variable 2 (enthalpy). + * \param[in] name_prog - string name of controlling variable 1 (progress variable). + * \param[in] name_enth - string name of controlling variable 2 (enthalpy). + * \returns 1 if the lookup and subsequent interpolation was a success, 0 if not. + */ + unsigned long LookUp_ProgEnth(const std::vector& val_names_var, std::vector& val_vars, su2double val_prog, + su2double val_enth, const std::string& name_prog, const std::string& name_enth); + + /*! + * \brief determine the minimum and maximum value of the enthalpy (controlling variable 2). + * \returns pair of minimum and maximum value of controlling variable 2. + */ + inline std::pair GetTableLimitsEnth() const { + return std::make_pair(limits_table_enth[0], limits_table_enth[1]); + } + + /*! + * \brief determine the minimum and maximum value of the progress variable (controlling variable 1). + * \returns pair of minimum and maximum value of controlling variable 1. + */ + inline std::pair GetTableLimitsProg() const { + return std::make_pair(limits_table_prog[0], limits_table_prog[1]); + } +}; diff --git a/Common/include/containers/CTrapezoidalMap.hpp b/Common/include/containers/CTrapezoidalMap.hpp new file mode 100644 index 00000000000..f3de71d832e --- /dev/null +++ b/Common/include/containers/CTrapezoidalMap.hpp @@ -0,0 +1,106 @@ +/*! + * \file TrapezoidalMap.hpp + * \brief Implementation of the trapezoidal map for tabulation and lookup of fluid properties + * \author D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#pragma once + +#include +#include + +#include "../../Common/include/linear_algebra/blas_structure.hpp" +#include "../../Common/include/toolboxes/CSquareMatrixCM.hpp" + +/*! + * \class CTrapezoidalMap + * \brief Construction of trapezoidal map for tabulated lookup + * \author: D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + */ + class CTrapezoidalMap { + protected: + + /* The unique values of x which exist in the data */ + std::vector unique_bands_x; + + su2activematrix edge_limits_x; + su2activematrix edge_limits_y; + + std::vector > edge_to_triangle; + + /* The value that each edge which intersects the band takes within that + * same band. Used to sort the edges */ + std::vector > > y_edge_at_band_mid; + + public: + + CTrapezoidalMap() = default; + + CTrapezoidalMap(const su2double* samples_x, + const su2double* samples_y, + const unsigned long size, + const std::vector >& edges, + const std::vector >& edge_to_triangle); + + /*! + * \brief return the index to the triangle that contains the coordinates (val_x,val_y) + * \param[in] val_x - x-coordinate or first independent variable + * \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); + + + /*! + * \brief get the indices of the vertical coordinate band (xmin,xmax) in the 2D search space + * that contains the coordinate val_x + * \param[in] val_x - x-coordinate or first independent variable + * \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); + + + /*! + * \brief for a given coordinate (val_x,value), known to be in the band (xmin,xmax) with band index (i_low,i_up), + * find the edges in the band (these edges come from the triangulation) that enclose the coordinate + * \param[in] val_band - pair i_low,i_up + * \param[in] val_x - x-coordinate or first independent variable + * \param[in] val_y - y-coordinate or first independent variable + * \param[out] pair (edge_low,edge_up) - lower edge and upper edge of a triangle that encloses the coordinate + */ + std::pair GetEdges(std::pair val_band, + su2double val_x, + su2double val_y) const; + + /*! + * \brief determine if the x-coordinate falls within the bounds xmin,xmax of the table + * \param[in] val_x - x-coordinate or first independent variable + * \param[out] bool - true if val_x is within (xmin,xmax) + */ + inline bool IsInsideHullX(su2double val_x) { + return (val_x >= unique_bands_x.front()) && (val_x <= unique_bands_x.back()); + } +}; diff --git a/Common/lib/Makefile.am b/Common/lib/Makefile.am index c15191bac45..dd5cba7fb81 100644 --- a/Common/lib/Makefile.am +++ b/Common/lib/Makefile.am @@ -132,7 +132,10 @@ lib_sources = \ ../src/linear_algebra/CSysMatrix.cpp \ ../src/linear_algebra/CSysSolve.cpp \ ../src/linear_algebra/CSysSolve_b.cpp \ - ../src/linear_algebra/CPastixWrapper.cpp + ../src/linear_algebra/CPastixWrapper.cpp \ + ../src/containers/CLookUpTable.cpp \ + ../src/containers/CTrapezoidalMap.cpp \ + ../src/containers/CFileReaderLUT.cpp lib_cxxflags = -fPIC -std=c++11 lib_ldadd = diff --git a/Common/src/containers/CFileReaderLUT.cpp b/Common/src/containers/CFileReaderLUT.cpp new file mode 100644 index 00000000000..0c4b2e79652 --- /dev/null +++ b/Common/src/containers/CFileReaderLUT.cpp @@ -0,0 +1,252 @@ +/*! + * \file CFileReaderLUT.hpp + * \brief reading lookup table for tabulated fluid properties + * \author D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../Common/include/containers/CFileReaderLUT.hpp" +#include "../../Common/include/option_structure.hpp" +#include "../../Common/include/parallelization/mpi_structure.hpp" +#include "../../../Common/include/linear_algebra/blas_structure.hpp" +#include "../../../Common/include/toolboxes/CSquareMatrixCM.hpp" + +using namespace std; + +void CFileReaderLUT::ReadRawDRG(const string& file_name) { + version_reader = "1.0.0"; + + /*--- Store MPI rank. ---*/ + rank = SU2_MPI::GetRank(); + + string line; + string word; + + istringstream stream_names_var; + + ifstream file_stream; + + int ixColon; + + bool eoHeader = false; + bool eoData = false; + bool eoConnectivity = false; + bool eoHull = false; + + file_stream.open(file_name.c_str(), ifstream::in); + + if (!file_stream.is_open()) { + SU2_MPI::Error(string("There is no look-up-table file called ") + file_name, CURRENT_FUNCTION); + } + + /* Read header */ + line = SkipToFlag(&file_stream, "
"); + + while (getline(file_stream, line) && !eoHeader) { + /* number of points in LUT */ + if (line.compare("[version]") == 0) { + getline(file_stream, line); + version_lut = line; + } + + /* number of points in LUT */ + if (line.compare("[number of points]") == 0) { + getline(file_stream, line); + n_points = stoi(line); + } + + /* number of triangles in LUT */ + if (line.compare("[number of triangles]") == 0) { + getline(file_stream, line); + n_triangles = stoi(line); + } + + /* number of points on the hull */ + if (line.compare("[number of hull points]") == 0) { + getline(file_stream, line); + n_hull_points = stoi(line); + } + + /* number of variables in LUT */ + if (line.compare("[number of variables]") == 0) { + getline(file_stream, line); + n_variables = stoi(line); + } + + /* variable names */ + if (line.compare("[variable names]") == 0) { + getline(file_stream, line); + stream_names_var.str(line); + while (stream_names_var) { + stream_names_var >> word; + ixColon = (int)word.find(":"); + + names_var.push_back(word.substr(ixColon + 1, word.size() - 1)); + } + names_var.pop_back(); // removes last redundant element + } + + // check if end of header is reached + if (line.compare("
") == 0) eoHeader = true; + } + + // check version_lut + if (version_lut.compare(version_reader) != 0) + SU2_MPI::Error("Version conflict between Dragon reader and Dragon library file.", CURRENT_FUNCTION); + + // check header quantities + if (n_points == 0 || n_triangles == 0 || n_variables == 0 || n_hull_points == 0) + SU2_MPI::Error( + "Number of points, triangles, hull points, or variables in Dragon " + "library header is zero.", + CURRENT_FUNCTION); + + // check if number of variables is consistent + if (n_variables != names_var.size()) + SU2_MPI::Error( + "Number of read variables does not match number of " + "variables specified in Dragon " + "library header.", + CURRENT_FUNCTION); + + /* now that n_variables, n_points, n_hull_points and n_variables is available, + * allocate memory */ + if (rank == MASTER_NODE) cout << "allocating memory for the data" << endl; + table_data.resize(GetNVariables(), GetNPoints()); + + + if (rank == MASTER_NODE) cout << "allocating memory for the triangles" << endl; + triangles.resize(GetNTriangles(), 3); + + if (rank == MASTER_NODE) cout << "allocating memory for the hull points" << endl; + hull.resize(GetNHullPoints()); + + /* flush any cout */ + if (rank == MASTER_NODE) cout << endl; + + // read data block + if (rank == MASTER_NODE) cout << "loading data block" << endl; + + line = SkipToFlag(&file_stream, ""); + + unsigned long pointCounter = 0; + while (getline(file_stream, line) && !eoData) { + // check if end of data is reached + if (line.compare("") == 0) eoData = true; + + if (!eoData) { + // one line contains values for one point for all variables + istringstream streamDataLine(line); + + // add next line to table array + for (unsigned long iVar = 0; iVar < n_variables; iVar++) { + streamDataLine >> word; + passivedouble tmp = stod(word); + table_data[iVar][pointCounter] = (su2double)tmp; + } + } + pointCounter++; + } + + if (n_points != pointCounter - 1) + SU2_MPI::Error( + "Number of read points does not match number of points " + "specified in Dragon library header.", + CURRENT_FUNCTION); + + // read connectivity + if (rank == MASTER_NODE) cout << "loading connectivity block" << endl; + + line = SkipToFlag(&file_stream, ""); + + unsigned long triCounter = 0; + while (getline(file_stream, line) && !eoConnectivity) { + // check if end of data is reached + if (line.compare("") == 0) eoConnectivity = true; + + if (!eoConnectivity) { + // one line contains values for one triangle (3 points) + istringstream streamTriLine(line); + + // add next line to triangles + for (int iPoint = 0; iPoint < 3; iPoint++) { + streamTriLine >> word; + // Dragon table index starts with 1, convert to c++ indexing starting + // with 0: + triangles[triCounter][iPoint] = stol(word) - 1; + } + } + triCounter++; + } + + if (n_triangles != triCounter - 1) + SU2_MPI::Error( + "Number of read triangles does not match number of points " + "specified in Dragon library header.", + CURRENT_FUNCTION); + + // read hull points + if (rank == MASTER_NODE) cout << "loading hull block" << endl; + + line = SkipToFlag(&file_stream, ""); + + unsigned long hullCounter = 0; + while (getline(file_stream, line) && !eoHull) { + // check if end of data is reached + if (line.compare("") == 0) eoHull = true; + + if (!eoHull) { + // one line contains one point ID for one point on the hull + istringstream streamHullLine(line); + + streamHullLine >> word; + + // Dragon table indices start with 1, convert to c++ indexing starting + // with 0: + hull[hullCounter] = stol(word) - 1; + } + hullCounter++; + } + + if (n_hull_points != hullCounter - 1) + SU2_MPI::Error( + "Number of read hull points does not match number of points " + "specified in Dragon library header.", + CURRENT_FUNCTION); + + file_stream.close(); + +} + +string CFileReaderLUT::SkipToFlag(ifstream* file_stream, const string& flag) { + string line; + getline(*file_stream, line); + + while (line.compare(flag) != 0 && !(*file_stream).eof()) { + getline(*file_stream, line); + } + + if ((*file_stream).eof()) SU2_MPI::Error("Flag not found in file", CURRENT_FUNCTION); + + return line; +} diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp new file mode 100644 index 00000000000..66908a5158f --- /dev/null +++ b/Common/src/containers/CLookUpTable.cpp @@ -0,0 +1,484 @@ +/*! + * \file CLookupTable.cpp + * \brief tabulation of fluid properties + * \author D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../../Common/include/containers/CLookUpTable.hpp" +#include "../../../Common/include/linear_algebra/blas_structure.hpp" +#include "../../../Common/include/toolboxes/CSquareMatrixCM.hpp" + +using namespace std; + +CLookUpTable::CLookUpTable(const string& var_file_name_lut, const string& name_prog, const string& name_enth) { + file_name_lut = var_file_name_lut; + + rank = SU2_MPI::GetRank(); + + LoadTableRaw(var_file_name_lut); + + FindTableLimits(name_prog, name_enth); + + if (rank == MASTER_NODE) + cout << "Detecting all unique edges and setting edge to triangle connectivity " + "..." + << endl; + + IdentifyUniqueEdges(); + + if (rank == MASTER_NODE) + cout << " done." << endl; + + PrintTableInfo(); + + if (rank == MASTER_NODE) + cout << "Building a trapezoidal map for the (progress variable, enthalpy) " + "space ..." + << endl; + + trap_map_prog_enth = CTrapezoidalMap(GetDataP(name_prog), GetDataP(name_enth), table_data.cols(), edges, edge_to_triangle); + + if (rank == MASTER_NODE) + cout << " done." << endl; + + if (rank == MASTER_NODE) + cout << "Precomputing interpolation coefficients..." << endl; + + + ComputeInterpCoeffs(name_prog, name_enth); + + if (rank == MASTER_NODE) + cout << "LUT fluid model ready for use" << endl; +} + +void CLookUpTable::LoadTableRaw(const string& var_file_name_lut) { + CFileReaderLUT file_reader; + + if (rank == MASTER_NODE) + cout << "Loading lookup table, filename = " << var_file_name_lut << " ..." << endl; + + file_reader.ReadRawDRG(var_file_name_lut); + + n_points = file_reader.GetNPoints(); + n_triangles = file_reader.GetNTriangles(); + n_variables = file_reader.GetNVariables(); + n_hull_points = file_reader.GetNHullPoints(); + version_lut = file_reader.GetVersionLUT(); + version_reader = file_reader.GetVersionReader(); + names_var = file_reader.GetNamesVar(); + table_data = file_reader.GetTableData(); + triangles = file_reader.GetTriangles(); + hull = file_reader.GetHull(); + + if (rank == MASTER_NODE) + cout << " done." << endl; +} + +void CLookUpTable::FindTableLimits(const string& name_prog, const string& name_enth) { + int ixEnth = GetIndexOfVar(name_enth); + int ixProg = GetIndexOfVar(name_prog); + + /* we find the lowest and highest value of enthalpy and progress variable in the table */ + limits_table_enth[0] = *min_element(&table_data[ixEnth][0], &table_data[ixEnth][0]+table_data.cols()); + limits_table_enth[1] = *max_element(&table_data[ixEnth][0], &table_data[ixEnth][0]+table_data.cols()); + limits_table_prog[0] = *min_element(&table_data[ixProg][0], &table_data[ixProg][0]+table_data.cols()); + limits_table_prog[1] = *max_element(&table_data[ixProg][0], &table_data[ixProg][0]+table_data.cols()); + +} + +void CLookUpTable::PrintTableInfo() { + if (rank == MASTER_NODE) { + cout << setfill(' '); + cout << endl; + cout << "+------------------------------------------------------------------+\n"; + cout << "| Look-Up-Table (LUT) info |\n"; + cout << "+------------------------------------------------------------------+" << endl; + cout << "| File name:" << setw(54) << right << file_name_lut << " |" << endl; + cout << "| Table version:" << setw(50) << right << version_lut << " |" << endl; + cout << "| Table reader version:" << setw(43) << right << version_reader << " |" << endl; + cout << "| Number of variables:" << setw(44) << right << n_variables << " |" << endl; + cout << "| Number of points:" << setw(47) << right << n_points << " |" << endl; + cout << "| Number of triangles:" << setw(44) << right << n_triangles << " |" << endl; + cout << "| Number of edges:" << setw(48) << right << edges.size() << " |" << endl; + cout << "+------------------------------------------------------------------+" << endl; + cout << "| Minimum enthalpy:" << setw(47) << right << limits_table_enth[0] << " |" << endl; + cout << "| Maximum enthalpy:" << setw(47) << right << limits_table_enth[1] << " |" << endl; + cout << "| Minimum progress variable:" << setw(38) << right << limits_table_prog[0] << " |" << endl; + cout << "| Maximum progress variable:" << setw(38) << right << limits_table_prog[1] << " |" << endl; + cout << "+------------------------------------------------------------------+" << endl; + cout << "| Variable names: |" << endl; + cout << "| |" << endl; + + for (unsigned long i_var = 0; i_var < names_var.size(); i_var++) + cout << "| " << right << setw(3) << i_var << ": " << left << setw(59) << names_var[i_var] << " |" << endl; + + cout << "+------------------------------------------------------------------+" << endl; + + cout << endl; + } +} + +void CLookUpTable::IdentifyUniqueEdges() { + + /* loop through elements and store the vector of element IDs (neighbors) + for each of the points in the table */ + + vector > neighborElemsOfPoint; + + neighborElemsOfPoint.resize(n_points); + for (unsigned long iElem = 0; iElem < n_triangles; iElem++) { + /* loop over 3 points per triangle */ + for (unsigned long iPoint = 0; iPoint < N_POINTS_TRIANGLE; iPoint++) { + /* get the global ID of the current point */ + const unsigned long GlobalIndex = triangles[iElem][iPoint]; + + /* add the current element ID to the neighbor list for this point */ + neighborElemsOfPoint[GlobalIndex].push_back(iElem); + } + } + + /* remove duplicates from the neighboring element lists*/ + vector::iterator vecIt; + for (unsigned long iPoint = 0; iPoint < n_points; iPoint++) { + /* sort neighboring elements for each point */ + sort(neighborElemsOfPoint[iPoint].begin(), neighborElemsOfPoint[iPoint].end()); + + /* uniquify list of neighboring elements */ + vecIt = unique(neighborElemsOfPoint[iPoint].begin(), neighborElemsOfPoint[iPoint].end()); + + /* adjust size of vector */ + neighborElemsOfPoint[iPoint].resize(vecIt - neighborElemsOfPoint[iPoint].begin()); + } + + /* loop through all neighboring elements of each point and store + the point IDs that are neighboring points */ + vector > neighborPointsOfPoint; + neighborPointsOfPoint.resize(n_points); + for (unsigned long iPoint = 0; iPoint < n_points; iPoint++) { + for (unsigned long iElem = 0; iElem < neighborElemsOfPoint[iPoint].size(); iElem++) { + /* loop over element points */ + for (unsigned long jPoint = 0; jPoint < N_POINTS_TRIANGLE; jPoint++) { + /* get the global ID of the current point */ + const unsigned long GlobalIndex = triangles[neighborElemsOfPoint[iPoint][iElem]][jPoint]; + + /* add the current element ID to the neighbor list for this point */ + if (GlobalIndex != iPoint) neighborPointsOfPoint[iPoint].push_back(GlobalIndex); + } + } + } + + /* remove duplicates from the neighboring points lists */ + for (unsigned long iPoint = 0; iPoint < n_points; iPoint++) { + /* sort neighboring points for each point */ + sort(neighborPointsOfPoint[iPoint].begin(), neighborPointsOfPoint[iPoint].end()); + + /* uniquify list of neighboring elements */ + vecIt = unique(neighborPointsOfPoint[iPoint].begin(), neighborPointsOfPoint[iPoint].end()); + + /* adjust size of vector */ + neighborPointsOfPoint[iPoint].resize(vecIt - neighborPointsOfPoint[iPoint].begin()); + } + + /* loop through our point neighbors and fill the vector of the unique + point pairs making up each edge in the grid. We impose a requirement + 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; iPoint++) { + for (unsigned long jPoint = 0; jPoint < neighborPointsOfPoint[iPoint].size(); jPoint++) { + /* Store the neighbor index more clearly. */ + const unsigned long GlobalIndex = neighborPointsOfPoint[iPoint][jPoint]; + + /* Store the edge so that the lower index of the pair is always first. */ + if (iPoint < GlobalIndex) { + vector edge(2); + edge[0] = iPoint; + edge[1] = GlobalIndex; + edges.push_back(edge); + } + } + } + + /* Loop over our edges data structure. For the first point in each + 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.resize(edges.size()); + for (unsigned long iEdge = 0; iEdge < edges.size(); iEdge++) { + /* Store the two points of the edge more clearly. */ + const unsigned long iPoint = edges[iEdge][0]; + const unsigned long jPoint = edges[iEdge][1]; + + /* Loop over all neighboring elements to iPoint. */ + for (unsigned long iElem = 0; iElem < neighborElemsOfPoint[iPoint].size(); iElem++) { + /* loop over 3 points per triangle */ + for (unsigned long kPoint = 0; kPoint < N_POINTS_TRIANGLE; kPoint++) { + /* Get the global ID of the current point. */ + const unsigned long GlobalIndex = triangles[neighborElemsOfPoint[iPoint][iElem]][kPoint]; + + /* Add the current element ID to the neighbor list for this point. */ + if (GlobalIndex == jPoint) edge_to_triangle[iEdge].push_back(neighborElemsOfPoint[iPoint][iElem]); + } + } + } +} + +void CLookUpTable::ComputeInterpCoeffs(const string& name_prog, const string& name_enth) { + /* build KD tree for enthalpy, progress variable space */ + + std::array next_triangle; + + const su2double* prog = GetDataP(name_prog); + const su2double* enth = GetDataP(name_enth); + + /* calculate weights for each triangle (basically a distance function) and + * build inverse interpolation matrices */ + for (unsigned long i_triangle = 0; i_triangle < n_triangles; i_triangle++) { + + for (int p = 0; p < 3; p++) { + next_triangle[p] = triangles[i_triangle][p]; + } + + su2activematrix prog_interp_mat_inv(3, 3); + GetInterpMatInv(prog, enth, next_triangle, prog_interp_mat_inv); + interp_mat_inv_prog_enth.push_back(prog_interp_mat_inv); + } +} + +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); + + /* 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]]; + + 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.Invert(); + global_M.Transpose(); + + for (unsigned int i=0; i= limits_table_prog[0] && val_prog <= limits_table_prog[1]) && + (val_enth >= limits_table_enth[0] && val_enth <= limits_table_enth[1])){ + + /* find the triangle that holds the (prog, enth) point */ + unsigned long id_triangle = trap_map_prog_enth.GetTriangle(val_prog, val_enth); + + if (IsInTriangle(val_prog, val_enth, id_triangle, name_prog, name_enth)) { + /* get interpolation coefficients for point on triangle */ + std::array interp_coeffs{0}; + GetInterpCoeffs(val_prog, val_enth, interp_mat_inv_prog_enth[id_triangle], interp_coeffs); + + /* first, copy the single triangle from the large triangle list*/ + std::array triangle{0}; + for (int p = 0; p < 3; p++) + triangle[p] = triangles[id_triangle][p]; + + *val_var = Interpolate(GetDataP(val_name_var), triangle, interp_coeffs); + exit_code = 0; + } else { + /* in bounding box but outside of table */ + unsigned long nearest_neighbor = FindNearestNeighborOnHull(val_prog, val_enth, name_prog, name_enth); + *val_var = GetDataP(val_name_var)[nearest_neighbor]; + exit_code = 1; + } + } else { + /* lookup is outside of table bounding box */ + unsigned long nearest_neighbor = FindNearestNeighborOnHull(val_prog, val_enth, name_prog, name_enth); + *val_var = GetDataP(val_name_var)[nearest_neighbor]; + exit_code = 1; + } + return exit_code; +} + +unsigned long CLookUpTable::LookUp_ProgEnth(const vector& val_names_var, vector& val_vars, + su2double val_prog, su2double val_enth, const string& name_prog, + const string& name_enth) { + vector look_up_data; + + for (long unsigned int i_var = 0; i_var < val_vars.size(); ++i_var) { + look_up_data.push_back(&val_vars[i_var]); + } + + unsigned long exit_code = LookUp_ProgEnth(val_names_var, look_up_data, val_prog, val_enth, name_prog, name_enth); + + return exit_code; +} + +unsigned long CLookUpTable::LookUp_ProgEnth(const vector& val_names_var, vector& val_vars, + su2double val_prog, su2double val_enth, const string& name_prog, + const string& name_enth) { + unsigned long exit_code = 0; + unsigned long nearest_neighbor = 0; + unsigned long id_triangle; + std::array interp_coeffs{0}; + std::array triangle{0}; + + /* check if progress variable value is in progress variable table range + * and if enthalpy is in enthalpy table range */ + if (val_prog >= limits_table_prog[0] && val_prog <= limits_table_prog[1] && val_enth >= limits_table_enth[0] && + val_enth <= limits_table_enth[1]) { + /* if so, try to find the triangle that holds the (prog, enth) point */ + id_triangle = trap_map_prog_enth.GetTriangle(val_prog, val_enth); + + /* 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_prog, val_enth, id_triangle, name_prog, name_enth)) { + /* if so, get interpolation coefficients for point in the triangle */ + GetInterpCoeffs(val_prog, val_enth, interp_mat_inv_prog_enth[id_triangle], interp_coeffs); + + /* exit_code 0 means point was in triangle */ + exit_code = 0; + + } else { + /* if point is not inside a triangle (outside table domain) search nearest neighbor */ + nearest_neighbor = FindNearestNeighborOnHull(val_prog, val_enth, name_prog, name_enth); + exit_code = 1; + } + + } else { + /* if point is outside of table ranges, find nearest neighbor */ + nearest_neighbor = FindNearestNeighborOnHull(val_prog, val_enth, name_prog, name_enth); + exit_code = 1; + } + + /* loop over variable names and interpolate / get values */ + for (long unsigned int i_var = 0; i_var < val_names_var.size(); ++i_var) { + if (val_names_var[i_var].compare("NULL") == 0) { + *val_vars[i_var] = 0.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[id_triangle][p]; + *val_vars[i_var] = Interpolate(GetDataP(val_names_var[i_var]), triangle, interp_coeffs); + } + else + *val_vars[i_var] = GetDataP(val_names_var[i_var])[nearest_neighbor]; + } + } + return exit_code; +} + +void CLookUpTable::GetInterpCoeffs(su2double val_x, su2double val_y, su2activematrix& interp_mat_inv, + std::array& interp_coeffs) { + + std::array query_vector = {1,val_x,val_y}; + + su2double d; + for (int i = 0; i < 3; i++) { + d = 0; + for (int j = 0; j < 3; j++) { + d = d + interp_mat_inv[i][j] * query_vector[j]; + } + interp_coeffs[i] = d; + } +} + +su2double CLookUpTable::Interpolate(const su2double* val_samples, std::array& val_triangle, + std::array& val_interp_coeffs) { + su2double result = 0; + su2double z = 0; + + for (int i_point = 0; i_point < N_POINTS_TRIANGLE; i_point++) { + z = val_samples[val_triangle[i_point]]; + result += val_interp_coeffs[i_point] * z; + } + + return result; +} + +unsigned long CLookUpTable::FindNearestNeighborOnHull(su2double val_prog, su2double val_enth, const string& name_prog, + const string& name_enth) { + su2double min_distance = 1.e99; + su2double next_distance = 1.e99; + su2double next_prog_norm; + su2double next_enth_norm; + unsigned long neighbor_id = 0; + + const su2double* prog_table = GetDataP(name_prog); + const su2double* enth_table = GetDataP(name_enth); + + su2double norm_coeff_prog = 1. / (limits_table_prog[1] - limits_table_prog[0]); + su2double norm_coeff_enth = 1. / (limits_table_enth[1] - limits_table_enth[0]); + su2double val_prog_norm = val_prog / (limits_table_prog[1] - limits_table_prog[0]); + su2double val_enth_norm = val_enth / (limits_table_enth[1] - limits_table_enth[0]); + + for (unsigned long i_point = 0; i_point < n_hull_points; ++i_point) { + next_prog_norm = prog_table[hull[i_point]] * norm_coeff_prog; + next_enth_norm = enth_table[hull[i_point]] * norm_coeff_enth; + + next_distance = sqrt(pow(val_prog_norm - next_prog_norm, 2) + pow(val_enth_norm - next_enth_norm, 2)); + + if (next_distance < min_distance) { + min_distance = next_distance; + neighbor_id = hull[i_point]; + } + } + return neighbor_id; +} + +bool CLookUpTable::IsInTriangle(su2double val_prog, su2double val_enth, unsigned long val_id_triangle, const string& name_prog, + const string& name_enth) { + su2double tri_prog_0 = GetDataP(name_prog)[triangles[val_id_triangle][0]]; + su2double tri_enth_0 = GetDataP(name_enth)[triangles[val_id_triangle][0]]; + + su2double tri_prog_1 = GetDataP(name_prog)[triangles[val_id_triangle][1]]; + su2double tri_enth_1 = GetDataP(name_enth)[triangles[val_id_triangle][1]]; + + su2double tri_prog_2 = GetDataP(name_prog)[triangles[val_id_triangle][2]]; + su2double tri_enth_2 = GetDataP(name_enth)[triangles[val_id_triangle][2]]; + + su2double area_tri = TriArea(tri_prog_0, tri_enth_0, tri_prog_1, tri_enth_1, tri_prog_2, tri_enth_2); + + su2double area_0 = TriArea(val_prog, val_enth, tri_prog_1, tri_enth_1, tri_prog_2, tri_enth_2); + su2double area_1 = TriArea(tri_prog_0, tri_enth_0, val_prog, val_enth, tri_prog_2, tri_enth_2); + su2double area_2 = TriArea(tri_prog_0, tri_enth_0, tri_prog_1, tri_enth_1, val_prog, val_enth); + + return (abs(area_tri - (area_0 + area_1 + area_2)) < area_tri * 1e-10); +} diff --git a/Common/src/containers/CTrapezoidalMap.cpp b/Common/src/containers/CTrapezoidalMap.cpp new file mode 100644 index 00000000000..a8f884b9c97 --- /dev/null +++ b/Common/src/containers/CTrapezoidalMap.cpp @@ -0,0 +1,233 @@ +/*! + * \file TrapezoidalMap.cpp + * \brief Implementation of the trapezoidal map for tabulation and lookup of fluid properties + * \author D. Mayer, T. Economon + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "../../Common/include/option_structure.hpp" +#include "../../Common/include/containers/CTrapezoidalMap.hpp" + +using namespace std; + +/* Trapezoidal map implementation. Reference: + * M. de Berg, O. Cheong M. van Kreveld, M. Overmars, + * Computational Geometry, Algorithms and Applications pp. 121-146 (2008) + */ +CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size, + vector > const& edges, + vector > const& val_edge_to_triangle) { + int rank = SU2_MPI::GetRank(); + su2double startTime = SU2_MPI::Wtime(); + + edge_to_triangle = vector >(val_edge_to_triangle); + + unique_bands_x.assign(samples_x,samples_x+size); + + /* sort x_bands and make them unique */ + sort(unique_bands_x.begin(), unique_bands_x.end()); + + auto iter = unique(unique_bands_x.begin(), unique_bands_x.end()); + + unique_bands_x.resize(distance(unique_bands_x.begin(), iter)); + + edge_limits_x.resize(edges.size(), 2); + edge_limits_y.resize(edges.size(), 2); + + /* store x and y values of each edge in a vector for a slight speed up + * as it prevents some uncoalesced accesses */ + for (unsigned long j = 0; j < edges.size(); j++) { + edge_limits_x[j][0] = samples_x[edges[j][0]]; + edge_limits_x[j][1] = samples_x[edges[j][1]]; + edge_limits_y[j][0] = samples_y[edges[j][0]]; + edge_limits_y[j][1] = samples_y[edges[j][1]]; + } + + /* number of bands */ + unsigned long n_bands_x = unique_bands_x.size() - 1; + /* band index */ + unsigned long i_band = 0; + /* number of edges */ + unsigned long n_edges = edges.size(); + /* edge index */ + unsigned long i_edge = 0; + /* counter for edges intersects */ + unsigned long n_intersects = 0; + /* lower and upper x value of each band */ + su2double band_lower_x = 0; + su2double band_upper_x = 0; + + su2double x_0; + su2double y_0; + su2double dy_edge; + su2double dx_edge; + su2double x_band_mid; + + /* y values of all intersecting edges for every band */ + y_edge_at_band_mid.resize(unique_bands_x.size() - 1); + + /* loop over bands */ + while (i_band < n_bands_x) { + band_lower_x = unique_bands_x[i_band]; + band_upper_x = unique_bands_x[i_band + 1]; + i_edge = 0; + n_intersects = 0; + + /* loop over edges and determine which edges appear in current band */ + while (i_edge < n_edges) { + /* check if edge intersects the band + * (vertical edges are automatically discarded) */ + if (((edge_limits_x[i_edge][0] <= band_lower_x) and (edge_limits_x[i_edge][1] >= band_upper_x)) or + ((edge_limits_x[i_edge][1] <= band_lower_x) and (edge_limits_x[i_edge][0] >= band_upper_x))) { + y_edge_at_band_mid[i_band].push_back(make_pair(0.0, 0)); + + x_0 = edge_limits_x[i_edge][0]; + y_0 = edge_limits_y[i_edge][0]; + + dy_edge = edge_limits_y[i_edge][1] - edge_limits_y[i_edge][0]; + dx_edge = edge_limits_x[i_edge][1] - edge_limits_x[i_edge][0]; + x_band_mid = (band_lower_x + band_upper_x) / 2.0; + + y_edge_at_band_mid[i_band][n_intersects].first = y_0 + dy_edge / dx_edge * (x_band_mid - x_0); + + /* save edge index so it can later be recalled when searching */ + y_edge_at_band_mid[i_band][n_intersects].second = i_edge; + + n_intersects++; + } + i_edge++; + } + + /* sort edges by their y values. + * note that these y values are unique (i.e. edges cannot + * intersect in a band) */ + sort(y_edge_at_band_mid[i_band].begin(), y_edge_at_band_mid[i_band].end()); + + i_band++; + } + + su2double stopTime = SU2_MPI::Wtime(); + + if (rank == MASTER_NODE) cout << "Construction of trapezoidal map took " << stopTime-startTime << " seconds\n" << endl; +} + + +unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) { +//unsigned long CTrapezoidalMap::GetTriangle(su2double val_x, su2double val_y) const { + /* find x band in which val_x sits */ + pair band = GetBand(val_x); + + /* within that band, find edges which enclose the (val_x, val_y) point */ + pair edges = GetEdges(band, val_x, val_y); + + /* identify the triangle using the two edges */ + std::array triangles_edge_low; + + for (int i=0;i<3;i++) + triangles_edge_low[i] = edge_to_triangle[edges.first][i]; + + std::array triangles_edge_up; + for (int i=0;i<3;i++) + triangles_edge_up[i] = edge_to_triangle[edges.second][i]; + + sort(triangles_edge_low.begin(), triangles_edge_low.end()); + sort(triangles_edge_up.begin(), triangles_edge_up.end()); + + // The intersection of the faces to which upper or lower belongs is + // the face that both belong to. + vector triangle(1); + set_intersection(triangles_edge_up.begin(), triangles_edge_up.end(), triangles_edge_low.begin(), + triangles_edge_low.end(), triangle.begin()); + + return triangle[0]; +} + +pair CTrapezoidalMap::GetBand(su2double val_x) { + unsigned long i_low = 0; + unsigned long i_up = 0; + + /* 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(); + + std::pair::iterator,std::vector::iterator> bounds; + bounds = std::equal_range (unique_bands_x.begin(), unique_bands_x.end(), val_x); + i_up = bounds.first - unique_bands_x.begin(); + i_low = i_up-1; + + return make_pair(i_low, i_up); +} + +pair CTrapezoidalMap::GetEdges(pair val_band, + su2double val_x, su2double val_y) const{ + su2double next_y; + su2double y_edge_low; + su2double y_edge_up; + su2double x_edge_low; + su2double x_edge_up; + + unsigned long i_band_low = val_band.first; + + unsigned long next_edge; + + unsigned long j_low = 0; + unsigned long j_mid = 0; + unsigned long j_up = 0; + + j_up = y_edge_at_band_mid[i_band_low].size() - 1; + j_low = 0; + + while (j_up - j_low > 1) { + j_mid = (j_up + j_low) / 2; + + // Select the edge associated with the x band (i_band_low) + // Search for the RunEdge in the y direction (second value is index of + // edge) + next_edge = y_edge_at_band_mid[i_band_low][j_mid].second; + + y_edge_low = edge_limits_y[next_edge][0]; + y_edge_up = edge_limits_y[next_edge][1]; + x_edge_low = edge_limits_x[next_edge][0]; + x_edge_up = edge_limits_x[next_edge][1]; + + // The search variable in j should be interpolated in i as well + next_y = y_edge_low + (y_edge_up - y_edge_low) / (x_edge_up - x_edge_low) * (val_x - x_edge_low); + + if (next_y > val_y) { + j_up = j_mid; + + } else if (next_y < val_y) { + j_low = j_mid; + + } else if (next_y == val_y) { + j_low = j_mid; + j_up = j_low + 1; + break; + } + } + + unsigned long edge_low = y_edge_at_band_mid[i_band_low][j_low].second; + unsigned long edge_up = y_edge_at_band_mid[i_band_low][j_up].second; + + return make_pair(edge_low, edge_up); +} diff --git a/Common/src/containers/meson.build b/Common/src/containers/meson.build new file mode 100644 index 00000000000..9a0fdd1a1ae --- /dev/null +++ b/Common/src/containers/meson.build @@ -0,0 +1,3 @@ +common_src += files(['CTrapezoidalMap.cpp', + 'CFileReaderLUT.cpp', + 'CLookUpTable.cpp']) diff --git a/Common/src/meson.build b/Common/src/meson.build index b3e0726e70c..3d84e6ea17a 100644 --- a/Common/src/meson.build +++ b/Common/src/meson.build @@ -13,6 +13,7 @@ subdir('geometry/elements') subdir('geometry/dual_grid') subdir('geometry/primal_grid') subdir('geometry/meshreader') +subdir('containers') subdir('interface_interpolation') subdir('fem') subdir('grid_movement') diff --git a/SU2_CFD/obj/Makefile.am b/SU2_CFD/obj/Makefile.am index 612d0483888..eef545d729f 100644 --- a/SU2_CFD/obj/Makefile.am +++ b/SU2_CFD/obj/Makefile.am @@ -86,6 +86,9 @@ libSU2Core_sources = ../src/definition_structure.cpp \ ../src/numerics/transition.cpp \ ../src/numerics/heat.cpp \ ../src/numerics/radiation.cpp \ + ../src/numerics/CLookUpTable.cpp \ + ../src/numerics/CTrapezoidalMap.cpp \ + ../src/numerics/CFileReaderLUT.cpp \ ../src/numerics/flow/convection/roe.cpp \ ../src/numerics/flow/convection/fds.cpp \ ../src/numerics/flow/convection/fvs.cpp \ diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 6d71f7377b0..0725aca2282 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -2339,11 +2339,16 @@ void CIncEulerSolver::BC_Inlet(CGeometry *geometry, CSolver **solver_container, } /*--- check if the inlet node is shared with a viscous wall ---*/ + if (geometry->nodes->GetViscousBoundary(iPoint)) { + /*--- match the velocity and pressure for the viscous wall---*/ + for (iDim = 0; iDim < nDim; iDim++) V_inlet[iDim+prim_idx.Velocity()] = nodes->GetVelocity(iPoint,iDim); - /* pressure obtained from interior */ + + /*--- pressure obtained from interior ---*/ + V_inlet[prim_idx.Pressure()] = nodes->GetPressure(iPoint); } diff --git a/UnitTests/Common/containers/CLookupTable_tests.cpp b/UnitTests/Common/containers/CLookupTable_tests.cpp new file mode 100644 index 00000000000..1d2b439af0d --- /dev/null +++ b/UnitTests/Common/containers/CLookupTable_tests.cpp @@ -0,0 +1,88 @@ +/*! + * \file CLookupTable_tests.cpp + * \brief Unit tests for the lookup table. + * \author N. Beishuizen + * \version 7.3.1 "Blackbird" + * + * SU2 Project Website: https://su2code.github.io + * + * The SU2 Project is maintained by the SU2 Foundation + * (http://su2foundation.org) + * + * Copyright 2012-2022, SU2 Contributors (cf. AUTHORS.md) + * + * SU2 is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + * SU2 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with SU2. If not, see . + */ + +#include "catch.hpp" +#include +#include + +#include +#include + +#include "../../../Common/include/CConfig.hpp" +#include "../../../Common/include/containers/CTrapezoidalMap.hpp" +#include "../../../Common/include/containers/CLookUpTable.hpp" +#include "../../../Common/include/containers/CFileReaderLUT.hpp" + + +TEST_CASE("LUTreader", "[tabulated chemistry]") { + + /*--- smaller and trivial lookup table ---*/ + + CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable.drg","PROGVAR","ENTHALPY"); + + /*--- string names of the controlling variables ---*/ + + string name_CV1 = "PROGVAR"; + string name_CV2 = "ENTHALPY"; + + /*--- look up a single value for density ---*/ + + su2double prog = 0.55; + su2double enth = -0.5; + string look_up_tag = "Density"; + su2double look_up_dat; + look_up_table.LookUp_ProgEnth(look_up_tag, &look_up_dat, prog, enth, name_CV1, name_CV2); + CHECK(look_up_dat == Approx(1.02)); + + /*--- look up a single value for viscosity ---*/ + + prog = 0.6; + enth = 0.9; + look_up_tag = "Viscosity"; + look_up_table.LookUp_ProgEnth(look_up_tag, &look_up_dat, prog, enth, name_CV1, name_CV2); + CHECK(look_up_dat == Approx(0.00007)); + + /* find the table limits */ + + auto limitsEnth = look_up_table.GetTableLimitsEnth(); + CHECK(limitsEnth.first == Approx(-1.0)); + CHECK(limitsEnth.second == Approx(1.0)); + + auto limitsProgvar = look_up_table.GetTableLimitsProg(); + CHECK(limitsProgvar.first == Approx(0.0)); + CHECK(limitsProgvar.second == Approx(1.0)); + + /* lookup value outside of lookup table */ + + prog = 1.10; + enth = 1.1; + look_up_tag = "Density"; + look_up_table.LookUp_ProgEnth(look_up_tag, &look_up_dat, prog, enth, name_CV1, name_CV2); + CHECK(look_up_dat == Approx(1.2)); + +} + diff --git a/UnitTests/Common/containers/lookuptable.drg b/UnitTests/Common/containers/lookuptable.drg new file mode 100644 index 00000000000..797e6ec8c61 --- /dev/null +++ b/UnitTests/Common/containers/lookuptable.drg @@ -0,0 +1,73 @@ +Dragon library + +
+[version] +1.0.0 + +[number of points] +12 + +[number of triangles] +14 + +[number of hull points] +8 + +progress variable definition +prog var = 1*Y-CH4 + +[progress variable range] +0.0 | 1.0 + +[enthalpy range] +-1.0 | 1.0 + +[number of variables] +4 + +[variable names] +1:PROGVAR 2:ENTHALPY 3:Density 4:Viscosity +
+ + +0.0 -1.0 1.0 1.0e-5 +0.3 -1.0 1.0 1.2e-5 +0.6 -1.0 1.0 1.4e-5 +1.0 -1.0 1.0 1.6e-5 +0.5 0.25 1.05 2.0e-5 +0.75 0.25 1.05 3.0e-5 +0.3 0.5 1.1 3.0e-5 +0.6 0.5 1.1 4.0e-5 +1.0 0.5 1.1 5.0e-5 +0.0 1.0 1.2 6.0e-5 +0.3 1.0 1.2 7.0e-5 +1.0 1.0 1.2 8.0e-5 + + + +1 2 7 +1 7 11 +1 11 10 +2 5 7 +2 3 5 +3 8 5 +5 8 7 +7 8 11 +3 6 8 +3 4 6 +4 9 6 +6 9 8 +8 9 12 +8 12 11 + + + +1 +2 +3 +4 +9 +12 +11 +10 + diff --git a/UnitTests/meson.build b/UnitTests/meson.build index b0be91c1356..3875cf1b576 100644 --- a/UnitTests/meson.build +++ b/UnitTests/meson.build @@ -11,6 +11,7 @@ su2_cfd_tests = files(['Common/geometry/primal_grid/CPrimalGrid_tests.cpp', 'Common/toolboxes/C1DInterpolation_tests.cpp', 'Common/vectorization.cpp', 'Common/toolboxes/ndflattener_tests.cpp', + 'Common/containers/CLookupTable_tests.cpp', 'SU2_CFD/numerics/CNumerics_tests.cpp', 'SU2_CFD/gradients.cpp', 'SU2_CFD/windowing.cpp']) diff --git a/UnitTests/test_driver.cpp b/UnitTests/test_driver.cpp index 0eb02a4a98d..2fb37dd1004 100644 --- a/UnitTests/test_driver.cpp +++ b/UnitTests/test_driver.cpp @@ -25,7 +25,7 @@ * License along with SU2. If not, see . */ -/*--- This marco tells Catch2 that we are defining a custom main (see +/*--- This macro tells Catch2 that we are defining a custom main (see * below) This should only be done in one file. It should not be done * in any unit test files. ---*/ #define CATCH_CONFIG_RUNNER