Hybrid Parallel AD (Part 3/?)

Common/include/basic_types/ad_structure.hpp

@@ -252,7 +252,7 @@ namespace AD{
 
   /*!
    * \brief Start a passive region, i.e. stop recording.
-   * \return True is tape was active.
+   * \return True if tape was active.
    */
   inline bool BeginPassive() { return false; }
 
@@ -262,6 +262,18 @@ namespace AD{
    */
   inline void EndPassive(bool wasActive) {}
 
+  /*!
+   * \brief Pause the use of preaccumulation.
+   * \return True if preaccumulation was active.
+   */
+  inline bool PausePreaccumulation() { return false; }
+
+  /*!
+   * \brief Resume the use of preaccumulation.
+   * \param[in] wasActive - Whether preaccumulation was active before pausing.
+   */
+  inline void ResumePreaccumulation(bool wasActive) {}
+
 #else
   using CheckpointHandler = codi::DataStore;
 
@@ -273,7 +285,7 @@ namespace AD{
 
   extern bool PreaccActive;
 #ifdef HAVE_OPDI
-  #pragma omp threadprivate(PreaccActive)
+  SU2_OMP(threadprivate(PreaccActive))
 #endif
 
   extern bool PreaccEnabled;
@@ -292,7 +304,7 @@ namespace AD{
 
   extern codi::PreaccumulationHelper<su2double> PreaccHelper;
 #ifdef HAVE_OPDI
-  #pragma omp threadprivate(PreaccHelper)
+  SU2_OMP(threadprivate(PreaccHelper))
 #endif
 
   /*--- Reference to the tape. ---*/
@@ -527,6 +539,26 @@ namespace AD{
 
   FORCEINLINE void EndPassive(bool wasActive) { if(wasActive) StartRecording(); }
 
+  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
+    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
+  }
+
 #endif // CODI_REVERSE_TYPE
 
 } // namespace AD

Common/include/toolboxes/graph_toolbox.hpp

@@ -527,7 +527,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 = 32, size_t MaxMB = 128>
+template<class T, typename Color_t = char, size_t MaxColors = 64, size_t MaxMB = 128>
 T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors = false,
                      std::vector<Color_t>* indexColor = nullptr)
 {

Common/src/basic_types/ad_structure.cpp

@@ -36,14 +36,14 @@ namespace AD {
 
   bool PreaccActive = false;
 #ifdef HAVE_OPDI
-  #pragma omp threadprivate(PreaccActive)
+  SU2_OMP(threadprivate(PreaccActive))
 #endif
 
   bool PreaccEnabled = true;
 
   codi::PreaccumulationHelper<su2double> PreaccHelper;
 #ifdef HAVE_OPDI
-  #pragma omp threadprivate(PreaccHelper)
+  SU2_OMP(threadprivate(PreaccHelper))
 #endif
 
   ExtFuncHelper* FuncHelper;

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -319,7 +319,11 @@ void CFVMFlowSolverBase<V, R>::HybridParallelInitialization(const CConfig& confi
       cout << "WARNING: On " << numRanksUsingReducer << " MPI ranks the coloring efficiency was less than "
            << COLORING_EFF_THRESH << " (min value was " << minEff << ").\n"
            << "         Those ranks will now use a fallback strategy, better performance may be possible\n"
-           << "         with a different value of config option EDGE_COLORING_GROUP_SIZE (default 512)." << endl;
+           << "         with a different value of config option EDGE_COLORING_GROUP_SIZE (default 512)."
+#ifdef HAVE_OPDI
+           << "\n         The memory usage of the discrete adjoint solver is higher when using the fallback."
+#endif
+           << endl;
     }
 
     if (config.GetUseVectorization() && (omp_get_max_threads() > 1) &&
@@ -1535,6 +1539,10 @@ void CFVMFlowSolverBase<V, R>::EdgeFluxResidual(const CGeometry *geometry,
     InstantiateEdgeNumerics(solvers, config);
   }
 
+#ifdef HAVE_OPDI
+  const auto preaccEnabled = ReducerStrategy && AD::PausePreaccumulation();
+#endif
+
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring) {
     /*--- Chunk size is at least OMP_MIN_SIZE and a multiple of the color group size. ---*/
@@ -1557,6 +1565,10 @@ void CFVMFlowSolverBase<V, R>::EdgeFluxResidual(const CGeometry *geometry,
     END_SU2_OMP_FOR
   }
 
+#ifdef HAVE_OPDI
+  AD::ResumePreaccumulation(preaccEnabled);
+#endif
+
   if (ReducerStrategy) {
     SumEdgeFluxes(geometry);
     if (config->GetKind_TimeIntScheme() == EULER_IMPLICIT) {

SU2_CFD/src/integration/CIntegration.cpp

@@ -77,9 +77,7 @@ void CIntegration::Space_Integration(CGeometry *geometry,
   CNumerics* visc_bound_numerics = numerics[VISC_BOUND_TERM + omp_get_thread_num()*MAX_TERMS];
 
 #ifdef HAVE_OPDI
-  /* disable preaccumulation for parallel boundary numerics */
-  bool preaccEnabled = AD::PreaccEnabled;
-  AD::PreaccEnabled = false;
+  const auto preaccEnabled = AD::PausePreaccumulation();
 #endif
 
   /*--- Boundary conditions that depend on other boundaries (they require MPI sincronization)---*/
@@ -185,7 +183,7 @@ void CIntegration::Space_Integration(CGeometry *geometry,
   }
 
 #ifdef HAVE_OPDI
-  AD::PreaccEnabled = preaccEnabled;
+  AD::ResumePreaccumulation(preaccEnabled);
 #endif
 
 }

SU2_CFD/src/solvers/CEulerSolver.cpp

@@ -1992,6 +1992,10 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
   su2double Primitive_i[MAXNVAR] = {0.0}, Primitive_j[MAXNVAR] = {0.0};
   su2double Secondary_i[MAXNVAR] = {0.0}, Secondary_j[MAXNVAR] = {0.0};
 
+#ifdef HAVE_OPDI
+  const auto preaccEnabled = ReducerStrategy && AD::PausePreaccumulation();
+#endif
+
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring)
   {
@@ -2176,6 +2180,10 @@ void CEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_contain
   END_SU2_OMP_FOR
   } // end color loop
 
+#ifdef HAVE_OPDI
+  AD::ResumePreaccumulation(preaccEnabled);
+#endif
+
   if (ReducerStrategy) {
     SumEdgeFluxes(geometry);
     if (implicit)

SU2_CFD/src/solvers/CIncEulerSolver.cpp

@@ -1018,6 +1018,10 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
   bool implicit    = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT);
   bool jst_scheme  = ((config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0));
 
+#ifdef HAVE_OPDI
+  const auto preaccEnabled = ReducerStrategy && AD::PausePreaccumulation();
+#endif
+
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring)
   {
@@ -1082,6 +1086,10 @@ void CIncEulerSolver::Centered_Residual(CGeometry *geometry, CSolver **solver_co
   END_SU2_OMP_FOR
   } // end color loop
 
+#ifdef HAVE_OPDI
+  AD::ResumePreaccumulation(preaccEnabled);
+#endif
+
   if (ReducerStrategy) {
     SumEdgeFluxes(geometry);
     if (implicit)
@@ -1110,6 +1118,10 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
   const bool limiter    = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER);
   const bool van_albada = (config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE);
 
+#ifdef HAVE_OPDI
+  const auto preaccEnabled = ReducerStrategy && AD::PausePreaccumulation();
+#endif
+
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring)
   {
@@ -1250,6 +1262,10 @@ void CIncEulerSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_cont
   END_SU2_OMP_FOR
   } // end color loop
 
+#ifdef HAVE_OPDI
+  AD::ResumePreaccumulation(preaccEnabled);
+#endif
+
   if (ReducerStrategy) {
     SumEdgeFluxes(geometry);
     if (implicit)

SU2_CFD/src/solvers/CTurbSolver.cpp

@@ -106,6 +106,10 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe
   su2double solution_i[MAXNVAR] = {0.0}, flowPrimVar_i[MAXNVARFLOW] = {0.0};
   su2double solution_j[MAXNVAR] = {0.0}, flowPrimVar_j[MAXNVARFLOW] = {0.0};
 
+#ifdef HAVE_OPDI
+  const auto preaccEnabled = ReducerStrategy && AD::PausePreaccumulation();
+#endif
+
   /*--- Loop over edge colors. ---*/
   for (auto color : EdgeColoring)
   {
@@ -232,6 +236,10 @@ void CTurbSolver::Upwind_Residual(CGeometry *geometry, CSolver **solver_containe
   END_SU2_OMP_FOR
   } // end color loop
 
+#ifdef HAVE_OPDI
+  AD::ResumePreaccumulation(preaccEnabled);
+#endif
+
   if (ReducerStrategy) {
     SumEdgeFluxes(geometry);
     if (implicit) Jacobian.SetDiagonalAsColumnSum();