Option to use level scheduling for OMP parallelization of ILU (#2793)

* use levels for FE connectivity

* config option and template etc.

* fix

* try to fix ARM 64 warning

* Apply suggestion from @pcarruscag

* regressions

Author: pcarruscag

Common/include/CConfig.hpp

@@ -645,6 +645,7 @@ class CConfig {
   unsigned long Linear_Solver_Restart_Deflation; /*!< \brief Number of vectors used for deflated restarts. */
   unsigned long Linear_Solver_Prec_Threads;      /*!< \brief Number of threads per rank for ILU and LU_SGS preconditioners. */
   unsigned short Linear_Solver_ILU_n;            /*!< \brief ILU fill=in level. */
+  bool Linear_Solver_ILU_levels;                 /*!< \brief Use level scheduling for OMP parallelization of ILU. */
   su2double SemiSpan;                   /*!< \brief Wing Semi span. */
   su2double MSW_Alpha;                  /*!< \brief Coefficient for blending states in the MSW scheme. */
   su2double Roe_Kappa;                  /*!< \brief Relaxation of the Roe scheme. */
@@ -4377,6 +4378,11 @@ class CConfig {
    */
   unsigned short GetLinear_Solver_ILU_n(void) const { return Linear_Solver_ILU_n; }
 
+  /*!
+   * \brief Get whether to use level scheduling for OMP parallelization of ILU.
+   */
+  bool GetLinear_Solver_ILU_levels(void) const { return Linear_Solver_ILU_levels; }
+
   /*!
    * \brief Get restart frequency of the linear solver for the implicit formulation.
    * \return Restart frequency of the linear solver for the implicit formulation.

Common/include/linear_algebra/CSysMatrix.hpp

@@ -28,9 +28,10 @@
 
 #pragma once
 
-#include "../../include/CConfig.hpp"
+#include "../CConfig.hpp"
 #include "CSysVector.hpp"
 #include "CPastixWrapper.hpp"
+#include "../toolboxes/graph_toolbox.hpp"
 
 #include <cstdlib>
 #include <vector>
@@ -120,10 +121,9 @@ class CSysMatrix {
   const int rank; /*!< \brief MPI Rank. */
   const int size; /*!< \brief MPI Size. */
 
-  enum : size_t {
-    MAXNVAR = 20
-  }; /*!< \brief Maximum number of variables the matrix can handle. The static
-                 size is needed for fast, per-thread, static memory allocation. */
+  /*!< \brief Maximum number of variables the matrix can handle. The static
+   * size is needed for fast, per-thread, static memory allocation. */
+  enum : size_t { MAXNVAR = 20 };
 
   enum { OMP_MAX_SIZE_L = 8192 }; /*!< \brief Max. chunk size used in light parallel for loops. */
   enum { OMP_MAX_SIZE_H = 512 };  /*!< \brief Max. chunk size used in heavy parallel for loops. */
@@ -158,7 +158,10 @@ class CSysMatrix {
   const unsigned long* col_ind_ilu; /*!< \brief Column index for each of the elements in val() (ILU). */
   unsigned short ilu_fill_in;       /*!< \brief Fill in level for the ILU preconditioner. */
 
-  ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner, or diagonal of ILU. */
+  /*!< \brief Level structure for alternative shared memory parallelization of ILU. */
+  CCompressedSparsePatternUL levels_ilu;
+
+  ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */
 
   /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/
   mutable vector<vector<const ScalarType*> >

Common/include/toolboxes/graph_toolbox.hpp

@@ -516,7 +516,7 @@ T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool includeOuterId
 
   {
     /*--- For each color keep track of the inner indices that are in it. ---*/
-    std::vector<std::vector<bool> > innerInColor;
+    std::vector<std::vector<bool>> innerInColor;
     innerInColor.emplace_back(nInner, false);
 
     /*--- Order in which we look for space in the colors to insert a new group. ---*/
@@ -672,4 +672,31 @@ su2double coloringEfficiency(const SparsePattern& coloring, int numThreads, int
   return ideal / real;
 }
 
+/*!
+ * \brief Compute the levels for the lower part of a sparse pattern.
+ * For example, corresponding to the dependencies of forward substitution.
+ */
+template <class T>
+T computeLevels(const T& pattern) {
+  using Index = typename T::IndexType;
+
+  std::vector<std::vector<Index>> levels;
+  {
+    const auto n = pattern.getOuterSize();
+    su2vector<int> level(n);
+    for (Index i = 0; i < n; ++i) {
+      level(i) = 0;
+      for (const auto j : pattern.getInnerIter(i)) {
+        if (j >= i) continue;
+        level(i) = std::max(level(i), level(j) + 1);
+      }
+      if (static_cast<std::size_t>(level(i) + 1) > levels.size()) {
+        levels.emplace_back();
+      }
+      levels[level(i)].push_back(i);
+    }
+  }
+  return T(levels);
+}
+
 /// @}

Common/src/CConfig.cpp

@@ -244,7 +244,7 @@ CConfig::CConfig(char case_filename[MAX_STRING_SIZE], CConfig *config) {
   /*--- Update original config file ---*/
 
   if (runtime_file) {
-    if (all_options.find("TIME_ITER") == all_options.end())
+    if (OptionIsSet("TIME_ITER"))
       config->SetnTime_Iter(nTimeIter);
   }
 }
@@ -1896,6 +1896,8 @@ void CConfig::SetConfig_Options() {
   addUnsignedLongOption("LINEAR_SOLVER_ITER", Linear_Solver_Iter, 10);
   /* DESCRIPTION: Fill in level for the ILU preconditioner */
   addUnsignedShortOption("LINEAR_SOLVER_ILU_FILL_IN", Linear_Solver_ILU_n, 0);
+  /* DESCRIPTION: Use level scheduling for OMP parallelization of the ILU preconditioner */
+  addBoolOption("LINEAR_SOLVER_ILU_LEVEL_SCHEDULING", Linear_Solver_ILU_levels, false);
   /* DESCRIPTION: Maximum number of iterations of the linear solver for the implicit formulation */
   addUnsignedLongOption("LINEAR_SOLVER_RESTART_FREQUENCY", Linear_Solver_Restart_Frequency, 10);
   /* DESCRIPTION: Number of vectors used for deflated restarts */
@@ -4047,7 +4049,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 
   /*--- Set the number of external iterations to 1 for the steady state problem ---*/
 
-  if (Kind_Solver == MAIN_SOLVER::FEM_ELASTICITY) nMGLevels = 0;
+  if (Kind_Solver == MAIN_SOLVER::FEM_ELASTICITY) {
+    nMGLevels = 0;
+    if (!OptionIsSet("LINEAR_SOLVER_ILU_LEVEL_SCHEDULING")) {
+      /*--- Different default behavior for this solver type. ---*/
+      Linear_Solver_ILU_levels = true;
+    }
+  }
 
   Radiation = (Kind_Radiation != RADIATION_MODEL::NONE);
 
@@ -5746,11 +5754,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     /*--- Check that spark ignition has required parameters defined ---*/
     if (flamelet_ParsedOptions.ignition_method == FLAMELET_INIT_TYPE::SPARK) {
       /*--- Check if SPARK_INIT was explicitly set in config file ---*/
-      if (all_options.find("SPARK_INIT") != all_options.end()) {
+      if (!OptionIsSet("SPARK_INIT")) {
         SU2_MPI::Error("FLAME_INIT_METHOD= SPARK requires SPARK_INIT to be defined in the config file.", CURRENT_FUNCTION);
       }
       /*--- Check if SPARK_REACTION_RATES was explicitly set in config file ---*/
-      if (all_options.find("SPARK_REACTION_RATES") != all_options.end()) {
+      if (!OptionIsSet("SPARK_REACTION_RATES")) {
         SU2_MPI::Error("FLAME_INIT_METHOD= SPARK requires SPARK_REACTION_RATES to be defined in the config file.", CURRENT_FUNCTION);
       }
       if (flamelet_ParsedOptions.nspark < flamelet_ParsedOptions.n_scalars) {

Common/src/geometry/CGeometry.cpp

@@ -4524,7 +4524,7 @@ su2double NearestNeighborDistance(CGeometry* geometry, const CConfig* config, co
    *    Compute an alternative distance based on volume and wall area. ---*/
 
   const auto nDim = geometry->GetnDim();
-  su2double Normal[3] = {};
+  su2double Normal[3] = {0, 0, 0};
   for (auto iMarker = 0u; iMarker < config->GetnMarker_All(); iMarker++) {
     if (!config->GetViscous_Wall(iMarker)) continue;