[WIP] Tag / debug tape recordings

Common/include/CConfig.hpp

@@ -1052,7 +1052,8 @@ class CConfig {
   long ParMETIS_pointWgt;           /*!< \brief Load balancing weight given to points. */
   long ParMETIS_edgeWgt;            /*!< \brief Load balancing weight given to edges. */
   unsigned short DirectDiff;        /*!< \brief Direct Differentation mode. */
-  bool DiscreteAdjoint;                /*!< \brief AD-based discrete adjoint mode. */
+  bool DiscreteAdjoint,                /*!< \brief AD-based discrete adjoint mode. */
+  DiscreteAdjointDebug;                /*!< \brief Discrete adjoint debug mode using tags. */
   su2double Const_DES;                 /*!< \brief Detached Eddy Simulation Constant. */
   WINDOW_FUNCTION Kind_WindowFct;      /*!< \brief Type of window (weight) function for objective functional. */
   unsigned short Kind_HybridRANSLES;   /*!< \brief Kind of Hybrid RANS/LES. */
@@ -8904,6 +8905,12 @@ class CConfig {
    */
   bool GetDiscrete_Adjoint(void) const { return DiscreteAdjoint; }
 
+  /*!
+   * \brief Get the indicator whether a debug run for the discrete adjoint solver will be started.
+   * \return the discrete adjoint debug indicator.
+   */
+  bool GetDiscrete_Adjoint_Debug(void) const { return DiscreteAdjointDebug; }
+
   /*!
    * \brief Get the number of subiterations while a ramp is applied.
    * \return Number of FSI subiters.

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 @@
  * \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,13 @@
  * \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 Sets the tag tape to a specific tag.
+ * \param[in] tag - the number to which the tag is set.
+ */
+inline void SetTag(int tag) {}
 
 /*!
  * \brief Pushes back the current tape position to the tape position's vector.
@@ -304,6 +319,7 @@
 using CheckpointHandler = codi::ExternalFunctionUserData;
 
 using Tape = su2double::Tape;
+using Identifier = su2double::Identifier;
 
 #ifdef HAVE_OPDI
 using ExtFuncHelper = codi::OpenMPExternalFunctionHelper<su2double>;
@@ -470,14 +486,16 @@
 
 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(); }
+
+FORCEINLINE void SetTag(int tag) { AD::getTape().setCurTag(tag); }
 
 // 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 +506,12 @@
 // 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());
 }
@@ -502,7 +522,8 @@
 template <class T, class... Ts, su2enable_if<std::is_same<T, su2double>::value> = 0>
 FORCEINLINE void SetPreaccIn(const T& data, Ts&&... moreData) {
   if (!PreaccActive) return;
-  if (IsIdentifierActive(data)) PreaccHelper.addInput(data);
+  // if (IsIdentifierActive(data))
+    PreaccHelper.addInput(data);
   SetPreaccIn(moreData...);
 }
 
@@ -515,9 +536,9 @@
 FORCEINLINE void SetPreaccIn(const T& data, const int size) {
   if (PreaccActive) {
     for (int i = 0; i < size; i++) {
-      if (IsIdentifierActive(data[i])) {
+      // if (IsIdentifierActive(data[i])) {
         PreaccHelper.addInput(data[i]);
-      }
+      // }
     }
   }
 }
@@ -527,9 +548,9 @@
   if (!PreaccActive) return;
   for (int i = 0; i < size_x; i++) {
     for (int j = 0; j < size_y; j++) {
-      if (IsIdentifierActive(data[i][j])) {
+      // if (IsIdentifierActive(data[i][j])) {
         PreaccHelper.addInput(data[i][j]);
-      }
+      // }
     }
   }
 }
@@ -547,17 +568,18 @@
 template <class T, class... Ts, su2enable_if<std::is_same<T, su2double>::value> = 0>
 FORCEINLINE void SetPreaccOut(T& data, Ts&&... moreData) {
   if (!PreaccActive) return;
-  if (IsIdentifierActive(data)) PreaccHelper.addOutput(data);
+  // if (IsIdentifierActive(data))
+    PreaccHelper.addOutput(data);
   SetPreaccOut(moreData...);
 }
 
 template <class T>
 FORCEINLINE void SetPreaccOut(T&& data, const int size) {
   if (PreaccActive) {
     for (int i = 0; i < size; i++) {
-      if (IsIdentifierActive(data[i])) {
+      // if (IsIdentifierActive(data[i])) {
         PreaccHelper.addOutput(data[i]);
-      }
+      // }
     }
   }
 }
@@ -567,9 +589,9 @@
   if (!PreaccActive) return;
   for (int i = 0; i < size_x; i++) {
     for (int j = 0; j < size_y; j++) {
-      if (IsIdentifierActive(data[i][j])) {
+      // if (IsIdentifierActive(data[i][j])) {
         PreaccHelper.addOutput(data[i][j]);
-      }
+      // }
     }
   }
 }

Common/include/code_config.hpp

@@ -110,6 +110,8 @@ using su2double = codi::RealReversePrimal;
 using su2double = codi::RealReversePrimalIndexGen<double, double, codi::ReuseIndexManager<int> >;
 #elif defined(CODI_PRIMAL_MULTIUSE_TAPE)
 using su2double = codi::RealReversePrimalIndex;
+#elif defined(CODI_TAG_TAPE)
+using su2double = codi::RealReverseTag;
 #else
 #error "Please define a CoDiPack tape."
 #endif

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/include/option_structure.hpp

@@ -2493,6 +2493,8 @@ enum class RECORDING {
   MESH_COORDS,
   MESH_DEFORM,
   SOLUTION_AND_MESH,
+  TAG_INIT_SOLUTION_VARIABLES,
+  TAG_CHECK_SOLUTION_VARIABLES
 };
 
 /*!

Common/include/parallelization/vectorization.hpp

@@ -90,11 +90,11 @@ class Array : public CVecExpr<Array<Scalar_t, N>, Scalar_t> {
  public:
   using Scalar = Scalar_t;
   enum : size_t { Size = N };
-  enum : size_t { Align = Size * sizeof(Scalar) };
+  enum : size_t { Align = Size * 32 };
   static constexpr bool StoreAsRef = true;
 
  private:
-  alignas(Size * sizeof(Scalar)) Scalar x_[N];
+   alignas(Size * 32) Scalar x_[N];
 
  public:
 #define ARRAY_BOILERPLATE                                                  \

Common/src/CConfig.cpp

@@ -1105,11 +1105,19 @@ void CConfig::SetConfig_Options() {
   addBoolOption("MULTIZONE", Multizone_Problem, NO);
   /*!\brief PHYSICAL_PROBLEM \n DESCRIPTION: Physical governing equations \n Options: see \link Solver_Map \endlink \n DEFAULT: NONE \ingroup Config*/
   addEnumOption("MULTIZONE_SOLVER", Kind_MZSolver, Multizone_Map, ENUM_MULTIZONE::MZ_BLOCK_GAUSS_SEIDEL);
-#ifdef CODI_REVERSE_TYPE
+#if defined (CODI_REVERSE_TYPE)
   const bool discAdjDefault = true;
+# if defined (CODI_TAG_TAPE)
+    const bool discAdjDebugDefault = true;
+# else
+    const bool discAdjDebugDefault = false;
+# endif
 #else
   const bool discAdjDefault = false;
+  const bool discAdjDebugDefault = false;
 #endif
+  // TODO Set the indicator through MATH_PROBLEM
+  DiscreteAdjointDebug = discAdjDebugDefault;
   /*!\brief MATH_PROBLEM  \n DESCRIPTION: Mathematical problem \n  Options: DIRECT, ADJOINT \ingroup Config*/
   addMathProblemOption("MATH_PROBLEM", ContinuousAdjoint, false, DiscreteAdjoint, discAdjDefault, Restart_Flow, discAdjDefault);
   /*!\brief KIND_TURB_MODEL \n DESCRIPTION: Specify turbulence model \n Options: see \link Turb_Model_Map \endlink \n DEFAULT: NONE \ingroup Config*/

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

@@ -69,7 +69,7 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
   };
 
   /*!
-   * \brief Kinds of recordings (three different ones).
+   * \brief Kinds of recordings.
    */
   enum class Kind_Tape {
     FULL_TAPE,                /*!< \brief Entire derivative information for a coupled adjoint
@@ -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. */
@@ -141,6 +141,11 @@ class CDiscAdjMultizoneDriver : public CMultizoneDriver {
    */
   void StartSolver() override;
 
+  /*!
+   * \brief [Overload] Launch the debug mode for the discrete adjoint multizone solver.
+   */
+  void DebugRun() override;
+
   /*!
    * \brief Preprocess the multizone iteration
    */

SU2_CFD/include/drivers/CDriver.hpp

@@ -425,6 +425,11 @@ class CDriver : public CDriverBase {
    */
   virtual void StartSolver() {}
 
+  /*!
+   * \brief Launch a debug run of the solver.
+   */
+  virtual void DebugRun() {}
+
   /*!
    * \brief Deallocation routine
    */

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -672,11 +672,11 @@ void CFVMFlowSolverBase<V, R>::ComputeVorticityAndStrainMag(const CConfig& confi
     StrainMag(iPoint) = sqrt(2.0*StrainMag(iPoint));
     AD::SetPreaccOut(StrainMag(iPoint));
 
+    AD::EndPreacc();
+
     /*--- Max is not differentiable, so we not register them for preacc. ---*/
     strainMax = max(strainMax, StrainMag(iPoint));
     omegaMax = max(omegaMax, GeometryToolbox::Norm(3, Vorticity));
-
-    AD::EndPreacc();
   }
   END_SU2_OMP_FOR
 

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/SU2_CFD.cpp

@@ -102,6 +102,7 @@ int main(int argc, char *argv[]) {
    and perform all the preprocessing. ---*/
 
   const bool disc_adj = config.GetDiscrete_Adjoint();
+  const bool disc_adj_debug = config.GetDiscrete_Adjoint_Debug();
   const bool multizone = config.GetMultizone_Problem();
   const bool harmonic_balance = (config.GetTime_Marching() == TIME_MARCHING::HARMONIC_BALANCE);