From 54801809083d15a3e6126da93a12425c500e59ee Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Wed, 10 Dec 2025 20:01:18 +0100 Subject: [PATCH 1/8] CTrapezoidalMapFAST implementation for FGM lookup tables --- Common/include/containers/CLookUpTable.hpp | 142 +--- .../containers/CTrapezoidalMapFAST.hpp | 755 ++++++++++++++++++ Common/include/option_structure.hpp | 4 +- Common/src/containers/CLookUpTable.cpp | 128 ++- SU2_CFD/src/fluid/CFluidFlamelet.cpp | 26 +- .../Common/containers/CLookupTable_tests.cpp | 115 ++- 6 files changed, 1030 insertions(+), 140 deletions(-) create mode 100644 Common/include/containers/CTrapezoidalMapFAST.hpp diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index e7d8fd17e6e5..5d83e19bff8c 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -35,6 +35,7 @@ #include "../../Common/include/option_structure.hpp" #include "CFileReaderLUT.hpp" #include "CTrapezoidalMap.hpp" +#include "CTrapezoidalMapFAST.hpp" /*! * \brief Look up table. @@ -63,6 +64,7 @@ class CLookUpTable { std::vector z_values_levels; /*!< \brief Constant z-values of each table level.*/ unsigned short table_dim = 2; /*!< \brief Table dimension.*/ + /*! * \brief The lower and upper limits of the z, y and x variable for each table level. */ @@ -100,10 +102,16 @@ class CLookUpTable { su2vector> hull; /*! \brief - * Trapezoidal map objects for the table levels. + * Trapezoidal map objects for the table levels (original SU2 implementation). */ su2vector trap_map_x_y; + /*! \brief + * Memory-efficient trapezoidal map objects for LUT_FAST. + */ + su2vector trap_map_x_y_FAST; + bool use_fast_trap_map_{false}; + /*! \brief * Vector of all the weight factors for the interpolation. */ @@ -186,24 +194,12 @@ class CLookUpTable { /*! * \brief Perform linear interpolation between two table levels for a single variable. - * \param[in] val_CV3 - Value of the third controlling variable at the query point. - * \param[in] lower_level - Table level index of the table level directly below the query point. - * \param[in] upper_level - Table level index of the table level directly above the query point. - * \param[in] lower_value - Result from x-y interpolation on the lower table level. - * \param[in] upper_value - Result from x-y interpolation on the upper table level. - * \param[in] val_vars - Pointer to interpolation result. */ void Linear_Interpolation(const su2double val_CV3, const unsigned long lower_level, const unsigned long upper_level, su2double& lower_value, su2double& upper_value, su2double*& var_vals) const; /*! * \brief Perform linear interpolation between two table levels for a vector of variables. - * \param[in] val_CV3 - Value of the third controlling variable at the query point. - * \param[in] lower_level - Table level index of the table level directly below the query point. - * \param[in] upper_level - Table level index of the table level directly above the query point. - * \param[in] lower_values - Results from x-y interpolation on the lower table level. - * \param[in] upper_values - Results from x-y interpolation on the upper table level. - * \param[in] val_vars - Pointer to interpolation results for all interpolation variables. */ void Linear_Interpolation(const su2double val_CV3, const unsigned long lower_level, const unsigned long upper_level, std::vector& lower_values, std::vector& upper_values, @@ -211,73 +207,34 @@ class CLookUpTable { /*! * \brief Find the point on the hull (boundary of the table) that is closest to the point P(val_CV1,val_CV2). - * \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. - * \returns Point id of the nearest neighbor on the hull. */ unsigned long FindNearestNeighborOnHull(su2double val_CV1, su2double val_CV2, 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] 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. */ void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, 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. - * \param[in] val_CV2 - Second coordinate of point P(val_CV1,val_CV2) to check. - * \param[in] name_var - Variable name to look up. - * \param[out] var_val - Pointer to the to be interpolated variable. - */ void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, const std::string& name_var, su2double* var_val, const unsigned long i_level = 0); /*! * \brief Determine if a point P(val_CV1,val_CV2) is inside the triangle val_id_triangle. - * \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] val_id_triangle - ID of the triangle to check. - * \returns True if the point is in the triangle, false if it is outside. */ bool IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, unsigned long i_level = 0); /*! * \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); @@ -285,11 +242,6 @@ class CLookUpTable { /*! * \brief Compute the values of the first and second controlling variable based on normalized query coordinates - * \param[in] inclusion_levels - Pair containing lower(first) and upper(second) table inclusion level indices. - * \param[in] val_CV1 - First controlling variable value. - * \param[in] val_CV2 - Second controlling variable value. - * \param[in] val_CV3 - Third controlling variable value. - * \returns Array with first and second controlling variable values for the upper and lower table inclusion levels. */ std::array, 2> ComputeNormalizedXY(std::pair& inclusion_levels, const su2double val_CV1, const su2double val_CV2, @@ -297,21 +249,12 @@ class CLookUpTable { /*! * \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); @@ -319,6 +262,24 @@ class CLookUpTable { public: CLookUpTable(const std::string& file_name_lut, std::string name_CV1_in, std::string name_CV2_in); + /*! + * \brief Build the original SU2 trapezoidal map (DAG-based). + * WARNING: This can use a lot of memory for large tables! + * Use EnableFastTrapMap() for large tables instead. + */ + void BuildOriginalTrapMap(); + + /*! + * \brief Enable Memory-efficient trapezoidal map (LUT_FAST mode). + * This builds band-based maps with O(n) memory instead of the original O(n log n) DAG. + */ + void EnableFastTrapMap(); + + /*! + * \brief Check if fast trapezoidal map is enabled. + */ + inline bool UsingFastTrapMap() const { return use_fast_trap_map_; } + /*! * \brief Print information to screen. */ @@ -326,76 +287,32 @@ class CLookUpTable { /*! * \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 whether query is inside (true) or outside (false) data set. */ 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 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 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, const su2double val_CV2, const unsigned long i_level = 0); - /*! - * \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. - * \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. - */ + 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 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. - */ 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. - * \param[in] val_CV3 - Value of controlling variable 3. - * \returns Pair of inclusion level indices (first = lower level index, second = upper level index). */ std::pair FindInclusionLevels(const su2double val_CV3); /*! * \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(const unsigned long i_level = 0) const { return limits_table_y[i_level]; @@ -403,7 +320,6 @@ 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(const unsigned long i_level = 0) const { return limits_table_x[i_level]; @@ -416,8 +332,6 @@ class CLookUpTable { /*! * \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(); diff --git a/Common/include/containers/CTrapezoidalMapFAST.hpp b/Common/include/containers/CTrapezoidalMapFAST.hpp new file mode 100644 index 000000000000..94e5053b3f56 --- /dev/null +++ b/Common/include/containers/CTrapezoidalMapFAST.hpp @@ -0,0 +1,755 @@ +/*! + * \file CTrapezoidalMapFAST.hpp + * \brief Fast trapezoidal map implementation for 2D LUT queries in SU2. + * Based on Pedro Gomes' (pcarruscag) LUT implementation: + * https://github.com/pcarruscag/LUT + * + * This is a faithful adaptation of Pedro's algorithm for SU2's FGM module. + * The algorithm uses band-based spatial indexing with O(log n) query time. + * + * NOTE: This implementation is C++11 compatible for SU2. + */ + +#ifndef CTRAPEZOIDALMAP_FAST_HPP +#define CTRAPEZOIDALMAP_FAST_HPP + +#pragma once + +#include +#include +#include +#include +#include +#include + +namespace su2_lut { + +/*--- Type aliases matching Pedro's implementation ---*/ +using IntT = int32_t; +using RealT = su2double; + +/*--- Simple Matrix container (row-major) ---*/ +template +struct Matrix { + std::vector data; + + void resize(size_t rows, size_t) { data.resize(rows * N); } + size_t rows() const { return data.size() / N; } + + const T& operator()(size_t i, size_t j) const { return data[i * N + j]; } + T& operator()(size_t i, size_t j) { return data[i * N + j]; } +}; + +using Matrix2i = Matrix; +using Matrix3i = Matrix; +using VectorInt = std::vector; +using VectorReal = std::vector; + +/*--- Trapezoidal Map structure ---*/ +struct TrapezoidalMap { + VectorInt offsets, edge_id; + VectorReal x_bands, edge_y; +}; + +/*--- C++11 compatible comparison functors ---*/ +struct PointComparator { + const VectorReal& x; + const VectorReal& y; + + PointComparator(const VectorReal& x_ref, const VectorReal& y_ref) + : x(x_ref), y(y_ref) {} + + bool operator()(const IntT i, const IntT j) const { + return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; + } +}; + +struct EdgeComparator { + bool operator()(const std::array& a, const std::array& b) const { + return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); + } +}; + +struct PairSecondComparator { + bool operator()(const std::pair& a, const std::pair& b) const { + return a.second < b.second; + } +}; + +/*! + * \brief Orders points by ascending x coordinates and updates triangle indices. + * This is critical - the algorithm assumes points are sorted by x. + */ +inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { + const IntT n_pts = static_cast(x.size()); + + // Create sort permutation based on (x, y) coordinates + std::vector perm(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + perm[i] = i; + } + std::sort(perm.begin(), perm.end(), PointComparator(x, y)); + + // Reorder coordinate arrays + VectorReal x_tmp(n_pts), y_tmp(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + x_tmp[i] = x[perm[i]]; + y_tmp[i] = y[perm[i]]; + } + x = x_tmp; + y = y_tmp; + + // Build inverse permutation to update triangle indices + std::vector inv_perm(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + inv_perm[perm[i]] = i; + } + + // Update all triangle vertex indices to use new point ordering + for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { + for (IntT j = 0; j < 3; ++j) { + triangles(i, j) = inv_perm[triangles(i, j)]; + } + } +} + +/*! + * \brief Helper to check if two edges are equal (same endpoints). + */ +inline bool EdgesEqual(const std::array& a, const std::array& b) { + return a[0] == b[0] && a[1] == b[1]; +} + +/*! + * \brief Extracts unique edges from triangles with face adjacency information. + * Each edge stores which triangles (up to 2) share it. + * Boundary edges have edge_faces(i,1) = -1. + */ +inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, + Matrix2i& edge_faces) { + // Extract all edges (3 per triangle), storing (min_pt, max_pt, triangle_id) + std::vector > edges; + edges.resize(3 * triangles.rows()); + + for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { + for (IntT i = 0; i < 3; ++i) { + const IntT j = (i + 1) % 3; + const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); + const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); + std::array edge_data = {{i_pt, j_pt, i_tri}}; + edges[3 * i_tri + i] = edge_data; + } + } + + // Sort to identify duplicate edges (shared between triangles) + std::sort(edges.begin(), edges.end(), EdgeComparator()); + + // Count unique edges + IntT n_edges = 1; + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + n_edges += static_cast(!EdgesEqual(edges[i], edges[i - 1])); + } + + // Build edge_pts and edge_faces arrays + edge_pts.resize(n_edges, 2); + edge_faces.resize(n_edges, 2); + + IntT pos = 0; + + // First edge + edge_pts(pos, 0) = edges[0][0]; + edge_pts(pos, 1) = edges[0][1]; + edge_faces(pos, 0) = edges[0][2]; + edge_faces(pos, 1) = -1; + ++pos; + + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + if (EdgesEqual(edges[i], edges[i - 1])) { + // Duplicate edge - record second adjacent triangle + edge_faces(pos - 1, 1) = edges[i][2]; + } else { + // New edge + edge_pts(pos, 0) = edges[i][0]; + edge_pts(pos, 1) = edges[i][1]; + edge_faces(pos, 0) = edges[i][2]; + edge_faces(pos, 1) = -1; + ++pos; + } + } +} + +/*! + * \brief Band detection result structure (for C++11 compatibility). + */ +struct BandResult { + IntT n_bands; + std::vector pt_to_band; + VectorReal x_bands; +}; + +/*! + * \brief Detects bands based on unique x-coordinates. + * Uses coarse banding to limit memory for large tables. + */ +inline BandResult DetectBands(const VectorReal& x, IntT max_bands = 0) { + BandResult result; + + // Safety check: empty input + if (x.empty()) { + result.n_bands = 0; + return result; + } + + // First pass: count unique x-coordinates + IntT n_unique = 1; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != x[i - 1]) n_unique++; + } + + // Determine if we need to use coarser bands + // Default max_bands: limit to sqrt(n_points) * 4, capped at 5000 + if (max_bands <= 0) { + max_bands = std::min(IntT(5000), static_cast(4.0 * std::sqrt(static_cast(x.size())))); + } + + const bool use_coarse_bands = (n_unique > max_bands); + + if (!use_coarse_bands) { + // Original algorithm: use all unique x-coordinates as band boundaries + result.pt_to_band.resize(x.size()); + result.pt_to_band[0] = 0; + result.n_bands = 0; + + for (IntT i = 1; i < static_cast(x.size()); ++i) { + result.n_bands += static_cast(x[i] != x[i - 1]); + result.pt_to_band[i] = result.n_bands; + } + + result.x_bands.resize(result.n_bands + 1); + IntT pos = 0; + result.x_bands[pos] = x[0]; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != result.x_bands[pos]) { + result.x_bands[++pos] = x[i]; + } + } + } else { + // Coarse banding: divide x-range into max_bands equal intervals + const RealT x_min = x.front(); // x is sorted + const RealT x_max = x.back(); + const RealT x_range = x_max - x_min; + + // Handle degenerate case where all x are the same + if (x_range < 1e-30) { + result.n_bands = 1; + result.pt_to_band.resize(x.size(), 0); + result.x_bands.resize(2); + result.x_bands[0] = x_min; + result.x_bands[1] = x_max; + return result; + } + + result.n_bands = max_bands; + const RealT band_width = x_range / max_bands; + + // Build band boundaries + result.x_bands.resize(max_bands + 1); + for (IntT i = 0; i <= max_bands; ++i) { + result.x_bands[i] = x_min + i * band_width; + } + // Ensure last band includes x_max exactly + result.x_bands[max_bands] = x_max; + + // Map each point to its band + result.pt_to_band.resize(x.size()); + for (IntT i = 0; i < static_cast(x.size()); ++i) { + // Calculate band index: floor((x - x_min) / band_width) + IntT band = static_cast((x[i] - x_min) / band_width); + // Clamp to valid range [0, n_bands-1] + band = std::max(IntT(0), std::min(band, max_bands - 1)); + result.pt_to_band[i] = band; + } + } + + return result; +} + +/*! + * \brief Builds the trapezoidal map for efficient 2D queries. + * For each band, stores edges crossing it, sorted by y-coordinate at band midpoint. + */ +inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, + const VectorReal& y, TrapezoidalMap& map) { + VectorReal& x_bands = map.x_bands; + VectorInt& offsets = map.offsets; + VectorInt& edge_id = map.edge_id; + VectorReal& edge_y = map.edge_y; + + BandResult band_result = DetectBands(x); + const IntT n_bands = band_result.n_bands; + x_bands = band_result.x_bands; + + // Safety check + if (n_bands <= 0 || x_bands.size() < 2) { + x_bands.clear(); + offsets.clear(); + edge_id.clear(); + edge_y.clear(); + return; + } + + // Helper function to find band index for an x-coordinate + // Uses binary search on x_bands + auto find_band = [&x_bands, n_bands](RealT x_val) -> IntT { + // Binary search for the band containing x_val + VectorReal::const_iterator it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); + IntT idx = static_cast(it - x_bands.begin()); + // lower_bound returns iterator to first element >= x_val + // We want the band to the left, so subtract 1 (but not below 0) + idx = std::max(IntT(0), idx - 1); + // Also cap at n_bands - 1 + idx = std::min(idx, n_bands - 1); + return idx; + }; + + // Count edges per band - use actual x-coordinates to determine band range + VectorInt& counts = offsets; + counts.clear(); + counts.resize(n_bands + 1, 0); + + for (IntT i = 0; i < static_cast(edge_pts.rows()); ++i) { + const IntT pt_0 = edge_pts(i, 0); + const IntT pt_1 = edge_pts(i, 1); + const RealT x_0 = x[pt_0]; + const RealT x_1 = x[pt_1]; + + // Find bands for each endpoint using actual x-coordinates + const IntT band_0 = find_band(x_0); + const IntT band_1 = find_band(x_1); + + // Edge spans from min_band to max_band + const IntT min_band = std::min(band_0, band_1); + const IntT max_band = std::max(band_0, band_1); + + // Count edge in each band it crosses + for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { + ++counts[j + 1]; + } + } + + // Convert counts to offsets (CSR format) + for (IntT i = 2; i < static_cast(offsets.size()); ++i) { + offsets[i] += offsets[i - 1]; + } + + // Check total memory requirement before allocation + const size_t total_entries = static_cast(offsets.back()); + const size_t memory_bytes = total_entries * (sizeof(IntT) + sizeof(RealT)); + const size_t memory_mb = memory_bytes / (1024 * 1024); + + // Memory limit: 2GB for edge storage (adjustable) + const size_t MAX_MEMORY_MB = 2048; + if (memory_mb > MAX_MEMORY_MB) { + // Too much memory - clear and return empty map + x_bands.clear(); + offsets.clear(); + edge_id.clear(); + edge_y.clear(); + return; + } + + // Fill edge data for each band + edge_id.resize(total_entries); + edge_y.resize(total_entries); + VectorInt pos = offsets; // Working copy of offsets + + for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { + const IntT pt_0 = edge_pts(i_edge, 0); + const IntT pt_1 = edge_pts(i_edge, 1); + const RealT x_0 = x[pt_0]; + const RealT x_1 = x[pt_1]; + const RealT y_0 = y[pt_0]; + const RealT y_1 = y[pt_1]; + + // Find bands for each endpoint + const IntT band_0 = find_band(x_0); + const IntT band_1 = find_band(x_1); + const IntT min_band = std::min(band_0, band_1); + const IntT max_band = std::max(band_0, band_1); + + // Calculate slope (handle vertical edges) + const RealT dx = x_1 - x_0; + const bool is_vertical = (std::abs(dx) < 1e-30); + const RealT dy_dx = is_vertical ? 0.0 : (y_1 - y_0) / dx; + + // Add edge to each band it crosses + for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { + edge_id[pos[j]] = i_edge; + const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2.0; + // Calculate y at band midpoint + if (is_vertical) { + edge_y[pos[j]] = (y_0 + y_1) / 2.0; // Use average y for vertical edges + } else { + edge_y[pos[j]] = y_0 + dy_dx * (x_mid - x_0); + } + ++pos[j]; + } + } + + // Sort edges within each band by y-coordinate + std::vector > tmp; + for (IntT i = 0; i < n_bands; ++i) { + const IntT begin = offsets[i]; + const IntT end = offsets[i + 1]; + tmp.resize(end - begin); + + for (IntT k = begin; k < end; ++k) { + tmp[k - begin] = std::make_pair(edge_id[k], edge_y[k]); + } + std::sort(tmp.begin(), tmp.end(), PairSecondComparator()); + + for (IntT k = begin; k < end; ++k) { + edge_id[k] = tmp[k - begin].first; + edge_y[k] = tmp[k - begin].second; + } + } +} + +/*! + * \brief Query result structure (for C++11 compatibility). + */ +struct QueryResult { + IntT edge_below; + IntT edge_above; +}; + +/*! + * \brief Query the trapezoidal map for edges bounding a point. + * Returns edges below and above the point. Either can be -1 if point is at boundary. + * With coarse banding, recomputes edge y-positions at the actual query x. + * Also searches adjacent bands to avoid missing candidates. + */ +inline QueryResult QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, + const VectorReal& x_coords, const VectorReal& y_coords, + const RealT& x, const RealT& y) { + QueryResult result; + result.edge_below = -1; + result.edge_above = -1; + + // Safety check: empty map + if (map.x_bands.size() < 2 || map.offsets.empty()) { + return result; + } + + const IntT n_bands = static_cast(map.x_bands.size()) - 1; + + // Find the band containing x + const VectorReal& x_bands = map.x_bands; + VectorReal::const_iterator it = std::lower_bound(x_bands.begin(), x_bands.end(), x); + IntT d = static_cast(it - x_bands.begin()); + IntT center_band = std::min(std::max(IntT(0), d - 1), n_bands - 1); + + // Search in current band AND adjacent bands to catch edges near boundaries + RealT best_y_below = -1e300; + RealT best_y_above = 1e300; + + for (IntT band_offset = -1; band_offset <= 1; ++band_offset) { + IntT band_idx = center_band + band_offset; + if (band_idx < 0 || band_idx >= n_bands) continue; + if (band_idx + 1 >= static_cast(map.offsets.size())) continue; + + const IntT begin = map.offsets[band_idx]; + const IntT end = map.offsets[band_idx + 1]; + + if (begin < 0 || end > static_cast(map.edge_id.size()) || begin >= end) { + continue; + } + + for (IntT k = begin; k < end; ++k) { + const IntT e_id = map.edge_id[k]; + + // Get edge endpoints + const IntT p0 = edge_pts(e_id, 0); + const IntT p1 = edge_pts(e_id, 1); + const RealT x0 = x_coords[p0], y0 = y_coords[p0]; + const RealT x1 = x_coords[p1], y1 = y_coords[p1]; + + // Check if query x is within the edge's x-extent + const RealT x_min_edge = std::min(x0, x1); + const RealT x_max_edge = std::max(x0, x1); + if (x < x_min_edge - 1e-10 || x > x_max_edge + 1e-10) { + continue; // Query x is outside this edge's range + } + + // Compute y-position of edge at query x + RealT edge_y_at_x; + const RealT dx = x1 - x0; + if (std::abs(dx) < 1e-30) { + // Vertical edge - use average y + edge_y_at_x = (y0 + y1) / 2.0; + } else { + // Linear interpolation + const RealT t = (x - x0) / dx; + edge_y_at_x = y0 + t * (y1 - y0); + } + + // Check if this edge is below or above the query point + if (edge_y_at_x <= y + 1e-10) { + if (edge_y_at_x > best_y_below) { + best_y_below = edge_y_at_x; + result.edge_below = e_id; + } + } + if (edge_y_at_x >= y - 1e-10) { + if (edge_y_at_x < best_y_above) { + best_y_above = edge_y_at_x; + result.edge_above = e_id; + } + } + } + } + + return result; +} + +/*! + * \brief Get triangles adjacent to two query edges (up to 3 triangles). + */ +inline std::array AdjacentTriangles(const IntT edge_0, const IntT edge_1, + const Matrix2i& edge_faces) { + std::array tris = {{-1, -1, -1}}; + IntT pos = 0; + + // Helper to insert unique triangle + if (edge_0 >= 0) { + IntT t0 = edge_faces(edge_0, 0); + IntT t1 = edge_faces(edge_0, 1); + if (t0 >= 0) { + tris[pos++] = t0; + } + if (t1 >= 0) { + bool found = false; + for (IntT k = 0; k < pos; ++k) { + if (t1 == tris[k]) { found = true; break; } + } + if (!found) tris[pos++] = t1; + } + } + + if (edge_1 >= 0) { + IntT t0 = edge_faces(edge_1, 0); + IntT t1 = edge_faces(edge_1, 1); + if (t0 >= 0) { + bool found = false; + for (IntT k = 0; k < pos; ++k) { + if (t0 == tris[k]) { found = true; break; } + } + if (!found && pos < 3) tris[pos++] = t0; + } + if (t1 >= 0) { + bool found = false; + for (IntT k = 0; k < pos; ++k) { + if (t1 == tris[k]) { found = true; break; } + } + if (!found && pos < 3) tris[pos++] = t1; + } + } + + return tris; +} + +/*! + * \brief Compute barycentric coordinates and check if point is inside triangle. + */ +inline std::array TriangleCoords(const IntT i_tri, const Matrix3i& triangles, + const VectorReal& x, const VectorReal& y, + const RealT x_q, const RealT y_q) { + const IntT i0 = triangles(i_tri, 0); + const IntT i1 = triangles(i_tri, 1); + const IntT i2 = triangles(i_tri, 2); + + const RealT x0 = x[i0], y0 = y[i0]; + const RealT x1 = x[i1], y1 = y[i1]; + const RealT x2 = x[i2], y2 = y[i2]; + + const RealT det = (y1 - y2) * (x0 - x2) + (x2 - x1) * (y0 - y2); + const RealT l0 = ((y1 - y2) * (x_q - x2) + (x2 - x1) * (y_q - y2)) / det; + const RealT l1 = ((y2 - y0) * (x_q - x2) + (x0 - x2) * (y_q - y2)) / det; + const RealT l2 = 1.0 - l0 - l1; + + std::array coords = {{l0, l1, l2}}; + return coords; +} + +/*! + * \brief Check if point is inside triangle (barycentric coords all >= 0). + */ +inline bool InTriangle(const std::array& coords, const RealT tol = 0.0) { + return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; +} + +/*! + * \brief Find triangle containing point (x_q, y_q) using the trapezoidal map. + * Returns triangle index and barycentric coordinates. + * Returns -1 if point is outside the mesh. + */ +inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, + const Matrix2i& edge_pts, const Matrix2i& edge_faces, + const VectorReal& x, const VectorReal& y, + const RealT x_q, const RealT y_q, std::array& bary_out) { + + QueryResult qr = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); + std::array candidates = AdjacentTriangles(qr.edge_below, qr.edge_above, edge_faces); + + const RealT tol = 1e-12; + for (IntT i = 0; i < 3; ++i) { + if (candidates[i] < 0) continue; + + std::array coords = TriangleCoords(candidates[i], triangles, x, y, x_q, y_q); + if (InTriangle(coords, tol)) { + bary_out = coords; + return candidates[i]; + } + } + + // Not found - return -1 + bary_out[0] = bary_out[1] = bary_out[2] = 0.0; + return -1; +} + +} // namespace su2_lut + + +/*! + * \class CTrapezoidalMapFAST + * \brief Wrapper class for SU2 integration of Pedro's trapezoidal map algorithm. + */ +class CTrapezoidalMapFAST { +private: + su2_lut::Matrix3i triangles; + su2_lut::Matrix2i edge_pts, edge_faces; + su2_lut::VectorReal x_coords, y_coords; + su2_lut::TrapezoidalMap map; + + unsigned long n_points; + unsigned long n_triangles; + bool build_successful; + +public: + CTrapezoidalMapFAST() : n_points(0), n_triangles(0), build_successful(false) {} + + /*! + * \brief Build interface matching CLookUpTable.cpp's expected signature. + * \param[in] num_points - Number of mesh points + * \param[in] num_triangles - Number of triangles + * \param[in] x - X coordinates of all points + * \param[in] y - Y coordinates of all points + * \param[in] connectivity - Triangle connectivity (flat array: tri0_p0, tri0_p1, tri0_p2, ...) + * \return true if map was built successfully, false if memory limit exceeded + */ + bool Build(unsigned long num_points, unsigned long num_triangles, + const su2double* x, const su2double* y, + const unsigned long* connectivity) { + + n_points = num_points; + n_triangles = num_triangles; + build_successful = false; + + // Safety check: handle empty input + if (num_points == 0 || num_triangles == 0) { + x_coords.clear(); + y_coords.clear(); + triangles.resize(0, 3); + edge_pts.resize(0, 2); + edge_faces.resize(0, 2); + map.x_bands.clear(); + map.offsets.clear(); + map.edge_id.clear(); + map.edge_y.clear(); + return false; + } + + // Copy coordinates + x_coords.assign(x, x + num_points); + y_coords.assign(y, y + num_points); + + // Copy triangle connectivity with bounds checking + triangles.resize(num_triangles, 3); + for (size_t i = 0; i < num_triangles; ++i) { + unsigned long v0 = connectivity[3*i + 0]; + unsigned long v1 = connectivity[3*i + 1]; + unsigned long v2 = connectivity[3*i + 2]; + + // Bounds check + if (v0 >= num_points || v1 >= num_points || v2 >= num_points) { + // Invalid triangle - skip by setting to valid indices + triangles(i, 0) = 0; + triangles(i, 1) = 0; + triangles(i, 2) = 0; + } else { + triangles(i, 0) = static_cast(v0); + triangles(i, 1) = static_cast(v1); + triangles(i, 2) = static_cast(v2); + } + } + + // Build the map using Pedro's algorithm + su2_lut::ReorderPoints(triangles, x_coords, y_coords); + su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); + su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); + + // Check if map was successfully built (not empty due to memory limit) + if (map.x_bands.empty() || map.offsets.empty()) { + build_successful = false; + return false; + } + + build_successful = true; + return true; + } + + /*! + * \brief Check if the map was successfully built. + */ + bool IsBuilt() const { return build_successful; } + + /*! + * \brief Find the triangle containing a query point. + * \param[in] val_x - X coordinate of query point + * \param[in] val_y - Y coordinate of query point + * \param[out] triangle_id - ID of the containing triangle + * \param[out] bary_coords - Barycentric coordinates of query point in triangle + * \return true if point is inside the mesh, false otherwise + */ + bool FindTriangle(su2double val_x, su2double val_y, + unsigned long& triangle_id, + std::array& bary_coords) const { + + // Safety check: empty map + if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { + bary_coords[0] = bary_coords[1] = bary_coords[2] = 0.0; + return false; + } + + std::array bary; + su2_lut::IntT tri_id = su2_lut::FindTriangle( + map, triangles, edge_pts, edge_faces, x_coords, y_coords, val_x, val_y, bary); + + if (tri_id >= 0) { + triangle_id = static_cast(tri_id); + bary_coords[0] = bary[0]; + bary_coords[1] = bary[1]; + bary_coords[2] = bary[2]; + return true; + } + return false; + } + + /*--- Accessors ---*/ + unsigned long GetNTriangles() const { return n_triangles; } + unsigned long GetNPoints() const { return n_points; } + unsigned long GetNEdges() const { return edge_pts.rows(); } + unsigned long GetNBands() const { return map.x_bands.size() > 0 ? map.x_bands.size() - 1 : 0; } +}; + +#endif // CTRAPEZOIDALMAP_FAST_HPP diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index bf99257b1387..ad7b2ac9f13e 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -599,12 +599,14 @@ MakePair("ONESPECIES", ONESPECIES) */ enum class ENUM_DATADRIVEN_METHOD { LUT = 0, - MLP = 1 + MLP = 1, + LUT_FAST = 2 /* New LUT implementation based on Pedro's code */ }; static const MapType DataDrivenMethod_Map = { MakePair("LUT", ENUM_DATADRIVEN_METHOD::LUT) MakePair("MLP", ENUM_DATADRIVEN_METHOD::MLP) + MakePair("LUT_FAST", ENUM_DATADRIVEN_METHOD::LUT_FAST) }; /*! diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 98722a6b5b28..05344a5484c9 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -47,9 +47,7 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, FindTableLimits(name_CV1, name_CV2); if (rank == MASTER_NODE) - cout << "Detecting all unique edges and setting edge to triangle connectivity " - "..." - << endl; + cout << "Detecting all unique edges and setting edge to triangle connectivity ..." << endl; IdentifyUniqueEdges(); @@ -60,18 +58,27 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, /* Add additional variable index which will always result in zero when looked up. */ idx_null = names_var.size(); + /*--- NOTE: Trapezoidal map is NOT built here anymore! ---*/ + /*--- Call EnableFastTrapMap() for LUT_FAST (efficient O(n) memory) ---*/ + /*--- Call BuildOriginalTrapMap() for standard LUT (original SU2 implementation) ---*/ + + ComputeInterpCoeffs(); + + if (rank == MASTER_NODE) cout << "LUT fluid model ready for use" << endl; +} + +void CLookUpTable::BuildOriginalTrapMap() { + /*--- Build the original SU2 trapezoidal map (DAG-based, O(n log n) memory) ---*/ + /*--- WARNING: This can use a lot of memory for large tables! ---*/ + if (rank == MASTER_NODE) switch (table_dim) { case 2: cout << "Building a trapezoidal map for the (" + name_CV1 + ", " + name_CV2 + - ") " - "space ..." - << endl; + ") space ..." << endl; break; case 3: cout << "Building trapezoidal map stack for the (" + name_CV1 + ", " + name_CV2 + - ") " - "space ..." - << endl; + ") space ..." << endl; break; default: break; @@ -82,11 +89,13 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, unsigned short barwidth = 65; bool display_map_info = (n_table_levels < 2); double tmap_memory_footprint = 0; + for (auto i_level = 0ul; i_level < n_table_levels; i_level++) { trap_map_x_y[i_level] = CTrapezoidalMap(GetDataP(name_CV1, i_level), GetDataP(name_CV2, i_level), table_data[i_level].cols(), edges[i_level], edge_to_triangle[i_level], display_map_info); tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint(); + /* Display a progress bar to monitor table generation process */ if (rank == MASTER_NODE) { su2double progress = su2double(i_level) / n_table_levels; @@ -114,10 +123,72 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, cout << "Trapezoidal map memory footprint: " << tmap_memory_footprint << " MB\n"; cout << "Table data memory footprint: " << memory_footprint_data << " MB\n" << endl; } +} - ComputeInterpCoeffs(); - - if (rank == MASTER_NODE) cout << "LUT fluid model ready for use" << endl; +void CLookUpTable::EnableFastTrapMap() { + /*--- Build Memory-efficient trapezoidal map (band-based, O(n) memory) ---*/ + + use_fast_trap_map_ = true; + + trap_map_x_y_FAST.resize(n_table_levels); + + if (rank == MASTER_NODE) { + cout << endl; + cout << "+--------------------------------------------------+" << endl; + cout << "| Building Memory-efficient trapezoidal maps |" << endl; + cout << "| (O(n) memory, O(log n) queries) |" << endl; + cout << "+--------------------------------------------------+" << endl; + } + + su2double startTime = SU2_MPI::Wtime(); + + for (unsigned long i_level = 0; i_level < n_table_levels; ++i_level) { + + const auto n_pts = n_points[i_level]; + const auto n_tris = n_triangles[i_level]; + + if (rank == MASTER_NODE) { + cout << " Level " << i_level << ": " << n_pts << " points, " + << n_tris << " triangles ... " << flush; + } + + // Get point coordinates + std::vector x_coords(n_pts); + std::vector y_coords(n_pts); + + for (unsigned long i_point = 0; i_point < n_pts; ++i_point) { + x_coords[i_point] = table_data[i_level][idx_CV1][i_point]; + y_coords[i_point] = table_data[i_level][idx_CV2][i_point]; + } + + // Get triangle connectivity + std::vector tri_conn(3 * n_tris); + for (unsigned long i_tri = 0; i_tri < n_tris; ++i_tri) { + tri_conn[3 * i_tri + 0] = triangles[i_level][i_tri][0]; + tri_conn[3 * i_tri + 1] = triangles[i_level][i_tri][1]; + tri_conn[3 * i_tri + 2] = triangles[i_level][i_tri][2]; + } + + // Build FAST trapezoidal map for this level + trap_map_x_y_FAST[i_level].Build(n_pts, n_tris, + x_coords.data(), y_coords.data(), + tri_conn.data()); + + if (rank == MASTER_NODE) { + cout << "done." << endl; + } + } + + su2double stopTime = SU2_MPI::Wtime(); + + if (rank == MASTER_NODE) { + cout << "+--------------------------------------------------+" << endl; + cout << "| FAST trapezoidal maps built successfully! |" << endl; + cout << "| Build time: " << setw(10) << setprecision(3) << (stopTime - startTime) + << " seconds |" << endl; + cout << "+--------------------------------------------------+" << endl; + cout << endl; + } } void CLookUpTable::LoadTableRaw(const string& var_file_name_lut) { @@ -617,17 +688,30 @@ bool CLookUpTable::LookUp_XY(const vector& idx_var, vector= *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[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) */ - return IsInTriangle(val_CV1, val_CV2, id_triangle, iLevel); + /*--- Use Memory-efficient trapezoidal map if enabled (LUT_FAST) ---*/ + if (use_fast_trap_map_ && iLevel < trap_map_x_y_FAST.size()) { + std::array bary_coords; + if (trap_map_x_y_FAST[iLevel].FindTriangle(val_CV1, val_CV2, id_triangle, bary_coords)) { + return true; + } + // Fall through to bounds check if not found (point may be outside domain) + } + + /*--- Use original trapezoidal map if available ---*/ + if (!use_fast_trap_map_ && iLevel < trap_map_x_y.size() && trap_map_x_y.size() > 0) { + /* 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[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[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) */ + return IsInTriangle(val_CV1, val_CV2, id_triangle, iLevel); + } } + return false; } diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index fd1417b3f56d..7fb363dcedd5 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -78,7 +78,23 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati } look_up_table = new CLookUpTable(datadriven_fluid_options.datadriven_filenames[0], table_scalar_names[I_PROGVAR], table_scalar_names[I_ENTH]); + /*--- Build the original trapezoidal map (can use lots of memory!) ---*/ + look_up_table->BuildOriginalTrapMap(); break; + + case ENUM_DATADRIVEN_METHOD::LUT_FAST: + if (rank == MASTER_NODE) { + cout << "**************************************************" << endl; + cout << "*** initializing the FAST lookup table ***" << endl; + cout << "*** (Memory-efficient trapezoidal map) ***" << endl; + cout << "**************************************************" << endl; + } + look_up_table = new CLookUpTable(config->GetDataDriven_FileNames()[0], table_scalar_names[I_PROGVAR], + table_scalar_names[I_ENTH]); + /*--- Build Memory-efficient trapezoidal map (O(n) memory) ---*/ + look_up_table->EnableFastTrapMap(); + break; + default: if (rank == MASTER_NODE) { cout << "***********************************************" << endl; @@ -104,8 +120,10 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati } CFluidFlamelet::~CFluidFlamelet() { - if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT) + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT || + Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT_FAST) delete look_up_table; + #ifdef USE_MLPCPP if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { delete iomap_TD; @@ -113,7 +131,7 @@ CFluidFlamelet::~CFluidFlamelet() { delete iomap_LookUp; delete lookup_mlp; if (preferential_diffusion) delete iomap_PD; - } + } #endif } @@ -214,7 +232,10 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { preferential_diffusion = flamelet_options.preferential_diffusion; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: + case ENUM_DATADRIVEN_METHOD::LUT_FAST: + if (preferential_diffusion) { preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); + } break; case ENUM_DATADRIVEN_METHOD::MLP: #ifdef USE_MLPCPP @@ -320,6 +341,7 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca bool inside; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: + case ENUM_DATADRIVEN_METHOD::LUT_FAST: if (output_refs.size() != LUT_idx.size()) SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); if (include_mixture_fraction) { diff --git a/UnitTests/Common/containers/CLookupTable_tests.cpp b/UnitTests/Common/containers/CLookupTable_tests.cpp index 4b9f06eb9469..c2b1d7abe9b0 100644 --- a/UnitTests/Common/containers/CLookupTable_tests.cpp +++ b/UnitTests/Common/containers/CLookupTable_tests.cpp @@ -35,11 +35,12 @@ #include "../../../Common/include/containers/CLookUpTable.hpp" #include "../../../Common/include/containers/CFileReaderLUT.hpp" +/*--- Original 2D LUT test (using default trapezoidal map) ---*/ TEST_CASE("LUTreader", "[tabulated chemistry]") { /*--- smaller and trivial lookup table ---*/ CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable.drg", "ProgressVariable", "EnthalpyTot"); - + look_up_table.BuildOriginalTrapMap(); /*--- string names of the controlling variables ---*/ string name_CV1 = "ProgressVariable"; @@ -84,11 +85,122 @@ TEST_CASE("LUTreader", "[tabulated chemistry]") { CHECK(look_up_dat == Approx(1.1738796125)); } +/*--- Original 3D LUT test (using default trapezoidal map) ---*/ TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { /*--- smaller and trivial lookup table ---*/ CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable_3D.drg", "ProgressVariable", "EnthalpyTot"); + look_up_table.BuildOriginalTrapMap(); + /*--- string names of the controlling variables ---*/ + + string name_CV1 = "ProgressVariable"; + string name_CV2 = "EnthalpyTot"; + + /*--- look up a single value for density ---*/ + + su2double prog = 0.55; + 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(idx_tag, &look_up_dat, prog, enth, mfrac); + CHECK(look_up_dat == Approx(1.02)); + + /*--- look up a single value for viscosity ---*/ + + prog = 0.6; + enth = 0.9; + mfrac = 0.8; + look_up_tag = "Viscosity"; + 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 */ + + auto limitsEnth = look_up_table.GetTableLimitsY(); + CHECK(SU2_TYPE::GetValue(*limitsEnth.first) == Approx(-1.0)); + CHECK(SU2_TYPE::GetValue(*limitsEnth.second) == Approx(1.0)); + + auto limitsProgvar = look_up_table.GetTableLimitsX(); + CHECK(SU2_TYPE::GetValue(*limitsProgvar.first) == Approx(0.0)); + CHECK(SU2_TYPE::GetValue(*limitsProgvar.second) == Approx(1.0)); + + /* lookup value outside of lookup table */ + + prog = 1.10; + enth = 1.1; + mfrac = 2.0; + look_up_tag = "Density"; + 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)); +} + +/*--- 2D LUT test using Pedro's FAST trapezoidal map (LUT_FAST) ---*/ +TEST_CASE("LUTreader_FAST", "[tabulated chemistry]") { + /*--- smaller and trivial lookup table ---*/ + + CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable.drg", "ProgressVariable", "EnthalpyTot"); + + /*--- Enable Pedro's FAST trapezoidal map ---*/ + look_up_table.EnableFastTrapMap(); + + /*--- string names of the controlling variables ---*/ + + string name_CV1 = "ProgressVariable"; + string name_CV2 = "EnthalpyTot"; + + /*--- look up a single value for density ---*/ + + 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(idx_tag, &look_up_dat, prog, enth); + CHECK(look_up_dat == Approx(1.02)); + + /*--- look up a single value for viscosity ---*/ + + prog = 0.6; + enth = 0.9; + look_up_tag = "Viscosity"; + 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 */ + + auto limitsEnth = look_up_table.GetTableLimitsY(); + CHECK(SU2_TYPE::GetValue(*limitsEnth.first) == Approx(-1.0)); + CHECK(SU2_TYPE::GetValue(*limitsEnth.second) == Approx(1.0)); + + auto limitsProgvar = look_up_table.GetTableLimitsX(); + CHECK(SU2_TYPE::GetValue(*limitsProgvar.first) == Approx(0.0)); + CHECK(SU2_TYPE::GetValue(*limitsProgvar.second) == Approx(1.0)); + + /* lookup value outside of lookup table */ + + prog = 1.10; + enth = 1.1; + look_up_tag = "Density"; + 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)); +} + +/*--- 3D LUT test using Pedro's FAST trapezoidal map (LUT_FAST) ---*/ +TEST_CASE("LUTreader_3D_FAST", "[tabulated chemistry]") { + /*--- smaller and trivial lookup table ---*/ + + CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable_3D.drg", "ProgressVariable", + "EnthalpyTot"); + + /*--- Enable Pedro's FAST trapezoidal map ---*/ + look_up_table.EnableFastTrapMap(); /*--- string names of the controlling variables ---*/ @@ -136,3 +248,4 @@ TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { look_up_table.LookUp_XYZ(idx_tag, &look_up_dat, prog, enth, mfrac); CHECK(look_up_dat == Approx(1.1738796125)); } + From 4b65f7835bfee39930b88e70df104bfe66e6bd53 Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Thu, 11 Dec 2025 14:16:54 +0100 Subject: [PATCH 2/8] LUT_FAST with C++11 compatibility --- Common/include/containers/CLookUpTable.hpp | 1 - Common/include/containers/CTrapezoidalMapFAST.hpp | 13 ++++--------- 2 files changed, 4 insertions(+), 10 deletions(-) diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index 5d83e19bff8c..fb7fdc18c3dd 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -264,7 +264,6 @@ class CLookUpTable { /*! * \brief Build the original SU2 trapezoidal map (DAG-based). - * WARNING: This can use a lot of memory for large tables! * Use EnableFastTrapMap() for large tables instead. */ void BuildOriginalTrapMap(); diff --git a/Common/include/containers/CTrapezoidalMapFAST.hpp b/Common/include/containers/CTrapezoidalMapFAST.hpp index 94e5053b3f56..80864c589a58 100644 --- a/Common/include/containers/CTrapezoidalMapFAST.hpp +++ b/Common/include/containers/CTrapezoidalMapFAST.hpp @@ -4,10 +4,7 @@ * Based on Pedro Gomes' (pcarruscag) LUT implementation: * https://github.com/pcarruscag/LUT * - * This is a faithful adaptation of Pedro's algorithm for SU2's FGM module. - * The algorithm uses band-based spatial indexing with O(log n) query time. - * - * NOTE: This implementation is C++11 compatible for SU2. + * */ #ifndef CTRAPEZOIDALMAP_FAST_HPP @@ -24,7 +21,6 @@ namespace su2_lut { -/*--- Type aliases matching Pedro's implementation ---*/ using IntT = int32_t; using RealT = su2double; @@ -51,7 +47,7 @@ struct TrapezoidalMap { VectorReal x_bands, edge_y; }; -/*--- C++11 compatible comparison functors ---*/ +/*--- Comparison functors ---*/ struct PointComparator { const VectorReal& x; const VectorReal& y; @@ -78,7 +74,6 @@ struct PairSecondComparator { /*! * \brief Orders points by ascending x coordinates and updates triangle indices. - * This is critical - the algorithm assumes points are sorted by x. */ inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { const IntT n_pts = static_cast(x.size()); @@ -179,7 +174,7 @@ inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, } /*! - * \brief Band detection result structure (for C++11 compatibility). + * \brief Band detection result structure. */ struct BandResult { IntT n_bands; @@ -416,7 +411,7 @@ inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, } /*! - * \brief Query result structure (for C++11 compatibility). + * \brief Query result structure. */ struct QueryResult { IntT edge_below; From a6ed7264cedd3e852e956b793317a03b64529b14 Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Thu, 11 Dec 2025 16:18:58 +0100 Subject: [PATCH 3/8] LUT_FAST: memory-efficient trapezoidal map (C++17) - All unit tests passing (4 test cases, 28 assertions) --- .../containers/CTrapezoidalMapFAST.hpp | 497 +++++++----------- 1 file changed, 189 insertions(+), 308 deletions(-) diff --git a/Common/include/containers/CTrapezoidalMapFAST.hpp b/Common/include/containers/CTrapezoidalMapFAST.hpp index 80864c589a58..1c57abecc7bb 100644 --- a/Common/include/containers/CTrapezoidalMapFAST.hpp +++ b/Common/include/containers/CTrapezoidalMapFAST.hpp @@ -3,8 +3,6 @@ * \brief Fast trapezoidal map implementation for 2D LUT queries in SU2. * Based on Pedro Gomes' (pcarruscag) LUT implementation: * https://github.com/pcarruscag/LUT - * - * */ #ifndef CTRAPEZOIDALMAP_FAST_HPP @@ -16,6 +14,7 @@ #include #include #include +#include #include #include @@ -47,31 +46,6 @@ struct TrapezoidalMap { VectorReal x_bands, edge_y; }; -/*--- Comparison functors ---*/ -struct PointComparator { - const VectorReal& x; - const VectorReal& y; - - PointComparator(const VectorReal& x_ref, const VectorReal& y_ref) - : x(x_ref), y(y_ref) {} - - bool operator()(const IntT i, const IntT j) const { - return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; - } -}; - -struct EdgeComparator { - bool operator()(const std::array& a, const std::array& b) const { - return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); - } -}; - -struct PairSecondComparator { - bool operator()(const std::pair& a, const std::pair& b) const { - return a.second < b.second; - } -}; - /*! * \brief Orders points by ascending x coordinates and updates triangle indices. */ @@ -80,19 +54,21 @@ inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { // Create sort permutation based on (x, y) coordinates std::vector perm(n_pts); - for (IntT i = 0; i < n_pts; ++i) { - perm[i] = i; - } - std::sort(perm.begin(), perm.end(), PointComparator(x, y)); + std::iota(perm.begin(), perm.end(), 0); + std::sort(perm.begin(), perm.end(), [&x, &y](const auto i, const auto j) { + return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; + }); // Reorder coordinate arrays - VectorReal x_tmp(n_pts), y_tmp(n_pts); - for (IntT i = 0; i < n_pts; ++i) { - x_tmp[i] = x[perm[i]]; - y_tmp[i] = y[perm[i]]; - } - x = x_tmp; - y = y_tmp; + auto reorder = [n_pts, &perm](const auto& v) { + VectorReal tmp(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + tmp[i] = v[perm[i]]; + } + return tmp; + }; + x = reorder(x); + y = reorder(y); // Build inverse permutation to update triangle indices std::vector inv_perm(n_pts); @@ -100,7 +76,7 @@ inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { inv_perm[perm[i]] = i; } - // Update all triangle vertex indices to use new point ordering + // Update triangle connectivity for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { for (IntT j = 0; j < 3; ++j) { triangles(i, j) = inv_perm[triangles(i, j)]; @@ -109,21 +85,19 @@ inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { } /*! - * \brief Helper to check if two edges are equal (same endpoints). + * \brief Check if two edges are equal (same point IDs). */ inline bool EdgesEqual(const std::array& a, const std::array& b) { return a[0] == b[0] && a[1] == b[1]; } /*! - * \brief Extracts unique edges from triangles with face adjacency information. - * Each edge stores which triangles (up to 2) share it. - * Boundary edges have edge_faces(i,1) = -1. + * \brief Extract unique edges from triangles. */ inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, Matrix2i& edge_faces) { - // Extract all edges (3 per triangle), storing (min_pt, max_pt, triangle_id) - std::vector > edges; + // Extract edges from triangles + std::vector> edges; edges.resize(3 * triangles.rows()); for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { @@ -131,13 +105,14 @@ inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, const IntT j = (i + 1) % 3; const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); - std::array edge_data = {{i_pt, j_pt, i_tri}}; - edges[3 * i_tri + i] = edge_data; + edges[3 * i_tri + i] = {i_pt, j_pt, i_tri}; } } - // Sort to identify duplicate edges (shared between triangles) - std::sort(edges.begin(), edges.end(), EdgeComparator()); + // Sort to identify duplicates + std::sort(edges.begin(), edges.end(), [](const auto& a, const auto& b) { + return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); + }); // Count unique edges IntT n_edges = 1; @@ -145,54 +120,38 @@ inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, n_edges += static_cast(!EdgesEqual(edges[i], edges[i - 1])); } - // Build edge_pts and edge_faces arrays + // Map edge points and edge faces edge_pts.resize(n_edges, 2); edge_faces.resize(n_edges, 2); - IntT pos = 0; - - // First edge - edge_pts(pos, 0) = edges[0][0]; - edge_pts(pos, 1) = edges[0][1]; - edge_faces(pos, 0) = edges[0][2]; - edge_faces(pos, 1) = -1; - ++pos; - + + auto new_edge = [&](const auto& edge) { + edge_pts(pos, 0) = edge[0]; + edge_pts(pos, 1) = edge[1]; + edge_faces(pos, 0) = edge[2]; + edge_faces(pos, 1) = -1; + ++pos; + }; + + new_edge(edges[0]); for (IntT i = 1; i < static_cast(edges.size()); ++i) { if (EdgesEqual(edges[i], edges[i - 1])) { - // Duplicate edge - record second adjacent triangle edge_faces(pos - 1, 1) = edges[i][2]; } else { - // New edge - edge_pts(pos, 0) = edges[i][0]; - edge_pts(pos, 1) = edges[i][1]; - edge_faces(pos, 0) = edges[i][2]; - edge_faces(pos, 1) = -1; - ++pos; + new_edge(edges[i]); } } } -/*! - * \brief Band detection result structure. - */ -struct BandResult { - IntT n_bands; - std::vector pt_to_band; - VectorReal x_bands; -}; - /*! * \brief Detects bands based on unique x-coordinates. * Uses coarse banding to limit memory for large tables. + * Returns tuple of (n_bands, pt_to_band, x_bands). */ -inline BandResult DetectBands(const VectorReal& x, IntT max_bands = 0) { - BandResult result; - +inline auto DetectBands(const VectorReal& x, IntT max_bands = 0) { // Safety check: empty input if (x.empty()) { - result.n_bands = 0; - return result; + return std::make_tuple(IntT{0}, std::vector{}, VectorReal{}); } // First pass: count unique x-coordinates @@ -202,87 +161,80 @@ inline BandResult DetectBands(const VectorReal& x, IntT max_bands = 0) { } // Determine if we need to use coarser bands - // Default max_bands: limit to sqrt(n_points) * 4, capped at 5000 if (max_bands <= 0) { - max_bands = std::min(IntT(5000), static_cast(4.0 * std::sqrt(static_cast(x.size())))); + max_bands = std::min(IntT{5000}, static_cast(4.0 * std::sqrt(static_cast(x.size())))); } const bool use_coarse_bands = (n_unique > max_bands); if (!use_coarse_bands) { // Original algorithm: use all unique x-coordinates as band boundaries - result.pt_to_band.resize(x.size()); - result.pt_to_band[0] = 0; - result.n_bands = 0; + std::vector pt_to_band(x.size()); + pt_to_band[0] = 0; + IntT n_bands = 0; for (IntT i = 1; i < static_cast(x.size()); ++i) { - result.n_bands += static_cast(x[i] != x[i - 1]); - result.pt_to_band[i] = result.n_bands; + n_bands += static_cast(x[i] != x[i - 1]); + pt_to_band[i] = n_bands; } - result.x_bands.resize(result.n_bands + 1); + VectorReal x_bands(n_bands + 1); IntT pos = 0; - result.x_bands[pos] = x[0]; + x_bands[pos] = x[0]; for (IntT i = 1; i < static_cast(x.size()); ++i) { - if (x[i] != result.x_bands[pos]) { - result.x_bands[++pos] = x[i]; + if (x[i] != x_bands[pos]) { + x_bands[++pos] = x[i]; } } + + return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); } else { // Coarse banding: divide x-range into max_bands equal intervals - const RealT x_min = x.front(); // x is sorted + const RealT x_min = x.front(); const RealT x_max = x.back(); const RealT x_range = x_max - x_min; // Handle degenerate case where all x are the same if (x_range < 1e-30) { - result.n_bands = 1; - result.pt_to_band.resize(x.size(), 0); - result.x_bands.resize(2); - result.x_bands[0] = x_min; - result.x_bands[1] = x_max; - return result; + std::vector pt_to_band(x.size(), 0); + VectorReal x_bands = {x_min, x_max}; + return std::make_tuple(IntT{1}, std::move(pt_to_band), std::move(x_bands)); } - result.n_bands = max_bands; + const IntT n_bands = max_bands; const RealT band_width = x_range / max_bands; // Build band boundaries - result.x_bands.resize(max_bands + 1); + VectorReal x_bands(max_bands + 1); for (IntT i = 0; i <= max_bands; ++i) { - result.x_bands[i] = x_min + i * band_width; + x_bands[i] = x_min + i * band_width; } - // Ensure last band includes x_max exactly - result.x_bands[max_bands] = x_max; + x_bands[max_bands] = x_max; // Map each point to its band - result.pt_to_band.resize(x.size()); + std::vector pt_to_band(x.size()); for (IntT i = 0; i < static_cast(x.size()); ++i) { - // Calculate band index: floor((x - x_min) / band_width) IntT band = static_cast((x[i] - x_min) / band_width); - // Clamp to valid range [0, n_bands-1] - band = std::max(IntT(0), std::min(band, max_bands - 1)); - result.pt_to_band[i] = band; + band = std::max(IntT{0}, std::min(band, max_bands - 1)); + pt_to_band[i] = band; } + + return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); } - - return result; } /*! * \brief Builds the trapezoidal map for efficient 2D queries. - * For each band, stores edges crossing it, sorted by y-coordinate at band midpoint. */ inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, const VectorReal& y, TrapezoidalMap& map) { - VectorReal& x_bands = map.x_bands; - VectorInt& offsets = map.offsets; - VectorInt& edge_id = map.edge_id; - VectorReal& edge_y = map.edge_y; - - BandResult band_result = DetectBands(x); - const IntT n_bands = band_result.n_bands; - x_bands = band_result.x_bands; + auto& x_bands = map.x_bands; + auto& offsets = map.offsets; + auto& edge_id = map.edge_id; + auto& edge_y = map.edge_y; + + const auto [n_bands, pt_to_band, bands] = DetectBands(x); + x_bands = std::move(bands); // Safety check if (n_bands <= 0 || x_bands.size() < 2) { @@ -294,21 +246,16 @@ inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, } // Helper function to find band index for an x-coordinate - // Uses binary search on x_bands auto find_band = [&x_bands, n_bands](RealT x_val) -> IntT { - // Binary search for the band containing x_val - VectorReal::const_iterator it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); IntT idx = static_cast(it - x_bands.begin()); - // lower_bound returns iterator to first element >= x_val - // We want the band to the left, so subtract 1 (but not below 0) - idx = std::max(IntT(0), idx - 1); - // Also cap at n_bands - 1 + idx = std::max(IntT{0}, idx - 1); idx = std::min(idx, n_bands - 1); return idx; }; - // Count edges per band - use actual x-coordinates to determine band range - VectorInt& counts = offsets; + // Count edges per band + auto& counts = offsets; counts.clear(); counts.resize(n_bands + 1, 0); @@ -318,15 +265,12 @@ inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, const RealT x_0 = x[pt_0]; const RealT x_1 = x[pt_1]; - // Find bands for each endpoint using actual x-coordinates const IntT band_0 = find_band(x_0); const IntT band_1 = find_band(x_1); - // Edge spans from min_band to max_band const IntT min_band = std::min(band_0, band_1); const IntT max_band = std::max(band_0, band_1); - // Count edge in each band it crosses for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { ++counts[j + 1]; } @@ -337,15 +281,11 @@ inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, offsets[i] += offsets[i - 1]; } - // Check total memory requirement before allocation + // Check memory requirement const size_t total_entries = static_cast(offsets.back()); - const size_t memory_bytes = total_entries * (sizeof(IntT) + sizeof(RealT)); - const size_t memory_mb = memory_bytes / (1024 * 1024); + const size_t memory_mb = total_entries * (sizeof(IntT) + sizeof(RealT)) / (1024 * 1024); - // Memory limit: 2GB for edge storage (adjustable) - const size_t MAX_MEMORY_MB = 2048; - if (memory_mb > MAX_MEMORY_MB) { - // Too much memory - clear and return empty map + if (memory_mb > 2048) { x_bands.clear(); offsets.clear(); edge_id.clear(); @@ -353,100 +293,83 @@ inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, return; } - // Fill edge data for each band - edge_id.resize(total_entries); - edge_y.resize(total_entries); - VectorInt pos = offsets; // Working copy of offsets + // Allocate storage + edge_id.resize(offsets.back()); + edge_y.resize(offsets.back()); + auto pos = offsets; + // Store edges in each band for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { const IntT pt_0 = edge_pts(i_edge, 0); const IntT pt_1 = edge_pts(i_edge, 1); - const RealT x_0 = x[pt_0]; - const RealT x_1 = x[pt_1]; - const RealT y_0 = y[pt_0]; - const RealT y_1 = y[pt_1]; + const RealT x_0 = x[pt_0], y_0 = y[pt_0]; + const RealT x_1 = x[pt_1], y_1 = y[pt_1]; - // Find bands for each endpoint const IntT band_0 = find_band(x_0); const IntT band_1 = find_band(x_1); + + if (band_0 == band_1) continue; + const IntT min_band = std::min(band_0, band_1); const IntT max_band = std::max(band_0, band_1); + const RealT dy_dx = (y_1 - y_0) / (x_1 - x_0); - // Calculate slope (handle vertical edges) - const RealT dx = x_1 - x_0; - const bool is_vertical = (std::abs(dx) < 1e-30); - const RealT dy_dx = is_vertical ? 0.0 : (y_1 - y_0) / dx; - - // Add edge to each band it crosses - for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { - edge_id[pos[j]] = i_edge; - const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2.0; - // Calculate y at band midpoint - if (is_vertical) { - edge_y[pos[j]] = (y_0 + y_1) / 2.0; // Use average y for vertical edges - } else { + for (IntT j = min_band; j < max_band && j < n_bands; ++j) { + if (pos[j] < static_cast(edge_id.size())) { + edge_id[pos[j]] = i_edge; + const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2; edge_y[pos[j]] = y_0 + dy_dx * (x_mid - x_0); + ++pos[j]; } - ++pos[j]; } } - // Sort edges within each band by y-coordinate - std::vector > tmp; + // Sort edges in each band by y coordinate + std::vector> tmp; for (IntT i = 0; i < n_bands; ++i) { const IntT begin = offsets[i]; const IntT end = offsets[i + 1]; - tmp.resize(end - begin); + if (begin >= end) continue; - for (IntT k = begin; k < end; ++k) { - tmp[k - begin] = std::make_pair(edge_id[k], edge_y[k]); + tmp.resize(end - begin); + for (auto k = begin; k < end; ++k) { + tmp[k - begin] = {edge_id[k], edge_y[k]}; } - std::sort(tmp.begin(), tmp.end(), PairSecondComparator()); - - for (IntT k = begin; k < end; ++k) { + std::sort(tmp.begin(), tmp.end(), [](const auto& a, const auto& b) { + return a.second < b.second; + }); + for (auto k = begin; k < end; ++k) { edge_id[k] = tmp[k - begin].first; edge_y[k] = tmp[k - begin].second; } } } -/*! - * \brief Query result structure. - */ -struct QueryResult { - IntT edge_below; - IntT edge_above; -}; - /*! * \brief Query the trapezoidal map for edges bounding a point. - * Returns edges below and above the point. Either can be -1 if point is at boundary. - * With coarse banding, recomputes edge y-positions at the actual query x. - * Also searches adjacent bands to avoid missing candidates. + * Returns pair of (edge_below, edge_above). Either can be -1 if at boundary. */ -inline QueryResult QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, - const VectorReal& x_coords, const VectorReal& y_coords, - const RealT& x, const RealT& y) { - QueryResult result; - result.edge_below = -1; - result.edge_above = -1; - - // Safety check: empty map +inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, + const VectorReal& x_coords, const VectorReal& y_coords, + const RealT& x, const RealT& y) { + // Safety check if (map.x_bands.size() < 2 || map.offsets.empty()) { - return result; + return std::make_pair(IntT{-1}, IntT{-1}); } const IntT n_bands = static_cast(map.x_bands.size()) - 1; // Find the band containing x - const VectorReal& x_bands = map.x_bands; - VectorReal::const_iterator it = std::lower_bound(x_bands.begin(), x_bands.end(), x); + const auto& x_bands = map.x_bands; + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x); IntT d = static_cast(it - x_bands.begin()); - IntT center_band = std::min(std::max(IntT(0), d - 1), n_bands - 1); + IntT center_band = std::min(std::max(IntT{0}, d - 1), n_bands - 1); - // Search in current band AND adjacent bands to catch edges near boundaries + // Search in current band AND adjacent bands RealT best_y_below = -1e300; RealT best_y_above = 1e300; + IntT edge_below = -1; + IntT edge_above = -1; for (IntT band_offset = -1; band_offset <= 1; ++band_offset) { IntT band_idx = center_band + band_offset; @@ -463,152 +386,132 @@ inline QueryResult QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i for (IntT k = begin; k < end; ++k) { const IntT e_id = map.edge_id[k]; - // Get edge endpoints const IntT p0 = edge_pts(e_id, 0); const IntT p1 = edge_pts(e_id, 1); const RealT x0 = x_coords[p0], y0 = y_coords[p0]; const RealT x1 = x_coords[p1], y1 = y_coords[p1]; - // Check if query x is within the edge's x-extent - const RealT x_min_edge = std::min(x0, x1); - const RealT x_max_edge = std::max(x0, x1); - if (x < x_min_edge - 1e-10 || x > x_max_edge + 1e-10) { - continue; // Query x is outside this edge's range + // Check if query x is within edge's x-extent + if (x < std::min(x0, x1) - 1e-10 || x > std::max(x0, x1) + 1e-10) { + continue; } // Compute y-position of edge at query x RealT edge_y_at_x; const RealT dx = x1 - x0; if (std::abs(dx) < 1e-30) { - // Vertical edge - use average y edge_y_at_x = (y0 + y1) / 2.0; } else { - // Linear interpolation const RealT t = (x - x0) / dx; edge_y_at_x = y0 + t * (y1 - y0); } - // Check if this edge is below or above the query point - if (edge_y_at_x <= y + 1e-10) { - if (edge_y_at_x > best_y_below) { - best_y_below = edge_y_at_x; - result.edge_below = e_id; - } + if (edge_y_at_x <= y + 1e-10 && edge_y_at_x > best_y_below) { + best_y_below = edge_y_at_x; + edge_below = e_id; } - if (edge_y_at_x >= y - 1e-10) { - if (edge_y_at_x < best_y_above) { - best_y_above = edge_y_at_x; - result.edge_above = e_id; - } + if (edge_y_at_x >= y - 1e-10 && edge_y_at_x < best_y_above) { + best_y_above = edge_y_at_x; + edge_above = e_id; } } } - return result; + return std::make_pair(edge_below, edge_above); } /*! * \brief Get triangles adjacent to two query edges (up to 3 triangles). */ -inline std::array AdjacentTriangles(const IntT edge_0, const IntT edge_1, - const Matrix2i& edge_faces) { - std::array tris = {{-1, -1, -1}}; +inline auto AdjacentTriangles(const IntT edge_0, const IntT edge_1, + const Matrix2i& edge_faces) { + std::array tris = {-1, -1, -1}; IntT pos = 0; - // Helper to insert unique triangle - if (edge_0 >= 0) { - IntT t0 = edge_faces(edge_0, 0); - IntT t1 = edge_faces(edge_0, 1); - if (t0 >= 0) { - tris[pos++] = t0; - } - if (t1 >= 0) { - bool found = false; - for (IntT k = 0; k < pos; ++k) { - if (t1 == tris[k]) { found = true; break; } - } - if (!found) tris[pos++] = t1; - } - } - - if (edge_1 >= 0) { - IntT t0 = edge_faces(edge_1, 0); - IntT t1 = edge_faces(edge_1, 1); - if (t0 >= 0) { - bool found = false; - for (IntT k = 0; k < pos; ++k) { - if (t0 == tris[k]) { found = true; break; } - } - if (!found && pos < 3) tris[pos++] = t0; + auto insert = [&tris, &pos](const IntT t) { + if (t < 0) return; + for (IntT i = 0; i < pos; ++i) { + if (t == tris[i]) return; } - if (t1 >= 0) { - bool found = false; - for (IntT k = 0; k < pos; ++k) { - if (t1 == tris[k]) { found = true; break; } - } - if (!found && pos < 3) tris[pos++] = t1; + tris[pos++] = t; + }; + + auto get_tris = [&edge_faces](const IntT e) { + if (e < 0) return std::array{IntT{-1}, IntT{-1}}; + return std::array{edge_faces(e, 0), edge_faces(e, 1)}; + }; + + for (const auto e : {edge_0, edge_1}) { + for (const auto t : get_tris(e)) { + insert(t); } } - return tris; } /*! - * \brief Compute barycentric coordinates and check if point is inside triangle. + * \brief Compute barycentric coordinates of point (x, y) in triangle. */ -inline std::array TriangleCoords(const IntT i_tri, const Matrix3i& triangles, - const VectorReal& x, const VectorReal& y, - const RealT x_q, const RealT y_q) { - const IntT i0 = triangles(i_tri, 0); - const IntT i1 = triangles(i_tri, 1); - const IntT i2 = triangles(i_tri, 2); - - const RealT x0 = x[i0], y0 = y[i0]; - const RealT x1 = x[i1], y1 = y[i1]; - const RealT x2 = x[i2], y2 = y[i2]; - - const RealT det = (y1 - y2) * (x0 - x2) + (x2 - x1) * (y0 - y2); - const RealT l0 = ((y1 - y2) * (x_q - x2) + (x2 - x1) * (y_q - y2)) / det; - const RealT l1 = ((y2 - y0) * (x_q - x2) + (x0 - x2) * (y_q - y2)) / det; - const RealT l2 = 1.0 - l0 - l1; - - std::array coords = {{l0, l1, l2}}; - return coords; +inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, + const VectorReal& x, const VectorReal& y, + const RealT x_q, const RealT y_q) { + const IntT p0 = triangles(i_tri, 0); + const IntT p1 = triangles(i_tri, 1); + const IntT p2 = triangles(i_tri, 2); + + const RealT x0 = x[p0], y0 = y[p0]; + const RealT x1 = x[p1], y1 = y[p1]; + const RealT x2 = x[p2], y2 = y[p2]; + + const RealT dx1 = x1 - x0, dy1 = y1 - y0; + const RealT dx2 = x2 - x0, dy2 = y2 - y0; + + auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { + return ux * vy - uy * vx; + }; + + const RealT det = cross(dx1, dy1, dx2, dy2); + if (std::abs(det) < 1e-30) { + return std::array{RealT{0}, RealT{0}, RealT{0}}; + } + + const RealT inv_det = 1.0 / det; + const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; + const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; + + return std::array{1 - a - b, a, b}; } /*! - * \brief Check if point is inside triangle (barycentric coords all >= 0). + * \brief Check if point is inside triangle. */ inline bool InTriangle(const std::array& coords, const RealT tol = 0.0) { return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; } /*! - * \brief Find triangle containing point (x_q, y_q) using the trapezoidal map. - * Returns triangle index and barycentric coordinates. - * Returns -1 if point is outside the mesh. + * \brief Find triangle containing point using the trapezoidal map. */ inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, const Matrix2i& edge_pts, const Matrix2i& edge_faces, const VectorReal& x, const VectorReal& y, const RealT x_q, const RealT y_q, std::array& bary_out) { - QueryResult qr = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); - std::array candidates = AdjacentTriangles(qr.edge_below, qr.edge_above, edge_faces); + const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); + const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); - const RealT tol = 1e-12; - for (IntT i = 0; i < 3; ++i) { - if (candidates[i] < 0) continue; + constexpr RealT tol = 1e-12; + for (const auto t : candidates) { + if (t < 0) continue; - std::array coords = TriangleCoords(candidates[i], triangles, x, y, x_q, y_q); + const auto coords = TriangleCoords(t, triangles, x, y, x_q, y_q); if (InTriangle(coords, tol)) { bary_out = coords; - return candidates[i]; + return t; } } - // Not found - return -1 - bary_out[0] = bary_out[1] = bary_out[2] = 0.0; + bary_out = {0.0, 0.0, 0.0}; return -1; } @@ -617,7 +520,7 @@ inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, /*! * \class CTrapezoidalMapFAST - * \brief Wrapper class for SU2 integration of Pedro's trapezoidal map algorithm. + * \brief Wrapper class for SU2 integration. */ class CTrapezoidalMapFAST { private: @@ -634,13 +537,7 @@ class CTrapezoidalMapFAST { CTrapezoidalMapFAST() : n_points(0), n_triangles(0), build_successful(false) {} /*! - * \brief Build interface matching CLookUpTable.cpp's expected signature. - * \param[in] num_points - Number of mesh points - * \param[in] num_triangles - Number of triangles - * \param[in] x - X coordinates of all points - * \param[in] y - Y coordinates of all points - * \param[in] connectivity - Triangle connectivity (flat array: tri0_p0, tri0_p1, tri0_p2, ...) - * \return true if map was built successfully, false if memory limit exceeded + * \brief Build the trapezoidal map from mesh data. */ bool Build(unsigned long num_points, unsigned long num_triangles, const su2double* x, const su2double* y, @@ -650,7 +547,6 @@ class CTrapezoidalMapFAST { n_triangles = num_triangles; build_successful = false; - // Safety check: handle empty input if (num_points == 0 || num_triangles == 0) { x_coords.clear(); y_coords.clear(); @@ -668,16 +564,14 @@ class CTrapezoidalMapFAST { x_coords.assign(x, x + num_points); y_coords.assign(y, y + num_points); - // Copy triangle connectivity with bounds checking + // Copy triangle connectivity triangles.resize(num_triangles, 3); for (size_t i = 0; i < num_triangles; ++i) { unsigned long v0 = connectivity[3*i + 0]; unsigned long v1 = connectivity[3*i + 1]; unsigned long v2 = connectivity[3*i + 2]; - // Bounds check if (v0 >= num_points || v1 >= num_points || v2 >= num_points) { - // Invalid triangle - skip by setting to valid indices triangles(i, 0) = 0; triangles(i, 1) = 0; triangles(i, 2) = 0; @@ -688,12 +582,11 @@ class CTrapezoidalMapFAST { } } - // Build the map using Pedro's algorithm + // Build the map su2_lut::ReorderPoints(triangles, x_coords, y_coords); su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); - // Check if map was successfully built (not empty due to memory limit) if (map.x_bands.empty() || map.offsets.empty()) { build_successful = false; return false; @@ -703,26 +596,17 @@ class CTrapezoidalMapFAST { return true; } - /*! - * \brief Check if the map was successfully built. - */ bool IsBuilt() const { return build_successful; } /*! * \brief Find the triangle containing a query point. - * \param[in] val_x - X coordinate of query point - * \param[in] val_y - Y coordinate of query point - * \param[out] triangle_id - ID of the containing triangle - * \param[out] bary_coords - Barycentric coordinates of query point in triangle - * \return true if point is inside the mesh, false otherwise */ bool FindTriangle(su2double val_x, su2double val_y, unsigned long& triangle_id, std::array& bary_coords) const { - // Safety check: empty map if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { - bary_coords[0] = bary_coords[1] = bary_coords[2] = 0.0; + bary_coords = {0.0, 0.0, 0.0}; return false; } @@ -732,15 +616,12 @@ class CTrapezoidalMapFAST { if (tri_id >= 0) { triangle_id = static_cast(tri_id); - bary_coords[0] = bary[0]; - bary_coords[1] = bary[1]; - bary_coords[2] = bary[2]; + bary_coords = {bary[0], bary[1], bary[2]}; return true; } return false; } - /*--- Accessors ---*/ unsigned long GetNTriangles() const { return n_triangles; } unsigned long GetNPoints() const { return n_points; } unsigned long GetNEdges() const { return edge_pts.rows(); } From 403d7911074e1910dc00ff4070c70615088128fa Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Tue, 6 Jan 2026 14:31:12 +0100 Subject: [PATCH 4/8] Merged with develop version --- SU2_CFD/src/fluid/CFluidFlamelet.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 7c582ed1b12a..901b4a80aa88 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -89,7 +89,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati cout << "*** (Memory-efficient trapezoidal map) ***" << endl; cout << "**************************************************" << endl; } - look_up_table = new CLookUpTable(config->GetDataDriven_FileNames()[0], table_scalar_names[I_PROGVAR], + look_up_table = new CLookUpTable(datadriven_fluid_options.datadriven_filenames[0], table_scalar_names[I_PROGVAR], table_scalar_names[I_ENTH]); /*--- Build Memory-efficient trapezoidal map (O(n) memory) ---*/ look_up_table->EnableFastTrapMap(); From a75c8bb67932e244f9fd03f2885621314c19a36e Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Tue, 6 Jan 2026 16:13:58 +0100 Subject: [PATCH 5/8] Restore Doxygen comments and add FAST trapezoidal map support --- Common/include/containers/CLookUpTable.hpp | 115 ++++++++++++++++++++- 1 file changed, 113 insertions(+), 2 deletions(-) diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index fb7fdc18c3dd..2458da32f424 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -64,7 +64,6 @@ class CLookUpTable { std::vector z_values_levels; /*!< \brief Constant z-values of each table level.*/ unsigned short table_dim = 2; /*!< \brief Table dimension.*/ - /*! * \brief The lower and upper limits of the z, y and x variable for each table level. */ @@ -194,12 +193,24 @@ class CLookUpTable { /*! * \brief Perform linear interpolation between two table levels for a single variable. + * \param[in] val_CV3 - Value of the third controlling variable at the query point. + * \param[in] lower_level - Table level index of the table level directly below the query point. + * \param[in] upper_level - Table level index of the table level directly above the query point. + * \param[in] lower_value - Result from x-y interpolation on the lower table level. + * \param[in] upper_value - Result from x-y interpolation on the upper table level. + * \param[in] val_vars - Pointer to interpolation result. */ void Linear_Interpolation(const su2double val_CV3, const unsigned long lower_level, const unsigned long upper_level, su2double& lower_value, su2double& upper_value, su2double*& var_vals) const; /*! * \brief Perform linear interpolation between two table levels for a vector of variables. + * \param[in] val_CV3 - Value of the third controlling variable at the query point. + * \param[in] lower_level - Table level index of the table level directly below the query point. + * \param[in] upper_level - Table level index of the table level directly above the query point. + * \param[in] lower_values - Results from x-y interpolation on the lower table level. + * \param[in] upper_values - Results from x-y interpolation on the upper table level. + * \param[in] val_vars - Pointer to interpolation results for all interpolation variables. */ void Linear_Interpolation(const su2double val_CV3, const unsigned long lower_level, const unsigned long upper_level, std::vector& lower_values, std::vector& upper_values, @@ -207,34 +218,73 @@ class CLookUpTable { /*! * \brief Find the point on the hull (boundary of the table) that is closest to the point P(val_CV1,val_CV2). + * \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. + * \returns Point id of the nearest neighbor on the hull. */ unsigned long FindNearestNeighborOnHull(su2double val_CV1, su2double val_CV2, 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] 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. */ void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, 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. + * \param[in] val_CV2 - Second coordinate of point P(val_CV1,val_CV2) to check. + * \param[in] name_var - Variable name to look up. + * \param[out] var_val - Pointer to the to be interpolated variable. + */ void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, const std::string& name_var, su2double* var_val, const unsigned long i_level = 0); /*! * \brief Determine if a point P(val_CV1,val_CV2) is inside the triangle val_id_triangle. + * \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] val_id_triangle - ID of the triangle to check. + * \returns True if the point is in the triangle, false if it is outside. */ bool IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, unsigned long i_level = 0); /*! * \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); @@ -242,6 +292,11 @@ class CLookUpTable { /*! * \brief Compute the values of the first and second controlling variable based on normalized query coordinates + * \param[in] inclusion_levels - Pair containing lower(first) and upper(second) table inclusion level indices. + * \param[in] val_CV1 - First controlling variable value. + * \param[in] val_CV2 - Second controlling variable value. + * \param[in] val_CV3 - Third controlling variable value. + * \returns Array with first and second controlling variable values for the upper and lower table inclusion levels. */ std::array, 2> ComputeNormalizedXY(std::pair& inclusion_levels, const su2double val_CV1, const su2double val_CV2, @@ -249,12 +304,21 @@ class CLookUpTable { /*! * \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); @@ -286,32 +350,76 @@ class CLookUpTable { /*! * \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 whether query is inside (true) or outside (false) data set. */ 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 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 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, const su2double val_CV2, const unsigned long i_level = 0); + /*! + * \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. + * \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. + */ 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 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. + */ 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. + * \param[in] val_CV3 - Value of controlling variable 3. + * \returns Pair of inclusion level indices (first = lower level index, second = upper level index). */ std::pair FindInclusionLevels(const su2double val_CV3); /*! * \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(const unsigned long i_level = 0) const { return limits_table_y[i_level]; @@ -319,6 +427,7 @@ 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(const unsigned long i_level = 0) const { return limits_table_x[i_level]; @@ -331,6 +440,8 @@ class CLookUpTable { /*! * \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(); From 6e14aa545f380c7095a16c3325e60e16877f21d4 Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Tue, 6 Jan 2026 16:36:15 +0100 Subject: [PATCH 6/8] SU2 style changes --- .../containers/CTrapezoidalMapFAST.hpp | 983 +++++++++--------- Common/src/containers/CLookUpTable.cpp | 57 +- SU2_CFD/src/fluid/CFluidFlamelet.cpp | 44 +- .../Common/containers/CLookupTable_tests.cpp | 1 - 4 files changed, 532 insertions(+), 553 deletions(-) diff --git a/Common/include/containers/CTrapezoidalMapFAST.hpp b/Common/include/containers/CTrapezoidalMapFAST.hpp index 1c57abecc7bb..55e26d4ac4eb 100644 --- a/Common/include/containers/CTrapezoidalMapFAST.hpp +++ b/Common/include/containers/CTrapezoidalMapFAST.hpp @@ -26,13 +26,13 @@ using RealT = su2double; /*--- Simple Matrix container (row-major) ---*/ template struct Matrix { - std::vector data; + std::vector data; - void resize(size_t rows, size_t) { data.resize(rows * N); } - size_t rows() const { return data.size() / N; } + void resize(size_t rows, size_t) { data.resize(rows * N); } + size_t rows() const { return data.size() / N; } - const T& operator()(size_t i, size_t j) const { return data[i * N + j]; } - T& operator()(size_t i, size_t j) { return data[i * N + j]; } + const T& operator()(size_t i, size_t j) const { return data[i * N + j]; } + T& operator()(size_t i, size_t j) { return data[i * N + j]; } }; using Matrix2i = Matrix; @@ -42,105 +42,102 @@ using VectorReal = std::vector; /*--- Trapezoidal Map structure ---*/ struct TrapezoidalMap { - VectorInt offsets, edge_id; - VectorReal x_bands, edge_y; + VectorInt offsets, edge_id; + VectorReal x_bands, edge_y; }; /*! * \brief Orders points by ascending x coordinates and updates triangle indices. */ inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { - const IntT n_pts = static_cast(x.size()); - - // Create sort permutation based on (x, y) coordinates - std::vector perm(n_pts); - std::iota(perm.begin(), perm.end(), 0); - std::sort(perm.begin(), perm.end(), [&x, &y](const auto i, const auto j) { - return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; - }); - - // Reorder coordinate arrays - auto reorder = [n_pts, &perm](const auto& v) { - VectorReal tmp(n_pts); - for (IntT i = 0; i < n_pts; ++i) { - tmp[i] = v[perm[i]]; - } - return tmp; - }; - x = reorder(x); - y = reorder(y); - - // Build inverse permutation to update triangle indices - std::vector inv_perm(n_pts); + const IntT n_pts = static_cast(x.size()); + + // Create sort permutation based on (x, y) coordinates + std::vector perm(n_pts); + std::iota(perm.begin(), perm.end(), 0); + std::sort(perm.begin(), perm.end(), + [&x, &y](const auto i, const auto j) { return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; }); + + // Reorder coordinate arrays + auto reorder = [n_pts, &perm](const auto& v) { + VectorReal tmp(n_pts); for (IntT i = 0; i < n_pts; ++i) { - inv_perm[perm[i]] = i; + tmp[i] = v[perm[i]]; } - - // Update triangle connectivity - for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { - for (IntT j = 0; j < 3; ++j) { - triangles(i, j) = inv_perm[triangles(i, j)]; - } + return tmp; + }; + x = reorder(x); + y = reorder(y); + + // Build inverse permutation to update triangle indices + std::vector inv_perm(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + inv_perm[perm[i]] = i; + } + + // Update triangle connectivity + for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { + for (IntT j = 0; j < 3; ++j) { + triangles(i, j) = inv_perm[triangles(i, j)]; } + } } /*! * \brief Check if two edges are equal (same point IDs). */ inline bool EdgesEqual(const std::array& a, const std::array& b) { - return a[0] == b[0] && a[1] == b[1]; + return a[0] == b[0] && a[1] == b[1]; } /*! * \brief Extract unique edges from triangles. */ -inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, - Matrix2i& edge_faces) { - // Extract edges from triangles - std::vector> edges; - edges.resize(3 * triangles.rows()); - - for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { - for (IntT i = 0; i < 3; ++i) { - const IntT j = (i + 1) % 3; - const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); - const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); - edges[3 * i_tri + i] = {i_pt, j_pt, i_tri}; - } +inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, Matrix2i& edge_faces) { + // Extract edges from triangles + std::vector> edges; + edges.resize(3 * triangles.rows()); + + for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { + for (IntT i = 0; i < 3; ++i) { + const IntT j = (i + 1) % 3; + const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); + const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); + edges[3 * i_tri + i] = {i_pt, j_pt, i_tri}; } - - // Sort to identify duplicates - std::sort(edges.begin(), edges.end(), [](const auto& a, const auto& b) { - return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); - }); - - // Count unique edges - IntT n_edges = 1; - for (IntT i = 1; i < static_cast(edges.size()); ++i) { - n_edges += static_cast(!EdgesEqual(edges[i], edges[i - 1])); - } - - // Map edge points and edge faces - edge_pts.resize(n_edges, 2); - edge_faces.resize(n_edges, 2); - IntT pos = 0; - - auto new_edge = [&](const auto& edge) { - edge_pts(pos, 0) = edge[0]; - edge_pts(pos, 1) = edge[1]; - edge_faces(pos, 0) = edge[2]; - edge_faces(pos, 1) = -1; - ++pos; - }; - - new_edge(edges[0]); - for (IntT i = 1; i < static_cast(edges.size()); ++i) { - if (EdgesEqual(edges[i], edges[i - 1])) { - edge_faces(pos - 1, 1) = edges[i][2]; - } else { - new_edge(edges[i]); - } + } + + // Sort to identify duplicates + std::sort(edges.begin(), edges.end(), + [](const auto& a, const auto& b) { return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); }); + + // Count unique edges + IntT n_edges = 1; + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + n_edges += static_cast(!EdgesEqual(edges[i], edges[i - 1])); + } + + // Map edge points and edge faces + edge_pts.resize(n_edges, 2); + edge_faces.resize(n_edges, 2); + IntT pos = 0; + + auto new_edge = [&](const auto& edge) { + edge_pts(pos, 0) = edge[0]; + edge_pts(pos, 1) = edge[1]; + edge_faces(pos, 0) = edge[2]; + edge_faces(pos, 1) = -1; + ++pos; + }; + + new_edge(edges[0]); + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + if (EdgesEqual(edges[i], edges[i - 1])) { + edge_faces(pos - 1, 1) = edges[i][2]; + } else { + new_edge(edges[i]); } + } } /*! @@ -149,483 +146,469 @@ inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, * Returns tuple of (n_bands, pt_to_band, x_bands). */ inline auto DetectBands(const VectorReal& x, IntT max_bands = 0) { - // Safety check: empty input - if (x.empty()) { - return std::make_tuple(IntT{0}, std::vector{}, VectorReal{}); + // Safety check: empty input + if (x.empty()) { + return std::make_tuple(IntT{0}, std::vector{}, VectorReal{}); + } + + // First pass: count unique x-coordinates + IntT n_unique = 1; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != x[i - 1]) n_unique++; + } + + // Determine if we need to use coarser bands + if (max_bands <= 0) { + max_bands = std::min(IntT{5000}, static_cast(4.0 * std::sqrt(static_cast(x.size())))); + } + + const bool use_coarse_bands = (n_unique > max_bands); + + if (!use_coarse_bands) { + // Original algorithm: use all unique x-coordinates as band boundaries + std::vector pt_to_band(x.size()); + pt_to_band[0] = 0; + IntT n_bands = 0; + + for (IntT i = 1; i < static_cast(x.size()); ++i) { + n_bands += static_cast(x[i] != x[i - 1]); + pt_to_band[i] = n_bands; } - - // First pass: count unique x-coordinates - IntT n_unique = 1; + + VectorReal x_bands(n_bands + 1); + IntT pos = 0; + x_bands[pos] = x[0]; for (IntT i = 1; i < static_cast(x.size()); ++i) { - if (x[i] != x[i - 1]) n_unique++; + if (x[i] != x_bands[pos]) { + x_bands[++pos] = x[i]; + } } - - // Determine if we need to use coarser bands - if (max_bands <= 0) { - max_bands = std::min(IntT{5000}, static_cast(4.0 * std::sqrt(static_cast(x.size())))); + + return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); + } else { + // Coarse banding: divide x-range into max_bands equal intervals + const RealT x_min = x.front(); + const RealT x_max = x.back(); + const RealT x_range = x_max - x_min; + + // Handle degenerate case where all x are the same + if (x_range < 1e-30) { + std::vector pt_to_band(x.size(), 0); + VectorReal x_bands = {x_min, x_max}; + return std::make_tuple(IntT{1}, std::move(pt_to_band), std::move(x_bands)); } - - const bool use_coarse_bands = (n_unique > max_bands); - - if (!use_coarse_bands) { - // Original algorithm: use all unique x-coordinates as band boundaries - std::vector pt_to_band(x.size()); - pt_to_band[0] = 0; - IntT n_bands = 0; - - for (IntT i = 1; i < static_cast(x.size()); ++i) { - n_bands += static_cast(x[i] != x[i - 1]); - pt_to_band[i] = n_bands; - } - - VectorReal x_bands(n_bands + 1); - IntT pos = 0; - x_bands[pos] = x[0]; - for (IntT i = 1; i < static_cast(x.size()); ++i) { - if (x[i] != x_bands[pos]) { - x_bands[++pos] = x[i]; - } - } - - return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); - } else { - // Coarse banding: divide x-range into max_bands equal intervals - const RealT x_min = x.front(); - const RealT x_max = x.back(); - const RealT x_range = x_max - x_min; - - // Handle degenerate case where all x are the same - if (x_range < 1e-30) { - std::vector pt_to_band(x.size(), 0); - VectorReal x_bands = {x_min, x_max}; - return std::make_tuple(IntT{1}, std::move(pt_to_band), std::move(x_bands)); - } - - const IntT n_bands = max_bands; - const RealT band_width = x_range / max_bands; - - // Build band boundaries - VectorReal x_bands(max_bands + 1); - for (IntT i = 0; i <= max_bands; ++i) { - x_bands[i] = x_min + i * band_width; - } - x_bands[max_bands] = x_max; - - // Map each point to its band - std::vector pt_to_band(x.size()); - for (IntT i = 0; i < static_cast(x.size()); ++i) { - IntT band = static_cast((x[i] - x_min) / band_width); - band = std::max(IntT{0}, std::min(band, max_bands - 1)); - pt_to_band[i] = band; - } - - return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); + + const IntT n_bands = max_bands; + const RealT band_width = x_range / max_bands; + + // Build band boundaries + VectorReal x_bands(max_bands + 1); + for (IntT i = 0; i <= max_bands; ++i) { + x_bands[i] = x_min + i * band_width; + } + x_bands[max_bands] = x_max; + + // Map each point to its band + std::vector pt_to_band(x.size()); + for (IntT i = 0; i < static_cast(x.size()); ++i) { + IntT band = static_cast((x[i] - x_min) / band_width); + band = std::max(IntT{0}, std::min(band, max_bands - 1)); + pt_to_band[i] = band; } + + return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); + } } /*! * \brief Builds the trapezoidal map for efficient 2D queries. */ -inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, - const VectorReal& y, TrapezoidalMap& map) { - auto& x_bands = map.x_bands; - auto& offsets = map.offsets; - auto& edge_id = map.edge_id; - auto& edge_y = map.edge_y; - - const auto [n_bands, pt_to_band, bands] = DetectBands(x); - x_bands = std::move(bands); - - // Safety check - if (n_bands <= 0 || x_bands.size() < 2) { - x_bands.clear(); - offsets.clear(); - edge_id.clear(); - edge_y.clear(); - return; +inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, const VectorReal& y, + TrapezoidalMap& map) { + auto& x_bands = map.x_bands; + auto& offsets = map.offsets; + auto& edge_id = map.edge_id; + auto& edge_y = map.edge_y; + + const auto [n_bands, pt_to_band, bands] = DetectBands(x); + x_bands = std::move(bands); + + // Safety check + if (n_bands <= 0 || x_bands.size() < 2) { + x_bands.clear(); + offsets.clear(); + edge_id.clear(); + edge_y.clear(); + return; + } + + // Helper function to find band index for an x-coordinate + auto find_band = [&x_bands, n_bands](RealT x_val) -> IntT { + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); + IntT idx = static_cast(it - x_bands.begin()); + idx = std::max(IntT{0}, idx - 1); + idx = std::min(idx, n_bands - 1); + return idx; + }; + + // Count edges per band + auto& counts = offsets; + counts.clear(); + counts.resize(n_bands + 1, 0); + + for (IntT i = 0; i < static_cast(edge_pts.rows()); ++i) { + const IntT pt_0 = edge_pts(i, 0); + const IntT pt_1 = edge_pts(i, 1); + const RealT x_0 = x[pt_0]; + const RealT x_1 = x[pt_1]; + + const IntT band_0 = find_band(x_0); + const IntT band_1 = find_band(x_1); + + const IntT min_band = std::min(band_0, band_1); + const IntT max_band = std::max(band_0, band_1); + + for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { + ++counts[j + 1]; } - - // Helper function to find band index for an x-coordinate - auto find_band = [&x_bands, n_bands](RealT x_val) -> IntT { - auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); - IntT idx = static_cast(it - x_bands.begin()); - idx = std::max(IntT{0}, idx - 1); - idx = std::min(idx, n_bands - 1); - return idx; - }; - - // Count edges per band - auto& counts = offsets; - counts.clear(); - counts.resize(n_bands + 1, 0); - - for (IntT i = 0; i < static_cast(edge_pts.rows()); ++i) { - const IntT pt_0 = edge_pts(i, 0); - const IntT pt_1 = edge_pts(i, 1); - const RealT x_0 = x[pt_0]; - const RealT x_1 = x[pt_1]; - - const IntT band_0 = find_band(x_0); - const IntT band_1 = find_band(x_1); - - const IntT min_band = std::min(band_0, band_1); - const IntT max_band = std::max(band_0, band_1); - - for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { - ++counts[j + 1]; - } - } - - // Convert counts to offsets (CSR format) - for (IntT i = 2; i < static_cast(offsets.size()); ++i) { - offsets[i] += offsets[i - 1]; + } + + // Convert counts to offsets (CSR format) + for (IntT i = 2; i < static_cast(offsets.size()); ++i) { + offsets[i] += offsets[i - 1]; + } + + // Check memory requirement + const size_t total_entries = static_cast(offsets.back()); + const size_t memory_mb = total_entries * (sizeof(IntT) + sizeof(RealT)) / (1024 * 1024); + + if (memory_mb > 2048) { + x_bands.clear(); + offsets.clear(); + edge_id.clear(); + edge_y.clear(); + return; + } + + // Allocate storage + edge_id.resize(offsets.back()); + edge_y.resize(offsets.back()); + auto pos = offsets; + + // Store edges in each band + for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { + const IntT pt_0 = edge_pts(i_edge, 0); + const IntT pt_1 = edge_pts(i_edge, 1); + const RealT x_0 = x[pt_0], y_0 = y[pt_0]; + const RealT x_1 = x[pt_1], y_1 = y[pt_1]; + + const IntT band_0 = find_band(x_0); + const IntT band_1 = find_band(x_1); + + if (band_0 == band_1) continue; + + const IntT min_band = std::min(band_0, band_1); + const IntT max_band = std::max(band_0, band_1); + const RealT dy_dx = (y_1 - y_0) / (x_1 - x_0); + + for (IntT j = min_band; j < max_band && j < n_bands; ++j) { + if (pos[j] < static_cast(edge_id.size())) { + edge_id[pos[j]] = i_edge; + const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2; + edge_y[pos[j]] = y_0 + dy_dx * (x_mid - x_0); + ++pos[j]; + } } - - // Check memory requirement - const size_t total_entries = static_cast(offsets.back()); - const size_t memory_mb = total_entries * (sizeof(IntT) + sizeof(RealT)) / (1024 * 1024); - - if (memory_mb > 2048) { - x_bands.clear(); - offsets.clear(); - edge_id.clear(); - edge_y.clear(); - return; - } - - // Allocate storage - edge_id.resize(offsets.back()); - edge_y.resize(offsets.back()); - auto pos = offsets; - - // Store edges in each band - for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { - const IntT pt_0 = edge_pts(i_edge, 0); - const IntT pt_1 = edge_pts(i_edge, 1); - const RealT x_0 = x[pt_0], y_0 = y[pt_0]; - const RealT x_1 = x[pt_1], y_1 = y[pt_1]; - - const IntT band_0 = find_band(x_0); - const IntT band_1 = find_band(x_1); - - if (band_0 == band_1) continue; - - const IntT min_band = std::min(band_0, band_1); - const IntT max_band = std::max(band_0, band_1); - const RealT dy_dx = (y_1 - y_0) / (x_1 - x_0); - - for (IntT j = min_band; j < max_band && j < n_bands; ++j) { - if (pos[j] < static_cast(edge_id.size())) { - edge_id[pos[j]] = i_edge; - const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2; - edge_y[pos[j]] = y_0 + dy_dx * (x_mid - x_0); - ++pos[j]; - } - } + } + + // Sort edges in each band by y coordinate + std::vector> tmp; + for (IntT i = 0; i < n_bands; ++i) { + const IntT begin = offsets[i]; + const IntT end = offsets[i + 1]; + if (begin >= end) continue; + + tmp.resize(end - begin); + for (auto k = begin; k < end; ++k) { + tmp[k - begin] = {edge_id[k], edge_y[k]}; } - - // Sort edges in each band by y coordinate - std::vector> tmp; - for (IntT i = 0; i < n_bands; ++i) { - const IntT begin = offsets[i]; - const IntT end = offsets[i + 1]; - if (begin >= end) continue; - - tmp.resize(end - begin); - for (auto k = begin; k < end; ++k) { - tmp[k - begin] = {edge_id[k], edge_y[k]}; - } - std::sort(tmp.begin(), tmp.end(), [](const auto& a, const auto& b) { - return a.second < b.second; - }); - for (auto k = begin; k < end; ++k) { - edge_id[k] = tmp[k - begin].first; - edge_y[k] = tmp[k - begin].second; - } + std::sort(tmp.begin(), tmp.end(), [](const auto& a, const auto& b) { return a.second < b.second; }); + for (auto k = begin; k < end; ++k) { + edge_id[k] = tmp[k - begin].first; + edge_y[k] = tmp[k - begin].second; } + } } /*! * \brief Query the trapezoidal map for edges bounding a point. * Returns pair of (edge_below, edge_above). Either can be -1 if at boundary. */ -inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, - const VectorReal& x_coords, const VectorReal& y_coords, - const RealT& x, const RealT& y) { - // Safety check - if (map.x_bands.size() < 2 || map.offsets.empty()) { - return std::make_pair(IntT{-1}, IntT{-1}); +inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, const VectorReal& x_coords, + const VectorReal& y_coords, const RealT& x, const RealT& y) { + // Safety check + if (map.x_bands.size() < 2 || map.offsets.empty()) { + return std::make_pair(IntT{-1}, IntT{-1}); + } + + const IntT n_bands = static_cast(map.x_bands.size()) - 1; + + // Find the band containing x + const auto& x_bands = map.x_bands; + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x); + IntT d = static_cast(it - x_bands.begin()); + IntT center_band = std::min(std::max(IntT{0}, d - 1), n_bands - 1); + + // Search in current band AND adjacent bands + RealT best_y_below = -1e300; + RealT best_y_above = 1e300; + IntT edge_below = -1; + IntT edge_above = -1; + + for (IntT band_offset = -1; band_offset <= 1; ++band_offset) { + IntT band_idx = center_band + band_offset; + if (band_idx < 0 || band_idx >= n_bands) continue; + if (band_idx + 1 >= static_cast(map.offsets.size())) continue; + + const IntT begin = map.offsets[band_idx]; + const IntT end = map.offsets[band_idx + 1]; + + if (begin < 0 || end > static_cast(map.edge_id.size()) || begin >= end) { + continue; } - - const IntT n_bands = static_cast(map.x_bands.size()) - 1; - - // Find the band containing x - const auto& x_bands = map.x_bands; - auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x); - IntT d = static_cast(it - x_bands.begin()); - IntT center_band = std::min(std::max(IntT{0}, d - 1), n_bands - 1); - - // Search in current band AND adjacent bands - RealT best_y_below = -1e300; - RealT best_y_above = 1e300; - IntT edge_below = -1; - IntT edge_above = -1; - - for (IntT band_offset = -1; band_offset <= 1; ++band_offset) { - IntT band_idx = center_band + band_offset; - if (band_idx < 0 || band_idx >= n_bands) continue; - if (band_idx + 1 >= static_cast(map.offsets.size())) continue; - - const IntT begin = map.offsets[band_idx]; - const IntT end = map.offsets[band_idx + 1]; - - if (begin < 0 || end > static_cast(map.edge_id.size()) || begin >= end) { - continue; - } - - for (IntT k = begin; k < end; ++k) { - const IntT e_id = map.edge_id[k]; - - const IntT p0 = edge_pts(e_id, 0); - const IntT p1 = edge_pts(e_id, 1); - const RealT x0 = x_coords[p0], y0 = y_coords[p0]; - const RealT x1 = x_coords[p1], y1 = y_coords[p1]; - - // Check if query x is within edge's x-extent - if (x < std::min(x0, x1) - 1e-10 || x > std::max(x0, x1) + 1e-10) { - continue; - } - - // Compute y-position of edge at query x - RealT edge_y_at_x; - const RealT dx = x1 - x0; - if (std::abs(dx) < 1e-30) { - edge_y_at_x = (y0 + y1) / 2.0; - } else { - const RealT t = (x - x0) / dx; - edge_y_at_x = y0 + t * (y1 - y0); - } - - if (edge_y_at_x <= y + 1e-10 && edge_y_at_x > best_y_below) { - best_y_below = edge_y_at_x; - edge_below = e_id; - } - if (edge_y_at_x >= y - 1e-10 && edge_y_at_x < best_y_above) { - best_y_above = edge_y_at_x; - edge_above = e_id; - } - } + + for (IntT k = begin; k < end; ++k) { + const IntT e_id = map.edge_id[k]; + + const IntT p0 = edge_pts(e_id, 0); + const IntT p1 = edge_pts(e_id, 1); + const RealT x0 = x_coords[p0], y0 = y_coords[p0]; + const RealT x1 = x_coords[p1], y1 = y_coords[p1]; + + // Check if query x is within edge's x-extent + if (x < std::min(x0, x1) - 1e-10 || x > std::max(x0, x1) + 1e-10) { + continue; + } + + // Compute y-position of edge at query x + RealT edge_y_at_x; + const RealT dx = x1 - x0; + if (std::abs(dx) < 1e-30) { + edge_y_at_x = (y0 + y1) / 2.0; + } else { + const RealT t = (x - x0) / dx; + edge_y_at_x = y0 + t * (y1 - y0); + } + + if (edge_y_at_x <= y + 1e-10 && edge_y_at_x > best_y_below) { + best_y_below = edge_y_at_x; + edge_below = e_id; + } + if (edge_y_at_x >= y - 1e-10 && edge_y_at_x < best_y_above) { + best_y_above = edge_y_at_x; + edge_above = e_id; + } } - - return std::make_pair(edge_below, edge_above); + } + + return std::make_pair(edge_below, edge_above); } /*! * \brief Get triangles adjacent to two query edges (up to 3 triangles). */ -inline auto AdjacentTriangles(const IntT edge_0, const IntT edge_1, - const Matrix2i& edge_faces) { - std::array tris = {-1, -1, -1}; - IntT pos = 0; +inline auto AdjacentTriangles(const IntT edge_0, const IntT edge_1, const Matrix2i& edge_faces) { + std::array tris = {-1, -1, -1}; + IntT pos = 0; + + auto insert = [&tris, &pos](const IntT t) { + if (t < 0) return; + for (IntT i = 0; i < pos; ++i) { + if (t == tris[i]) return; + } + tris[pos++] = t; + }; + + auto get_tris = [&edge_faces](const IntT e) { + if (e < 0) return std::array{IntT{-1}, IntT{-1}}; + return std::array{edge_faces(e, 0), edge_faces(e, 1)}; + }; - auto insert = [&tris, &pos](const IntT t) { - if (t < 0) return; - for (IntT i = 0; i < pos; ++i) { - if (t == tris[i]) return; - } - tris[pos++] = t; - }; - - auto get_tris = [&edge_faces](const IntT e) { - if (e < 0) return std::array{IntT{-1}, IntT{-1}}; - return std::array{edge_faces(e, 0), edge_faces(e, 1)}; - }; - - for (const auto e : {edge_0, edge_1}) { - for (const auto t : get_tris(e)) { - insert(t); - } + for (const auto e : {edge_0, edge_1}) { + for (const auto t : get_tris(e)) { + insert(t); } - return tris; + } + return tris; } /*! * \brief Compute barycentric coordinates of point (x, y) in triangle. */ -inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, - const VectorReal& x, const VectorReal& y, +inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, const VectorReal& x, const VectorReal& y, const RealT x_q, const RealT y_q) { - const IntT p0 = triangles(i_tri, 0); - const IntT p1 = triangles(i_tri, 1); - const IntT p2 = triangles(i_tri, 2); - - const RealT x0 = x[p0], y0 = y[p0]; - const RealT x1 = x[p1], y1 = y[p1]; - const RealT x2 = x[p2], y2 = y[p2]; - - const RealT dx1 = x1 - x0, dy1 = y1 - y0; - const RealT dx2 = x2 - x0, dy2 = y2 - y0; - - auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { - return ux * vy - uy * vx; - }; - - const RealT det = cross(dx1, dy1, dx2, dy2); - if (std::abs(det) < 1e-30) { - return std::array{RealT{0}, RealT{0}, RealT{0}}; - } - - const RealT inv_det = 1.0 / det; - const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; - const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; - - return std::array{1 - a - b, a, b}; + const IntT p0 = triangles(i_tri, 0); + const IntT p1 = triangles(i_tri, 1); + const IntT p2 = triangles(i_tri, 2); + + const RealT x0 = x[p0], y0 = y[p0]; + const RealT x1 = x[p1], y1 = y[p1]; + const RealT x2 = x[p2], y2 = y[p2]; + + const RealT dx1 = x1 - x0, dy1 = y1 - y0; + const RealT dx2 = x2 - x0, dy2 = y2 - y0; + + auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { return ux * vy - uy * vx; }; + + const RealT det = cross(dx1, dy1, dx2, dy2); + if (std::abs(det) < 1e-30) { + return std::array{RealT{0}, RealT{0}, RealT{0}}; + } + + const RealT inv_det = 1.0 / det; + const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; + const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; + + return std::array{1 - a - b, a, b}; } /*! * \brief Check if point is inside triangle. */ inline bool InTriangle(const std::array& coords, const RealT tol = 0.0) { - return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; + return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; } /*! * \brief Find triangle containing point using the trapezoidal map. */ -inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, - const Matrix2i& edge_pts, const Matrix2i& edge_faces, - const VectorReal& x, const VectorReal& y, - const RealT x_q, const RealT y_q, std::array& bary_out) { - - const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); - const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); - - constexpr RealT tol = 1e-12; - for (const auto t : candidates) { - if (t < 0) continue; - - const auto coords = TriangleCoords(t, triangles, x, y, x_q, y_q); - if (InTriangle(coords, tol)) { - bary_out = coords; - return t; - } +inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, const Matrix2i& edge_pts, + const Matrix2i& edge_faces, const VectorReal& x, const VectorReal& y, const RealT x_q, + const RealT y_q, std::array& bary_out) { + const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); + const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); + + constexpr RealT tol = 1e-12; + for (const auto t : candidates) { + if (t < 0) continue; + + const auto coords = TriangleCoords(t, triangles, x, y, x_q, y_q); + if (InTriangle(coords, tol)) { + bary_out = coords; + return t; } - - bary_out = {0.0, 0.0, 0.0}; - return -1; + } + + bary_out = {0.0, 0.0, 0.0}; + return -1; } } // namespace su2_lut - /*! * \class CTrapezoidalMapFAST * \brief Wrapper class for SU2 integration. */ class CTrapezoidalMapFAST { -private: - su2_lut::Matrix3i triangles; - su2_lut::Matrix2i edge_pts, edge_faces; - su2_lut::VectorReal x_coords, y_coords; - su2_lut::TrapezoidalMap map; - - unsigned long n_points; - unsigned long n_triangles; - bool build_successful; - -public: - CTrapezoidalMapFAST() : n_points(0), n_triangles(0), build_successful(false) {} - - /*! - * \brief Build the trapezoidal map from mesh data. - */ - bool Build(unsigned long num_points, unsigned long num_triangles, - const su2double* x, const su2double* y, - const unsigned long* connectivity) { - - n_points = num_points; - n_triangles = num_triangles; - build_successful = false; - - if (num_points == 0 || num_triangles == 0) { - x_coords.clear(); - y_coords.clear(); - triangles.resize(0, 3); - edge_pts.resize(0, 2); - edge_faces.resize(0, 2); - map.x_bands.clear(); - map.offsets.clear(); - map.edge_id.clear(); - map.edge_y.clear(); - return false; - } - - // Copy coordinates - x_coords.assign(x, x + num_points); - y_coords.assign(y, y + num_points); - - // Copy triangle connectivity - triangles.resize(num_triangles, 3); - for (size_t i = 0; i < num_triangles; ++i) { - unsigned long v0 = connectivity[3*i + 0]; - unsigned long v1 = connectivity[3*i + 1]; - unsigned long v2 = connectivity[3*i + 2]; - - if (v0 >= num_points || v1 >= num_points || v2 >= num_points) { - triangles(i, 0) = 0; - triangles(i, 1) = 0; - triangles(i, 2) = 0; - } else { - triangles(i, 0) = static_cast(v0); - triangles(i, 1) = static_cast(v1); - triangles(i, 2) = static_cast(v2); - } - } - - // Build the map - su2_lut::ReorderPoints(triangles, x_coords, y_coords); - su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); - su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); - - if (map.x_bands.empty() || map.offsets.empty()) { - build_successful = false; - return false; - } - - build_successful = true; - return true; + private: + su2_lut::Matrix3i triangles; + su2_lut::Matrix2i edge_pts, edge_faces; + su2_lut::VectorReal x_coords, y_coords; + su2_lut::TrapezoidalMap map; + + unsigned long n_points; + unsigned long n_triangles; + bool build_successful; + + public: + CTrapezoidalMapFAST() : n_points(0), n_triangles(0), build_successful(false) {} + + /*! + * \brief Build the trapezoidal map from mesh data. + */ + bool Build(unsigned long num_points, unsigned long num_triangles, const su2double* x, const su2double* y, + const unsigned long* connectivity) { + n_points = num_points; + n_triangles = num_triangles; + build_successful = false; + + if (num_points == 0 || num_triangles == 0) { + x_coords.clear(); + y_coords.clear(); + triangles.resize(0, 3); + edge_pts.resize(0, 2); + edge_faces.resize(0, 2); + map.x_bands.clear(); + map.offsets.clear(); + map.edge_id.clear(); + map.edge_y.clear(); + return false; + } + + // Copy coordinates + x_coords.assign(x, x + num_points); + y_coords.assign(y, y + num_points); + + // Copy triangle connectivity + triangles.resize(num_triangles, 3); + for (size_t i = 0; i < num_triangles; ++i) { + unsigned long v0 = connectivity[3 * i + 0]; + unsigned long v1 = connectivity[3 * i + 1]; + unsigned long v2 = connectivity[3 * i + 2]; + + if (v0 >= num_points || v1 >= num_points || v2 >= num_points) { + triangles(i, 0) = 0; + triangles(i, 1) = 0; + triangles(i, 2) = 0; + } else { + triangles(i, 0) = static_cast(v0); + triangles(i, 1) = static_cast(v1); + triangles(i, 2) = static_cast(v2); + } + } + + // Build the map + su2_lut::ReorderPoints(triangles, x_coords, y_coords); + su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); + su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); + + if (map.x_bands.empty() || map.offsets.empty()) { + build_successful = false; + return false; } - - bool IsBuilt() const { return build_successful; } - - /*! - * \brief Find the triangle containing a query point. - */ - bool FindTriangle(su2double val_x, su2double val_y, - unsigned long& triangle_id, - std::array& bary_coords) const { - - if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { - bary_coords = {0.0, 0.0, 0.0}; - return false; - } - - std::array bary; - su2_lut::IntT tri_id = su2_lut::FindTriangle( - map, triangles, edge_pts, edge_faces, x_coords, y_coords, val_x, val_y, bary); - - if (tri_id >= 0) { - triangle_id = static_cast(tri_id); - bary_coords = {bary[0], bary[1], bary[2]}; - return true; - } - return false; + + build_successful = true; + return true; + } + + bool IsBuilt() const { return build_successful; } + + /*! + * \brief Find the triangle containing a query point. + */ + bool FindTriangle(su2double val_x, su2double val_y, unsigned long& triangle_id, + std::array& bary_coords) const { + if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { + bary_coords = {0.0, 0.0, 0.0}; + return false; + } + + std::array bary; + su2_lut::IntT tri_id = + su2_lut::FindTriangle(map, triangles, edge_pts, edge_faces, x_coords, y_coords, val_x, val_y, bary); + + if (tri_id >= 0) { + triangle_id = static_cast(tri_id); + bary_coords = {bary[0], bary[1], bary[2]}; + return true; } + return false; + } - unsigned long GetNTriangles() const { return n_triangles; } - unsigned long GetNPoints() const { return n_points; } - unsigned long GetNEdges() const { return edge_pts.rows(); } - unsigned long GetNBands() const { return map.x_bands.size() > 0 ? map.x_bands.size() - 1 : 0; } + unsigned long GetNTriangles() const { return n_triangles; } + unsigned long GetNPoints() const { return n_points; } + unsigned long GetNEdges() const { return edge_pts.rows(); } + unsigned long GetNBands() const { return map.x_bands.size() > 0 ? map.x_bands.size() - 1 : 0; } }; #endif // CTRAPEZOIDALMAP_FAST_HPP diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 05344a5484c9..005ea10d17c6 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -46,8 +46,7 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, FindTableLimits(name_CV1, name_CV2); - if (rank == MASTER_NODE) - cout << "Detecting all unique edges and setting edge to triangle connectivity ..." << endl; + if (rank == MASTER_NODE) cout << "Detecting all unique edges and setting edge to triangle connectivity ..." << endl; IdentifyUniqueEdges(); @@ -61,7 +60,7 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, /*--- NOTE: Trapezoidal map is NOT built here anymore! ---*/ /*--- Call EnableFastTrapMap() for LUT_FAST (efficient O(n) memory) ---*/ /*--- Call BuildOriginalTrapMap() for standard LUT (original SU2 implementation) ---*/ - + ComputeInterpCoeffs(); if (rank == MASTER_NODE) cout << "LUT fluid model ready for use" << endl; @@ -70,15 +69,13 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, void CLookUpTable::BuildOriginalTrapMap() { /*--- Build the original SU2 trapezoidal map (DAG-based, O(n log n) memory) ---*/ /*--- WARNING: This can use a lot of memory for large tables! ---*/ - + if (rank == MASTER_NODE) switch (table_dim) { case 2: - cout << "Building a trapezoidal map for the (" + name_CV1 + ", " + name_CV2 + - ") space ..." << endl; + cout << "Building a trapezoidal map for the (" + name_CV1 + ", " + name_CV2 + ") space ..." << endl; break; case 3: - cout << "Building trapezoidal map stack for the (" + name_CV1 + ", " + name_CV2 + - ") space ..." << endl; + cout << "Building trapezoidal map stack for the (" + name_CV1 + ", " + name_CV2 + ") space ..." << endl; break; default: break; @@ -89,13 +86,13 @@ void CLookUpTable::BuildOriginalTrapMap() { unsigned short barwidth = 65; bool display_map_info = (n_table_levels < 2); double tmap_memory_footprint = 0; - + for (auto i_level = 0ul; i_level < n_table_levels; i_level++) { trap_map_x_y[i_level] = CTrapezoidalMap(GetDataP(name_CV1, i_level), GetDataP(name_CV2, i_level), table_data[i_level].cols(), edges[i_level], edge_to_triangle[i_level], display_map_info); tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint(); - + /* Display a progress bar to monitor table generation process */ if (rank == MASTER_NODE) { su2double progress = su2double(i_level) / n_table_levels; @@ -127,11 +124,11 @@ void CLookUpTable::BuildOriginalTrapMap() { void CLookUpTable::EnableFastTrapMap() { /*--- Build Memory-efficient trapezoidal map (band-based, O(n) memory) ---*/ - + use_fast_trap_map_ = true; - + trap_map_x_y_FAST.resize(n_table_levels); - + if (rank == MASTER_NODE) { cout << endl; cout << "+--------------------------------------------------+" << endl; @@ -139,28 +136,26 @@ void CLookUpTable::EnableFastTrapMap() { cout << "| (O(n) memory, O(log n) queries) |" << endl; cout << "+--------------------------------------------------+" << endl; } - + su2double startTime = SU2_MPI::Wtime(); - + for (unsigned long i_level = 0; i_level < n_table_levels; ++i_level) { - const auto n_pts = n_points[i_level]; const auto n_tris = n_triangles[i_level]; - + if (rank == MASTER_NODE) { - cout << " Level " << i_level << ": " << n_pts << " points, " - << n_tris << " triangles ... " << flush; + cout << " Level " << i_level << ": " << n_pts << " points, " << n_tris << " triangles ... " << flush; } - + // Get point coordinates std::vector x_coords(n_pts); std::vector y_coords(n_pts); - + for (unsigned long i_point = 0; i_point < n_pts; ++i_point) { x_coords[i_point] = table_data[i_level][idx_CV1][i_point]; y_coords[i_point] = table_data[i_level][idx_CV2][i_point]; } - + // Get triangle connectivity std::vector tri_conn(3 * n_tris); for (unsigned long i_tri = 0; i_tri < n_tris; ++i_tri) { @@ -168,24 +163,22 @@ void CLookUpTable::EnableFastTrapMap() { tri_conn[3 * i_tri + 1] = triangles[i_level][i_tri][1]; tri_conn[3 * i_tri + 2] = triangles[i_level][i_tri][2]; } - + // Build FAST trapezoidal map for this level - trap_map_x_y_FAST[i_level].Build(n_pts, n_tris, - x_coords.data(), y_coords.data(), - tri_conn.data()); - + trap_map_x_y_FAST[i_level].Build(n_pts, n_tris, x_coords.data(), y_coords.data(), tri_conn.data()); + if (rank == MASTER_NODE) { cout << "done." << endl; } } - + su2double stopTime = SU2_MPI::Wtime(); - + if (rank == MASTER_NODE) { cout << "+--------------------------------------------------+" << endl; cout << "| FAST trapezoidal maps built successfully! |" << endl; - cout << "| Build time: " << setw(10) << setprecision(3) << (stopTime - startTime) - << " seconds |" << endl; + cout << "| Build time: " << setw(10) << setprecision(3) << (stopTime - startTime) << " seconds |" + << endl; cout << "+--------------------------------------------------+" << endl; cout << endl; } @@ -711,7 +704,7 @@ bool CLookUpTable::FindInclusionTriangle(const su2double val_CV1, const su2doubl return IsInTriangle(val_CV1, val_CV2, id_triangle, iLevel); } } - + return false; } diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 901b4a80aa88..12ae1a9dc54b 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -55,7 +55,8 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati scalars_vector.resize(n_scalars); table_scalar_names.resize(n_scalars); - for (auto iCV = 0u; iCV < n_control_vars; iCV++) table_scalar_names[iCV] = flamelet_options.controlling_variable_names[iCV]; + for (auto iCV = 0u; iCV < n_control_vars; iCV++) + table_scalar_names[iCV] = flamelet_options.controlling_variable_names[iCV]; /*--- auxiliary species transport equations---*/ for (auto i_aux = 0u; i_aux < n_user_scalars; i_aux++) { @@ -64,10 +65,11 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati controlling_variable_names.resize(n_control_vars); for (auto iCV = 0u; iCV < n_control_vars; iCV++) - controlling_variable_names[iCV] =flamelet_options.controlling_variable_names[iCV]; + controlling_variable_names[iCV] = flamelet_options.controlling_variable_names[iCV]; passive_specie_names.resize(n_user_scalars); - for (auto i_aux = 0u; i_aux < n_user_scalars; i_aux++) passive_specie_names[i_aux] = flamelet_options.user_scalar_names[i_aux]; + for (auto i_aux = 0u; i_aux < n_user_scalars; i_aux++) + passive_specie_names[i_aux] = flamelet_options.user_scalar_names[i_aux]; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: @@ -81,7 +83,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati /*--- Build the original trapezoidal map (can use lots of memory!) ---*/ look_up_table->BuildOriginalTrapMap(); break; - + case ENUM_DATADRIVEN_METHOD::LUT_FAST: if (rank == MASTER_NODE) { cout << "**************************************************" << endl; @@ -94,7 +96,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati /*--- Build Memory-efficient trapezoidal map (O(n) memory) ---*/ look_up_table->EnableFastTrapMap(); break; - + default: if (rank == MASTER_NODE) { cout << "***********************************************" << endl; @@ -102,7 +104,8 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati cout << "***********************************************" << endl; } #ifdef USE_MLPCPP - lookup_mlp = new MLPToolbox::CLookUp_ANN(datadriven_fluid_options.n_filenames, datadriven_fluid_options.datadriven_filenames); + lookup_mlp = new MLPToolbox::CLookUp_ANN(datadriven_fluid_options.n_filenames, + datadriven_fluid_options.datadriven_filenames); if ((rank == MASTER_NODE)) lookup_mlp->DisplayNetworkInfo(); #else SU2_MPI::Error("SU2 was not compiled with MLPCpp enabled (-Denable-mlpcpp=true).", CURRENT_FUNCTION); @@ -120,7 +123,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati } CFluidFlamelet::~CFluidFlamelet() { - if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT || + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT || Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT_FAST) delete look_up_table; @@ -202,8 +205,7 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { size_t n_sources = n_control_vars + 2 * n_user_scalars; varnames_Sources.resize(n_sources); val_vars_Sources.resize(n_sources); - for (auto iCV = 0u; iCV < n_control_vars; iCV++) - varnames_Sources[iCV] = flamelet_options.cv_source_names[iCV]; + for (auto iCV = 0u; iCV < n_control_vars; iCV++) varnames_Sources[iCV] = flamelet_options.cv_source_names[iCV]; /*--- No source term for enthalpy ---*/ /*--- For the auxiliary equations, we use a positive (production) and a negative (consumption) term: @@ -224,7 +226,8 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { } else { varnames_LookUp.resize(n_lookups); val_vars_LookUp.resize(n_lookups); - for (auto iLookup = 0u; iLookup < n_lookups; iLookup++) varnames_LookUp[iLookup] = flamelet_options.lookup_names[iLookup]; + for (auto iLookup = 0u; iLookup < n_lookups; iLookup++) + varnames_LookUp[iLookup] = flamelet_options.lookup_names[iLookup]; } /*--- Preferential diffusion scalars ---*/ @@ -246,7 +249,7 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { case ENUM_DATADRIVEN_METHOD::LUT: case ENUM_DATADRIVEN_METHOD::LUT_FAST: if (preferential_diffusion) { - preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); + preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); } break; case ENUM_DATADRIVEN_METHOD::MLP: @@ -280,10 +283,10 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { } #endif } else { - for (auto iVar=0u; iVar < varnames_TD.size(); iVar++) { + 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++) { + for (auto iVar = 0u; iVar < varnames_Sources.size(); iVar++) { unsigned long LUT_idx; if (noSource(varnames_Sources[iVar])) { LUT_idx = look_up_table->GetNullIndex(); @@ -292,16 +295,16 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { } LUT_idx_Sources.push_back(LUT_idx); } - for (auto iVar=0u; iVar < varnames_LookUp.size(); iVar++) { + for (auto iVar = 0u; iVar < varnames_LookUp.size(); iVar++) { unsigned long LUT_idx; if (noSource(varnames_LookUp[iVar])) LUT_idx = look_up_table->GetNullIndex(); - else + else LUT_idx = look_up_table->GetIndexOfVar(varnames_LookUp[iVar]); LUT_idx_LookUp.push_back(LUT_idx); } if (preferential_diffusion) { - for (auto iVar=0u; iVar < varnames_PD.size(); iVar++) { + for (auto iVar = 0u; iVar < varnames_PD.size(); iVar++) { LUT_idx_PD.push_back(look_up_table->GetIndexOfVar(varnames_PD[iVar])); } } @@ -312,7 +315,7 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca vector& output_refs) { AD::StartPreacc(); for (auto iVar = 0u; iVar < input_scalar.size(); iVar++) AD::SetPreaccIn(input_scalar[iVar]); - + su2double val_enth = input_scalar[I_ENTH]; su2double val_prog = input_scalar[I_PROGVAR]; su2double val_mixfrac = include_mixture_fraction ? input_scalar[I_MIXFRAC] : 0.0; @@ -347,7 +350,6 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca default: break; } - /*--- Add all quantities and their names to the look up vectors. ---*/ bool inside; @@ -361,8 +363,10 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca } else { inside = look_up_table->LookUp_XY(LUT_idx, output_refs, val_prog, val_enth); } - if (inside) extrapolation = 0; - else extrapolation = 1; + 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 c2b1d7abe9b0..fa5668bdc6ac 100644 --- a/UnitTests/Common/containers/CLookupTable_tests.cpp +++ b/UnitTests/Common/containers/CLookupTable_tests.cpp @@ -248,4 +248,3 @@ TEST_CASE("LUTreader_3D_FAST", "[tabulated chemistry]") { look_up_table.LookUp_XYZ(idx_tag, &look_up_dat, prog, enth, mfrac); CHECK(look_up_dat == Approx(1.1738796125)); } - From 7e452d8c1db17fa19f962d3c53b5996c8c2e9f66 Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Tue, 6 Jan 2026 17:11:08 +0100 Subject: [PATCH 7/8] Fix AD compatibility in CTrapezoidalMapFAST --- Common/include/containers/CTrapezoidalMapFAST.hpp | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/Common/include/containers/CTrapezoidalMapFAST.hpp b/Common/include/containers/CTrapezoidalMapFAST.hpp index 55e26d4ac4eb..7b18f4919ac1 100644 --- a/Common/include/containers/CTrapezoidalMapFAST.hpp +++ b/Common/include/containers/CTrapezoidalMapFAST.hpp @@ -18,10 +18,13 @@ #include #include +#include "../basic_types/datatype_structure.hpp" + namespace su2_lut { using IntT = int32_t; using RealT = su2double; +using PassiveReal = passivedouble; /*--- Simple Matrix container (row-major) ---*/ template @@ -211,7 +214,7 @@ inline auto DetectBands(const VectorReal& x, IntT max_bands = 0) { // Map each point to its band std::vector pt_to_band(x.size()); for (IntT i = 0; i < static_cast(x.size()); ++i) { - IntT band = static_cast((x[i] - x_min) / band_width); + IntT band = static_cast(SU2_TYPE::GetValue((x[i] - x_min) / band_width)); band = std::max(IntT{0}, std::min(band, max_bands - 1)); pt_to_band[i] = band; } @@ -393,7 +396,7 @@ inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_ // Compute y-position of edge at query x RealT edge_y_at_x; const RealT dx = x1 - x0; - if (std::abs(dx) < 1e-30) { + if (std::abs(SU2_TYPE::GetValue(dx)) < 1e-30) { edge_y_at_x = (y0 + y1) / 2.0; } else { const RealT t = (x - x0) / dx; @@ -461,15 +464,15 @@ inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, const Ve auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { return ux * vy - uy * vx; }; const RealT det = cross(dx1, dy1, dx2, dy2); - if (std::abs(det) < 1e-30) { - return std::array{RealT{0}, RealT{0}, RealT{0}}; + if (std::abs(SU2_TYPE::GetValue(det)) < 1e-30) { + return std::array{RealT{0}, RealT{0}, RealT{0}}; } const RealT inv_det = 1.0 / det; const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; - return std::array{1 - a - b, a, b}; + return std::array{1 - a - b, a, b}; } /*! @@ -488,7 +491,7 @@ inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, c const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); - constexpr RealT tol = 1e-12; + const RealT tol = 1e-12; for (const auto t : candidates) { if (t < 0) continue; From 792ba49cc994884990252326a7b8a283e1ecb802 Mon Sep 17 00:00:00 2001 From: tkiymaz <79564236+tkiymaz@users.noreply.github.com> Date: Tue, 28 Jul 2026 13:26:44 +0300 Subject: [PATCH 8/8] replace the old trapezoidal map with the memory-efficient implementation --- Common/include/containers/CLookUpTable.hpp | 49 +- Common/include/containers/CTrapezoidalMap.hpp | 525 +++++++++++++-- .../containers/CTrapezoidalMapFAST.hpp | 617 ------------------ Common/include/option_structure.hpp | 4 +- Common/src/containers/CLookUpTable.cpp | 150 +---- Common/src/containers/CTrapezoidalMap.cpp | 283 -------- Common/src/containers/meson.build | 3 +- SU2_CFD/src/fluid/CFluidFlamelet.cpp | 26 +- .../Common/containers/CLookupTable_tests.cpp | 114 +--- 9 files changed, 491 insertions(+), 1280 deletions(-) delete mode 100644 Common/include/containers/CTrapezoidalMapFAST.hpp delete mode 100644 Common/src/containers/CTrapezoidalMap.cpp diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index 2458da32f424..8148631d4b46 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -35,7 +35,6 @@ #include "../../Common/include/option_structure.hpp" #include "CFileReaderLUT.hpp" #include "CTrapezoidalMap.hpp" -#include "CTrapezoidalMapFAST.hpp" /*! * \brief Look up table. @@ -101,16 +100,10 @@ class CLookUpTable { su2vector> hull; /*! \brief - * Trapezoidal map objects for the table levels (original SU2 implementation). + * Trapezoidal map objects for the table levels. */ su2vector trap_map_x_y; - /*! \brief - * Memory-efficient trapezoidal map objects for LUT_FAST. - */ - su2vector trap_map_x_y_FAST; - bool use_fast_trap_map_{false}; - /*! \brief * Vector of all the weight factors for the interpolation. */ @@ -267,29 +260,6 @@ class CLookUpTable { void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, const std::string& name_var, su2double* var_val, const unsigned long i_level = 0); - /*! - * \brief Determine if a point P(val_CV1,val_CV2) is inside the triangle val_id_triangle. - * \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] val_id_triangle - ID of the triangle to check. - * \returns True if the point is in the triangle, false if it is outside. - */ - bool IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, unsigned long i_level = 0); - - /*! - * \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); - } - /*! * \brief Compute the values of the first and second controlling variable based on normalized query coordinates * \param[in] inclusion_levels - Pair containing lower(first) and upper(second) table inclusion level indices. @@ -326,23 +296,6 @@ class CLookUpTable { public: CLookUpTable(const std::string& file_name_lut, std::string name_CV1_in, std::string name_CV2_in); - /*! - * \brief Build the original SU2 trapezoidal map (DAG-based). - * Use EnableFastTrapMap() for large tables instead. - */ - void BuildOriginalTrapMap(); - - /*! - * \brief Enable Memory-efficient trapezoidal map (LUT_FAST mode). - * This builds band-based maps with O(n) memory instead of the original O(n log n) DAG. - */ - void EnableFastTrapMap(); - - /*! - * \brief Check if fast trapezoidal map is enabled. - */ - inline bool UsingFastTrapMap() const { return use_fast_trap_map_; } - /*! * \brief Print information to screen. */ diff --git a/Common/include/containers/CTrapezoidalMap.hpp b/Common/include/containers/CTrapezoidalMap.hpp index afd9eebd3ed5..091bca917a33 100644 --- a/Common/include/containers/CTrapezoidalMap.hpp +++ b/Common/include/containers/CTrapezoidalMap.hpp @@ -1,7 +1,8 @@ /*! * \file CTrapezoidalMap.hpp - * \brief Implementation of the trapezoidal map for tabulation and lookup of fluid properties - * \author D. Mayer, T. Economon + * \brief Memory-efficient trapezoidal map for 2D lookup table queries, + * based on the LUT implementation of P. Gomes (https://github.com/pcarruscag/LUT). + * \author T. Kiymaz * \version 8.3.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -27,82 +28,494 @@ #pragma once -#include +#include +#include +#include +#include +#include +#include #include -#include "../../Common/include/linear_algebra/blas_structure.hpp" -#include "../../Common/include/toolboxes/CSquareMatrixCM.hpp" +#include "../basic_types/datatype_structure.hpp" + +namespace su2_lut { + +using IntT = int32_t; +using RealT = su2double; + +/*--- Simple row-major matrix with a fixed number of columns. ---*/ +template +struct Matrix { + std::vector data; + + void resize(size_t rows, size_t) { data.resize(rows * N); } + size_t rows() const { return data.size() / N; } + + const T& operator()(size_t i, size_t j) const { return data[i * N + j]; } + T& operator()(size_t i, size_t j) { return data[i * N + j]; } +}; + +using Matrix2i = Matrix; +using Matrix3i = Matrix; +using VectorInt = std::vector; +using VectorReal = std::vector; + +/*--- The map is defined by the limits of the bands in the x direction and a CSR of + * the edge IDs in each band, sorted by the edge y position at the band midpoint. ---*/ +struct TrapezoidalMap { + VectorInt offsets, edge_id; + VectorReal x_bands, edge_y; +}; /*! - * \class CTrapezoidalMap - * \ingroup LookUpInterp - * \brief Construction of trapezoidal map for tabulated lookup - * \author: D. Mayer, T. Economon - * \version 8.3.0 "Harrier" + * \brief Orders points by ascending x coordinates and updates triangle indices. */ -class CTrapezoidalMap { - protected: - /* The unique values of x which exist in the data */ - std::vector unique_bands_x; +inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { + const IntT n_pts = static_cast(x.size()); + + std::vector perm(n_pts); + std::iota(perm.begin(), perm.end(), 0); + std::sort(perm.begin(), perm.end(), + [&x, &y](const auto i, const auto j) { return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; }); + + auto reorder = [n_pts, &perm](const auto& v) { + VectorReal tmp(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + tmp[i] = v[perm[i]]; + } + return tmp; + }; + x = reorder(x); + y = reorder(y); + + std::vector inv_perm(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + inv_perm[perm[i]] = i; + } + for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { + for (IntT j = 0; j < 3; ++j) { + triangles(i, j) = inv_perm[triangles(i, j)]; + } + } +} + +/*! + * \brief Extracts unique edges from triangles. Edges are defined by two point IDs and + * up to two adjacent triangles (boundary edges have the second triangle ID < 0). + */ +inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, Matrix2i& edge_faces) { + std::vector> edges; + edges.resize(3 * triangles.rows()); + + for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { + for (IntT i = 0; i < 3; ++i) { + const IntT j = (i + 1) % 3; + const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); + const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); + edges[3 * i_tri + i] = {i_pt, j_pt, i_tri}; + } + } + + /*--- Sort to identify duplicates. ---*/ + std::sort(edges.begin(), edges.end(), + [](const auto& a, const auto& b) { return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); }); + + auto is_equal = [](const auto& a, const auto& b) { return a[0] == b[0] && a[1] == b[1]; }; + + IntT n_edges = 1; + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + n_edges += static_cast(!is_equal(edges[i], edges[i - 1])); + } + + edge_pts.resize(n_edges, 2); + edge_faces.resize(n_edges, 2); + IntT pos = 0; + + auto new_edge = [&](const auto& edge) { + edge_pts(pos, 0) = edge[0]; + edge_pts(pos, 1) = edge[1]; + edge_faces(pos, 0) = edge[2]; + edge_faces(pos, 1) = -1; + ++pos; + }; + + new_edge(edges[0]); + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + if (is_equal(edges[i], edges[i - 1])) { + edge_faces(pos - 1, 1) = edges[i][2]; + } else { + new_edge(edges[i]); + } + } +} + +/*! + * \brief Detects the x bands of the map. One band per unique x coordinate is used unless + * that would exceed max_bands, in which case equal-width bands are used to limit memory. + * \return Tuple of (n_bands, x_bands). + */ +inline auto DetectBands(const VectorReal& x, IntT max_bands = 0) { + if (max_bands <= 0) { + max_bands = std::min(IntT{5000}, static_cast(4.0 * std::sqrt(static_cast(x.size())))); + } + + IntT n_unique = 1; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != x[i - 1]) n_unique++; + } + + if (n_unique <= max_bands) { + const IntT n_bands = n_unique - 1; + VectorReal x_bands(n_unique); + IntT pos = 0; + x_bands[pos] = x[0]; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != x_bands[pos]) { + x_bands[++pos] = x[i]; + } + } + return std::make_tuple(n_bands, std::move(x_bands)); + } + + const RealT x_min = x.front(); + const RealT x_max = x.back(); + const RealT band_width = (x_max - x_min) / max_bands; + + VectorReal x_bands(max_bands + 1); + for (IntT i = 0; i <= max_bands; ++i) { + x_bands[i] = x_min + i * band_width; + } + x_bands[max_bands] = x_max; + + return std::make_tuple(max_bands, std::move(x_bands)); +} + +/*! + * \brief Builds the trapezoidal map for a set of edges (points must be ordered by x). + */ +inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, const VectorReal& y, + TrapezoidalMap& map) { + auto& x_bands = map.x_bands; + auto& offsets = map.offsets; + auto& edge_id = map.edge_id; + auto& edge_y = map.edge_y; + + auto clear_map = [&]() { + x_bands.clear(); + offsets.clear(); + edge_id.clear(); + edge_y.clear(); + }; + + const auto [n_bands, bands] = DetectBands(x); + x_bands = std::move(bands); + + if (n_bands <= 0) { + clear_map(); + return; + } + + auto find_band = [&x_bands, n_bands = n_bands](RealT x_val) -> IntT { + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); + const IntT idx = static_cast(it - x_bands.begin()); + return std::min(std::max(IntT{0}, idx - 1), n_bands - 1); + }; + + /*--- Count edges per band. Each edge is stored in every band between the bands of its + * two endpoints (inclusive), a superset of the bands it overlaps. ---*/ + auto& counts = offsets; + counts.clear(); + counts.resize(n_bands + 1, 0); + + for (IntT i = 0; i < static_cast(edge_pts.rows()); ++i) { + const IntT band_0 = find_band(x[edge_pts(i, 0)]); + const IntT band_1 = find_band(x[edge_pts(i, 1)]); + + for (IntT j = std::min(band_0, band_1); j <= std::max(band_0, band_1); ++j) { + ++counts[j + 1]; + } + } + + /*--- Convert counts to offsets (CSR format). ---*/ + for (IntT i = 2; i < static_cast(offsets.size()); ++i) { + offsets[i] += offsets[i - 1]; + } + + /*--- Give up (build failure) rather than allocating an excessive amount of memory. ---*/ + const size_t memory_mb = static_cast(offsets.back()) * (sizeof(IntT) + sizeof(RealT)) / (1024 * 1024); + if (memory_mb > 2048) { + clear_map(); + return; + } + + edge_id.resize(offsets.back()); + edge_y.resize(offsets.back()); + auto pos = offsets; + + for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { + const IntT pt_0 = edge_pts(i_edge, 0); + const IntT pt_1 = edge_pts(i_edge, 1); + const RealT x_0 = x[pt_0], y_0 = y[pt_0]; + const RealT x_1 = x[pt_1], y_1 = y[pt_1]; + + const IntT band_0 = find_band(x_0); + const IntT band_1 = find_band(x_1); + + const RealT dx = x_1 - x_0; + const bool vertical = std::abs(SU2_TYPE::GetValue(dx)) < 1e-30; + const RealT dy_dx = vertical ? RealT{0} : (y_1 - y_0) / dx; - su2activematrix edge_limits_x; - su2activematrix edge_limits_y; + for (IntT j = std::min(band_0, band_1); j <= std::max(band_0, band_1); ++j) { + edge_id[pos[j]] = i_edge; + const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2; + edge_y[pos[j]] = vertical ? (y_0 + y_1) / 2 : y_0 + dy_dx * (x_mid - x_0); + ++pos[j]; + } + } + + /*--- Sort the edges in each band by y coordinate. ---*/ + std::vector> tmp; + for (IntT i = 0; i < n_bands; ++i) { + const IntT begin = offsets[i]; + const IntT end = offsets[i + 1]; + if (begin >= end) continue; + + tmp.resize(end - begin); + for (auto k = begin; k < end; ++k) { + tmp[k - begin] = {edge_id[k], edge_y[k]}; + } + std::sort(tmp.begin(), tmp.end(), [](const auto& a, const auto& b) { return a.second < b.second; }); + for (auto k = begin; k < end; ++k) { + edge_id[k] = tmp[k - begin].first; + edge_y[k] = tmp[k - begin].second; + } + } +} + +/*! + * \brief Returns the IDs of the edges directly below and above a query point + * (either ID can be -1 if the point is at a boundary). + */ +inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, const VectorReal& x_coords, + const VectorReal& y_coords, const RealT& x, const RealT& y) { + if (map.x_bands.size() < 2 || map.offsets.empty()) { + return std::make_pair(IntT{-1}, IntT{-1}); + } + + const auto& x_bands = map.x_bands; + const IntT n_bands = static_cast(x_bands.size()) - 1; + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x); + const IntT d = static_cast(it - x_bands.begin()); + const IntT band_idx = std::min(std::max(IntT{0}, d - 1), n_bands - 1); + + RealT best_y_below = -1e300; + RealT best_y_above = 1e300; + IntT edge_below = -1; + IntT edge_above = -1; + + const IntT begin = map.offsets[band_idx]; + const IntT end = map.offsets[band_idx + 1]; + + for (IntT k = begin; k < end; ++k) { + const IntT e_id = map.edge_id[k]; + + const IntT p0 = edge_pts(e_id, 0); + const IntT p1 = edge_pts(e_id, 1); + const RealT x0 = x_coords[p0], y0 = y_coords[p0]; + const RealT x1 = x_coords[p1], y1 = y_coords[p1]; + + if (x < std::min(x0, x1) - 1e-10 || x > std::max(x0, x1) + 1e-10) { + continue; + } + + /*--- y position of the edge at the query x. ---*/ + RealT edge_y_at_x; + const RealT dx = x1 - x0; + if (std::abs(SU2_TYPE::GetValue(dx)) < 1e-30) { + edge_y_at_x = (y0 + y1) / 2.0; + } else { + edge_y_at_x = y0 + (x - x0) / dx * (y1 - y0); + } + + if (edge_y_at_x <= y + 1e-10 && edge_y_at_x > best_y_below) { + best_y_below = edge_y_at_x; + edge_below = e_id; + } + if (edge_y_at_x >= y - 1e-10 && edge_y_at_x < best_y_above) { + best_y_above = edge_y_at_x; + edge_above = e_id; + } + } + + return std::make_pair(edge_below, edge_above); +} + +/*! + * \brief Returns the IDs of the triangles adjacent to two query edges (up to 3 triangles). + */ +inline auto AdjacentTriangles(const IntT edge_0, const IntT edge_1, const Matrix2i& edge_faces) { + std::array tris = {-1, -1, -1}; + IntT pos = 0; + + auto insert = [&tris, &pos](const IntT t) { + if (t < 0) return; + for (IntT i = 0; i < pos; ++i) { + if (t == tris[i]) return; + } + tris[pos++] = t; + }; + + auto get_tris = [&edge_faces](const IntT e) { + if (e < 0) return std::array{IntT{-1}, IntT{-1}}; + return std::array{edge_faces(e, 0), edge_faces(e, 1)}; + }; + + for (const auto e : {edge_0, edge_1}) { + for (const auto t : get_tris(e)) { + insert(t); + } + } + return tris; +} + +/*! + * \brief Computes the barycentric coordinates of point (x_q, y_q) in a triangle. + */ +inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, const VectorReal& x, const VectorReal& y, + const RealT x_q, const RealT y_q) { + const IntT p0 = triangles(i_tri, 0); + const IntT p1 = triangles(i_tri, 1); + const IntT p2 = triangles(i_tri, 2); - su2vector > edge_to_triangle; + const RealT x0 = x[p0], y0 = y[p0]; + const RealT x1 = x[p1], y1 = y[p1]; + const RealT x2 = x[p2], y2 = y[p2]; + + const RealT dx1 = x1 - x0, dy1 = y1 - y0; + const RealT dx2 = x2 - x0, dy2 = y2 - y0; + + auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { return ux * vy - uy * vx; }; + + const RealT det = cross(dx1, dy1, dx2, dy2); + if (std::abs(SU2_TYPE::GetValue(det)) < 1e-30) { + return std::array{RealT{0}, RealT{0}, RealT{0}}; + } - /* The value that each edge which intersects the band takes within that - * same band. Used to sort the edges */ - su2vector > > y_edge_at_band_mid; + const RealT inv_det = 1.0 / det; + const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; + const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; - double memory_footprint = 0; + return std::array{1 - a - b, a, b}; +} + +/*! + * \brief Checks if a point is inside a triangle based on its barycentric coordinates. + */ +inline bool InTriangle(const std::array& coords, const RealT tol = 0.0) { + return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; +} + +/*! + * \brief Finds the triangle containing a point using the trapezoidal map. + */ +inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, const Matrix2i& edge_pts, + const Matrix2i& edge_faces, const VectorReal& x, const VectorReal& y, const RealT x_q, + const RealT y_q, std::array& bary_out) { + const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); + const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); + + const RealT tol = 1e-12; + for (const auto t : candidates) { + if (t < 0) continue; + + const auto coords = TriangleCoords(t, triangles, x, y, x_q, y_q); + if (InTriangle(coords, tol)) { + bary_out = coords; + return t; + } + } + + bary_out = {0.0, 0.0, 0.0}; + return -1; +} + +} // namespace su2_lut + +/*! + * \class CTrapezoidalMap + * \ingroup LookUpInterp + * \brief Trapezoidal map for finding the triangle containing a query point in a 2D triangulation. + */ +class CTrapezoidalMap { + private: + su2_lut::Matrix3i triangles; + su2_lut::Matrix2i edge_pts, edge_faces; + su2_lut::VectorReal x_coords, y_coords; + su2_lut::TrapezoidalMap map; + + unsigned long n_points = 0; + unsigned long n_triangles = 0; public: CTrapezoidalMap() = default; - CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size, - const std::vector >& edges, - const su2vector >& edge_to_triangle, bool display = false); - /*! - * \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 + * \brief Build the trapezoidal map from a triangulation. + * \return True on success. */ - unsigned long GetTriangle(const su2double val_x, const su2double val_y); + bool Build(unsigned long num_points, unsigned long num_triangles, const su2double* x, const su2double* y, + const unsigned long* connectivity) { + n_points = num_points; + n_triangles = num_triangles; - /*! - * \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(const su2double val_x); + if (num_points == 0 || num_triangles == 0) return false; - /*! - * \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; + x_coords.assign(x, x + num_points); + y_coords.assign(y, y + num_points); + + triangles.resize(num_triangles, 3); + for (size_t i = 0; i < 3 * num_triangles; ++i) { + triangles.data[i] = static_cast(connectivity[i]); + } + + su2_lut::ReorderPoints(triangles, x_coords, y_coords); + su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); + su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); + + return !map.x_bands.empty() && !map.offsets.empty(); + } /*! - * \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) + * \brief Find the triangle containing a query point. + * \return True if the point is inside the triangulation. */ - inline bool IsInsideHullX(su2double val_x) { - return (val_x >= unique_bands_x.front()) && (val_x <= unique_bands_x.back()); + bool FindTriangle(su2double val_x, su2double val_y, unsigned long& triangle_id, + std::array& bary_coords) const { + if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { + bary_coords = {0.0, 0.0, 0.0}; + return false; + } + + std::array bary; + const su2_lut::IntT tri_id = + su2_lut::FindTriangle(map, triangles, edge_pts, edge_faces, x_coords, y_coords, val_x, val_y, bary); + + if (tri_id < 0) return false; + + triangle_id = static_cast(tri_id); + bary_coords = {bary[0], bary[1], bary[2]}; + return true; } /*! - * \brief get memory footprint of trapezoidal map. - * \return - memory footprint in mega bytes. + * \brief Get the memory footprint of the map in MB. */ - double GetMemoryFootprint() const { return memory_footprint; } + double GetMemoryFootprint() const { + const size_t bytes = + (map.edge_id.size() + map.offsets.size() + edge_pts.data.size() + edge_faces.data.size() + + triangles.data.size()) * + sizeof(su2_lut::IntT) + + (map.edge_y.size() + map.x_bands.size() + x_coords.size() + y_coords.size()) * sizeof(su2_lut::RealT); + return double(bytes) / (1024.0 * 1024.0); + } }; diff --git a/Common/include/containers/CTrapezoidalMapFAST.hpp b/Common/include/containers/CTrapezoidalMapFAST.hpp deleted file mode 100644 index 7b18f4919ac1..000000000000 --- a/Common/include/containers/CTrapezoidalMapFAST.hpp +++ /dev/null @@ -1,617 +0,0 @@ -/*! - * \file CTrapezoidalMapFAST.hpp - * \brief Fast trapezoidal map implementation for 2D LUT queries in SU2. - * Based on Pedro Gomes' (pcarruscag) LUT implementation: - * https://github.com/pcarruscag/LUT - */ - -#ifndef CTRAPEZOIDALMAP_FAST_HPP -#define CTRAPEZOIDALMAP_FAST_HPP - -#pragma once - -#include -#include -#include -#include -#include -#include -#include - -#include "../basic_types/datatype_structure.hpp" - -namespace su2_lut { - -using IntT = int32_t; -using RealT = su2double; -using PassiveReal = passivedouble; - -/*--- Simple Matrix container (row-major) ---*/ -template -struct Matrix { - std::vector data; - - void resize(size_t rows, size_t) { data.resize(rows * N); } - size_t rows() const { return data.size() / N; } - - const T& operator()(size_t i, size_t j) const { return data[i * N + j]; } - T& operator()(size_t i, size_t j) { return data[i * N + j]; } -}; - -using Matrix2i = Matrix; -using Matrix3i = Matrix; -using VectorInt = std::vector; -using VectorReal = std::vector; - -/*--- Trapezoidal Map structure ---*/ -struct TrapezoidalMap { - VectorInt offsets, edge_id; - VectorReal x_bands, edge_y; -}; - -/*! - * \brief Orders points by ascending x coordinates and updates triangle indices. - */ -inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { - const IntT n_pts = static_cast(x.size()); - - // Create sort permutation based on (x, y) coordinates - std::vector perm(n_pts); - std::iota(perm.begin(), perm.end(), 0); - std::sort(perm.begin(), perm.end(), - [&x, &y](const auto i, const auto j) { return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; }); - - // Reorder coordinate arrays - auto reorder = [n_pts, &perm](const auto& v) { - VectorReal tmp(n_pts); - for (IntT i = 0; i < n_pts; ++i) { - tmp[i] = v[perm[i]]; - } - return tmp; - }; - x = reorder(x); - y = reorder(y); - - // Build inverse permutation to update triangle indices - std::vector inv_perm(n_pts); - for (IntT i = 0; i < n_pts; ++i) { - inv_perm[perm[i]] = i; - } - - // Update triangle connectivity - for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { - for (IntT j = 0; j < 3; ++j) { - triangles(i, j) = inv_perm[triangles(i, j)]; - } - } -} - -/*! - * \brief Check if two edges are equal (same point IDs). - */ -inline bool EdgesEqual(const std::array& a, const std::array& b) { - return a[0] == b[0] && a[1] == b[1]; -} - -/*! - * \brief Extract unique edges from triangles. - */ -inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, Matrix2i& edge_faces) { - // Extract edges from triangles - std::vector> edges; - edges.resize(3 * triangles.rows()); - - for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { - for (IntT i = 0; i < 3; ++i) { - const IntT j = (i + 1) % 3; - const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); - const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); - edges[3 * i_tri + i] = {i_pt, j_pt, i_tri}; - } - } - - // Sort to identify duplicates - std::sort(edges.begin(), edges.end(), - [](const auto& a, const auto& b) { return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); }); - - // Count unique edges - IntT n_edges = 1; - for (IntT i = 1; i < static_cast(edges.size()); ++i) { - n_edges += static_cast(!EdgesEqual(edges[i], edges[i - 1])); - } - - // Map edge points and edge faces - edge_pts.resize(n_edges, 2); - edge_faces.resize(n_edges, 2); - IntT pos = 0; - - auto new_edge = [&](const auto& edge) { - edge_pts(pos, 0) = edge[0]; - edge_pts(pos, 1) = edge[1]; - edge_faces(pos, 0) = edge[2]; - edge_faces(pos, 1) = -1; - ++pos; - }; - - new_edge(edges[0]); - for (IntT i = 1; i < static_cast(edges.size()); ++i) { - if (EdgesEqual(edges[i], edges[i - 1])) { - edge_faces(pos - 1, 1) = edges[i][2]; - } else { - new_edge(edges[i]); - } - } -} - -/*! - * \brief Detects bands based on unique x-coordinates. - * Uses coarse banding to limit memory for large tables. - * Returns tuple of (n_bands, pt_to_band, x_bands). - */ -inline auto DetectBands(const VectorReal& x, IntT max_bands = 0) { - // Safety check: empty input - if (x.empty()) { - return std::make_tuple(IntT{0}, std::vector{}, VectorReal{}); - } - - // First pass: count unique x-coordinates - IntT n_unique = 1; - for (IntT i = 1; i < static_cast(x.size()); ++i) { - if (x[i] != x[i - 1]) n_unique++; - } - - // Determine if we need to use coarser bands - if (max_bands <= 0) { - max_bands = std::min(IntT{5000}, static_cast(4.0 * std::sqrt(static_cast(x.size())))); - } - - const bool use_coarse_bands = (n_unique > max_bands); - - if (!use_coarse_bands) { - // Original algorithm: use all unique x-coordinates as band boundaries - std::vector pt_to_band(x.size()); - pt_to_band[0] = 0; - IntT n_bands = 0; - - for (IntT i = 1; i < static_cast(x.size()); ++i) { - n_bands += static_cast(x[i] != x[i - 1]); - pt_to_band[i] = n_bands; - } - - VectorReal x_bands(n_bands + 1); - IntT pos = 0; - x_bands[pos] = x[0]; - for (IntT i = 1; i < static_cast(x.size()); ++i) { - if (x[i] != x_bands[pos]) { - x_bands[++pos] = x[i]; - } - } - - return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); - } else { - // Coarse banding: divide x-range into max_bands equal intervals - const RealT x_min = x.front(); - const RealT x_max = x.back(); - const RealT x_range = x_max - x_min; - - // Handle degenerate case where all x are the same - if (x_range < 1e-30) { - std::vector pt_to_band(x.size(), 0); - VectorReal x_bands = {x_min, x_max}; - return std::make_tuple(IntT{1}, std::move(pt_to_band), std::move(x_bands)); - } - - const IntT n_bands = max_bands; - const RealT band_width = x_range / max_bands; - - // Build band boundaries - VectorReal x_bands(max_bands + 1); - for (IntT i = 0; i <= max_bands; ++i) { - x_bands[i] = x_min + i * band_width; - } - x_bands[max_bands] = x_max; - - // Map each point to its band - std::vector pt_to_band(x.size()); - for (IntT i = 0; i < static_cast(x.size()); ++i) { - IntT band = static_cast(SU2_TYPE::GetValue((x[i] - x_min) / band_width)); - band = std::max(IntT{0}, std::min(band, max_bands - 1)); - pt_to_band[i] = band; - } - - return std::make_tuple(n_bands, std::move(pt_to_band), std::move(x_bands)); - } -} - -/*! - * \brief Builds the trapezoidal map for efficient 2D queries. - */ -inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, const VectorReal& y, - TrapezoidalMap& map) { - auto& x_bands = map.x_bands; - auto& offsets = map.offsets; - auto& edge_id = map.edge_id; - auto& edge_y = map.edge_y; - - const auto [n_bands, pt_to_band, bands] = DetectBands(x); - x_bands = std::move(bands); - - // Safety check - if (n_bands <= 0 || x_bands.size() < 2) { - x_bands.clear(); - offsets.clear(); - edge_id.clear(); - edge_y.clear(); - return; - } - - // Helper function to find band index for an x-coordinate - auto find_band = [&x_bands, n_bands](RealT x_val) -> IntT { - auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); - IntT idx = static_cast(it - x_bands.begin()); - idx = std::max(IntT{0}, idx - 1); - idx = std::min(idx, n_bands - 1); - return idx; - }; - - // Count edges per band - auto& counts = offsets; - counts.clear(); - counts.resize(n_bands + 1, 0); - - for (IntT i = 0; i < static_cast(edge_pts.rows()); ++i) { - const IntT pt_0 = edge_pts(i, 0); - const IntT pt_1 = edge_pts(i, 1); - const RealT x_0 = x[pt_0]; - const RealT x_1 = x[pt_1]; - - const IntT band_0 = find_band(x_0); - const IntT band_1 = find_band(x_1); - - const IntT min_band = std::min(band_0, band_1); - const IntT max_band = std::max(band_0, band_1); - - for (IntT j = min_band; j <= max_band && j < n_bands; ++j) { - ++counts[j + 1]; - } - } - - // Convert counts to offsets (CSR format) - for (IntT i = 2; i < static_cast(offsets.size()); ++i) { - offsets[i] += offsets[i - 1]; - } - - // Check memory requirement - const size_t total_entries = static_cast(offsets.back()); - const size_t memory_mb = total_entries * (sizeof(IntT) + sizeof(RealT)) / (1024 * 1024); - - if (memory_mb > 2048) { - x_bands.clear(); - offsets.clear(); - edge_id.clear(); - edge_y.clear(); - return; - } - - // Allocate storage - edge_id.resize(offsets.back()); - edge_y.resize(offsets.back()); - auto pos = offsets; - - // Store edges in each band - for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { - const IntT pt_0 = edge_pts(i_edge, 0); - const IntT pt_1 = edge_pts(i_edge, 1); - const RealT x_0 = x[pt_0], y_0 = y[pt_0]; - const RealT x_1 = x[pt_1], y_1 = y[pt_1]; - - const IntT band_0 = find_band(x_0); - const IntT band_1 = find_band(x_1); - - if (band_0 == band_1) continue; - - const IntT min_band = std::min(band_0, band_1); - const IntT max_band = std::max(band_0, band_1); - const RealT dy_dx = (y_1 - y_0) / (x_1 - x_0); - - for (IntT j = min_band; j < max_band && j < n_bands; ++j) { - if (pos[j] < static_cast(edge_id.size())) { - edge_id[pos[j]] = i_edge; - const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2; - edge_y[pos[j]] = y_0 + dy_dx * (x_mid - x_0); - ++pos[j]; - } - } - } - - // Sort edges in each band by y coordinate - std::vector> tmp; - for (IntT i = 0; i < n_bands; ++i) { - const IntT begin = offsets[i]; - const IntT end = offsets[i + 1]; - if (begin >= end) continue; - - tmp.resize(end - begin); - for (auto k = begin; k < end; ++k) { - tmp[k - begin] = {edge_id[k], edge_y[k]}; - } - std::sort(tmp.begin(), tmp.end(), [](const auto& a, const auto& b) { return a.second < b.second; }); - for (auto k = begin; k < end; ++k) { - edge_id[k] = tmp[k - begin].first; - edge_y[k] = tmp[k - begin].second; - } - } -} - -/*! - * \brief Query the trapezoidal map for edges bounding a point. - * Returns pair of (edge_below, edge_above). Either can be -1 if at boundary. - */ -inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, const VectorReal& x_coords, - const VectorReal& y_coords, const RealT& x, const RealT& y) { - // Safety check - if (map.x_bands.size() < 2 || map.offsets.empty()) { - return std::make_pair(IntT{-1}, IntT{-1}); - } - - const IntT n_bands = static_cast(map.x_bands.size()) - 1; - - // Find the band containing x - const auto& x_bands = map.x_bands; - auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x); - IntT d = static_cast(it - x_bands.begin()); - IntT center_band = std::min(std::max(IntT{0}, d - 1), n_bands - 1); - - // Search in current band AND adjacent bands - RealT best_y_below = -1e300; - RealT best_y_above = 1e300; - IntT edge_below = -1; - IntT edge_above = -1; - - for (IntT band_offset = -1; band_offset <= 1; ++band_offset) { - IntT band_idx = center_band + band_offset; - if (band_idx < 0 || band_idx >= n_bands) continue; - if (band_idx + 1 >= static_cast(map.offsets.size())) continue; - - const IntT begin = map.offsets[band_idx]; - const IntT end = map.offsets[band_idx + 1]; - - if (begin < 0 || end > static_cast(map.edge_id.size()) || begin >= end) { - continue; - } - - for (IntT k = begin; k < end; ++k) { - const IntT e_id = map.edge_id[k]; - - const IntT p0 = edge_pts(e_id, 0); - const IntT p1 = edge_pts(e_id, 1); - const RealT x0 = x_coords[p0], y0 = y_coords[p0]; - const RealT x1 = x_coords[p1], y1 = y_coords[p1]; - - // Check if query x is within edge's x-extent - if (x < std::min(x0, x1) - 1e-10 || x > std::max(x0, x1) + 1e-10) { - continue; - } - - // Compute y-position of edge at query x - RealT edge_y_at_x; - const RealT dx = x1 - x0; - if (std::abs(SU2_TYPE::GetValue(dx)) < 1e-30) { - edge_y_at_x = (y0 + y1) / 2.0; - } else { - const RealT t = (x - x0) / dx; - edge_y_at_x = y0 + t * (y1 - y0); - } - - if (edge_y_at_x <= y + 1e-10 && edge_y_at_x > best_y_below) { - best_y_below = edge_y_at_x; - edge_below = e_id; - } - if (edge_y_at_x >= y - 1e-10 && edge_y_at_x < best_y_above) { - best_y_above = edge_y_at_x; - edge_above = e_id; - } - } - } - - return std::make_pair(edge_below, edge_above); -} - -/*! - * \brief Get triangles adjacent to two query edges (up to 3 triangles). - */ -inline auto AdjacentTriangles(const IntT edge_0, const IntT edge_1, const Matrix2i& edge_faces) { - std::array tris = {-1, -1, -1}; - IntT pos = 0; - - auto insert = [&tris, &pos](const IntT t) { - if (t < 0) return; - for (IntT i = 0; i < pos; ++i) { - if (t == tris[i]) return; - } - tris[pos++] = t; - }; - - auto get_tris = [&edge_faces](const IntT e) { - if (e < 0) return std::array{IntT{-1}, IntT{-1}}; - return std::array{edge_faces(e, 0), edge_faces(e, 1)}; - }; - - for (const auto e : {edge_0, edge_1}) { - for (const auto t : get_tris(e)) { - insert(t); - } - } - return tris; -} - -/*! - * \brief Compute barycentric coordinates of point (x, y) in triangle. - */ -inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, const VectorReal& x, const VectorReal& y, - const RealT x_q, const RealT y_q) { - const IntT p0 = triangles(i_tri, 0); - const IntT p1 = triangles(i_tri, 1); - const IntT p2 = triangles(i_tri, 2); - - const RealT x0 = x[p0], y0 = y[p0]; - const RealT x1 = x[p1], y1 = y[p1]; - const RealT x2 = x[p2], y2 = y[p2]; - - const RealT dx1 = x1 - x0, dy1 = y1 - y0; - const RealT dx2 = x2 - x0, dy2 = y2 - y0; - - auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { return ux * vy - uy * vx; }; - - const RealT det = cross(dx1, dy1, dx2, dy2); - if (std::abs(SU2_TYPE::GetValue(det)) < 1e-30) { - return std::array{RealT{0}, RealT{0}, RealT{0}}; - } - - const RealT inv_det = 1.0 / det; - const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; - const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; - - return std::array{1 - a - b, a, b}; -} - -/*! - * \brief Check if point is inside triangle. - */ -inline bool InTriangle(const std::array& coords, const RealT tol = 0.0) { - return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; -} - -/*! - * \brief Find triangle containing point using the trapezoidal map. - */ -inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, const Matrix2i& edge_pts, - const Matrix2i& edge_faces, const VectorReal& x, const VectorReal& y, const RealT x_q, - const RealT y_q, std::array& bary_out) { - const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); - const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); - - const RealT tol = 1e-12; - for (const auto t : candidates) { - if (t < 0) continue; - - const auto coords = TriangleCoords(t, triangles, x, y, x_q, y_q); - if (InTriangle(coords, tol)) { - bary_out = coords; - return t; - } - } - - bary_out = {0.0, 0.0, 0.0}; - return -1; -} - -} // namespace su2_lut - -/*! - * \class CTrapezoidalMapFAST - * \brief Wrapper class for SU2 integration. - */ -class CTrapezoidalMapFAST { - private: - su2_lut::Matrix3i triangles; - su2_lut::Matrix2i edge_pts, edge_faces; - su2_lut::VectorReal x_coords, y_coords; - su2_lut::TrapezoidalMap map; - - unsigned long n_points; - unsigned long n_triangles; - bool build_successful; - - public: - CTrapezoidalMapFAST() : n_points(0), n_triangles(0), build_successful(false) {} - - /*! - * \brief Build the trapezoidal map from mesh data. - */ - bool Build(unsigned long num_points, unsigned long num_triangles, const su2double* x, const su2double* y, - const unsigned long* connectivity) { - n_points = num_points; - n_triangles = num_triangles; - build_successful = false; - - if (num_points == 0 || num_triangles == 0) { - x_coords.clear(); - y_coords.clear(); - triangles.resize(0, 3); - edge_pts.resize(0, 2); - edge_faces.resize(0, 2); - map.x_bands.clear(); - map.offsets.clear(); - map.edge_id.clear(); - map.edge_y.clear(); - return false; - } - - // Copy coordinates - x_coords.assign(x, x + num_points); - y_coords.assign(y, y + num_points); - - // Copy triangle connectivity - triangles.resize(num_triangles, 3); - for (size_t i = 0; i < num_triangles; ++i) { - unsigned long v0 = connectivity[3 * i + 0]; - unsigned long v1 = connectivity[3 * i + 1]; - unsigned long v2 = connectivity[3 * i + 2]; - - if (v0 >= num_points || v1 >= num_points || v2 >= num_points) { - triangles(i, 0) = 0; - triangles(i, 1) = 0; - triangles(i, 2) = 0; - } else { - triangles(i, 0) = static_cast(v0); - triangles(i, 1) = static_cast(v1); - triangles(i, 2) = static_cast(v2); - } - } - - // Build the map - su2_lut::ReorderPoints(triangles, x_coords, y_coords); - su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); - su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); - - if (map.x_bands.empty() || map.offsets.empty()) { - build_successful = false; - return false; - } - - build_successful = true; - return true; - } - - bool IsBuilt() const { return build_successful; } - - /*! - * \brief Find the triangle containing a query point. - */ - bool FindTriangle(su2double val_x, su2double val_y, unsigned long& triangle_id, - std::array& bary_coords) const { - if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { - bary_coords = {0.0, 0.0, 0.0}; - return false; - } - - std::array bary; - su2_lut::IntT tri_id = - su2_lut::FindTriangle(map, triangles, edge_pts, edge_faces, x_coords, y_coords, val_x, val_y, bary); - - if (tri_id >= 0) { - triangle_id = static_cast(tri_id); - bary_coords = {bary[0], bary[1], bary[2]}; - return true; - } - return false; - } - - unsigned long GetNTriangles() const { return n_triangles; } - unsigned long GetNPoints() const { return n_points; } - unsigned long GetNEdges() const { return edge_pts.rows(); } - unsigned long GetNBands() const { return map.x_bands.size() > 0 ? map.x_bands.size() - 1 : 0; } -}; - -#endif // CTRAPEZOIDALMAP_FAST_HPP diff --git a/Common/include/option_structure.hpp b/Common/include/option_structure.hpp index 2e24be709da1..202062fd243d 100644 --- a/Common/include/option_structure.hpp +++ b/Common/include/option_structure.hpp @@ -599,14 +599,12 @@ MakePair("ONESPECIES", ONESPECIES) */ enum class ENUM_DATADRIVEN_METHOD { LUT = 0, - MLP = 1, - LUT_FAST = 2 /* New LUT implementation based on Pedro's code */ + MLP = 1 }; static const MapType DataDrivenMethod_Map = { MakePair("LUT", ENUM_DATADRIVEN_METHOD::LUT) MakePair("MLP", ENUM_DATADRIVEN_METHOD::MLP) - MakePair("LUT_FAST", ENUM_DATADRIVEN_METHOD::LUT_FAST) }; /*! diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 005ea10d17c6..2b4a6bb41dde 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -57,19 +57,6 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, /* Add additional variable index which will always result in zero when looked up. */ idx_null = names_var.size(); - /*--- NOTE: Trapezoidal map is NOT built here anymore! ---*/ - /*--- Call EnableFastTrapMap() for LUT_FAST (efficient O(n) memory) ---*/ - /*--- Call BuildOriginalTrapMap() for standard LUT (original SU2 implementation) ---*/ - - ComputeInterpCoeffs(); - - if (rank == MASTER_NODE) cout << "LUT fluid model ready for use" << endl; -} - -void CLookUpTable::BuildOriginalTrapMap() { - /*--- Build the original SU2 trapezoidal map (DAG-based, O(n log n) memory) ---*/ - /*--- WARNING: This can use a lot of memory for large tables! ---*/ - if (rank == MASTER_NODE) switch (table_dim) { case 2: cout << "Building a trapezoidal map for the (" + name_CV1 + ", " + name_CV2 + ") space ..." << endl; @@ -83,105 +70,43 @@ void CLookUpTable::BuildOriginalTrapMap() { trap_map_x_y.resize(n_table_levels); su2double startTime = SU2_MPI::Wtime(); - unsigned short barwidth = 65; - bool display_map_info = (n_table_levels < 2); double tmap_memory_footprint = 0; for (auto i_level = 0ul; i_level < n_table_levels; i_level++) { - trap_map_x_y[i_level] = - CTrapezoidalMap(GetDataP(name_CV1, i_level), GetDataP(name_CV2, i_level), table_data[i_level].cols(), - edges[i_level], edge_to_triangle[i_level], display_map_info); - tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint(); - - /* Display a progress bar to monitor table generation process */ - if (rank == MASTER_NODE) { - su2double progress = su2double(i_level) / n_table_levels; - auto completed = floor(progress * barwidth); - auto to_do = barwidth - completed; - cout << "[" << setfill('=') << setw(completed); - cout << '>'; - cout << setfill(' ') << setw(to_do) << std::right << "] " << 100 * progress << "%\r"; - cout.flush(); - } - } - su2double stopTime = SU2_MPI::Wtime(); - - if (rank == MASTER_NODE) { - switch (table_dim) { - case 2: - cout << "\nConstruction of trapezoidal map took " << stopTime - startTime << " seconds\n" << endl; - break; - case 3: - cout << "\nConstruction of trapezoidal map stack took " << stopTime - startTime << " seconds\n" << endl; - break; - default: - break; - } - cout << "Trapezoidal map memory footprint: " << tmap_memory_footprint << " MB\n"; - cout << "Table data memory footprint: " << memory_footprint_data << " MB\n" << endl; - } -} - -void CLookUpTable::EnableFastTrapMap() { - /*--- Build Memory-efficient trapezoidal map (band-based, O(n) memory) ---*/ - - use_fast_trap_map_ = true; - - trap_map_x_y_FAST.resize(n_table_levels); - - if (rank == MASTER_NODE) { - cout << endl; - cout << "+--------------------------------------------------+" << endl; - cout << "| Building Memory-efficient trapezoidal maps |" << endl; - cout << "| (O(n) memory, O(log n) queries) |" << endl; - cout << "+--------------------------------------------------+" << endl; - } - - su2double startTime = SU2_MPI::Wtime(); - - for (unsigned long i_level = 0; i_level < n_table_levels; ++i_level) { const auto n_pts = n_points[i_level]; const auto n_tris = n_triangles[i_level]; - if (rank == MASTER_NODE) { - cout << " Level " << i_level << ": " << n_pts << " points, " << n_tris << " triangles ... " << flush; - } - - // Get point coordinates std::vector x_coords(n_pts); std::vector y_coords(n_pts); - - for (unsigned long i_point = 0; i_point < n_pts; ++i_point) { + for (auto i_point = 0ul; i_point < n_pts; ++i_point) { x_coords[i_point] = table_data[i_level][idx_CV1][i_point]; y_coords[i_point] = table_data[i_level][idx_CV2][i_point]; } - // Get triangle connectivity std::vector tri_conn(3 * n_tris); - for (unsigned long i_tri = 0; i_tri < n_tris; ++i_tri) { + for (auto i_tri = 0ul; i_tri < n_tris; ++i_tri) { tri_conn[3 * i_tri + 0] = triangles[i_level][i_tri][0]; tri_conn[3 * i_tri + 1] = triangles[i_level][i_tri][1]; tri_conn[3 * i_tri + 2] = triangles[i_level][i_tri][2]; } - // Build FAST trapezoidal map for this level - trap_map_x_y_FAST[i_level].Build(n_pts, n_tris, x_coords.data(), y_coords.data(), tri_conn.data()); - - if (rank == MASTER_NODE) { - cout << "done." << endl; - } + if (!trap_map_x_y[i_level].Build(n_pts, n_tris, x_coords.data(), y_coords.data(), tri_conn.data())) + SU2_MPI::Error( + "Construction of trapezoidal map failed for level " + std::to_string(i_level) + " of table " + file_name_lut, + CURRENT_FUNCTION); + tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint(); } - su2double stopTime = SU2_MPI::Wtime(); if (rank == MASTER_NODE) { - cout << "+--------------------------------------------------+" << endl; - cout << "| FAST trapezoidal maps built successfully! |" << endl; - cout << "| Build time: " << setw(10) << setprecision(3) << (stopTime - startTime) << " seconds |" - << endl; - cout << "+--------------------------------------------------+" << endl; - cout << endl; + cout << "Construction of trapezoidal map took " << stopTime - startTime << " seconds\n"; + cout << "Trapezoidal map memory footprint: " << tmap_memory_footprint << " MB\n"; + cout << "Table data memory footprint: " << memory_footprint_data << " MB\n" << endl; } + + ComputeInterpCoeffs(); + + if (rank == MASTER_NODE) cout << "LUT fluid model ready for use" << endl; } void CLookUpTable::LoadTableRaw(const string& var_file_name_lut) { @@ -681,31 +606,8 @@ bool CLookUpTable::LookUp_XY(const vector& idx_var, vector bary_coords; - if (trap_map_x_y_FAST[iLevel].FindTriangle(val_CV1, val_CV2, id_triangle, bary_coords)) { - return true; - } - // Fall through to bounds check if not found (point may be outside domain) - } - - /*--- Use original trapezoidal map if available ---*/ - if (!use_fast_trap_map_ && iLevel < trap_map_x_y.size() && trap_map_x_y.size() > 0) { - /* 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[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[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) */ - return IsInTriangle(val_CV1, val_CV2, id_triangle, iLevel); - } - } - - return false; + std::array bary_coords; + return trap_map_x_y[iLevel].FindTriangle(val_CV1, val_CV2, id_triangle, bary_coords); } void CLookUpTable::GetInterpCoeffs(su2double val_CV1, su2double val_CV2, const su2activematrix& interp_mat_inv, @@ -860,26 +762,6 @@ void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const InterpolateToNearestNeighbors(val_CV1, val_CV2, names_var, val_names_var, i_level); } -bool CLookUpTable::IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, - unsigned long i_level) { - 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 = 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 = 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); - - su2double area_0 = TriArea(val_CV1, val_CV2, tri_x_1, tri_y_1, tri_x_2, tri_y_2); - su2double area_1 = TriArea(tri_x_0, tri_y_0, val_CV1, val_CV2, tri_x_2, tri_y_2); - su2double area_2 = TriArea(tri_x_0, tri_y_0, tri_x_1, tri_y_1, val_CV1, val_CV2); - - return (abs(area_tri - (area_0 + area_1 + area_2)) < area_tri * 1e-10); -} - bool CLookUpTable::CheckForVariables(const std::vector& vars_to_check) const { for (const string& var_to_check : vars_to_check) { if (!std::any_of(names_var.begin(), names_var.end(), diff --git a/Common/src/containers/CTrapezoidalMap.cpp b/Common/src/containers/CTrapezoidalMap.cpp deleted file mode 100644 index ea0472848455..000000000000 --- a/Common/src/containers/CTrapezoidalMap.cpp +++ /dev/null @@ -1,283 +0,0 @@ -/*! - * \file CTrapezoidalMap.cpp - * \brief Implementation of the trapezoidal map for tabulation and lookup of fluid properties - * \author D. Mayer, T. Economon, N. Beishuizen - * \version 8.3.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2025, 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 -#include - -#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) - * NOTE: the current implementation is actually the simpler 'slab' approach. - */ -CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size, - vector > const& edges, - su2vector > const& val_edge_to_triangle, bool display) { - int rank = SU2_MPI::GetRank(); - su2double startTime = SU2_MPI::Wtime(); - - edge_to_triangle = su2vector >(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; - unsigned long j_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].emplace_back(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(); - - /* calculate size of trapezoidal map components */ - double size_unique_bands = sizeof(su2double) * unique_bands_x.size() / 1e6; - double size_edge_limits_x = sizeof(su2double) * edge_limits_x.size() * 2 / 1e6; - double size_edge_limits_y = sizeof(su2double) * edge_limits_y.size() * 2 / 1e6; - - double size_edge_to_triangle = 0; - for (i_edge = 0; i_edge < edge_to_triangle.size(); i_edge++) - for (j_edge = 0; j_edge < edge_to_triangle[i_edge].size(); j_edge++) - size_edge_to_triangle += sizeof(unsigned long) / 1e6; - - double size_y_edge_at_band_mid = 0; - for (unsigned long i_y = 0; i_y < y_edge_at_band_mid.size(); i_y++) - for (unsigned long j_y = 0; j_y < y_edge_at_band_mid[i_y].size(); j_y++) - size_y_edge_at_band_mid += sizeof(su2double) / 1e6 + sizeof(unsigned long) / 1e6; - - memory_footprint = - size_unique_bands + size_edge_limits_x + size_edge_limits_y + size_edge_to_triangle + size_y_edge_at_band_mid; - - /* print size of trapezoidal map components to screen */ - if ((rank == MASTER_NODE) && display) { - cout << setfill(' '); - cout << "\n" << endl; - cout << "+------------------------------------------------------------------+\n"; - cout << "| Trapezoidal map info |\n"; - cout << "+------------------------------------------------------------------+" << endl; - - cout << "| Time to construct trapezoidal map: " << setw(22) << right << stopTime - startTime << " sec" - << " |" << endl; - cout << "| Size of unique_bands in memory: " << setw(22) << size_unique_bands << " MB " - << " |" << endl; - cout << "| Size of edge_limits_x in memory: " << setw(22) << size_edge_limits_x << " MB " - << " |" << endl; - cout << "| Size of edge_limits_y in memory: " << setw(22) << size_edge_limits_y << " MB " - << " |" << endl; - cout << "| Size of edge_to_triangle in memory: " << setw(22) << size_edge_to_triangle << " MB " - << " |" << endl; - cout << "| Size of y_edge_at_band_mid in memory: " << setw(22) << size_y_edge_at_band_mid << " MB " - << " |" << endl; - cout << "| Total: " << setw(22) << memory_footprint << " MB " - << " |" << endl; - cout << "+------------------------------------------------------------------+" << endl; - cout << "\n" << endl; - } -} - -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); - - /* within that band, find edges which enclose the (val_x, val_y) point */ - pair edges = GetEdges(band, val_x, val_y); - - /* identify the adjacent triangles using the two edges */ - std::array triangles_edge_low; - for (unsigned long i = 0; i < edge_to_triangle[edges.first].size(); i++) - triangles_edge_low[i] = edge_to_triangle[edges.first][i]; - - std::array triangles_edge_up; - for (unsigned long i = 0; i < edge_to_triangle[edges.second].size(); 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; - set_intersection(triangles_edge_up.begin(), triangles_edge_up.end(), triangles_edge_low.begin(), - triangles_edge_low.end(), std::back_inserter(triangle)); - - /*--- We failed to find an intersection, so take the lower triangle inside the band enclosing the point---*/ - if (triangle.size() < 1) { - triangle.resize(1, triangles_edge_low[0]); - } - - return triangle[0]; -} - -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_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_sample); - - /*--- if upper bound = 0, then use the range [0,1] ---*/ - i_up = max(1, 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 index 4c8d4fe618a4..0743fd6b60c0 100644 --- a/Common/src/containers/meson.build +++ b/Common/src/containers/meson.build @@ -1,3 +1,2 @@ -common_src += files(['CTrapezoidalMap.cpp', - 'CFileReaderLUT.cpp', +common_src += files(['CFileReaderLUT.cpp', 'CLookUpTable.cpp']) diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index 12ae1a9dc54b..76322a641c2b 100644 --- a/SU2_CFD/src/fluid/CFluidFlamelet.cpp +++ b/SU2_CFD/src/fluid/CFluidFlamelet.cpp @@ -80,21 +80,6 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati } look_up_table = new CLookUpTable(datadriven_fluid_options.datadriven_filenames[0], table_scalar_names[I_PROGVAR], table_scalar_names[I_ENTH]); - /*--- Build the original trapezoidal map (can use lots of memory!) ---*/ - look_up_table->BuildOriginalTrapMap(); - break; - - case ENUM_DATADRIVEN_METHOD::LUT_FAST: - if (rank == MASTER_NODE) { - cout << "**************************************************" << endl; - cout << "*** initializing the FAST lookup table ***" << endl; - cout << "*** (Memory-efficient trapezoidal map) ***" << endl; - cout << "**************************************************" << endl; - } - look_up_table = new CLookUpTable(datadriven_fluid_options.datadriven_filenames[0], table_scalar_names[I_PROGVAR], - table_scalar_names[I_ENTH]); - /*--- Build Memory-efficient trapezoidal map (O(n) memory) ---*/ - look_up_table->EnableFastTrapMap(); break; default: @@ -123,10 +108,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati } CFluidFlamelet::~CFluidFlamelet() { - if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT || - Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT_FAST) - delete look_up_table; - + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT) delete look_up_table; #ifdef USE_MLPCPP if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { delete iomap_TD; @@ -247,10 +229,7 @@ void CFluidFlamelet::PreprocessLookUp(CConfig* config) { preferential_diffusion = flamelet_options.preferential_diffusion; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: - case ENUM_DATADRIVEN_METHOD::LUT_FAST: - if (preferential_diffusion) { - preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); - } + preferential_diffusion = look_up_table->CheckForVariables(varnames_PD); break; case ENUM_DATADRIVEN_METHOD::MLP: #ifdef USE_MLPCPP @@ -355,7 +334,6 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca bool inside; switch (Kind_DataDriven_Method) { case ENUM_DATADRIVEN_METHOD::LUT: - case ENUM_DATADRIVEN_METHOD::LUT_FAST: if (output_refs.size() != LUT_idx.size()) SU2_MPI::Error(string("Output vector size incompatible with manifold lookup operation."), CURRENT_FUNCTION); if (include_mixture_fraction) { diff --git a/UnitTests/Common/containers/CLookupTable_tests.cpp b/UnitTests/Common/containers/CLookupTable_tests.cpp index fa5668bdc6ac..4b9f06eb9469 100644 --- a/UnitTests/Common/containers/CLookupTable_tests.cpp +++ b/UnitTests/Common/containers/CLookupTable_tests.cpp @@ -35,118 +35,10 @@ #include "../../../Common/include/containers/CLookUpTable.hpp" #include "../../../Common/include/containers/CFileReaderLUT.hpp" -/*--- Original 2D LUT test (using default trapezoidal map) ---*/ TEST_CASE("LUTreader", "[tabulated chemistry]") { /*--- smaller and trivial lookup table ---*/ CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable.drg", "ProgressVariable", "EnthalpyTot"); - look_up_table.BuildOriginalTrapMap(); - /*--- string names of the controlling variables ---*/ - - string name_CV1 = "ProgressVariable"; - string name_CV2 = "EnthalpyTot"; - - /*--- look up a single value for density ---*/ - - 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(idx_tag, &look_up_dat, prog, enth); - CHECK(look_up_dat == Approx(1.02)); - - /*--- look up a single value for viscosity ---*/ - - prog = 0.6; - enth = 0.9; - look_up_tag = "Viscosity"; - 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 */ - - auto limitsEnth = look_up_table.GetTableLimitsY(); - CHECK(SU2_TYPE::GetValue(*limitsEnth.first) == Approx(-1.0)); - CHECK(SU2_TYPE::GetValue(*limitsEnth.second) == Approx(1.0)); - - auto limitsProgvar = look_up_table.GetTableLimitsX(); - CHECK(SU2_TYPE::GetValue(*limitsProgvar.first) == Approx(0.0)); - CHECK(SU2_TYPE::GetValue(*limitsProgvar.second) == Approx(1.0)); - - /* lookup value outside of lookup table */ - - prog = 1.10; - enth = 1.1; - look_up_tag = "Density"; - 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)); -} - -/*--- Original 3D LUT test (using default trapezoidal map) ---*/ -TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { - /*--- smaller and trivial lookup table ---*/ - - CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable_3D.drg", "ProgressVariable", - "EnthalpyTot"); - look_up_table.BuildOriginalTrapMap(); - /*--- string names of the controlling variables ---*/ - - string name_CV1 = "ProgressVariable"; - string name_CV2 = "EnthalpyTot"; - - /*--- look up a single value for density ---*/ - - su2double prog = 0.55; - 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(idx_tag, &look_up_dat, prog, enth, mfrac); - CHECK(look_up_dat == Approx(1.02)); - - /*--- look up a single value for viscosity ---*/ - - prog = 0.6; - enth = 0.9; - mfrac = 0.8; - look_up_tag = "Viscosity"; - 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 */ - - auto limitsEnth = look_up_table.GetTableLimitsY(); - CHECK(SU2_TYPE::GetValue(*limitsEnth.first) == Approx(-1.0)); - CHECK(SU2_TYPE::GetValue(*limitsEnth.second) == Approx(1.0)); - - auto limitsProgvar = look_up_table.GetTableLimitsX(); - CHECK(SU2_TYPE::GetValue(*limitsProgvar.first) == Approx(0.0)); - CHECK(SU2_TYPE::GetValue(*limitsProgvar.second) == Approx(1.0)); - - /* lookup value outside of lookup table */ - - prog = 1.10; - enth = 1.1; - mfrac = 2.0; - look_up_tag = "Density"; - 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)); -} - -/*--- 2D LUT test using Pedro's FAST trapezoidal map (LUT_FAST) ---*/ -TEST_CASE("LUTreader_FAST", "[tabulated chemistry]") { - /*--- smaller and trivial lookup table ---*/ - - CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable.drg", "ProgressVariable", "EnthalpyTot"); - - /*--- Enable Pedro's FAST trapezoidal map ---*/ - look_up_table.EnableFastTrapMap(); /*--- string names of the controlling variables ---*/ @@ -192,16 +84,12 @@ TEST_CASE("LUTreader_FAST", "[tabulated chemistry]") { CHECK(look_up_dat == Approx(1.1738796125)); } -/*--- 3D LUT test using Pedro's FAST trapezoidal map (LUT_FAST) ---*/ -TEST_CASE("LUTreader_3D_FAST", "[tabulated chemistry]") { +TEST_CASE("LUTreader_3D", "[tabulated chemistry]") { /*--- smaller and trivial lookup table ---*/ CLookUpTable look_up_table("src/SU2/UnitTests/Common/containers/lookuptable_3D.drg", "ProgressVariable", "EnthalpyTot"); - /*--- Enable Pedro's FAST trapezoidal map ---*/ - look_up_table.EnableFastTrapMap(); - /*--- string names of the controlling variables ---*/ string name_CV1 = "ProgressVariable";