[WIP] Tag / debug tape recordings

Common/include/basic_types/ad_structure.hpp

@@ -38,6 +38,9 @@
  */
 namespace AD {
 #ifndef CODI_REVERSE_TYPE
+
+using Identifier = int;
+
 /*!
  * \brief Start the recording of the operations and involved variables.
  * If called, the computational graph of all operations occuring after the call will be stored,
@@ -101,14 +104,20 @@ inline void EndUseAdjoints() {}
  * \param[in] index - Position in the adjoint vector.
  * \param[in] val - adjoint value to be set.
  */
-inline void SetDerivative(int index, const double val) {}
+inline void SetDerivative(Identifier index, const double val) {}
 
 /*!
  * \brief Extracts the adjoint value at index
  * \param[in] index - position in the adjoint vector where the derivative will be extracted.
  * \return Derivative value.
  */
-inline double GetDerivative(int index) { return 0.0; }
+inline double GetDerivative(Identifier index) { return 0.0; }
+
+/*!
+ * \brief Returns the identifier that represents an inactive variable.
+ * \return Passive index.
+ */
+inline Identifier GetPassiveIndex() { return 0; }
 
 /*!
  * \brief Clears the currently stored adjoints but keeps the computational graph.
@@ -259,7 +268,7 @@ inline void SetExtFuncOut(T&& data, const int size_x, const int size_y) {}
  * \param[in] data - variable whose gradient information will be extracted.
  * \param[in] index - where obtained gradient information will be stored.
  */
-inline void SetIndex(int& index, const su2double& data) {}
+inline void SetIndex(Identifier& index, const su2double& data) {}
 
 /*!
  * \brief Pushes back the current tape position to the tape position's vector.
@@ -304,6 +313,7 @@ inline void EndNoSharedReading() {}
 using CheckpointHandler = codi::ExternalFunctionUserData;
 
 using Tape = su2double::Tape;
+using Identifier = su2double::Identifier;
 
 #ifdef HAVE_OPDI
 using ExtFuncHelper = codi::OpenMPExternalFunctionHelper<su2double>;
@@ -470,14 +480,14 @@ FORCEINLINE void BeginUseAdjoints() { AD::getTape().beginUseAdjointVector(); }
 
 FORCEINLINE void EndUseAdjoints() { AD::getTape().endUseAdjointVector(); }
 
-FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); }
+FORCEINLINE void SetIndex(Identifier& index, const su2double& data) { index = data.getIdentifier(); }
 
 // 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.
+FORCEINLINE void SetDerivative(Identifier index, const double val) {
+  if (!AD::getTape().isIdentifierActive(index))  // Allow multiple threads to "set the derivative" of passive variables without causing data races.
     return;
 
   AD::getTape().setGradient(index, val, codi::AdjointsManagement::Manual);
@@ -488,10 +498,12 @@ FORCEINLINE void SetDerivative(int index, const double val) {
 // 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) {
+FORCEINLINE double GetDerivative(Identifier index) {
   return AD::getTape().getGradient(index, codi::AdjointsManagement::Manual);
 }
 
+FORCEINLINE Identifier GetPassiveIndex() { return AD::getTape().getPassiveIndex(); }
+
 FORCEINLINE bool IsIdentifierActive(su2double const& value) {
   return getTape().isIdentifierActive(value.getIdentifier());
 }

Common/include/geometry/dual_grid/CPoint.hpp

@@ -113,9 +113,8 @@ class CPoint {
   su2activevector MaxLength;          /*!< \brief The maximum cell-center to cell-center length. */
   su2activevector RoughnessHeight;    /*!< \brief Roughness of the nearest wall. */
 
-  su2matrix<int> AD_InputIndex; /*!< \brief Indices of Coord variables in the adjoint vector. */
-  su2matrix<int>
-      AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */
+  su2matrix<AD::Identifier> AD_InputIndex;  /*!< \brief Indices of Coord variables in the adjoint vector. */
+  su2matrix<AD::Identifier> AD_OutputIndex; /*!< \brief Indices of Coord variables in the adjoint vector after having been updated. */
 
   /*!
    * \brief Allocate fields required by the minimal constructor.

Common/src/geometry/dual_grid/CPoint.cpp

@@ -67,8 +67,8 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig* config) {
   }
 
   if (config->GetDiscrete_Adjoint()) {
-    AD_InputIndex.resize(npoint, nDim) = 0;
-    AD_OutputIndex.resize(npoint, nDim) = 0;
+    AD_InputIndex.resize(npoint, nDim) = AD::GetPassiveIndex();
+    AD_OutputIndex.resize(npoint, nDim) = AD::GetPassiveIndex();
   }
 
   /*--- Multigrid structures. ---*/

SU2_CFD/include/drivers/CDiscAdjMultizoneDriver.hpp

@@ -96,7 +96,7 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
 
   bool eval_transfer = false;     /*!< \brief Evaluate the transfer section of the tape. */
   su2double ObjFunc;              /*!< \brief Value of the objective function. */
-  int ObjFunc_Index;              /*!< \brief Index of the value of the objective function. */
+  AD::Identifier ObjFunc_Index;   /*!< \brief Index of the value of the objective function. */
 
   CIteration*** direct_iteration;       /*!< \brief Array of pointers to the direct iterations. */
   COutput** direct_output;              /*!< \brief Array of pointers to the direct outputs. */

SU2_CFD/include/variables/CVariable.hpp

@@ -92,8 +92,8 @@ class CVariable {
 
   MatrixType Solution_BGS_k;     /*!< \brief Old solution container for BGS iterations. */
 
-  su2matrix<int> AD_InputIndex;    /*!< \brief Indices of Solution variables in the adjoint vector. */
-  su2matrix<int> AD_OutputIndex;   /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */
+  su2matrix<AD::Identifier> AD_InputIndex;    /*!< \brief Indices of Solution variables in the adjoint vector. */
+  su2matrix<AD::Identifier> AD_OutputIndex;   /*!< \brief Indices of Solution variables in the adjoint vector after having been updated. */
 
   VectorType SolutionExtra; /*!< \brief Stores adjoint solution for extra solution variables.
                                         Currently only streamwise periodic pressure-drop for massflow prescribed flows. */
@@ -118,7 +118,7 @@ class CVariable {
     assert(false && "A base method of CVariable was used, but it should have been overridden by the derived class.");
   }
 
-  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<int>* ad_index = nullptr) {
+  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<AD::Identifier>* ad_index = nullptr) {
     const auto nPoint = variable.rows();
     SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads()))
     for (unsigned long iPoint = 0; iPoint < nPoint; ++iPoint) {
@@ -133,7 +133,7 @@ class CVariable {
     END_SU2_OMP_FOR
   }
 
-  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<int>& ad_index) {
+  void RegisterContainer(bool input, su2activematrix& variable, su2matrix<AD::Identifier>& ad_index) {
     RegisterContainer(input, variable, &ad_index);
   }
 
@@ -2180,15 +2180,15 @@ class CVariable {
   }
 
   inline void GetAdjointSolution_time_n(unsigned long iPoint, su2double *adj_sol) const {
-    int index = 0;
+    AD::Identifier index = AD::GetPassiveIndex();
     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 {
-    int index = 0;
+    AD::Identifier index = AD::GetPassiveIndex();
     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/CDiscAdjMultizoneDriver.cpp

@@ -744,7 +744,7 @@
     if (kind_recording == RECORDING::SOLUTION_VARIABLES) {
       cout << " Objective function                   : " << ObjFunc;
       if (driver_config->GetWrt_AD_Statistics()){
-        cout << " (" << ObjFunc_Index << ")\n";
+        // cout << " (" << ObjFunc_Index << ")\n";
       }
       cout << endl;
     }

SU2_CFD/src/solvers/CSolver.cpp

@@ -4071,7 +4071,7 @@ void CSolver::SetVertexTractionsAdjoint(CGeometry *geometry, const CConfig *conf
 
   unsigned short iMarker, iDim;
   unsigned long iVertex, iPoint;
-  int index;
+  AD::Identifier index;
 
   /*--- Loop over all the markers ---*/
   for (iMarker = 0; iMarker < config->GetnMarker_All(); iMarker++) {

SU2_CFD/src/variables/CVariable.cpp

@@ -69,8 +69,8 @@ CVariable::CVariable(unsigned long npoint, unsigned long ndim, unsigned long nva
       External.resize(nPoint,nVar) = su2double(0.0);
 
     if (!adjoint) {
-      AD_InputIndex.resize(nPoint,nVar) = -1;
-      AD_OutputIndex.resize(nPoint,nVar) = -1;
+      AD_InputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex();
+      AD_OutputIndex.resize(nPoint,nVar) = AD::GetPassiveIndex();
     }
   }
 

meson.build

@@ -95,6 +95,8 @@ if get_option('enable-autodiff') and not omp
     codi_rev_args += '-DCODI_PRIMAL_REUSE_TAPE'
   elif get_option('codi-tape') == 'PrimalMultiUse'
     codi_rev_args += '-DCODI_PRIMAL_MULTIUSE_TAPE'
+  elif get_option('codi-tape') == 'Tag'
+    codi_rev_args += '-DCODI_TAG_TAPE'
   else
     error('Invalid CoDiPack tape choice @0@'.format(get_option('codi-tape')))
   endif

meson_options.txt

@@ -23,7 +23,7 @@ option('enable-coolprop',  type : 'boolean', value : false, description: 'enable
 option('enable-mlpcpp', type : 'boolean', value : false, description: 'enable profiling through gprof')
 option('enable-gprof', type : 'boolean', value : false, description: 'enable MLPCpp support')
 option('opdi-backend', type : 'combo', choices : ['auto', 'macro', 'ompt'], value : 'auto', description: 'OpDiLib backend choice')
-option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse'], value : 'JacobianLinear', description: 'CoDiPack tape choice')
+option('codi-tape', type : 'combo', choices : ['JacobianLinear', 'JacobianReuse', 'JacobianMultiUse', 'PrimalLinear', 'PrimalReuse', 'PrimalMultiUse', 'Tag'], value : 'JacobianLinear', description: 'CoDiPack tape choice')
 option('opdi-shared-read-opt', type : 'boolean', value : true, description : 'OpDiLib shared reading optimization')
 option('librom_root', type : 'string', value : '', description: 'libROM base directory')
 option('enable-librom', type : 'boolean', value : false, description: 'enable LLNL libROM support')