Fix develop tsan failure due to cantera build

.github/workflows/regression.yml

@@ -122,6 +122,8 @@ jobs:
         with:
           entrypoint: /bin/rm
           args: -rf install install_bin.tgz src ccache ${{ matrix.config_set }}
+      - name: Lower ASLR entropy for TSan
+        run: sudo sysctl vm.mmap_rnd_bits=28
       - name: Build
         uses: docker://ghcr.io/su2code/su2/build-su2-tsan:261011-1107
         with:
@@ -324,6 +326,8 @@ jobs:
           done
           chmod a+x $BIN_FOLDER/*
           ls -lahR $BIN_FOLDER
+      - name: Lower ASLR entropy for TSan
+        run: sudo sysctl vm.mmap_rnd_bits=28
       - name: Run Tests in Container
         uses: docker://ghcr.io/su2code/su2/test-su2-tsan:261011-1107
         env:

Common/include/CConfig.hpp

@@ -601,7 +601,8 @@ class CConfig {
   MUSCL_Heat,              /*!< \brief MUSCL scheme for the (fvm) heat equation.*/
   MUSCL_AdjFlow,           /*!< \brief MUSCL scheme for the adj flow equations.*/
   MUSCL_AdjTurb,           /*!< \brief MUSCL scheme for the adj turbulence equations.*/
-  MUSCL_Species;           /*!< \brief MUSCL scheme for the species equations.*/
+  MUSCL_Species,           /*!< \brief MUSCL scheme for the species equations.*/
+  MUSCL_FirstOrderWallEnds; /*!< \brief First-order reconstruction at no-slip wall ends and sharp wall edges.*/
   su2double MUSCL_Kappa,   /*!< \brief Blending coefficient for U-MUSCL scheme (for the runtime eq. system). */
   MUSCL_Kappa_Flow,        /*!< \brief Blending coefficient for U-MUSCL scheme for the flow equations.*/
   MUSCL_Kappa_Turb,        /*!< \brief Blending coefficient for U-MUSCL scheme for the turbulence equations.*/
@@ -1117,6 +1118,7 @@ class CConfig {
   su2double ParMETIS_tolerance;     /*!< \brief Load balancing tolerance for ParMETIS. */
   long ParMETIS_pointWgt;           /*!< \brief Load balancing weight given to points. */
   long ParMETIS_edgeWgt;            /*!< \brief Load balancing weight given to edges. */
+  bool ParMETIS_anisotropyWgt;      /*!< \brief Weight graph edges by cell anisotropy and keep anisotropic lines on one rank. */
   unsigned short DirectDiff;        /*!< \brief Direct Differentation mode. */
   bool DiscreteAdjoint,                /*!< \brief AD-based discrete adjoint mode. */
   DiscreteAdjointDebug;                /*!< \brief Discrete adjoint debug mode using tags. */
@@ -1177,6 +1179,8 @@ class CConfig {
   unsigned short *MG_PreSmooth_p{nullptr}, *MG_PostSmooth_p{nullptr}, *MG_CorrecSmooth_p{nullptr};
   unsigned short nMG_CflScaling_p{0};
   su2double *MG_CflScaling_p{nullptr};
+  unsigned short nMG_LinSolIter_p{0};
+  unsigned short *MG_LinSolIter_p{nullptr};
 
   ENUM_STREAMWISE_PERIODIC Kind_Streamwise_Periodic; /*!< \brief Kind of Streamwise periodic flow (pressure drop or massflow) */
   bool Streamwise_Periodic_Temperature;              /*!< \brief Use real periodicity for Energy equation or otherwise outlet source term. */
@@ -1752,6 +1756,10 @@ class CConfig {
    */
   bool GetCFL_Adapt(void) const { return CFL_Adapt; }
 
+  /*!
+   * \brief Whether the adaptive CFL follows the nonlinear residual trend instead of the linear solve.
+   */
+
   /*!
    * \brief Get the outlier mitigation parameters.
    */
@@ -4859,6 +4867,11 @@ class CConfig {
    */
   bool GetMUSCL_Flow(void) const { return MUSCL_Flow; }
 
+  /*!
+   * \brief Get if edges at no-slip wall ends and sharp wall edges are reconstructed to first order.
+   */
+  bool GetMUSCL_FirstOrderWallEnds(void) const { return MUSCL_FirstOrderWallEnds; }
+
   /*!
    * \brief Get if the upwind scheme used MUSCL or not.
    * \return MUSCL scheme.
@@ -10334,6 +10347,11 @@ class CConfig {
    */
   long GetParMETIS_EdgeWeight() const { return ParMETIS_edgeWgt; }
 
+  /*!
+   * \brief Whether ParMETIS weights graph edges by cell anisotropy and partitions with coordinates.
+   */
+  bool GetParMETIS_AnisotropyWeight() const { return ParMETIS_anisotropyWgt; }
+
   /*!
    * \brief Find the marker index (if any) that is part of a given interface pair.
    * \param[in] iInterface - Number of the interface pair being tested, starting at 0.

Common/include/containers/container_decorators.hpp

@@ -126,7 +126,8 @@ class C3DContainerDecorator {
 
  private:
   Storage m_storage;
-  Index m_innerSz;
+  /*--- One, not zero, so rows() on a container that was never resized does not divide by it. ---*/
+  Index m_innerSz = 1;
 
  public:
   C3DContainerDecorator() = default;

Common/include/geometry/CMultiGridGeometry.hpp

@@ -110,7 +110,8 @@ class CMultiGridGeometry final : public CGeometry {
     vector<std::array<su2double, MAXNDIM>> normal; /*!< \brief Unit normal there, pointing into the domain. */
     vector<su2double> strength;                    /*!< \brief Local anisotropy, the ordering key. */
     vector<char> tier;                             /*!< \brief Paving order: 0 viscous wall, 1 Euler
-                                                        wall, 2 anything else. */
+                                                        wall, 2 anything else; one more once the
+                                                        mixed nodes are put first. */
     unsigned long nRefusedCurvature = 0;           /*!< \brief Euler wall nodes the curvature limit kept out. */
   };
 

Common/include/geometry/CPhysicalGeometry.hpp

@@ -337,6 +337,42 @@ class CPhysicalGeometry final : public CGeometry {
    */
   void PrepareAdjacency(const CConfig* config);
 
+#if defined(HAVE_MPI) && defined(HAVE_PARMETIS)
+  /*!
+   * \brief ParMETIS edge weights that make cutting across the short edges of anisotropic cells expensive.
+   * \param[in] maxWeight - Cap on the weight of any edge.
+   * \return One weight per entry of the adjacency array.
+   */
+  vector<idx_t> ComputeAnisotropyEdgeWeights(long maxWeight) const;
+
+  /*!
+   * \brief Reads, for any global points of the linear partition, a value each rank holds for its own points.
+   * \param[in] wanted - Global indices to read.
+   * \param[in] localValue - Value of every local point.
+   * \return The value of each wanted point.
+   */
+  vector<unsigned long> FetchFromOwners(const vector<unsigned long>& wanted,
+                                        const vector<unsigned long>& localValue) const;
+
+  /*!
+   * \brief Joins points into lines along the strong anisotropic edges and names each line by its lowest global index.
+   * \param[in] adjwgt - Anisotropy weight of every entry of the adjacency array.
+   * \param[out] lineEdge - Whether each entry of the adjacency array is a line edge.
+   * \return Global index of the root of the line of every local point, the point itself when it is on no line.
+   */
+  vector<unsigned long> LinkAnisotropicLines(const vector<idx_t>& adjwgt, vector<char>& lineEdge) const;
+
+  /*!
+   * \brief Moves each line to the partition that holds most of its points.
+   * \param[in] root - Line roots from LinkAnisotropicLines.
+   * \param[in,out] part - Partition of every local point.
+   * \param[out] longestLine - Number of points of the longest line.
+   * \return Number of points that changed partition, summed over all ranks.
+   */
+  unsigned long KeepLinesTogether(const vector<unsigned long>& root, vector<idx_t>& part,
+                                  unsigned long& longestLine) const;
+#endif
+
   /*!
    * \brief Find repeated nodes between two elements to identify the common face.
    * \param[in] first_elem - Identification of the first element.

Common/include/linear_algebra/CSysMatrix.hpp

@@ -34,6 +34,7 @@
 #include "CPastixWrapper.hpp"
 #include "../toolboxes/graph_toolbox.hpp"
 
+#include <atomic>
 #include <cstdint>
 #include <cstdlib>
 #include <vector>
@@ -363,6 +364,19 @@ class CSysMatrix {
    * into ilu_level_ptr / d_precond_level_idx below for the GPU triangular solves. */
   CCompressedSparsePatternUL levels_ilu;
 
+  /*!< \brief Per-thread progress counter, padded to a cache line. */
+  struct alignas(64) LevelProgress {
+    std::atomic<uint64_t> steps{0};
+    LevelProgress() = default;
+    LevelProgress(const LevelProgress& other) : steps(other.steps.load()) {}
+  };
+
+  /*!< \brief Steps each thread must see completed by every other thread, per (step, thread, other thread),
+   * for the forward and backward level scheduled sweeps. Empty when barriers separate the levels. */
+  std::vector<uint32_t> level_waits_fwd, level_waits_bwd;
+  mutable std::vector<LevelProgress> level_progress;
+  unsigned long level_num_threads = 0;
+
   /*!< \brief Coloring of the (domain-only) ILU dependency graph, used only by the GPU iterative
    * factorization (see IluFactorColorKernel and ilu_color_ptr / d_ilu_color_idx below). The
    * host/OMP path and the GPU triangular solves use levels_ilu instead. */
@@ -666,6 +680,25 @@ class CSysMatrix {
   void ComputeJacobiPreconditionerGPU(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod,
                                       CGeometry* geometry, const CConfig* config) const;
 
+  /*!
+   * \brief Compute the wait tables of the point-to-point synchronized level scheduled sweeps.
+   */
+  void BuildLevelWaits();
+
+  /*!
+   * \brief Apply a function to the rows of all levels in increasing or decreasing order, waiting only for
+   * the threads owning the rows depended on. Must be called by all threads of the team.
+   */
+  template <class F>
+  void LevelSweep(bool forward, const F& row_func) const;
+
+  /*!
+   * \brief Whether the calling thread team matches the one the wait tables were built for.
+   */
+  inline bool UseLevelWaits() const {
+    return !level_waits_fwd.empty() && static_cast<unsigned long>(omp_get_num_threads()) == level_num_threads;
+  }
+
   /*!
    * \brief Build the ILU preconditioner on the device, from the device copy of the matrix.
    * \note Requires the device matrix to be up to date, see HtDTransfer.

Common/include/linear_algebra/CSysSolve.hpp

@@ -455,9 +455,10 @@ class CSysSolve {
    * \param[in,out] LinSysSol - Linear system solution
    * \param[in] geometry -  Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
+   * \param[in] maxIter - Iteration cap replacing LINEAR_SOLVER_ITER in standard mode, 0 keeps the config value.
    */
   unsigned long Solve(MatrixType& Jacobian, const CSysVector<su2double>& LinSysRes, CSysVector<su2double>& LinSysSol,
-                      CGeometry* geometry, const CConfig* config);
+                      CGeometry* geometry, const CConfig* config, unsigned long maxIter = 0);
 
   /*!
    * \brief Solve the adjoint linear system using a Krylov subspace method

Common/include/option_structure.hpp

@@ -80,6 +80,7 @@ const unsigned int MAX_NUMBER_FFD = 15;       /*!< \brief Maximum number of FFDB
 enum: unsigned int{MAX_SOLS = 14};            /*!< \brief Maximum number of solutions at the same time (dimension of solution container array). */
 const unsigned int MAX_TERMS = 7;             /*!< \brief Maximum number of terms in the numerical equations (dimension of solver container array). */
 const unsigned int MAX_ZONES = 3;             /*!< \brief Maximum number of zones. */
+const unsigned short MAX_MGLEVELS = 10;       /*!< \brief Maximum number of coarse multigrid levels, which bounds the per-level arrays of the multigrid integration. */
 const unsigned int MAX_FE_KINDS = 7;          /*!< \brief Maximum number of Finite Elements. */
 const unsigned int NO_RK_ITER = 0;            /*!< \brief No Runge-Kutta iteration. */
 
@@ -1246,10 +1247,13 @@ struct CMGOptions {
   std::vector<unsigned short> MG_PostSmooth;   /*!< \brief Multigrid post-smoothing iterations per level. */
   std::vector<unsigned short> MG_CorrecSmooth; /*!< \brief Multigrid Jacobi correction-smoothing per level. */
   std::vector<su2double> MG_CflScaling;        /*!< \brief Per-level CFL scaling factors relative to the previous (finer) level. Entry [i] scales level i+1 from level i. Size = nMGLevels. */
+  std::vector<unsigned long> MG_LinSolIter;    /*!< \brief Linear solver iteration cap per MG level (entry 0 = finest). */
   bool MG_Smooth_EarlyExit{false};        /*!< \brief Enable early exit for MG smoothing iterations. */
   bool MG_Smooth_Output{false};           /*!< \brief Output compact per-cycle smoothing summary. */
   su2double MG_Smooth_StagnationTol{0.0}; /*!< \brief Stagnation early exit: stop if current_rms >= prev_rms * tol. 0 = disabled. */
   bool MG_Implicit_Lines{false};          /*!< \brief Enable implicit-lines agglomeration from walls. */
+  bool MG_Linear_Prolongation{false};     /*!< \brief Prolong the correction with a limited least-squares
+                                                 gradient instead of piecewise-constant injection. */
   unsigned long MG_Startup_Iter{100};     /*!< \brief Iterations per mesh during FMG startup, and the length of each level's CFL ramp. 0 = no iteration budget. */
   su2double MG_Startup_Convergence{-2.0}; /*!< \brief FMG: orders of magnitude (log10) that CONV_FIELD must drop on the
                                                  active level before promoting to the next finer one. Negative is a
@@ -2940,6 +2944,7 @@ enum class MPI_QUANTITIES {
   COORDINATES          ,  /*!< \brief Vertex coordinates communication. */
   COORDINATES_OLD      ,  /*!< \brief Old vertex coordinates communication. */
   MAX_LENGTH           ,  /*!< \brief Maximum length communication. */
+  WALL_DISTANCE        ,  /*!< \brief Wall distance and roughness of the nearest wall communication. */
   GRID_VELOCITY        ,  /*!< \brief Grid velocity communication. */
   SOLUTION_EDDY        ,  /*!< \brief Turbulent solution plus eddy viscosity communication. */
   STOCH_SOURCE_LANG    ,  /*!< \brief Stochastic source term for Langevin equations communication. */

Common/src/CConfig.cpp

@@ -1930,6 +1930,7 @@ void CConfig::SetConfig_Options() {
   addULongArrayOption("OUTLIER_MITIGATION_PARAM", 4, true, outlierMitigationParam);
   /* DESCRIPTION: Activate The adaptive CFL number. */
   addBoolOption("CFL_ADAPT", CFL_Adapt, false);
+  /* DESCRIPTION: Adapt the CFL on the finest-grid nonlinear residuals instead of the linear solve, with a cut on divergence. */
   /* !\brief CFL_ADAPT_PARAM
    * DESCRIPTION: Parameters of the adaptive CFL number (factor down, factor up, CFL limit (min and max)[, acceptable linear residual][, starting iteration]).
    * Parameters in square brackets are optional, parameter "starting iteration" only valid with parameter "acceptable linear residual".
@@ -2129,6 +2130,10 @@ void CConfig::SetConfig_Options() {
   /*!\brief MG_IMPLICIT_LINES\n DESCRIPTION: Pave the coarse grid with advancing fronts raised from boundaries
    * that carry a stretched layer normal to themselves. DEFAULT: NO \ingroup Config*/
   addBoolOption("MG_IMPLICIT_LINES", MGOptions.MG_Implicit_Lines, false);
+  /*!\brief MG_LINEAR_PROLONGATION\n DESCRIPTION: Prolong the multigrid correction with a limited least-squares
+   * gradient over the coarse control volume instead of injecting the parent value into every child. Coarse CVs on
+   * walls and symmetry planes keep the constant operator. DEFAULT: NO \ingroup Config*/
+  addBoolOption("MG_LINEAR_PROLONGATION", MGOptions.MG_Linear_Prolongation, false);
   /*!\brief MG_STARTUP_ITER\n DESCRIPTION: Max number of iterations spent on each mesh during the Full
    * Multigrid (FMG) startup phase. DEFAULT: 100 \ingroup Config*/
   addUnsignedLongOption("MG_STARTUP_ITER", MGOptions.MG_Startup_Iter, 100);
@@ -2137,15 +2142,19 @@ void CConfig::SetConfig_Options() {
    * DEFAULT: -2 \ingroup Config*/
   addDoubleOption("MG_STARTUP_CONVERGENCE", MGOptions.MG_Startup_Convergence, -2.0);
   /*!\brief MG_STARTUP_STAGNATION\n DESCRIPTION: Full-MG promotion on stagnation. If the active level's residual ratio
-   * between successive iterations exceeds this value for MG_STARTUP_STAGNATION_ITER consecutive iterations, promote to
-   * the next finer level without waiting out MG_STARTUP_ITER. 0 disables it. DEFAULT: 0.99 \ingroup Config*/
+   * between successive iterations stays between this value and its inverse for MG_STARTUP_STAGNATION_ITER consecutive
+   * iterations, promote to the next finer level without waiting out MG_STARTUP_ITER. 0 disables it.
+   * DEFAULT: 0.99 \ingroup Config*/
   addDoubleOption("MG_STARTUP_STAGNATION", MGOptions.MG_Startup_Stagnation, 0.99);
   /*!\brief MG_STARTUP_STAGNATION_ITER\n DESCRIPTION: Consecutive stalled iterations required before Full-MG promotes
    * on stagnation. 0 disables it, as MG_STARTUP_STAGNATION= 0 does. DEFAULT: 5 \ingroup Config*/
   addUnsignedLongOption("MG_STARTUP_STAGNATION_ITER", MGOptions.MG_Startup_Stagnation_Iter, 5);
   /*!\brief MG_CFL_SCALING\n DESCRIPTION: Per-level CFL scaling factors for coarse MG levels. Entry i is the ratio CFL(i+1)/CFL(i).
    * If fewer values than nMGLevels are given, the last value is repeated. DEFAULT: 0.25 (i.e., 1/4 per level) \ingroup Config*/
   addDoubleListOption("MG_CFL_SCALING", nMG_CflScaling_p, MG_CflScaling_p);
+  /*!\brief MG_LINEAR_SOLVER_ITER\n DESCRIPTION: Maximum linear solver iterations on coarse MG levels. Entry i applies to level i+1;
+   * the last value is repeated, and the cap never exceeds LINEAR_SOLVER_ITER. DEFAULT: LINEAR_SOLVER_ITER \ingroup Config*/
+  addUShortListOption("MG_LINEAR_SOLVER_ITER", nMG_LinSolIter_p, MG_LinSolIter_p);
 
   /*!\par CONFIG_CATEGORY: Spatial Discretization \ingroup Config*/
   /*--- Options related to the spatial discretization ---*/
@@ -2176,6 +2185,10 @@ void CConfig::SetConfig_Options() {
 
   /*!\brief MUSCL_FLOW \n DESCRIPTION: Check if the MUSCL scheme should be used \ingroup Config*/
   addBoolOption("MUSCL_FLOW", MUSCL_Flow, true);
+  /*!\brief MUSCL_FIRST_ORDER_WALL_ENDS \n DESCRIPTION: Reconstruct the flow to first order on the edges of points
+   * where a no-slip wall ends on a boundary continuing its surface, or has a sharp edge such as a trailing edge, and
+   * skip the multigrid correction of the flow there \ingroup Config*/
+  addBoolOption("MUSCL_FIRST_ORDER_WALL_ENDS", MUSCL_FirstOrderWallEnds, false);
   /*!\brief MUSCL_KAPPA_FLOW \n DESCRIPTION: Blending coefficient for the U-MUSCL scheme \ingroup Config*/
   addDoubleOption("MUSCL_KAPPA_FLOW", MUSCL_Kappa_Flow, 0.0);
   /*!\brief RAMP_MUSCL \n DESCRIPTION: Enable ramping of the MUSCL scheme from 1st to 2nd order using specified method*/
@@ -3115,6 +3128,9 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: ParMETIS load balancing weight for edges (equiv. to neighbors) */
   addLongOption("PARMETIS_EDGE_WEIGHT", ParMETIS_edgeWgt, 1);
 
+  /* DESCRIPTION: Weight ParMETIS graph edges by cell anisotropy and keep anisotropic lines on one rank */
+  addBoolOption("PARMETIS_ANISOTROPY_WEIGHT", ParMETIS_anisotropyWgt, false);
+
   /*--- options that are used in the Hybrid RANS/LES Simulations  ---*/
   /*!\par CONFIG_CATEGORY:Hybrid_RANSLES Options\ingroup Config*/
 
@@ -4882,6 +4898,11 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     }
   }
 
+  if (RampMUSCL && !(MUSCL_Flow || MUSCL_Turb || MUSCL_Species || MUSCL_Heat)) {
+    if (rank == MASTER_NODE) cout << "WARNING: RAMP_MUSCL is ignored because no MUSCL reconstruction is active." << endl;
+    RampMUSCL = false;
+  }
+
   if(RampMUSCL && !DiscreteAdjoint){
     if (RampMUSCLParam.RampMUSCLPower <= 0.0) SU2_MPI::Error("RAMP_MUSCL_POWER cannot be less than or equal to zero!", CURRENT_FUNCTION);
     rampMUSCLValue = 0.0;
@@ -5084,6 +5105,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
   Kappa_2nd_AdjFlow = jst_adj_coeff[0];
   Kappa_4th_AdjFlow = jst_adj_coeff[1];
 
+  /*--- The multigrid integration carries per-level arrays of this size. ---*/
+
+  if (nMGLevels > MAX_MGLEVELS) {
+    SU2_MPI::Error("MGLEVEL is larger than the supported maximum of " + std::to_string(MAX_MGLEVELS) + ".",
+                   CURRENT_FUNCTION);
+  }
+
   /*--- Fill MG smooth vectors to size nMGLevels+1.
    Use parsed values (truncating or extending by repeat) or defaults if not set. ---*/
 
@@ -5114,6 +5142,13 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
       for (unsigned short i = 0; i < nMGLevels; ++i)
         MGOptions.MG_CflScaling[i] = 0.25;
     }
+
+    /*--- Linear solver iteration cap per level; level 0 always uses LINEAR_SOLVER_ITER. ---*/
+    MGOptions.MG_LinSolIter.assign(nMGLevels + 1, Linear_Solver_Iter);
+    for (unsigned short i = 1; i <= nMGLevels && nMG_LinSolIter_p != 0; ++i) {
+      const auto cap = MG_LinSolIter_p[min<unsigned short>(i - 1, nMG_LinSolIter_p - 1)];
+      MGOptions.MG_LinSolIter[i] = min<unsigned long>(cap, Linear_Solver_Iter);
+    }
   }
 
   /*--- Override MG Smooth parameters ---*/