Hybrid Parallel AD (Part 2/?)

Common/include/CConfig.hpp

@@ -5235,7 +5235,7 @@ class CConfig {
    * \param[in] ext - the extension to be added.
    * \return The new filename
    */
-  string GetFilename(string filename, string ext, unsigned long Iter) const;
+  string GetFilename(string filename, string ext, int Iter) const;
 
   /*!
    * \brief Append the zone index to the restart or the solution files.
@@ -5444,7 +5444,8 @@ class CConfig {
    * \param[in] val_val - Value of the design variable that we want to read.
    * \return Design variable step.
    */
-  su2double GetDV_Value(unsigned short val_dv, unsigned short val_val = 0) const { return DV_Value[val_dv][val_val]; }
+  su2double& GetDV_Value(unsigned short val_dv, unsigned short val_val = 0) { return DV_Value[val_dv][val_val]; }
+  const su2double& GetDV_Value(unsigned short val_dv, unsigned short val_val = 0) const { return DV_Value[val_dv][val_val]; }
 
   /*!
    * \brief Set the value of the design variable step, we use this value in design problems.

Common/include/basic_types/ad_structure.hpp

@@ -271,14 +271,6 @@ namespace AD{
 
   extern ExtFuncHelper* FuncHelper;
 
-  /*--- Stores the indices of the input variables (they might be overwritten) ---*/
-
-  extern std::vector<su2double::GradientData> inputValues;
-
-  /*--- Current position inside the adjoint vector ---*/
-
-  extern int adjointVectorPosition;
-
   extern bool Status;
 
   extern bool PreaccActive;
@@ -297,28 +289,17 @@ namespace AD{
 
   extern std::vector<TapePosition> TapePositions;
 
-  extern std::vector<su2double::GradientData> localInputValues;
-
-  extern std::vector<su2double*> localOutputValues;
-
   extern codi::PreaccumulationHelper<su2double> PreaccHelper;
 
   /*--- Reference to the tape. ---*/
 
-  FORCEINLINE su2double::TapeType& getGlobalTape() {
-    return su2double::getGlobalTape();
-  }
+  FORCEINLINE su2double::TapeType& getGlobalTape() {return su2double::getGlobalTape();}
 
-  FORCEINLINE void RegisterInput(su2double &data, bool push_index = true) {
-    AD::getGlobalTape().registerInput(data);
-    if (push_index) {
-      inputValues.push_back(data.getGradientData());
-    }
-  }
+  FORCEINLINE void RegisterInput(su2double &data) {AD::getGlobalTape().registerInput(data);}
 
   FORCEINLINE void RegisterOutput(su2double& data) {AD::getGlobalTape().registerOutput(data);}
 
-  FORCEINLINE void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();}
+  FORCEINLINE void ResetInput(su2double &data) {data = data.getValue();}
 
   FORCEINLINE void StartRecording() {AD::getGlobalTape().setActive();}
 
@@ -335,7 +316,6 @@ namespace AD{
     opdi::logic->prepareEvaluate();
   #endif
     AD::getGlobalTape().evaluate();
-    adjointVectorPosition = 0;
   }
 
   FORCEINLINE void ComputeAdjoint(unsigned short enter, unsigned short leave) {
@@ -346,19 +326,13 @@ namespace AD{
   #else
     AD::getGlobalTape().evaluate(TapePositions[enter], TapePositions[leave]);
   #endif
-    if (leave == 0)
-      adjointVectorPosition = 0;
   }
 
   FORCEINLINE void Reset() {
     AD::getGlobalTape().reset();
   #if defined(HAVE_OPDI)
     opdi::logic->reset();
   #endif
-    if (inputValues.size() != 0) {
-      adjointVectorPosition = 0;
-      inputValues.clear();
-    }
     if (TapePositions.size() != 0) {
     #if defined(HAVE_OPDI)
       for (TapePosition& pos : TapePositions) {

Common/include/basic_types/datatype_structure.hpp

@@ -95,19 +95,9 @@ namespace SU2_TYPE {
 
   FORCEINLINE void SetDerivative(su2double& data, const passivedouble &val) {data.setGradient(val);}
 
-#ifdef CODI_REVERSE_TYPE
-  FORCEINLINE passivedouble GetSecondary(const su2double& data) {
-    return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]);
-  }
-
-  FORCEINLINE passivedouble GetDerivative(const su2double& data) {
-    return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]);
-  }
-#else // forward
   FORCEINLINE passivedouble GetSecondary(const su2double& data) {return data.getGradient();}
 
   FORCEINLINE passivedouble GetDerivative(const su2double& data) {return data.getGradient();}
-#endif
 
 #else // passive type, no AD
 

Common/include/code_config.hpp

@@ -123,9 +123,13 @@ using su2mixedfloat = passivedouble;
 
 /*--- Detect if OpDiLib has to be used. ---*/
 #if defined(HAVE_OMP) && defined(CODI_REVERSE_TYPE)
+#ifndef __INTEL_COMPILER
 #define HAVE_OPDI
+#else
+#warning Hybrid parallel reverse mode AD cannot be used with Intel compilers.
 #endif
 
 #if (_OPENMP >= 201811 && !defined(FORCE_OPDI_MACRO_BACKEND)) || defined(FORCE_OPDI_OMPT_BACKEND)
 #define HAVE_OMPT
 #endif
+#endif

Common/include/geometry/CGeometry.hpp

@@ -1234,14 +1234,14 @@ class CGeometry {
    * \brief Ray Intersects Triangle (Moller and Trumbore algorithm)
    */
   bool RayIntersectsTriangle(const su2double orig[3], const su2double dir[3],
-  const su2double vert0[3], const su2double vert1[3], const su2double vert2[3],
-  su2double *intersect);
+                             const su2double vert0[3], const su2double vert1[3], const su2double vert2[3],
+                             su2double *intersect);
 
   /*!
    * \brief Segment Intersects Triangle
    */
   bool SegmentIntersectsTriangle(su2double point0[3], const su2double point1[3],
-  su2double vert0[3], su2double vert1[3], su2double vert2[3]);
+                                 su2double vert0[3], su2double vert1[3], su2double vert2[3]);
 
   /*!
    * \brief Segment Intersects Line (for 2D FFD Intersection)
@@ -1250,16 +1250,15 @@ class CGeometry {
 
   /*!
    * \brief Register the coordinates of the mesh nodes.
-   * \param[in] config
    */
-  void RegisterCoordinates(const CConfig *config) const;
+  void RegisterCoordinates() const;
 
   /*!
    * \brief Update the multi-grid structure and the wall-distance.
    * \param geometry_container - Geometrical definition.
    * \param config - Config
    */
-  void UpdateGeometry(CGeometry **geometry_container, CConfig *config);
+  static void UpdateGeometry(CGeometry **geometry_container, CConfig *config);
 
   /*!
    * \brief Update the multi-grid structure for the customized boundary conditions

Common/include/geometry/dual_grid/CPoint.hpp

@@ -563,7 +563,7 @@ class CPoint {
   void SetVolume_n();
 
   /*!
-   * \brief Set the volume of the control volume at time n+1.
+   * \brief Set the volume of the control volume at time n-1.
    */
   void SetVolume_nM1();
 
@@ -831,52 +831,23 @@ class CPoint {
   }
 
   /*!
-   * \brief Set the adjoint values of the coordinates.
+   * \brief Register coordinates of a point.
    * \param[in] iPoint - Index of the point.
-   * \param[in] adj_sol - The adjoint values of the coordinates.
+   * \param[in] input - Register as input or output.
    */
-  inline void SetAdjointCoord(unsigned long iPoint, const su2double *adj_coor) {
-    for (unsigned long iDim = 0; iDim < nDim; iDim++)
-      SU2_TYPE::SetDerivative(Coord(iPoint,iDim), SU2_TYPE::GetValue(adj_coor[iDim]));
-  }
-
-  /*!
-   * \brief Set the adjoint values of the coordinates.
-   * \param[in] iPoint - Index of the point.
-   * \param[in] adj_sol - The adjoint values of the coordinates.
-   */
-  inline void SetAdjointCoord_LocalIndex(unsigned long iPoint, const su2double *adj_coor) {
-    for (unsigned long iDim = 0; iDim < nDim; iDim++)
-      AD::SetDerivative(AD_OutputIndex(iPoint,iDim), SU2_TYPE::GetValue(adj_coor[iDim]));
-  }
-
-  /*!
-   * \brief Get the adjoint values of the coordinates.
-   * \param[in] iPoint - Index of the point.
-   * \param[in] adj_sol - The adjoint values of the coordinates.
-   */
-  inline void GetAdjointCoord(unsigned long iPoint, su2double *adj_coor) const {
-    for (unsigned long iDim = 0; iDim < nDim; iDim++)
-      adj_coor[iDim] = SU2_TYPE::GetDerivative(Coord(iPoint,iDim));
-  }
-
-  /*!
-   * \brief Get the adjoint values of the coordinates.
-   * \param[in] iPoint - Index of the point.
-   * \param[in] adj_sol - The adjoint values of the coordinates.
-   */
-  inline void GetAdjointCoord_LocalIndex(unsigned long iPoint, su2double *adj_coor) const {
-    for (unsigned long iDim = 0; iDim < nDim; iDim++)
-      adj_coor[iDim] = AD::GetDerivative(AD_InputIndex(iPoint,iDim));
+  inline void RegisterCoordinates(unsigned long iPoint, bool input) {
+    for (unsigned long iDim = 0; iDim < nDim; iDim++) {
+      if(input) {
+        AD::RegisterInput(Coord(iPoint,iDim));
+        AD::SetIndex(AD_InputIndex(iPoint,iDim), Coord(iPoint,iDim));
+      }
+      else {
+        AD::RegisterOutput(Coord(iPoint,iDim));
+        AD::SetIndex(AD_OutputIndex(iPoint,iDim), Coord(iPoint,iDim));
+      }
+    }
   }
 
-  /*!
-   * \brief Set the adjoint vector indices of Coord vector.
-   * \param[in] iPoint - Index of the point.
-   * \param[in] input - Save them to the input or output indices vector.
-   */
-  void SetIndex(unsigned long iPoint, bool input);
-
   /*!
    * \brief Set wall roughnesses according to stored closest wall information.
    * \param[in] roughness - Mapping [rank][zone][marker] -> roughness
@@ -892,4 +863,5 @@ class CPoint {
       }
     }
   }
+
 };

Common/src/CConfig.cpp

@@ -7820,7 +7820,7 @@ CConfig::~CConfig(void) {
 
 }
 
-string CConfig::GetFilename(string filename, string ext, unsigned long Iter) const {
+string CConfig::GetFilename(string filename, string ext, int Iter) const {
 
   /*--- Remove any extension --- */
 
@@ -7840,7 +7840,7 @@ string CConfig::GetFilename(string filename, string ext, unsigned long Iter) con
     filename = GetMultiInstance_FileName(filename, GetiInst(), ext);
 
   if (GetTime_Domain()){
-    filename = GetUnsteady_FileName(filename, (int)Iter, ext);
+    filename = GetUnsteady_FileName(filename, Iter, ext);
   }
 
   return filename;

Common/src/basic_types/ad_structure.cpp

@@ -31,12 +31,6 @@ namespace AD {
 #ifdef CODI_REVERSE_TYPE
   /*--- Initialization of the global variables ---*/
 
-  int adjointVectorPosition = 0;
-
-  std::vector<su2double::GradientData> inputValues;
-  std::vector<su2double::GradientData> localInputValues;
-  std::vector<su2double*> localOutputValues;
-
   TapePosition StartPosition, EndPosition;
   std::vector<TapePosition> TapePositions;
 

Common/src/geometry/CGeometry.cpp

@@ -2476,18 +2476,12 @@ void CGeometry::ComputeAirfoil_Section(su2double *Plane_P0, su2double *Plane_Nor
 
 }
 
-void CGeometry::RegisterCoordinates(const CConfig *config) const {
+void CGeometry::RegisterCoordinates() const {
   const bool input = true;
-  const bool push_index = config->GetMultizone_Problem()? false : true;
 
   SU2_OMP_FOR_STAT(roundUpDiv(nPoint,omp_get_num_threads()))
   for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {
-    for (auto iDim = 0u; iDim < nDim; iDim++) {
-      AD::RegisterInput(nodes->GetCoord(iPoint)[iDim], push_index);
-    }
-    if(!push_index) {
-      nodes->SetIndex(iPoint, input);
-    }
+    nodes->RegisterCoordinates(iPoint, input);
   }
   END_SU2_OMP_FOR
 }
@@ -2496,10 +2490,6 @@ void CGeometry::UpdateGeometry(CGeometry **geometry_container, CConfig *config)
 
   geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, COORDINATES);
   geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, COORDINATES);
-  if (config->GetDynamic_Grid()){
-    geometry_container[MESH_0]->InitiateComms(geometry_container[MESH_0], config, GRID_VELOCITY);
-    geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, GRID_VELOCITY);
-  }
 
   geometry_container[MESH_0]->SetControlVolume(config, UPDATE);
   geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE);
@@ -2732,8 +2722,6 @@ void CGeometry::ComputeSurf_Curvature(CConfig *config) {
   su2double U[3] = {0.0}, V[3] = {0.0}, W[3] = {0.0}, Length_U, Length_V, Length_W, CosValue, Angle_Value, *K, *Angle_Defect, *Area_Vertex, *Angle_Alpha, *Angle_Beta, **NormalMeanK, MeanK, GaussK, MaxPrinK, cot_alpha, cot_beta, delta, X1, X2, X3, Y1, Y2, Y3, radius, *Buffer_Send_Coord, *Buffer_Receive_Coord, *Coord, Dist, MinDist, MaxK, MinK, SigmaK;
   bool *Check_Edge;
 
-  const bool fea = config->GetStructuralProblem();
-
   /*--- Allocate surface curvature ---*/
   K = new su2double [nPoint];
   for (iPoint = 0; iPoint < nPoint; iPoint++) K[iPoint] = 0.0;
@@ -2997,7 +2985,7 @@ void CGeometry::ComputeSurf_Curvature(CConfig *config) {
 
   SigmaK = sqrt(SigmaK/su2double(TotalnPointDomain));
 
-  if ((rank == MASTER_NODE) && (!fea))
+  if (rank == MASTER_NODE)
     cout << "Max K: " << MaxK << ". Mean K: " << MeanK << ". Standard deviation K: " << SigmaK << "." << endl;
 
   Point_Critical.clear();

Common/src/geometry/dual_grid/CPoint.cpp

@@ -73,7 +73,7 @@ void CPoint::FullAllocation(unsigned short imesh, const CConfig *config) {
     }
   }
 
-  if(config->GetAD_Mode() && config->GetMultizone_Problem()) {
+  if (config->GetDiscrete_Adjoint()) {
     AD_InputIndex.resize(npoint,nDim) = 0;
     AD_OutputIndex.resize(npoint,nDim) = 0;
   }
@@ -153,17 +153,6 @@ void CPoint::SetPoints(const vector<vector<unsigned long> >& pointsMatrix) {
   Edge = CCompressedSparsePatternL(Point.outerPtr(), Point.outerPtr()+Point.getOuterSize()+1, long(-1));
 }
 
-void CPoint::SetIndex(unsigned long iPoint, bool input) {
-  for (unsigned long iDim = 0; iDim < nDim; iDim++) {
-    if(input) {
-      AD::SetIndex(AD_InputIndex(iPoint,iDim), Coord(iPoint,iDim));
-    }
-    else {
-      AD::SetIndex(AD_OutputIndex(iPoint,iDim), Coord(iPoint,iDim));
-    }
-  }
-}
-
 void CPoint::SetVolume_n() {
   assert(Volume_n.size() == Volume.size());
   parallelCopy(Volume.size(), Volume.data(), Volume_n.data());