Further Explicit Adjoints Locking and Lock-Free Adjoints Access

Common/include/basic_types/ad_structure.hpp

@@ -396,6 +396,8 @@ FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getI
 
 // 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().
+// This method does not perform locking either.
+// It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints().
 FORCEINLINE void SetDerivative(int index, const double val) {
   if (index == 0)  // Allow multiple threads to "set the derivative" of passive variables without causing data races.
     return;
@@ -406,6 +408,8 @@ FORCEINLINE void SetDerivative(int index, const double val) {
 // 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().
+// This method does not perform locking either.
+// It should be safeguarded by calls to AD::BeginUseAdjoints() and AD::EndUseAdjoints().
 FORCEINLINE double GetDerivative(int index) {
   return AD::getTape().getGradient(index, codi::AdjointsManagement::Manual);
 }

SU2_CFD/include/variables/CVariable.hpp

@@ -2170,13 +2170,19 @@ class CVariable {
   }
 
   inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const {
-    for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++)
-      adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n(iPoint,iVar));
+    int index = 0;
+    for (unsigned long iVar = 0; iVar < Solution_time_n.cols(); iVar++) {
+      AD::SetIndex(index, Solution_time_n(iPoint, iVar));
+      adj_sol[iVar] = AD::GetDerivative(index);
+    }
   }
 
   inline void GetAdjointSolution_time_n1(unsigned long iPoint, su2double *adj_sol) const {
-    for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); iVar++)
-      adj_sol[iVar] = SU2_TYPE::GetDerivative(Solution_time_n1(iPoint,iVar));
+    int index = 0;
+    for (unsigned long iVar = 0; iVar < Solution_time_n1.cols(); iVar++) {
+      AD::SetIndex(index, Solution_time_n1(iPoint, iVar));
+      adj_sol[iVar] = AD::GetDerivative(index);
+    }
   }
 
   /*!

SU2_CFD/src/drivers/CDriver.cpp

@@ -248,6 +248,8 @@ CDriverBase(confFile, val_nZone, MPICommunicator), StopCalc(false), fsi(false),
       cout << endl <<"---------------------- Turbomachinery Preprocessing ---------------------" << endl;
 
     PreprocessTurbomachinery(config_container, geometry_container, solver_container, interface_container);
+  } else {
+    mixingplane = false;
   }
 
 

SU2_CFD/src/solvers/CDiscAdjFEASolver.cpp

@@ -232,22 +232,30 @@ void CDiscAdjFEASolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *co
 
   /*--- Extract and store the adjoint solution ---*/
 
+  AD::BeginUseAdjoints();
+
   for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
     su2double Solution[MAXNVAR] = {0.0};
     direct_solver->GetNodes()->GetAdjointSolution(iPoint,Solution);
     nodes->SetSolution(iPoint,Solution);
   }
 
+  AD::EndUseAdjoints();
+
   if (CrossTerm) return;
 
   /*--- Extract and store the adjoint solution at time n (including accel. and velocity) ---*/
 
   if (config->GetTime_Domain()) {
+    AD::BeginUseAdjoints();
+
     for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
       su2double Solution[MAXNVAR] = {0.0};
       direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution);
       nodes->Set_Solution_time_n(iPoint,Solution);
     }
+
+    AD::EndUseAdjoints();
   }
 
   /*--- Set the residuals ---*/
@@ -358,6 +366,8 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo
 
   /*--- Extract the geometric sensitivities ---*/
 
+  AD::BeginUseAdjoints();
+
   for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) {
 
     auto Coord = geometry->nodes->GetCoord(iPoint);
@@ -375,6 +385,8 @@ void CDiscAdjFEASolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSo
     }
   }
 
+  AD::EndUseAdjoints();
+
 }
 
 void CDiscAdjFEASolver::ReadDV(const CConfig *config) {

SU2_CFD/src/solvers/CDiscAdjSolver.cpp

@@ -317,6 +317,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
 
   if (!config->GetMultizone_Problem()) nodes->Set_OldSolution();
 
+  AD::BeginUseAdjoints();
+
   SU2_OMP_FOR_STAT(omp_chunk_size)
   for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
 
@@ -339,6 +341,8 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
   }
   END_SU2_OMP_FOR
 
+  AD::EndUseAdjoints();
+
   direct_solver->ExtractAdjoint_SolutionExtra(nodes->GetSolutionExtra(), config);
 
   /*--- Residuals and time_n terms are not needed when evaluating multizone cross terms. ---*/
@@ -355,24 +359,32 @@ void CDiscAdjSolver::ExtractAdjoint_Solution(CGeometry *geometry, CConfig *confi
 
   /*--- Extract and store the adjoint of the primal solution at time n ---*/
   if (time_n_needed) {
+    AD::BeginUseAdjoints();
+
     SU2_OMP_FOR_STAT(omp_chunk_size)
     for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
       su2double Solution[MAXNVAR] = {0.0};
       direct_solver->GetNodes()->GetAdjointSolution_time_n(iPoint,Solution);
       nodes->Set_Solution_time_n(iPoint,Solution);
     }
     END_SU2_OMP_FOR
+
+    AD::EndUseAdjoints();
   }
 
   /*--- Extract and store the adjoint of the primal solution at time n-1 ---*/
   if (time_n1_needed) {
+    AD::BeginUseAdjoints();
+
     SU2_OMP_FOR_STAT(omp_chunk_size)
     for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
       su2double Solution[MAXNVAR] = {0.0};
       direct_solver->GetNodes()->GetAdjointSolution_time_n1(iPoint,Solution);
       nodes->Set_Solution_time_n1(iPoint,Solution);
     }
     END_SU2_OMP_FOR
+
+    AD::EndUseAdjoints();
   }
 
 }
@@ -477,6 +489,8 @@ void CDiscAdjSolver::SetAdjoint_Output(CGeometry *geometry, CConfig *config) {
 
 void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolver*) {
 
+  AD::BeginUseAdjoints();
+
   SU2_OMP_PARALLEL {
 
   const bool time_stepping = (config->GetTime_Marching() != TIME_MARCHING::STEADY);
@@ -510,6 +524,8 @@ void CDiscAdjSolver::SetSensitivity(CGeometry *geometry, CConfig *config, CSolve
 
   }
   END_SU2_OMP_PARALLEL
+
+  AD::EndUseAdjoints();
 }
 
 void CDiscAdjSolver::SetSurface_Sensitivity(CGeometry *geometry, CConfig *config) {

TestCases/hybrid_regression_AD.py

@@ -242,7 +242,7 @@ def main():
     pywrapper_FEA_AD_FlowLoad.test_vals_aarch64 = [-0.131745, -0.553214, -0.000364, -0.003101]
     pywrapper_FEA_AD_FlowLoad.command       = TestCase.Command(exec = "python", param = "run_adjoint.py --parallel -f")
     pywrapper_FEA_AD_FlowLoad.timeout       = 1600
-    pywrapper_FEA_AD_FlowLoad.tol           = 1e-4
+    pywrapper_FEA_AD_FlowLoad.tol           = 1e-3
     pywrapper_FEA_AD_FlowLoad.new_output    = False
     pywrapper_FEA_AD_FlowLoad.enabled_with_tsan = False
     test_list.append(pywrapper_FEA_AD_FlowLoad)
@@ -257,7 +257,7 @@ def main():
     pywrapper_CFD_AD_MeshDisp.test_vals_aarch64 = [30.000000, -2.516536, 1.386443, 0.000000]
     pywrapper_CFD_AD_MeshDisp.command       = TestCase.Command(exec = "python", param = "run_adjoint.py --parallel -f")
     pywrapper_CFD_AD_MeshDisp.timeout       = 1600
-    pywrapper_CFD_AD_MeshDisp.tol           = 1e-4
+    pywrapper_CFD_AD_MeshDisp.tol           = 1e-3
     pywrapper_CFD_AD_MeshDisp.new_output    = False
     pywrapper_CFD_AD_MeshDisp.enabled_with_tsan = False
     test_list.append(pywrapper_CFD_AD_MeshDisp)

TestCases/serial_regression.py

@@ -1066,7 +1066,7 @@ def main():
     airfoilRBF.cfg_dir   = "fea_fsi/Airfoil_RBF"
     airfoilRBF.cfg_file  = "config.cfg"
     airfoilRBF.test_iter = 1
-    airfoilRBF.test_vals = [1.000000, -2.786183, -4.977959]
+    airfoilRBF.test_vals = [1.000000, -2.786186, -4.977944]
     airfoilRBF.multizone = True
     test_list.append(airfoilRBF)