Merge pull request #1856 from su2code/cleanup_and_visualize_linelets

Cleanup Linelets and create output to visualize them

Common/include/geometry/CGeometry.hpp

@@ -28,6 +28,7 @@
 
 #pragma once
 
+#include <limits>
 #include "../parallelization/mpi_structure.hpp"
 
 #ifdef HAVE_METIS
@@ -72,7 +73,7 @@ using namespace std;
  * \author F. Palacios
  */
 class CGeometry {
-protected:
+ protected:
   enum : size_t {OMP_MIN_SIZE = 32}; /*!< \brief Chunk size for small loops. */
   enum : size_t {MAXNDIM = 3};
 
@@ -186,7 +187,29 @@ class CGeometry {
 
   ColMajorMatrix<uint8_t> CoarseGridColor_;  /*!< \brief Coarse grid levels, colorized. */
 
-public:
+ public:
+  /*!< \brief Linelets (mesh lines perpendicular to stretching direction). */
+  struct CLineletInfo {
+    static constexpr passivedouble ALPHA_ISOTROPIC = 0.8; /*!< \brief Detect isotropic mesh region. */
+    enum : unsigned long {MAX_LINELET_POINTS = 32}; /*!< \brief Maximum points per linelet. */
+
+    std::vector<std::vector<unsigned long>> linelets; /*!< \brief Point indices for each linelet. */
+
+    /*!< \brief Index of the linelet of each point ("linelets" transfered to points). */
+    std::vector<unsigned> lineletIdx;
+
+    /*!< \brief Signals that a point is not on a linelet. */
+    enum : unsigned {NO_LINELET = std::numeric_limits<unsigned>::max()};
+
+    /*!< \brief Coloring for OpenMP parallelization, "linelets" is sorted by color. */
+    std::vector<unsigned long> colorOffsets;
+
+    std::vector<uint8_t> lineletColor;  /*!< \brief Coloring transfered to points, for visualization. */
+  };
+ protected:
+  mutable CLineletInfo lineletInfo;
+
+ public:
   /*--- Main geometric elements of the grid. ---*/
 
   CPrimalGrid** elem{nullptr};           /*!< \brief Element vector (primal grid information). */
@@ -1658,6 +1681,11 @@ class CGeometry {
    */
   inline unsigned long GetElementColorGroupSize() const { return elemColorGroupSize; }
 
+  /*!
+   * \brief Get the linelet definition, this function computes the linelets if that has not been done yet.
+   */
+  const CLineletInfo& GetLineletInfo(const CConfig* config) const;
+
   /*!
    * \brief Compute an ADT including the coordinates of all viscous markers
    * \param[in] config - Definition of the particular problem.

Common/include/linear_algebra/CPreconditioner.hpp

@@ -270,7 +270,7 @@ class CLineletPreconditioner final : public CPreconditioner<ScalarType> {
    * \note Request the associated matrix to build the preconditioner.
    */
   inline void Build() override {
-    sparse_matrix.BuildJacobiPreconditioner();
+    sparse_matrix.BuildLineletPreconditioner(geometry, config);
   }
 };
 

Common/include/linear_algebra/CSysMatrix.hpp

@@ -148,10 +148,6 @@ class CSysMatrix {
 
   ScalarType *invM;                 /*!< \brief Inverse of (Jacobi) preconditioner, or diagonal of ILU. */
 
-  unsigned long nLinelet;                      /*!< \brief Number of Linelets in the system. */
-  vector<bool> LineletBool;                    /*!< \brief Identify if a point belong to a Linelet. */
-  vector<vector<unsigned long> > LineletPoint; /*!< \brief Linelet structure. */
-
   /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/
   mutable vector<vector<const ScalarType*> > LineletUpper; /*!< \brief Pointers to the upper blocks of the tri-diag system (working memory). */
   mutable vector<vector<ScalarType> > LineletInvDiag;      /*!< \brief Inverse of the diagonal blocks of the tri-diag system (working memory). */
@@ -886,9 +882,8 @@ class CSysMatrix {
    * \brief Build the Linelet preconditioner.
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
-   * \return Average number of points per linelet.
    */
-  unsigned long BuildLineletPreconditioner(CGeometry *geometry, const CConfig *config);
+  void BuildLineletPreconditioner(const CGeometry *geometry, const CConfig *config);
 
   /*!
    * \brief Multiply CSysVector by the preconditioner

Common/include/linear_algebra/CSysMatrix.inl

@@ -192,7 +192,8 @@ FORCEINLINE void CSysMatrix<ScalarType>::UpperProduct(const CSysVector<ScalarTyp
 
   for (auto index = dia_ptr[row_i]+1; index < row_ptr[row_i+1]; index++) {
     auto col_j = col_ind[index];
-    if (col_j < col_ub)
+    /*--- Always include halos. ---*/
+    if (col_j < col_ub || col_j >= nPointDomain)
       MatrixVectorProductAdd(&matrix[index*nVar*nEqn], &vec[col_j*nEqn], prod);
   }
 }

Common/include/toolboxes/graph_toolbox.hpp

@@ -530,7 +530,7 @@ T createNaturalColoring(Index_t numInnerIndexes)
  * \param[out] indexColor - Optional, vector with colors given to the outer indices.
  * \return Coloring in the same type of the input pattern.
  */
-template<class T, typename Color_t = char, size_t MaxColors = 64, size_t MaxMB = 128>
+template<typename Color_t = char, size_t MaxColors = 64, size_t MaxMB = 128, class T>
 T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors = false,
                      std::vector<Color_t>* indexColor = nullptr)
 {

Common/src/geometry/CGeometry.cpp

@@ -25,6 +25,8 @@
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
 
+#include <unordered_set>
+
 #include "../../include/geometry/CGeometry.hpp"
 #include "../../include/geometry/elements/CElement.hpp"
 #include "../../include/parallelization/omp_structure.hpp"
@@ -3843,7 +3845,7 @@ void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const*
     /*--- Color the coarse points. ---*/
     vector<tColor> color;
     const auto& adjacency = geometry[iMesh]->nodes->GetPoints();
-    if (colorSparsePattern<CCompressedSparsePatternUL, tColor, nColor>(adjacency, 1, false, &color).empty())
+    if (colorSparsePattern<tColor, nColor>(adjacency, 1, false, &color).empty())
       continue;
 
     /*--- Propagate colors to fine mesh. ---*/
@@ -3859,6 +3861,178 @@ void CGeometry::ColorMGLevels(unsigned short nMGLevels, const CGeometry* const*
   }
 }
 
+const CGeometry::CLineletInfo& CGeometry::GetLineletInfo(const CConfig* config) const {
+  auto& li = lineletInfo;
+  if (!li.linelets.empty() || nPoint == 0) return li;
+
+  li.lineletIdx.resize(nPoint, CLineletInfo::NO_LINELET);
+
+  /*--- Estimate number of linelets. ---*/
+
+  unsigned long nLinelet = 0;
+  for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+    if (config->GetSolid_Wall(iMarker) || config->GetMarker_All_KindBC(iMarker) == DISPLACEMENT_BOUNDARY) {
+      nLinelet += nVertex[iMarker];
+    }
+  }
+
+  /*--- If the domain contains well defined Linelets ---*/
+
+  unsigned long maxNPoints = 0, sumNPoints = 0;
+
+  if (nLinelet != 0) {
+    /*--- Define the basic linelets, starting from each vertex, preventing duplication of points. ---*/
+
+    li.linelets.resize(nLinelet);
+    nLinelet = 0;
+    for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
+      if (config->GetSolid_Wall(iMarker) || config->GetMarker_All_KindBC(iMarker) == DISPLACEMENT_BOUNDARY) {
+        for (auto iVertex = 0ul; iVertex < nVertex[iMarker]; iVertex++) {
+          const auto iPoint = vertex[iMarker][iVertex]->GetNode();
+          if (li.lineletIdx[iPoint] == CLineletInfo::NO_LINELET && nodes->GetDomain(iPoint)) {
+            li.linelets[nLinelet].push_back(iPoint);
+            li.lineletIdx[iPoint] = nLinelet;
+            ++nLinelet;
+          }
+        }
+      }
+    }
+    li.linelets.resize(nLinelet);
+
+    /*--- Create the linelet structure. ---*/
+
+    nLinelet = 0;
+    for (auto& linelet : li.linelets) {
+      while (linelet.size() < CLineletInfo::MAX_LINELET_POINTS) {
+        const auto iPoint = linelet.back();
+
+        /*--- Compute the value of the max and min weights to detect if this region is isotropic. ---*/
+
+        su2double max_weight = 0.0, min_weight = std::numeric_limits<su2double>::max();
+        for (auto iNode = 0u; iNode < nodes->GetnPoint(iPoint); iNode++) {
+          const auto jPoint = nodes->GetPoint(iPoint, iNode);
+          const auto iEdge = nodes->GetEdge(iPoint, iNode);
+          const auto* normal = edges->GetNormal(iEdge);
+          const su2double area = GeometryToolbox::Norm(nDim, normal);
+          const su2double volume_iPoint = nodes->GetVolume(iPoint);
+          const su2double volume_jPoint = nodes->GetVolume(jPoint);
+          const su2double weight = 0.5 * area * (1.0 / volume_iPoint + 1.0 / volume_jPoint);
+          max_weight = max(max_weight, weight);
+          min_weight = min(min_weight, weight);
+        }
+
+        /*--- Isotropic, stop this linelet. ---*/
+        if (min_weight / max_weight > CLineletInfo::ALPHA_ISOTROPIC) break;
+
+        /*--- Otherwise, add the closest valid neighbor. ---*/
+
+        su2double min_dist2 = std::numeric_limits<su2double>::max();
+        auto next_Point = iPoint;
+        const auto* iCoord = nodes->GetCoord(iPoint);
+
+        for (const auto jPoint : nodes->GetPoints(iPoint)) {
+          if (li.lineletIdx[jPoint] == CLineletInfo::NO_LINELET && nodes->GetDomain(jPoint)) {
+            const auto* jCoord = nodes->GetCoord(jPoint);
+            const su2double d2 = GeometryToolbox::SquaredDistance(nDim, iCoord, jCoord);
+            su2double cosTheta = 1;
+            if (linelet.size() > 1) {
+              const auto* kCoord = nodes->GetCoord(linelet[linelet.size() - 2]);
+              su2double dij[3] = {0.0}, dki[3] = {0.0};
+              GeometryToolbox::Distance(nDim, iCoord, kCoord, dki);
+              GeometryToolbox::Distance(nDim, jCoord, iCoord, dij);
+              cosTheta = GeometryToolbox::DotProduct(3, dki, dij) / sqrt(d2 * GeometryToolbox::SquaredNorm(nDim, dki));
+            }
+            if (d2 < min_dist2 && cosTheta > 0.7071) {
+              next_Point = jPoint;
+              min_dist2 = d2;
+            }
+          }
+        }
+
+        /*--- Did not find a suitable point. ---*/
+        if (next_Point == iPoint) break;
+
+        linelet.push_back(next_Point);
+        li.lineletIdx[next_Point] = nLinelet;
+      }
+      ++nLinelet;
+
+      maxNPoints = max(maxNPoints, linelet.size());
+      sumNPoints += linelet.size();
+    }
+  }
+
+  /*--- Average linelet size over all ranks. ---*/
+
+  unsigned long globalNPoints, globalNLineLets;
+  SU2_MPI::Allreduce(&sumNPoints, &globalNPoints, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+  SU2_MPI::Allreduce(&nLinelet, &globalNLineLets, 1, MPI_UNSIGNED_LONG, MPI_SUM, SU2_MPI::GetComm());
+
+  if (rank == MASTER_NODE) {
+    std::cout << "Computed linelet structure, "
+              << static_cast<unsigned long>(passivedouble(globalNPoints) / globalNLineLets)
+              << " points in each line (average)." << std::endl;
+  }
+
+  /*--- Color the linelets for OpenMP parallelization and visualization. ---*/
+
+  /*--- Adjacency between lines, computed from point neighbors. ---*/
+  std::vector<std::vector<unsigned long>> adjacency(nLinelet);
+  for (auto iLine = 0ul; iLine < nLinelet; ++iLine) {
+    std::unordered_set<unsigned long> neighbors;
+    for (const auto iPoint : li.linelets[iLine]) {
+      neighbors.insert(iPoint);
+      for (const auto jPoint : nodes->GetPoints(iPoint)) {
+        neighbors.insert(jPoint);
+      }
+    }
+    adjacency[iLine].reserve(neighbors.size());
+    for (const auto iPoint : neighbors) {
+      adjacency[iLine].push_back(iPoint);
+    }
+  }
+
+  const auto coloring = colorSparsePattern<uint8_t, std::numeric_limits<uint8_t>::max()>(
+      CCompressedSparsePatternUL(adjacency), 1, true);
+  const auto nColors = coloring.getOuterSize();
+
+  /*--- Sort linelets by color. ---*/
+  std::vector<std::vector<unsigned long>> sortedLinelets;
+  sortedLinelets.reserve(nLinelet);
+  li.colorOffsets.reserve(nColors + 1);
+  li.colorOffsets.push_back(0);
+
+  for (auto iColor = 0ul; iColor < nColors; ++iColor) {
+    for (const auto iLine : coloring.getInnerIter(iColor)) {
+      /*--- Store the new linelet index for its points. ---*/
+      for (const auto iPoint : li.linelets[iLine]) {
+        li.lineletIdx[iPoint] = sortedLinelets.size();
+      }
+      sortedLinelets.push_back(std::move(li.linelets[iLine]));
+    }
+    li.colorOffsets.push_back(sortedLinelets.size());
+  }
+  li.linelets = std::move(sortedLinelets);
+
+  /*--- For visualization, offset colors to avoid coloring across ranks. ---*/
+  std::vector<unsigned long> allNColors(size);
+  SU2_MPI::Allgather(&nColors, 1, MPI_UNSIGNED_LONG, allNColors.data(), 1, MPI_UNSIGNED_LONG, SU2_MPI::GetComm());
+  unsigned long offset = 0;
+  for (int i = 0; i < rank; ++i) offset += allNColors[i];
+
+  /*--- Finally, transfer colors to points, using "0" as "no linelet". ---*/
+  li.lineletColor.resize(nPoint, 0);
+  for (auto iColor = 0ul; iColor < nColors; ++iColor) {
+    for (auto iLine = li.colorOffsets[iColor]; iLine < li.colorOffsets[iColor + 1]; ++iLine) {
+      for (const auto iPoint : li.linelets[iLine]) {
+        li.lineletColor[iPoint] = 1 + offset + iColor;
+      }
+    }
+  }
+
+  return li;
+}
+
 void CGeometry::ComputeWallDistance(const CConfig* const* config_container, CGeometry ****geometry_container){
 
   int nZone = config_container[ZONE_0]->GetnZone();