From 54801809083d15a3e6126da93a12425c500e59ee Mon Sep 17 00:00:00 2001 From: tkiymaz Date: Wed, 10 Dec 2025 20:01:18 +0100 Subject: [PATCH 1/7] 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/7] 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/7] 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/7] 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/7] 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/7] 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/7] 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;