su2code · jblueh · May 23, 2023 · May 15, 2023 · May 15, 2023 · May 15, 2023
diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp
@@ -74,6 +74,11 @@ inline void RegisterInput(su2double& data, bool push_index = true) {}
  */
 inline void RegisterOutput(su2double& data) {}
 
+/*!
+ * \brief Resize the adjoint vector, for subsequent access without bounds checking.
+ */
+inline void ResizeAdjoints() {}
+
 /*!
  * \brief Sets the adjoint value at index to val
  * \param[in] index - Position in the adjoint vector.
@@ -369,11 +374,24 @@ FORCEINLINE void Reset() {
   }
 }
 
+FORCEINLINE void ResizeAdjoints() { AD::getTape().resizeAdjointVector(); }
+
 FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); }
 
-FORCEINLINE void SetDerivative(int index, const double val) { AD::getTape().setGradient(index, val); }
+// WARNING: For performance reasons, this method does not perform bounds checking.
+// When using it, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
+FORCEINLINE void SetDerivative(int index, const double val) {
+  using BoundsChecking = codi::GradientAccessTapeInterface<su2double::Gradient, su2double::Identifier>::BoundsChecking;
+  AD::getTape().setGradient(index, val, BoundsChecking::False);
+}
 
-FORCEINLINE double GetDerivative(int index) { return AD::getTape().getGradient(index); }
+// WARNING: For performance reasons, this method does not perform bounds checking.
+// If called after tape evaluations, the adjoints should exist.
+// Otherwise, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
+FORCEINLINE double GetDerivative(int index) {
+  using BoundsChecking = codi::GradientAccessTapeInterface<su2double::Gradient, su2double::Identifier>::BoundsChecking;
+  return AD::getTape().getGradient(index, BoundsChecking::False);
+}
 
 FORCEINLINE bool IsIdentifierActive(su2double const& value) {
   return getTape().isIdentifierActive(value.getIdentifier());
@@ -523,26 +541,14 @@ FORCEINLINE void delete_handler(void* handler) {
 
 FORCEINLINE bool BeginPassive() {
   if (AD::getTape().isActive()) {
-    StopRecording();
+    AD::getTape().setPassive();
     return true;
   }
   return false;
 }
 
 FORCEINLINE void EndPassive(bool wasActive) {
-  if (wasActive) StartRecording();
-}
-
-FORCEINLINE bool PausePreaccumulation() {
-  const auto current = PreaccEnabled;
-  if (!current) return false;
-  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
-  return true;
-}
-
-FORCEINLINE void ResumePreaccumulation(bool wasActive) {
-  if (!wasActive) return;
-  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
+  if (wasActive) AD::getTape().setActive();
 }
 
 FORCEINLINE void StartNoSharedReading() {
@@ -558,6 +564,19 @@ FORCEINLINE void EndNoSharedReading() {
   opdi::logic->addReverseBarrier();
 #endif
 }
+
+FORCEINLINE bool PausePreaccumulation() {
+  const auto current = PreaccEnabled;
+  if (!current) return false;
+  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
+  return true;
+}
+
+FORCEINLINE void ResumePreaccumulation(bool wasActive) {
+  if (!wasActive) return;
+  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
+}
+
 #endif  // CODI_REVERSE_TYPE
 
 void Initialize();

diff --git a/Common/include/code_config.hpp b/Common/include/code_config.hpp
@@ -96,7 +96,7 @@ FORCEINLINE Out su2staticcast_p(In ptr) {
 #include "codi/tools/data/externalFunctionUserData.hpp"
 
 #if defined(HAVE_OMP)
-using su2double = codi::RealReverseIndexOpenMP;
+using su2double = codi::RealReverseIndexOpenMPGen<double, double>;
 #else
 #if defined(CODI_INDEX_TAPE)
 using su2double = codi::RealReverseIndex;

diff --git a/Common/include/parallelization/omp_structure.hpp b/Common/include/parallelization/omp_structure.hpp
@@ -188,12 +188,14 @@ void omp_finalize();
  * thread, with all threads and memory views synchronized both beforehand and afterwards.
  */
 
-#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS \
-  SU2_OMP_BARRIER                        \
+#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS             \
+  SU2_OMP_BARRIER                                    \
+  if (omp_in_parallel()) AD::StartNoSharedReading(); \
   SU2_OMP_MASTER
 
-#define END_SU2_OMP_SAFE_GLOBAL_ACCESS \
-  END_SU2_OMP_MASTER                   \
+#define END_SU2_OMP_SAFE_GLOBAL_ACCESS             \
+  END_SU2_OMP_MASTER                               \
+  if (omp_in_parallel()) AD::EndNoSharedReading(); \
   SU2_OMP_BARRIER
 
 #define SU2_OMP_SAFE_GLOBAL_ACCESS(...) BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS{__VA_ARGS__} END_SU2_OMP_SAFE_GLOBAL_ACCESS

diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp
@@ -192,7 +192,7 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
                                                   const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                   unsigned long m, ScalarType& residual, bool monitoring,
                                                   const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   ScalarType norm_r = 0.0, norm0 = 0.0;
   unsigned long i = 0;
 
@@ -241,16 +241,22 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
     if (tol_type == LinearToleranceType::RELATIVE) norm0 = norm_r;
 
     if ((norm_r < tol * norm0) || (norm_r < eps)) {
-      if (master && (lin_sol_mode != LINEAR_SOLVER_MODE::MESH_DEFORM))
+      if (masterRank && (lin_sol_mode != LINEAR_SOLVER_MODE::MESH_DEFORM)) {
+        SU2_OMP_MASTER
         cout << "CSysSolve::ConjugateGradient(): system solved by initial guess." << endl;
+        END_SU2_OMP_MASTER
+      }
       return 0;
     }
 
     /*--- Output header information including initial residual ---*/
 
-    if (monitoring && master) {
-      WriteHeader("CG", tol, norm_r);
-      WriteHistory(i, norm_r / norm0);
+    if (monitoring && masterRank) {
+      SU2_OMP_MASTER {
+        WriteHeader("CG", tol, norm_r);
+        WriteHistory(i, norm_r / norm0);
+      }
+      END_SU2_OMP_MASTER
     }
   }
 
@@ -281,7 +287,11 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
 
       norm_r = r.norm();
       if (norm_r < tol * norm0) break;
-      if (((monitoring) && (master)) && ((i + 1) % monitorFreq == 0)) WriteHistory(i + 1, norm_r / norm0);
+      if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
+        SU2_OMP_MASTER
+        WriteHistory(i + 1, norm_r / norm0);
+        END_SU2_OMP_MASTER
+      }
     }
 
     precond(r, z);
@@ -300,16 +310,22 @@ unsigned long CSysSolve<ScalarType>::CG_LinSolver(const CSysVector<ScalarType>&
   /*--- Recalculate final residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {
-    if (master) WriteFinalResidual("CG", i, norm_r / norm0);
+    if (masterRank) {
+      SU2_OMP_MASTER
+      WriteFinalResidual("CG", i, norm_r / norm0);
+      END_SU2_OMP_MASTER
+    }
 
     if (recomputeRes) {
       mat_vec(x, A_x);
       r = b - A_x;
       ScalarType true_res = r.norm();
 
       if (fabs(true_res - norm_r) > tol * 10.0) {
-        if (master) {
+        if (masterRank) {
+          SU2_OMP_MASTER
           WriteWarning(norm_r, true_res, tol);
+          END_SU2_OMP_MASTER
         }
       }
     }
@@ -325,7 +341,7 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
                                                       const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                       unsigned long m, ScalarType& residual, bool monitoring,
                                                       const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   const bool flexible = !precond.IsIdentity();
 
   /*---  Check the subspace size ---*/
@@ -386,7 +402,11 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   if ((beta < tol * norm0) || (beta < eps)) {
     /*--- System is already solved ---*/
 
-    if (master) cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl;
+    if (masterRank) {
+      SU2_OMP_MASTER
+      cout << "CSysSolve::FGMRES(): system solved by initial guess." << endl;
+      END_SU2_OMP_MASTER
+    }
     residual = beta;
     return 0;
   }
@@ -403,9 +423,12 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   /*--- Output header information including initial residual ---*/
 
   unsigned long i = 0;
-  if ((monitoring) && (master)) {
-    WriteHeader("FGMRES", tol, beta);
-    WriteHistory(i, beta / norm0);
+  if ((monitoring) && (masterRank)) {
+    SU2_OMP_MASTER {
+      WriteHeader("FGMRES", tol, beta);
+      WriteHistory(i, beta / norm0);
+    }
+    END_SU2_OMP_MASTER
   }
 
   /*---  Loop over all search directions ---*/
@@ -444,7 +467,11 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
 
     /*---  Output the relative residual if necessary ---*/
 
-    if ((((monitoring) && (master)) && ((i + 1) % monitorFreq == 0)) && (master)) WriteHistory(i + 1, beta / norm0);
+    if ((((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0))) {
+      SU2_OMP_MASTER
+      WriteHistory(i + 1, beta / norm0);
+      END_SU2_OMP_MASTER
+    }
   }
 
   /*---  Solve the least-squares system and update solution ---*/
@@ -460,16 +487,22 @@ unsigned long CSysSolve<ScalarType>::FGMRES_LinSolver(const CSysVector<ScalarTyp
   /*---  Recalculate final (neg.) residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {
-    if (master) WriteFinalResidual("FGMRES", i, beta / norm0);
+    if (masterRank) {
+      SU2_OMP_MASTER
+      WriteFinalResidual("FGMRES", i, beta / norm0);
+      END_SU2_OMP_MASTER
+    }
 
     if (recomputeRes) {
       mat_vec(x, W[0]);
       W[0] -= b;
       ScalarType res = W[0].norm();
 
       if (fabs(res - beta) > tol * 10) {
-        if (master) {
+        if (masterRank) {
+          SU2_OMP_MASTER
           WriteWarning(beta, res, tol);
+          END_SU2_OMP_MASTER
         }
       }
     }
@@ -511,7 +544,7 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
                                                        const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                        unsigned long m, ScalarType& residual, bool monitoring,
                                                        const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   ScalarType norm_r = 0.0, norm0 = 0.0;
   unsigned long i = 0;
 
@@ -561,15 +594,22 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
     if (tol_type == LinearToleranceType::RELATIVE) norm0 = norm_r;
 
     if ((norm_r < tol * norm0) || (norm_r < eps)) {
-      if (master) cout << "CSysSolve::BCGSTAB(): system solved by initial guess." << endl;
+      if (masterRank) {
+        SU2_OMP_MASTER
+        cout << "CSysSolve::BCGSTAB(): system solved by initial guess." << endl;
+        END_SU2_OMP_MASTER
+      }
       return 0;
     }
 
     /*--- Output header information including initial residual ---*/
 
-    if ((monitoring) && (master)) {
-      WriteHeader("BCGSTAB", tol, norm_r);
-      WriteHistory(i, norm_r / norm0);
+    if ((monitoring) && (masterRank)) {
+      SU2_OMP_MASTER {
+        WriteHeader("BCGSTAB", tol, norm_r);
+        WriteHistory(i, norm_r / norm0);
+      }
+      END_SU2_OMP_MASTER
     }
   }
 
@@ -637,22 +677,32 @@ unsigned long CSysSolve<ScalarType>::BCGSTAB_LinSolver(const CSysVector<ScalarTy
 
       norm_r = r.norm();
       if (norm_r < tol * norm0) break;
-      if (((monitoring) && (master)) && ((i + 1) % monitorFreq == 0) && (master)) WriteHistory(i + 1, norm_r / norm0);
+      if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
+        SU2_OMP_MASTER
+        WriteHistory(i + 1, norm_r / norm0);
+        END_SU2_OMP_MASTER
+      }
     }
   }
 
   /*--- Recalculate final residual (this should be optional) ---*/
 
   if ((monitoring) && (config->GetComm_Level() == COMM_FULL)) {
-    if (master) WriteFinalResidual("BCGSTAB", i, norm_r / norm0);
+    if (masterRank) {
+      SU2_OMP_MASTER
+      WriteFinalResidual("BCGSTAB", i, norm_r / norm0);
+      END_SU2_OMP_MASTER
+    }
 
     if (recomputeRes) {
       mat_vec(x, A_x);
       r = b - A_x;
       ScalarType true_res = r.norm();
 
-      if ((fabs(true_res - norm_r) > tol * 10.0) && (master)) {
+      if ((fabs(true_res - norm_r) > tol * 10.0) && (masterRank)) {
+        SU2_OMP_MASTER
         WriteWarning(norm_r, true_res, tol);
+        END_SU2_OMP_MASTER
       }
     }
   }
@@ -667,7 +717,7 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
                                                         const CPreconditioner<ScalarType>& precond, ScalarType tol,
                                                         unsigned long m, ScalarType& residual, bool monitoring,
                                                         const CConfig* config) const {
-  const bool master = (SU2_MPI::GetRank() == MASTER_NODE) && (omp_get_thread_num() == 0);
+  const bool masterRank = (SU2_MPI::GetRank() == MASTER_NODE);
   const bool fix_iter_mode = tol < eps;
   ScalarType norm_r = 0.0, norm0 = 0.0;
   unsigned long i = 0;
@@ -717,15 +767,22 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
     if (tol_type == LinearToleranceType::RELATIVE) norm0 = norm_r;
 
     if ((norm_r < tol * norm0) || (norm_r < eps)) {
-      if (master) cout << "CSysSolve::Smoother_LinSolver(): system solved by initial guess." << endl;
+      if (masterRank) {
+        SU2_OMP_MASTER
+        cout << "CSysSolve::Smoother_LinSolver(): system solved by initial guess." << endl;
+        END_SU2_OMP_MASTER
+      }
       return 0;
     }
 
     /*--- Output header information including initial residual. ---*/
 
-    if ((monitoring) && (master)) {
-      WriteHeader("Smoother", tol, norm_r);
-      WriteHistory(i, norm_r / norm0);
+    if ((monitoring) && (masterRank)) {
+      SU2_OMP_MASTER {
+        WriteHeader("Smoother", tol, norm_r);
+        WriteHistory(i, norm_r / norm0);
+      }
+      END_SU2_OMP_MASTER
     }
   }
 
@@ -759,14 +816,20 @@ unsigned long CSysSolve<ScalarType>::Smoother_LinSolver(const CSysVector<ScalarT
     if (!fix_iter_mode && config->GetComm_Level() == COMM_FULL) {
       norm_r = r.norm();
       if (norm_r < tol * norm0) break;
-      if (((monitoring) && (master)) && ((i + 1) % monitorFreq == 0)) WriteHistory(i + 1, norm_r / norm0);
+      if (((monitoring) && (masterRank)) && ((i + 1) % monitorFreq == 0)) {
+        SU2_OMP_MASTER
+        WriteHistory(i + 1, norm_r / norm0);
+        END_SU2_OMP_MASTER
+      }
     }
   }
 
   if (fix_iter_mode) norm_r = r.norm();
 
-  if ((monitoring) && (master) && (config->GetComm_Level() == COMM_FULL)) {
+  if ((monitoring) && (masterRank) && (config->GetComm_Level() == COMM_FULL)) {
+    SU2_OMP_MASTER
     WriteFinalResidual("Smoother", i, norm_r / norm0);
+    END_SU2_OMP_MASTER
   }
 
   residual = norm_r / norm0;