diff --git a/Common/include/containers/CLookUpTable.hpp b/Common/include/containers/CLookUpTable.hpp index e7d8fd17e6e5..8148631d4b46 100644 --- a/Common/include/containers/CLookUpTable.hpp +++ b/Common/include/containers/CLookUpTable.hpp @@ -260,29 +260,6 @@ class CLookUpTable { void InterpolateToNearestNeighbors(const su2double val_CV1, const su2double val_CV2, const std::string& name_var, su2double* var_val, const unsigned long i_level = 0); - /*! - * \brief Determine if a point P(val_CV1,val_CV2) is inside the triangle val_id_triangle. - * \param[in] val_CV1 - First coordinate of point P(val_CV1,val_CV2) to check. - * \param[in] val_CV2 - Second coordinate of point P(val_CV1,val_CV2) to check. - * \param[in] val_id_triangle - ID of the triangle to check. - * \returns True if the point is in the triangle, false if it is outside. - */ - bool IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, unsigned long i_level = 0); - - /*! - * \brief Compute the area of a triangle given the 3 points of the triangle. - * \param[in] x1 - The coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). - * \param[in] y1 - The coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). - * \param[in] x2 - The coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). - * \param[in] y2 - The coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). - * \param[in] x3 - The coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). - * \param[in] y3 - The coordinates of the points P1(x1,y1), P2(x2,y2) and P3(x3,y3). - * \returns The absolute value of the area of the triangle. - */ - inline su2double TriArea(su2double x1, su2double y1, su2double x2, su2double y2, su2double x3, su2double y3) { - return abs((x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) * 0.5); - } - /*! * \brief Compute the values of the first and second controlling variable based on normalized query coordinates * \param[in] inclusion_levels - Pair containing lower(first) and upper(second) table inclusion level indices. diff --git a/Common/include/containers/CTrapezoidalMap.hpp b/Common/include/containers/CTrapezoidalMap.hpp index afd9eebd3ed5..091bca917a33 100644 --- a/Common/include/containers/CTrapezoidalMap.hpp +++ b/Common/include/containers/CTrapezoidalMap.hpp @@ -1,7 +1,8 @@ /*! * \file CTrapezoidalMap.hpp - * \brief Implementation of the trapezoidal map for tabulation and lookup of fluid properties - * \author D. Mayer, T. Economon + * \brief Memory-efficient trapezoidal map for 2D lookup table queries, + * based on the LUT implementation of P. Gomes (https://github.com/pcarruscag/LUT). + * \author T. Kiymaz * \version 8.3.0 "Harrier" * * SU2 Project Website: https://su2code.github.io @@ -27,82 +28,494 @@ #pragma once -#include +#include +#include +#include +#include +#include +#include #include -#include "../../Common/include/linear_algebra/blas_structure.hpp" -#include "../../Common/include/toolboxes/CSquareMatrixCM.hpp" +#include "../basic_types/datatype_structure.hpp" + +namespace su2_lut { + +using IntT = int32_t; +using RealT = su2double; + +/*--- Simple row-major matrix with a fixed number of columns. ---*/ +template +struct Matrix { + std::vector data; + + void resize(size_t rows, size_t) { data.resize(rows * N); } + size_t rows() const { return data.size() / N; } + + const T& operator()(size_t i, size_t j) const { return data[i * N + j]; } + T& operator()(size_t i, size_t j) { return data[i * N + j]; } +}; + +using Matrix2i = Matrix; +using Matrix3i = Matrix; +using VectorInt = std::vector; +using VectorReal = std::vector; + +/*--- The map is defined by the limits of the bands in the x direction and a CSR of + * the edge IDs in each band, sorted by the edge y position at the band midpoint. ---*/ +struct TrapezoidalMap { + VectorInt offsets, edge_id; + VectorReal x_bands, edge_y; +}; /*! - * \class CTrapezoidalMap - * \ingroup LookUpInterp - * \brief Construction of trapezoidal map for tabulated lookup - * \author: D. Mayer, T. Economon - * \version 8.3.0 "Harrier" + * \brief Orders points by ascending x coordinates and updates triangle indices. */ -class CTrapezoidalMap { - protected: - /* The unique values of x which exist in the data */ - std::vector unique_bands_x; +inline void ReorderPoints(Matrix3i& triangles, VectorReal& x, VectorReal& y) { + const IntT n_pts = static_cast(x.size()); + + std::vector perm(n_pts); + std::iota(perm.begin(), perm.end(), 0); + std::sort(perm.begin(), perm.end(), + [&x, &y](const auto i, const auto j) { return x[i] != x[j] ? x[i] < x[j] : y[i] < y[j]; }); + + auto reorder = [n_pts, &perm](const auto& v) { + VectorReal tmp(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + tmp[i] = v[perm[i]]; + } + return tmp; + }; + x = reorder(x); + y = reorder(y); + + std::vector inv_perm(n_pts); + for (IntT i = 0; i < n_pts; ++i) { + inv_perm[perm[i]] = i; + } + for (IntT i = 0; i < static_cast(triangles.rows()); ++i) { + for (IntT j = 0; j < 3; ++j) { + triangles(i, j) = inv_perm[triangles(i, j)]; + } + } +} + +/*! + * \brief Extracts unique edges from triangles. Edges are defined by two point IDs and + * up to two adjacent triangles (boundary edges have the second triangle ID < 0). + */ +inline void ExtractEdges(const Matrix3i& triangles, Matrix2i& edge_pts, Matrix2i& edge_faces) { + std::vector> edges; + edges.resize(3 * triangles.rows()); + + for (IntT i_tri = 0; i_tri < static_cast(triangles.rows()); ++i_tri) { + for (IntT i = 0; i < 3; ++i) { + const IntT j = (i + 1) % 3; + const IntT i_pt = std::min(triangles(i_tri, i), triangles(i_tri, j)); + const IntT j_pt = std::max(triangles(i_tri, i), triangles(i_tri, j)); + edges[3 * i_tri + i] = {i_pt, j_pt, i_tri}; + } + } + + /*--- Sort to identify duplicates. ---*/ + std::sort(edges.begin(), edges.end(), + [](const auto& a, const auto& b) { return a[0] != b[0] ? (a[0] < b[0]) : (a[1] < b[1]); }); + + auto is_equal = [](const auto& a, const auto& b) { return a[0] == b[0] && a[1] == b[1]; }; + + IntT n_edges = 1; + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + n_edges += static_cast(!is_equal(edges[i], edges[i - 1])); + } + + edge_pts.resize(n_edges, 2); + edge_faces.resize(n_edges, 2); + IntT pos = 0; + + auto new_edge = [&](const auto& edge) { + edge_pts(pos, 0) = edge[0]; + edge_pts(pos, 1) = edge[1]; + edge_faces(pos, 0) = edge[2]; + edge_faces(pos, 1) = -1; + ++pos; + }; + + new_edge(edges[0]); + for (IntT i = 1; i < static_cast(edges.size()); ++i) { + if (is_equal(edges[i], edges[i - 1])) { + edge_faces(pos - 1, 1) = edges[i][2]; + } else { + new_edge(edges[i]); + } + } +} + +/*! + * \brief Detects the x bands of the map. One band per unique x coordinate is used unless + * that would exceed max_bands, in which case equal-width bands are used to limit memory. + * \return Tuple of (n_bands, x_bands). + */ +inline auto DetectBands(const VectorReal& x, IntT max_bands = 0) { + if (max_bands <= 0) { + max_bands = std::min(IntT{5000}, static_cast(4.0 * std::sqrt(static_cast(x.size())))); + } + + IntT n_unique = 1; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != x[i - 1]) n_unique++; + } + + if (n_unique <= max_bands) { + const IntT n_bands = n_unique - 1; + VectorReal x_bands(n_unique); + IntT pos = 0; + x_bands[pos] = x[0]; + for (IntT i = 1; i < static_cast(x.size()); ++i) { + if (x[i] != x_bands[pos]) { + x_bands[++pos] = x[i]; + } + } + return std::make_tuple(n_bands, std::move(x_bands)); + } + + const RealT x_min = x.front(); + const RealT x_max = x.back(); + const RealT band_width = (x_max - x_min) / max_bands; + + VectorReal x_bands(max_bands + 1); + for (IntT i = 0; i <= max_bands; ++i) { + x_bands[i] = x_min + i * band_width; + } + x_bands[max_bands] = x_max; + + return std::make_tuple(max_bands, std::move(x_bands)); +} + +/*! + * \brief Builds the trapezoidal map for a set of edges (points must be ordered by x). + */ +inline void BuildTrapezoidalMap(const Matrix2i& edge_pts, const VectorReal& x, const VectorReal& y, + TrapezoidalMap& map) { + auto& x_bands = map.x_bands; + auto& offsets = map.offsets; + auto& edge_id = map.edge_id; + auto& edge_y = map.edge_y; + + auto clear_map = [&]() { + x_bands.clear(); + offsets.clear(); + edge_id.clear(); + edge_y.clear(); + }; + + const auto [n_bands, bands] = DetectBands(x); + x_bands = std::move(bands); + + if (n_bands <= 0) { + clear_map(); + return; + } + + auto find_band = [&x_bands, n_bands = n_bands](RealT x_val) -> IntT { + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x_val); + const IntT idx = static_cast(it - x_bands.begin()); + return std::min(std::max(IntT{0}, idx - 1), n_bands - 1); + }; + + /*--- Count edges per band. Each edge is stored in every band between the bands of its + * two endpoints (inclusive), a superset of the bands it overlaps. ---*/ + auto& counts = offsets; + counts.clear(); + counts.resize(n_bands + 1, 0); + + for (IntT i = 0; i < static_cast(edge_pts.rows()); ++i) { + const IntT band_0 = find_band(x[edge_pts(i, 0)]); + const IntT band_1 = find_band(x[edge_pts(i, 1)]); + + for (IntT j = std::min(band_0, band_1); j <= std::max(band_0, band_1); ++j) { + ++counts[j + 1]; + } + } + + /*--- Convert counts to offsets (CSR format). ---*/ + for (IntT i = 2; i < static_cast(offsets.size()); ++i) { + offsets[i] += offsets[i - 1]; + } + + /*--- Give up (build failure) rather than allocating an excessive amount of memory. ---*/ + const size_t memory_mb = static_cast(offsets.back()) * (sizeof(IntT) + sizeof(RealT)) / (1024 * 1024); + if (memory_mb > 2048) { + clear_map(); + return; + } + + edge_id.resize(offsets.back()); + edge_y.resize(offsets.back()); + auto pos = offsets; + + for (IntT i_edge = 0; i_edge < static_cast(edge_pts.rows()); ++i_edge) { + const IntT pt_0 = edge_pts(i_edge, 0); + const IntT pt_1 = edge_pts(i_edge, 1); + const RealT x_0 = x[pt_0], y_0 = y[pt_0]; + const RealT x_1 = x[pt_1], y_1 = y[pt_1]; + + const IntT band_0 = find_band(x_0); + const IntT band_1 = find_band(x_1); + + const RealT dx = x_1 - x_0; + const bool vertical = std::abs(SU2_TYPE::GetValue(dx)) < 1e-30; + const RealT dy_dx = vertical ? RealT{0} : (y_1 - y_0) / dx; - su2activematrix edge_limits_x; - su2activematrix edge_limits_y; + for (IntT j = std::min(band_0, band_1); j <= std::max(band_0, band_1); ++j) { + edge_id[pos[j]] = i_edge; + const RealT x_mid = (x_bands[j] + x_bands[j + 1]) / 2; + edge_y[pos[j]] = vertical ? (y_0 + y_1) / 2 : y_0 + dy_dx * (x_mid - x_0); + ++pos[j]; + } + } + + /*--- Sort the edges in each band by y coordinate. ---*/ + std::vector> tmp; + for (IntT i = 0; i < n_bands; ++i) { + const IntT begin = offsets[i]; + const IntT end = offsets[i + 1]; + if (begin >= end) continue; + + tmp.resize(end - begin); + for (auto k = begin; k < end; ++k) { + tmp[k - begin] = {edge_id[k], edge_y[k]}; + } + std::sort(tmp.begin(), tmp.end(), [](const auto& a, const auto& b) { return a.second < b.second; }); + for (auto k = begin; k < end; ++k) { + edge_id[k] = tmp[k - begin].first; + edge_y[k] = tmp[k - begin].second; + } + } +} + +/*! + * \brief Returns the IDs of the edges directly below and above a query point + * (either ID can be -1 if the point is at a boundary). + */ +inline auto QueryTrapezoidalMap(const TrapezoidalMap& map, const Matrix2i& edge_pts, const VectorReal& x_coords, + const VectorReal& y_coords, const RealT& x, const RealT& y) { + if (map.x_bands.size() < 2 || map.offsets.empty()) { + return std::make_pair(IntT{-1}, IntT{-1}); + } + + const auto& x_bands = map.x_bands; + const IntT n_bands = static_cast(x_bands.size()) - 1; + auto it = std::lower_bound(x_bands.begin(), x_bands.end(), x); + const IntT d = static_cast(it - x_bands.begin()); + const IntT band_idx = std::min(std::max(IntT{0}, d - 1), n_bands - 1); + + RealT best_y_below = -1e300; + RealT best_y_above = 1e300; + IntT edge_below = -1; + IntT edge_above = -1; + + const IntT begin = map.offsets[band_idx]; + const IntT end = map.offsets[band_idx + 1]; + + for (IntT k = begin; k < end; ++k) { + const IntT e_id = map.edge_id[k]; + + const IntT p0 = edge_pts(e_id, 0); + const IntT p1 = edge_pts(e_id, 1); + const RealT x0 = x_coords[p0], y0 = y_coords[p0]; + const RealT x1 = x_coords[p1], y1 = y_coords[p1]; + + if (x < std::min(x0, x1) - 1e-10 || x > std::max(x0, x1) + 1e-10) { + continue; + } + + /*--- y position of the edge at the query x. ---*/ + RealT edge_y_at_x; + const RealT dx = x1 - x0; + if (std::abs(SU2_TYPE::GetValue(dx)) < 1e-30) { + edge_y_at_x = (y0 + y1) / 2.0; + } else { + edge_y_at_x = y0 + (x - x0) / dx * (y1 - y0); + } + + if (edge_y_at_x <= y + 1e-10 && edge_y_at_x > best_y_below) { + best_y_below = edge_y_at_x; + edge_below = e_id; + } + if (edge_y_at_x >= y - 1e-10 && edge_y_at_x < best_y_above) { + best_y_above = edge_y_at_x; + edge_above = e_id; + } + } + + return std::make_pair(edge_below, edge_above); +} + +/*! + * \brief Returns the IDs of the triangles adjacent to two query edges (up to 3 triangles). + */ +inline auto AdjacentTriangles(const IntT edge_0, const IntT edge_1, const Matrix2i& edge_faces) { + std::array tris = {-1, -1, -1}; + IntT pos = 0; + + auto insert = [&tris, &pos](const IntT t) { + if (t < 0) return; + for (IntT i = 0; i < pos; ++i) { + if (t == tris[i]) return; + } + tris[pos++] = t; + }; + + auto get_tris = [&edge_faces](const IntT e) { + if (e < 0) return std::array{IntT{-1}, IntT{-1}}; + return std::array{edge_faces(e, 0), edge_faces(e, 1)}; + }; + + for (const auto e : {edge_0, edge_1}) { + for (const auto t : get_tris(e)) { + insert(t); + } + } + return tris; +} + +/*! + * \brief Computes the barycentric coordinates of point (x_q, y_q) in a triangle. + */ +inline auto TriangleCoords(const IntT i_tri, const Matrix3i& triangles, const VectorReal& x, const VectorReal& y, + const RealT x_q, const RealT y_q) { + const IntT p0 = triangles(i_tri, 0); + const IntT p1 = triangles(i_tri, 1); + const IntT p2 = triangles(i_tri, 2); - su2vector > edge_to_triangle; + const RealT x0 = x[p0], y0 = y[p0]; + const RealT x1 = x[p1], y1 = y[p1]; + const RealT x2 = x[p2], y2 = y[p2]; + + const RealT dx1 = x1 - x0, dy1 = y1 - y0; + const RealT dx2 = x2 - x0, dy2 = y2 - y0; + + auto cross = [](const RealT ux, const RealT uy, const RealT vx, const RealT vy) { return ux * vy - uy * vx; }; + + const RealT det = cross(dx1, dy1, dx2, dy2); + if (std::abs(SU2_TYPE::GetValue(det)) < 1e-30) { + return std::array{RealT{0}, RealT{0}, RealT{0}}; + } - /* The value that each edge which intersects the band takes within that - * same band. Used to sort the edges */ - su2vector > > y_edge_at_band_mid; + const RealT inv_det = 1.0 / det; + const RealT a = (cross(x_q, y_q, dx2, dy2) - cross(x0, y0, dx2, dy2)) * inv_det; + const RealT b = (cross(x0, y0, dx1, dy1) - cross(x_q, y_q, dx1, dy1)) * inv_det; - double memory_footprint = 0; + return std::array{1 - a - b, a, b}; +} + +/*! + * \brief Checks if a point is inside a triangle based on its barycentric coordinates. + */ +inline bool InTriangle(const std::array& coords, const RealT tol = 0.0) { + return coords[0] >= -tol && coords[1] >= -tol && coords[2] >= -tol; +} + +/*! + * \brief Finds the triangle containing a point using the trapezoidal map. + */ +inline IntT FindTriangle(const TrapezoidalMap& map, const Matrix3i& triangles, const Matrix2i& edge_pts, + const Matrix2i& edge_faces, const VectorReal& x, const VectorReal& y, const RealT x_q, + const RealT y_q, std::array& bary_out) { + const auto [e_below, e_above] = QueryTrapezoidalMap(map, edge_pts, x, y, x_q, y_q); + const auto candidates = AdjacentTriangles(e_below, e_above, edge_faces); + + const RealT tol = 1e-12; + for (const auto t : candidates) { + if (t < 0) continue; + + const auto coords = TriangleCoords(t, triangles, x, y, x_q, y_q); + if (InTriangle(coords, tol)) { + bary_out = coords; + return t; + } + } + + bary_out = {0.0, 0.0, 0.0}; + return -1; +} + +} // namespace su2_lut + +/*! + * \class CTrapezoidalMap + * \ingroup LookUpInterp + * \brief Trapezoidal map for finding the triangle containing a query point in a 2D triangulation. + */ +class CTrapezoidalMap { + private: + su2_lut::Matrix3i triangles; + su2_lut::Matrix2i edge_pts, edge_faces; + su2_lut::VectorReal x_coords, y_coords; + su2_lut::TrapezoidalMap map; + + unsigned long n_points = 0; + unsigned long n_triangles = 0; public: CTrapezoidalMap() = default; - CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size, - const std::vector >& edges, - const su2vector >& edge_to_triangle, bool display = false); - /*! - * \brief return the index to the triangle that contains the coordinates (val_x,val_y) - * \param[in] val_x - x-coordinate or first independent variable - * \param[in] val_y - y-coordinate or second independent variable - * \param[out] val_index - index to the triangle + * \brief Build the trapezoidal map from a triangulation. + * \return True on success. */ - unsigned long GetTriangle(const su2double val_x, const su2double val_y); + bool Build(unsigned long num_points, unsigned long num_triangles, const su2double* x, const su2double* y, + const unsigned long* connectivity) { + n_points = num_points; + n_triangles = num_triangles; - /*! - * \brief get the indices of the vertical coordinate band (xmin,xmax) in the 2D search space - * that contains the coordinate val_x - * \param[in] val_x - x-coordinate or first independent variable - * \param[out] val_band - a pair(i_low,i_up) , the lower index and upper index between which the value val_x - * can be found - */ - std::pair GetBand(const su2double val_x); + if (num_points == 0 || num_triangles == 0) return false; - /*! - * \brief for a given coordinate (val_x,value), known to be in the band (xmin,xmax) with band index (i_low,i_up), - * find the edges in the band (these edges come from the triangulation) that enclose the coordinate - * \param[in] val_band - pair i_low,i_up - * \param[in] val_x - x-coordinate or first independent variable - * \param[in] val_y - y-coordinate or first independent variable - * \param[out] pair (edge_low,edge_up) - lower edge and upper edge of a triangle that encloses the coordinate - */ - std::pair GetEdges(std::pair val_band, su2double val_x, - su2double val_y) const; + x_coords.assign(x, x + num_points); + y_coords.assign(y, y + num_points); + + triangles.resize(num_triangles, 3); + for (size_t i = 0; i < 3 * num_triangles; ++i) { + triangles.data[i] = static_cast(connectivity[i]); + } + + su2_lut::ReorderPoints(triangles, x_coords, y_coords); + su2_lut::ExtractEdges(triangles, edge_pts, edge_faces); + su2_lut::BuildTrapezoidalMap(edge_pts, x_coords, y_coords, map); + + return !map.x_bands.empty() && !map.offsets.empty(); + } /*! - * \brief determine if the x-coordinate falls within the bounds xmin,xmax of the table - * \param[in] val_x - x-coordinate or first independent variable - * \param[out] bool - true if val_x is within (xmin,xmax) + * \brief Find the triangle containing a query point. + * \return True if the point is inside the triangulation. */ - inline bool IsInsideHullX(su2double val_x) { - return (val_x >= unique_bands_x.front()) && (val_x <= unique_bands_x.back()); + bool FindTriangle(su2double val_x, su2double val_y, unsigned long& triangle_id, + std::array& bary_coords) const { + if (n_triangles == 0 || n_points == 0 || map.x_bands.empty()) { + bary_coords = {0.0, 0.0, 0.0}; + return false; + } + + std::array bary; + const su2_lut::IntT tri_id = + su2_lut::FindTriangle(map, triangles, edge_pts, edge_faces, x_coords, y_coords, val_x, val_y, bary); + + if (tri_id < 0) return false; + + triangle_id = static_cast(tri_id); + bary_coords = {bary[0], bary[1], bary[2]}; + return true; } /*! - * \brief get memory footprint of trapezoidal map. - * \return - memory footprint in mega bytes. + * \brief Get the memory footprint of the map in MB. */ - double GetMemoryFootprint() const { return memory_footprint; } + double GetMemoryFootprint() const { + const size_t bytes = + (map.edge_id.size() + map.offsets.size() + edge_pts.data.size() + edge_faces.data.size() + + triangles.data.size()) * + sizeof(su2_lut::IntT) + + (map.edge_y.size() + map.x_bands.size() + x_coords.size() + y_coords.size()) * sizeof(su2_lut::RealT); + return double(bytes) / (1024.0 * 1024.0); + } }; diff --git a/Common/src/containers/CLookUpTable.cpp b/Common/src/containers/CLookUpTable.cpp index 98722a6b5b28..2b4a6bb41dde 100644 --- a/Common/src/containers/CLookUpTable.cpp +++ b/Common/src/containers/CLookUpTable.cpp @@ -46,10 +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(); @@ -62,16 +59,10 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, 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; @@ -79,38 +70,36 @@ CLookUpTable::CLookUpTable(const string& var_file_name_lut, string name_CV1_in, trap_map_x_y.resize(n_table_levels); su2double startTime = SU2_MPI::Wtime(); - unsigned short barwidth = 65; - bool display_map_info = (n_table_levels < 2); double tmap_memory_footprint = 0; + for (auto i_level = 0ul; i_level < n_table_levels; i_level++) { - trap_map_x_y[i_level] = - CTrapezoidalMap(GetDataP(name_CV1, i_level), GetDataP(name_CV2, i_level), table_data[i_level].cols(), - edges[i_level], edge_to_triangle[i_level], display_map_info); - tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint(); - /* Display a progress bar to monitor table generation process */ - if (rank == MASTER_NODE) { - su2double progress = su2double(i_level) / n_table_levels; - auto completed = floor(progress * barwidth); - auto to_do = barwidth - completed; - cout << "[" << setfill('=') << setw(completed); - cout << '>'; - cout << setfill(' ') << setw(to_do) << std::right << "] " << 100 * progress << "%\r"; - cout.flush(); + const auto n_pts = n_points[i_level]; + const auto n_tris = n_triangles[i_level]; + + std::vector x_coords(n_pts); + std::vector y_coords(n_pts); + for (auto i_point = 0ul; i_point < n_pts; ++i_point) { + x_coords[i_point] = table_data[i_level][idx_CV1][i_point]; + y_coords[i_point] = table_data[i_level][idx_CV2][i_point]; + } + + std::vector tri_conn(3 * n_tris); + for (auto i_tri = 0ul; i_tri < n_tris; ++i_tri) { + tri_conn[3 * i_tri + 0] = triangles[i_level][i_tri][0]; + tri_conn[3 * i_tri + 1] = triangles[i_level][i_tri][1]; + tri_conn[3 * i_tri + 2] = triangles[i_level][i_tri][2]; } + + if (!trap_map_x_y[i_level].Build(n_pts, n_tris, x_coords.data(), y_coords.data(), tri_conn.data())) + SU2_MPI::Error( + "Construction of trapezoidal map failed for level " + std::to_string(i_level) + " of table " + file_name_lut, + CURRENT_FUNCTION); + tmap_memory_footprint += trap_map_x_y[i_level].GetMemoryFootprint(); } su2double stopTime = SU2_MPI::Wtime(); if (rank == MASTER_NODE) { - switch (table_dim) { - case 2: - cout << "\nConstruction of trapezoidal map took " << stopTime - startTime << " seconds\n" << endl; - break; - case 3: - cout << "\nConstruction of trapezoidal map stack took " << stopTime - startTime << " seconds\n" << endl; - break; - default: - break; - } + cout << "Construction of trapezoidal map took " << stopTime - startTime << " seconds\n"; cout << "Trapezoidal map memory footprint: " << tmap_memory_footprint << " MB\n"; cout << "Table data memory footprint: " << memory_footprint_data << " MB\n" << endl; } @@ -617,18 +606,8 @@ 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); - } - return false; + std::array bary_coords; + return trap_map_x_y[iLevel].FindTriangle(val_CV1, val_CV2, id_triangle, bary_coords); } void CLookUpTable::GetInterpCoeffs(su2double val_CV1, su2double val_CV2, const su2activematrix& interp_mat_inv, @@ -783,26 +762,6 @@ void CLookUpTable::InterpolateToNearestNeighbors(const su2double val_CV1, const InterpolateToNearestNeighbors(val_CV1, val_CV2, names_var, val_names_var, i_level); } -bool CLookUpTable::IsInTriangle(su2double val_CV1, su2double val_CV2, unsigned long val_id_triangle, - unsigned long i_level) { - su2double tri_x_0 = table_data[i_level][idx_CV1][triangles[i_level][val_id_triangle][0]]; - su2double tri_y_0 = table_data[i_level][idx_CV2][triangles[i_level][val_id_triangle][0]]; - - su2double tri_x_1 = table_data[i_level][idx_CV1][triangles[i_level][val_id_triangle][1]]; - su2double tri_y_1 = table_data[i_level][idx_CV2][triangles[i_level][val_id_triangle][1]]; - - su2double tri_x_2 = table_data[i_level][idx_CV1][triangles[i_level][val_id_triangle][2]]; - su2double tri_y_2 = table_data[i_level][idx_CV2][triangles[i_level][val_id_triangle][2]]; - - su2double area_tri = TriArea(tri_x_0, tri_y_0, tri_x_1, tri_y_1, tri_x_2, tri_y_2); - - su2double area_0 = TriArea(val_CV1, val_CV2, tri_x_1, tri_y_1, tri_x_2, tri_y_2); - su2double area_1 = TriArea(tri_x_0, tri_y_0, val_CV1, val_CV2, tri_x_2, tri_y_2); - su2double area_2 = TriArea(tri_x_0, tri_y_0, tri_x_1, tri_y_1, val_CV1, val_CV2); - - return (abs(area_tri - (area_0 + area_1 + area_2)) < area_tri * 1e-10); -} - bool CLookUpTable::CheckForVariables(const std::vector& vars_to_check) const { for (const string& var_to_check : vars_to_check) { if (!std::any_of(names_var.begin(), names_var.end(), diff --git a/Common/src/containers/CTrapezoidalMap.cpp b/Common/src/containers/CTrapezoidalMap.cpp deleted file mode 100644 index ea0472848455..000000000000 --- a/Common/src/containers/CTrapezoidalMap.cpp +++ /dev/null @@ -1,283 +0,0 @@ -/*! - * \file CTrapezoidalMap.cpp - * \brief Implementation of the trapezoidal map for tabulation and lookup of fluid properties - * \author D. Mayer, T. Economon, N. Beishuizen - * \version 8.3.0 "Harrier" - * - * SU2 Project Website: https://su2code.github.io - * - * The SU2 Project is maintained by the SU2 Foundation - * (http://su2foundation.org) - * - * Copyright 2012-2025, SU2 Contributors (cf. AUTHORS.md) - * - * SU2 is free software; you can redistribute it and/or - * modify it under the terms of the GNU Lesser General Public - * License as published by the Free Software Foundation; either - * version 2.1 of the License, or (at your option) any later version. - * - * SU2 is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * Lesser General Public License for more details. - * - * You should have received a copy of the GNU Lesser General Public - * License along with SU2. If not, see . - */ - -#include -#include - -#include "../../Common/include/option_structure.hpp" -#include "../../Common/include/containers/CTrapezoidalMap.hpp" - -using namespace std; - -/* Trapezoidal map implementation. Reference: - * M. de Berg, O. Cheong M. van Kreveld, M. Overmars, - * Computational Geometry, Algorithms and Applications pp. 121-146 (2008) - * NOTE: the current implementation is actually the simpler 'slab' approach. - */ -CTrapezoidalMap::CTrapezoidalMap(const su2double* samples_x, const su2double* samples_y, const unsigned long size, - vector > const& edges, - su2vector > const& val_edge_to_triangle, bool display) { - int rank = SU2_MPI::GetRank(); - su2double startTime = SU2_MPI::Wtime(); - - edge_to_triangle = su2vector >(val_edge_to_triangle); - - unique_bands_x.assign(samples_x, samples_x + size); - - /* sort x_bands and make them unique */ - sort(unique_bands_x.begin(), unique_bands_x.end()); - - auto iter = unique(unique_bands_x.begin(), unique_bands_x.end()); - - unique_bands_x.resize(distance(unique_bands_x.begin(), iter)); - - edge_limits_x.resize(edges.size(), 2); - edge_limits_y.resize(edges.size(), 2); - - /* store x and y values of each edge in a vector for a slight speed up - * as it prevents some uncoalesced accesses */ - for (unsigned long j = 0; j < edges.size(); j++) { - edge_limits_x[j][0] = samples_x[edges[j][0]]; - edge_limits_x[j][1] = samples_x[edges[j][1]]; - edge_limits_y[j][0] = samples_y[edges[j][0]]; - edge_limits_y[j][1] = samples_y[edges[j][1]]; - } - - /* number of bands */ - unsigned long n_bands_x = unique_bands_x.size() - 1; - /* band index */ - unsigned long i_band = 0; - /* number of edges */ - unsigned long n_edges = edges.size(); - /* edge index */ - unsigned long i_edge = 0; - unsigned long j_edge = 0; - /* counter for edges intersects */ - unsigned long n_intersects = 0; - /* lower and upper x value of each band */ - su2double band_lower_x = 0; - su2double band_upper_x = 0; - - su2double x_0; - su2double y_0; - su2double dy_edge; - su2double dx_edge; - su2double x_band_mid; - - /* y values of all intersecting edges for every band */ - y_edge_at_band_mid.resize(unique_bands_x.size() - 1); - - /* loop over bands */ - while (i_band < n_bands_x) { - band_lower_x = unique_bands_x[i_band]; - band_upper_x = unique_bands_x[i_band + 1]; - i_edge = 0; - n_intersects = 0; - - /* loop over edges and determine which edges appear in current band */ - while (i_edge < n_edges) { - /* check if edge intersects the band - * (vertical edges are automatically discarded) */ - if (((edge_limits_x[i_edge][0] <= band_lower_x) and (edge_limits_x[i_edge][1] >= band_upper_x)) or - ((edge_limits_x[i_edge][1] <= band_lower_x) and (edge_limits_x[i_edge][0] >= band_upper_x))) { - y_edge_at_band_mid[i_band].emplace_back(0.0, 0); - - x_0 = edge_limits_x[i_edge][0]; - y_0 = edge_limits_y[i_edge][0]; - - dy_edge = edge_limits_y[i_edge][1] - edge_limits_y[i_edge][0]; - dx_edge = edge_limits_x[i_edge][1] - edge_limits_x[i_edge][0]; - x_band_mid = (band_lower_x + band_upper_x) / 2.0; - - y_edge_at_band_mid[i_band][n_intersects].first = y_0 + dy_edge / dx_edge * (x_band_mid - x_0); - - /* save edge index so it can later be recalled when searching */ - y_edge_at_band_mid[i_band][n_intersects].second = i_edge; - - n_intersects++; - } - i_edge++; - } - - /* sort edges by their y values. - * note that these y values are unique (i.e. edges cannot - * intersect in a band) */ - sort(y_edge_at_band_mid[i_band].begin(), y_edge_at_band_mid[i_band].end()); - - i_band++; - } - - su2double stopTime = SU2_MPI::Wtime(); - - /* calculate size of trapezoidal map components */ - double size_unique_bands = sizeof(su2double) * unique_bands_x.size() / 1e6; - double size_edge_limits_x = sizeof(su2double) * edge_limits_x.size() * 2 / 1e6; - double size_edge_limits_y = sizeof(su2double) * edge_limits_y.size() * 2 / 1e6; - - double size_edge_to_triangle = 0; - for (i_edge = 0; i_edge < edge_to_triangle.size(); i_edge++) - for (j_edge = 0; j_edge < edge_to_triangle[i_edge].size(); j_edge++) - size_edge_to_triangle += sizeof(unsigned long) / 1e6; - - double size_y_edge_at_band_mid = 0; - for (unsigned long i_y = 0; i_y < y_edge_at_band_mid.size(); i_y++) - for (unsigned long j_y = 0; j_y < y_edge_at_band_mid[i_y].size(); j_y++) - size_y_edge_at_band_mid += sizeof(su2double) / 1e6 + sizeof(unsigned long) / 1e6; - - memory_footprint = - size_unique_bands + size_edge_limits_x + size_edge_limits_y + size_edge_to_triangle + size_y_edge_at_band_mid; - - /* print size of trapezoidal map components to screen */ - if ((rank == MASTER_NODE) && display) { - cout << setfill(' '); - cout << "\n" << endl; - cout << "+------------------------------------------------------------------+\n"; - cout << "| Trapezoidal map info |\n"; - cout << "+------------------------------------------------------------------+" << endl; - - cout << "| Time to construct trapezoidal map: " << setw(22) << right << stopTime - startTime << " sec" - << " |" << endl; - cout << "| Size of unique_bands in memory: " << setw(22) << size_unique_bands << " MB " - << " |" << endl; - cout << "| Size of edge_limits_x in memory: " << setw(22) << size_edge_limits_x << " MB " - << " |" << endl; - cout << "| Size of edge_limits_y in memory: " << setw(22) << size_edge_limits_y << " MB " - << " |" << endl; - cout << "| Size of edge_to_triangle in memory: " << setw(22) << size_edge_to_triangle << " MB " - << " |" << endl; - cout << "| Size of y_edge_at_band_mid in memory: " << setw(22) << size_y_edge_at_band_mid << " MB " - << " |" << endl; - cout << "| Total: " << setw(22) << memory_footprint << " MB " - << " |" << endl; - cout << "+------------------------------------------------------------------+" << endl; - cout << "\n" << endl; - } -} - -unsigned long CTrapezoidalMap::GetTriangle(const su2double val_x, const su2double val_y) { - /* find x band in which val_x sits */ - pair band = GetBand(val_x); - - /* within that band, find edges which enclose the (val_x, val_y) point */ - pair edges = GetEdges(band, val_x, val_y); - - /* identify the adjacent triangles using the two edges */ - std::array triangles_edge_low; - for (unsigned long i = 0; i < edge_to_triangle[edges.first].size(); i++) - triangles_edge_low[i] = edge_to_triangle[edges.first][i]; - - std::array triangles_edge_up; - for (unsigned long i = 0; i < edge_to_triangle[edges.second].size(); i++) - triangles_edge_up[i] = edge_to_triangle[edges.second][i]; - - sort(triangles_edge_low.begin(), triangles_edge_low.end()); - sort(triangles_edge_up.begin(), triangles_edge_up.end()); - - /* The intersection of the faces to which upper or lower belongs is the face that both belong to. */ - vector triangle; - set_intersection(triangles_edge_up.begin(), triangles_edge_up.end(), triangles_edge_low.begin(), - triangles_edge_low.end(), std::back_inserter(triangle)); - - /*--- We failed to find an intersection, so take the lower triangle inside the band enclosing the point---*/ - if (triangle.size() < 1) { - triangle.resize(1, triangles_edge_low[0]); - } - - return triangle[0]; -} - -pair CTrapezoidalMap::GetBand(const su2double val_x) { - unsigned long i_low = 0; - unsigned long i_up = 0; - su2double val_x_sample = val_x; - /* check if val_x is in x-bounds of the table, if not then project val_x to either x-min or x-max */ - if (val_x_sample < unique_bands_x.front()) val_x_sample = unique_bands_x.front(); - if (val_x_sample > unique_bands_x.back()) val_x_sample = unique_bands_x.back(); - - std::pair::iterator, std::vector::iterator> bounds; - bounds = std::equal_range(unique_bands_x.begin(), unique_bands_x.end(), val_x_sample); - - /*--- if upper bound = 0, then use the range [0,1] ---*/ - i_up = max(1, bounds.first - unique_bands_x.begin()); - i_low = i_up - 1; - - return make_pair(i_low, i_up); -} - -pair CTrapezoidalMap::GetEdges(pair val_band, - su2double val_x, su2double val_y) const { - su2double next_y; - su2double y_edge_low; - su2double y_edge_up; - su2double x_edge_low; - su2double x_edge_up; - - unsigned long i_band_low = val_band.first; - - unsigned long next_edge; - - unsigned long j_low = 0; - unsigned long j_mid = 0; - unsigned long j_up = 0; - - j_up = y_edge_at_band_mid[i_band_low].size() - 1; - j_low = 0; - - while (j_up - j_low > 1) { - j_mid = (j_up + j_low) / 2; - - // Select the edge associated with the x band (i_band_low) - // Search for the RunEdge in the y direction (second value is index of - // edge) - next_edge = y_edge_at_band_mid[i_band_low][j_mid].second; - - y_edge_low = edge_limits_y[next_edge][0]; - y_edge_up = edge_limits_y[next_edge][1]; - x_edge_low = edge_limits_x[next_edge][0]; - x_edge_up = edge_limits_x[next_edge][1]; - - // The search variable in j should be interpolated in i as well - next_y = y_edge_low + (y_edge_up - y_edge_low) / (x_edge_up - x_edge_low) * (val_x - x_edge_low); - - if (next_y > val_y) { - j_up = j_mid; - - } else if (next_y < val_y) { - j_low = j_mid; - - } else if (next_y == val_y) { - j_low = j_mid; - j_up = j_low + 1; - break; - } - } - - unsigned long edge_low = y_edge_at_band_mid[i_band_low][j_low].second; - unsigned long edge_up = y_edge_at_band_mid[i_band_low][j_up].second; - - return make_pair(edge_low, edge_up); -} diff --git a/Common/src/containers/meson.build b/Common/src/containers/meson.build index 4c8d4fe618a4..0743fd6b60c0 100644 --- a/Common/src/containers/meson.build +++ b/Common/src/containers/meson.build @@ -1,3 +1,2 @@ -common_src += files(['CTrapezoidalMap.cpp', - 'CFileReaderLUT.cpp', +common_src += files(['CFileReaderLUT.cpp', 'CLookUpTable.cpp']) diff --git a/SU2_CFD/src/fluid/CFluidFlamelet.cpp b/SU2_CFD/src/fluid/CFluidFlamelet.cpp index f6fd529df832..76322a641c2b 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: @@ -79,6 +81,7 @@ 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]); break; + default: if (rank == MASTER_NODE) { cout << "***********************************************" << endl; @@ -86,7 +89,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); @@ -104,8 +108,7 @@ CFluidFlamelet::CFluidFlamelet(CConfig* config, su2double value_pressure_operati } CFluidFlamelet::~CFluidFlamelet() { - if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT) - delete look_up_table; + if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::LUT) delete look_up_table; #ifdef USE_MLPCPP if (Kind_DataDriven_Method == ENUM_DATADRIVEN_METHOD::MLP) { delete iomap_TD; @@ -113,7 +116,7 @@ CFluidFlamelet::~CFluidFlamelet() { delete iomap_LookUp; delete lookup_mlp; if (preferential_diffusion) delete iomap_PD; - } + } #endif } @@ -184,8 +187,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: @@ -206,7 +208,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 ---*/ @@ -259,10 +262,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(); @@ -271,16 +274,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])); } } @@ -291,7 +294,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; @@ -326,7 +329,6 @@ unsigned long CFluidFlamelet::EvaluateDataSet(const vector& input_sca default: break; } - /*--- Add all quantities and their names to the look up vectors. ---*/ bool inside; @@ -339,8 +341,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());