Merge branch 'develop' into Dependabot-GitHub-Actions

Author: pcarruscag

.github/workflows/release-management.yml

@@ -35,7 +35,7 @@ jobs:
           key: ${{ matrix.os_bin }}-${{ github.sha }}
           restore-keys: ${{ matrix.os_bin }}
       - name: Build
-        uses: docker://su2code/build-su2-cross:latest
+        uses: docker://ghcr.io/su2code/su2/build-su2-cross:220716-1459
         with:
           args: -b ${{ github.sha }} -f "${{matrix.flags}}"
       - name: Create Archive

Common/include/basic_types/ad_structure.hpp

@@ -559,21 +559,13 @@ namespace AD{
   FORCEINLINE bool PausePreaccumulation() {
     const auto current = PreaccEnabled;
     if (!current) return false;
-    SU2_OMP_BARRIER
-    SU2_OMP_MASTER
-    PreaccEnabled = false;
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
     return true;
   }
 
   FORCEINLINE void ResumePreaccumulation(bool wasActive) {
     if (!wasActive) return;
-    SU2_OMP_BARRIER
-    SU2_OMP_MASTER
-    PreaccEnabled = true;
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
   }
 
   FORCEINLINE void StartNoSharedReading() {

Common/include/geometry/dual_grid/CEdge.hpp

@@ -43,6 +43,7 @@ class CEdge {
   using NodeArray = C2DContainer<Index, Index, StorageType::ColumnMajor, 64, DynamicSize, 2>;
   NodeArray Nodes;           /*!< \brief Vector to store the node indices of the edge. */
   su2activematrix Normal;    /*!< \brief Normal (area) of the edge. */
+  const Index nEdge, nEdgeSIMD;
 
   friend class CPhysicalGeometry;
 
@@ -70,13 +71,27 @@ class CEdge {
   inline unsigned long GetNode(unsigned long iEdge, unsigned long iNode) const { return Nodes(iEdge,iNode); }
 
   /*!
-   * \brief SIMD version of GetNode, iNode returned for multiple contiguous iEdges
+   * \brief SIMD version of GetNode, iNode returned for contiguous iEdges.
    */
   template<class T, size_t N>
   FORCEINLINE simd::Array<T,N> GetNode(simd::Array<T,N> iEdge, unsigned long iNode) const {
     return simd::Array<T,N>(&Nodes(iEdge[0],iNode));
   }
 
+  /*!
+   * \brief Sets the tail of "Nodes" to repeat one of the last edges.
+   * \note This is needed when using the SIMD version of GetNode and
+   *       the number of edges is not a multiple of the simd width.
+   */
+  void SetPaddingNodes() {
+    for (auto iEdge = nEdge; iEdge < nEdgeSIMD; ++iEdge) {
+      /*--- Pad nodes by repeating the first edge in the last SIMD group. ---*/
+      const auto iEdge0 = nEdgeSIMD - simd::preferredLen<su2double>();
+      Nodes(iEdge, LEFT) = Nodes(iEdge0, LEFT);
+      Nodes(iEdge, RIGHT) = Nodes(iEdge0, RIGHT);
+    }
+  }
+
   /*!
    * \brief Set the node indices of an edge.
    * \param[in] iEdge - Edge index.

Common/include/linear_algebra/CSysSolve.hpp

@@ -219,12 +219,11 @@ class CSysSolve {
   void HandleTemporariesIn(const CSysVector<OtherType>& LinSysRes, CSysVector<OtherType>& LinSysSol) {
 
     /*--- Set the pointers. ---*/
-    SU2_OMP_MASTER {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       LinSysRes_ptr = &LinSysRes;
       LinSysSol_ptr = &LinSysSol;
     }
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
   /*!
@@ -241,12 +240,11 @@ class CSysSolve {
     LinSysSol_tmp.PassiveCopy(LinSysSol);
 
     /*--- Set the pointers. ---*/
-    SU2_OMP_MASTER {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       LinSysRes_ptr = &LinSysRes_tmp;
       LinSysSol_ptr = &LinSysSol_tmp;
     }
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
   /*!
@@ -258,13 +256,11 @@ class CSysSolve {
   void HandleTemporariesOut(CSysVector<OtherType>& LinSysSol) {
 
     /*--- Reset the pointers. ---*/
-    SU2_OMP_BARRIER
-    SU2_OMP_MASTER {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       LinSysRes_ptr = nullptr;
       LinSysSol_ptr = nullptr;
     }
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
   /*!
@@ -279,13 +275,11 @@ class CSysSolve {
     LinSysSol.PassiveCopy(LinSysSol_tmp);
 
     /*--- Reset the pointers. ---*/
-    SU2_OMP_BARRIER
-    SU2_OMP_MASTER {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       LinSysRes_ptr = nullptr;
       LinSysSol_ptr = nullptr;
     }
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
   }
 
 public:

Common/include/linear_algebra/CSysSolve_b.hpp

@@ -32,8 +32,8 @@
 #ifdef CODI_REVERSE_TYPE
 template<class ScalarType>
 struct CSysSolve_b {
-  static void Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::Real* x_b, size_t m,
-                      const codi::RealReverse::Real* y, const codi::RealReverse::Real* y_b, size_t n,
+  static void Solve_b(const su2double::Real* x, su2double::Real* x_b, size_t m,
+                      const su2double::Real* y, const su2double::Real* y_b, size_t n,
                       codi::DataStore* d);
 };
 #endif

Common/include/linear_algebra/CSysVector.hpp

@@ -186,10 +186,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     /*--- check if self-assignment, otherwise perform deep copy ---*/
     if ((const void*)this == (const void*)&other) return;
 
-    SU2_OMP_MASTER
-    Initialize(other.GetNBlk(), other.GetNBlkDomain(), other.GetNVar(), nullptr, true, false);
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    SU2_OMP_SAFE_GLOBAL_ACCESS(Initialize(other.GetNBlk(), other.GetNBlkDomain(), other.GetNVar(), nullptr, true, false);)
 
     CSYSVEC_PARFOR
     for (auto i = 0ul; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]);
@@ -297,11 +294,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   ScalarType dot(const VecExpr::CVecExpr<T, ScalarType>& expr) const {
     static ScalarType dotRes;
     /*--- All threads get the same "view" of the vectors and shared variable. ---*/
-    SU2_OMP_BARRIER
-    SU2_OMP_MASTER
-    dotRes = 0.0;
-    END_SU2_OMP_MASTER
-    SU2_OMP_BARRIER
+    SU2_OMP_SAFE_GLOBAL_ACCESS(dotRes = 0.0;)
 
     /*--- Local dot product for each thread. ---*/
     ScalarType sum = 0.0;
@@ -317,16 +310,16 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
 
 #ifdef HAVE_MPI
     /*--- Reduce across all mpi ranks, only master thread communicates. ---*/
-    SU2_OMP_BARRIER
-    SU2_OMP_MASTER {
+    BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
       sum = dotRes;
       const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE;
       SelectMPIWrapper<ScalarType>::W::Allreduce(&sum, &dotRes, 1, mpi_type, MPI_SUM, SU2_MPI::GetComm());
     }
-    END_SU2_OMP_MASTER
-#endif
+    END_SU2_OMP_SAFE_GLOBAL_ACCESS
+#else
     /*--- Make view of result consistent across threads. ---*/
     SU2_OMP_BARRIER
+#endif
 
     return dotRes;
   }

Common/include/option_structure.hpp

@@ -613,11 +613,13 @@ MakePair("ONESPECIES", ONESPECIES)
  * \brief types of coefficient transport model
  */
 enum class TRANSCOEFFMODEL {
+  SUTHERLAND,
   WILKE,
   GUPTAYOS,
   CHAPMANN_ENSKOG
 };
 static const MapType<std::string, TRANSCOEFFMODEL> TransCoeffModel_Map = {
+MakePair("SUTHERLAND", TRANSCOEFFMODEL::SUTHERLAND)
 MakePair("WILKE", TRANSCOEFFMODEL::WILKE)
 MakePair("GUPTA-YOS", TRANSCOEFFMODEL::GUPTAYOS)
 MakePair("CHAPMANN-ENSKOG", TRANSCOEFFMODEL::CHAPMANN_ENSKOG)

Common/include/parallelization/omp_structure.hpp

@@ -185,6 +185,25 @@ void omp_finalize();
 
 #endif
 
+/* The SU2_OMP_SAFE_GLOBAL_ACCESS constructs are used to safeguard code that should only be executed by the master
+ * thread, with all threads and memory views synchronized both beforehand and afterwards.
+ */
+
+#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS \
+  SU2_OMP_BARRIER \
+  SU2_OMP_MASTER
+
+#define END_SU2_OMP_SAFE_GLOBAL_ACCESS \
+  END_SU2_OMP_MASTER \
+  SU2_OMP_BARRIER
+
+#define SU2_OMP_SAFE_GLOBAL_ACCESS(...) \
+  BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS \
+  { \
+    __VA_ARGS__ \
+  } \
+  END_SU2_OMP_SAFE_GLOBAL_ACCESS
+
 /*--- Convenience functions (e.g. to compute chunk sizes). ---*/
 
 /*!

Common/include/parallelization/vectorization.hpp

@@ -56,6 +56,7 @@ constexpr size_t PREFERRED_SIZE = 8;
  */
 template<class T>
 constexpr size_t preferredLen() { return PREFERRED_SIZE / sizeof(T); }
+
 template<>
 constexpr size_t preferredLen<su2double>() { return PREFERRED_SIZE / sizeof(passivedouble); }
 

Common/src/CConfig.cpp

@@ -3792,8 +3792,8 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
     SU2_MPI::Error("Only STANDARD_AIR fluid model can be used with US Measurement System", CURRENT_FUNCTION);
   }
 
-  if (Kind_FluidModel == SU2_NONEQ && Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE ) {
-    SU2_MPI::Error("Only WILKE transport model is stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", CURRENT_FUNCTION);
+  if (Kind_FluidModel == SU2_NONEQ && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::SUTHERLAND) ) {
+    SU2_MPI::Error("Only WILKE and SUTHERLAND transport models are stable for the NEMO solver using SU2TClib. Use Mutation++ instead.", CURRENT_FUNCTION);
   }
 
   if (Kind_FluidModel == MUTATIONPP && (Kind_TransCoeffModel != TRANSCOEFFMODEL::WILKE && Kind_TransCoeffModel != TRANSCOEFFMODEL::CHAPMANN_ENSKOG)) {