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.

Common/include/code_config.hpp

@@ -116,9 +116,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

@@ -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

@@ -543,7 +543,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();
 
@@ -811,50 +811,21 @@ 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]));
+  inline void RegisterCoordinates(unsigned long iPoint, bool input) {
+    for (unsigned long iDim = 0; iDim < nDim; iDim++) {
+      if(input) {
+        AD::RegisterInput(Coord(iPoint,iDim),false);
+        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 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));
-  }
-
-  /*!
-   * \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);
-
 };

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/geometry/CGeometry.cpp

@@ -2475,18 +2475,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
 }
@@ -2495,10 +2489,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);

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;
   }
@@ -147,17 +147,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());

SU2_CFD/include/solvers/CDiscAdjFEASolver.hpp

@@ -49,7 +49,6 @@ class CDiscAdjFEASolver final : public CSolver {
     unsigned short size = 0;
     su2double* val = nullptr;         /*!< \brief Value of the variable. */
     int* AD_Idx = nullptr;            /*!< \brief Derivative index in the AD tape. */
-    bool localIdx = false;
     su2double* LocalSens = nullptr;   /*!< \brief Local sensitivity (domain). */
     su2double* GlobalSens = nullptr;  /*!< \brief Global sensitivity (mpi). */
     su2double* TotalSens = nullptr;   /*!< \brief Total sensitivity (time domain). */
@@ -69,28 +68,22 @@ class CDiscAdjFEASolver final : public CSolver {
 
     void clear() {
       size = 0;
-      localIdx = false;
       delete [] val;
       delete [] AD_Idx;
       delete [] LocalSens;
       delete [] GlobalSens;
       delete [] TotalSens;
     }
 
-    void Register(bool push_index) {
-      for (auto i = 0u; i < size; ++i) AD::RegisterInput(val[i], push_index);
-    }
-
-    void SetIndex() {
-      for (auto i = 0u; i < size; ++i) AD::SetIndex(AD_Idx[i], val[i]);
-      localIdx = true;
+    void Register() {
+      for (auto i = 0u; i < size; ++i) {
+        AD::RegisterInput(val[i], false);
+        AD::SetIndex(AD_Idx[i], val[i]);
+      }
     }
 
     void GetDerivative() {
-      if (localIdx)
-        for (auto i = 0u; i < size; ++i) LocalSens[i] = AD::GetDerivative(AD_Idx[i]);
-      else
-        for (auto i = 0u; i < size; ++i) LocalSens[i] = SU2_TYPE::GetDerivative(val[i]);
+      for (auto i = 0u; i < size; ++i) LocalSens[i] = AD::GetDerivative(AD_Idx[i]);
 
       SU2_MPI::Allreduce(LocalSens, GlobalSens, size, MPI_DOUBLE, MPI_SUM, SU2_MPI::GetComm());
     }

SU2_CFD/include/solvers/CFVMFlowSolverBase.inl

@@ -707,7 +707,7 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
 
   /*--- Restart the solution from file information ---*/
 
-  unsigned short iDim, iVar, iMesh, iMeshFine;
+  unsigned short iDim, iVar, iMesh;
   unsigned long iPoint, index, iChildren, Point_Fine;
   unsigned short turb_model = config->GetKind_Turb_Model();
   su2double Area_Children, Area_Parent;
@@ -746,6 +746,15 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
     Read_SU2_Restart_ASCII(geometry[MESH_0], config, restart_filename);
   }
 
+  if (update_geo && dynamic_grid) {
+    auto notFound = fields.end();
+    if (find(fields.begin(), notFound, string("\"Grid_Velocity_x\"")) == notFound) {
+      if (rank == MASTER_NODE)
+        cout << "\nWARNING: The restart file does not contain grid velocities, these will be set to zero.\n" << endl;
+      steady_restart = true;
+    }
+  }
+
   /*--- Load data from the restart into correct containers. ---*/
 
   unsigned long counter = 0, iPoint_Global = 0;
@@ -829,39 +838,22 @@ void CFVMFlowSolverBase<V, R>::LoadRestart_impl(CGeometry **geometry, CSolver **
   END_SU2_OMP_MASTER
   SU2_OMP_BARRIER
 
-  /*--- Update the geometry for flows on deforming meshes ---*/
+  /*--- Update the geometry for flows on deforming meshes. ---*/
 
   if ((dynamic_grid || static_fsi) && update_geo) {
 
-    /*--- Communicate the new coordinates and grid velocities at the halos ---*/
-
-    geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, COORDINATES);
-    geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, COORDINATES);
+    CGeometry::UpdateGeometry(geometry, config);
 
     if (dynamic_grid) {
-      geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY);
-      geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY);
-    }
+      for (iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
 
-    /*--- Recompute the edges and dual mesh control volumes in the
-     domain and on the boundaries. ---*/
-
-    geometry[MESH_0]->SetControlVolume(config, UPDATE);
-    geometry[MESH_0]->SetBoundControlVolume(config, UPDATE);
-    geometry[MESH_0]->SetMaxLength(config);
-
-    /*--- Update the multigrid structure after setting up the finest grid,
-     including computing the grid velocities on the coarser levels. ---*/
-
-    for (iMesh = 1; iMesh <= config->GetnMGLevels(); iMesh++) {
-      iMeshFine = iMesh-1;
-      geometry[iMesh]->SetControlVolume(config, geometry[iMeshFine], UPDATE);
-      geometry[iMesh]->SetBoundControlVolume(config, geometry[iMeshFine],UPDATE);
-      geometry[iMesh]->SetCoord(geometry[iMeshFine]);
-      if (dynamic_grid) {
-        geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMeshFine], config);
+        /*--- Compute the grid velocities on the coarser levels. ---*/
+        if (iMesh) geometry[iMesh]->SetRestricted_GridVelocity(geometry[iMesh-1], config);
+        else {
+          geometry[MESH_0]->InitiateComms(geometry[MESH_0], config, GRID_VELOCITY);
+          geometry[MESH_0]->CompleteComms(geometry[MESH_0], config, GRID_VELOCITY);
+        }
       }
-      geometry[iMesh]->SetMaxLength(config);
     }
   }
 
@@ -973,7 +965,7 @@ void CFVMFlowSolverBase<V, R>::SetInitialCondition(CGeometry **geometry, CSolver
 
   /*--- The value of the solution for the first iteration of the dual time ---*/
 
-  if (dual_time && (TimeIter == 0 || (restart && TimeIter == config->GetRestart_Iter()))) {
+  if (dual_time && TimeIter == config->GetRestart_Iter()) {
     PushSolutionBackInTime(TimeIter, restart, rans, solver_container, geometry, config);
   }
 
@@ -996,26 +988,36 @@ void CFVMFlowSolverBase<V, R>::PushSolutionBackInTime(unsigned long TimeIter, bo
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
       solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n1();
     }
+
+    if (dynamic_grid) {
+      geometry[iMesh]->nodes->SetVolume_n();
+      geometry[iMesh]->nodes->SetVolume_nM1();
+    }
+
+    if (config->GetGrid_Movement()) {
+      geometry[iMesh]->nodes->SetCoord_n();
+      geometry[iMesh]->nodes->SetCoord_n1();
+    }
   }
 
-  if (restart && (TimeIter == config->GetRestart_Iter()) && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) {
-    /*--- Load an additional restart file for a 2nd-order restart ---*/
+  if (restart && (config->GetTime_Marching() == TIME_MARCHING::DT_STEPPING_2ND)) {
+
+    /*--- Load an additional restart file for a 2nd-order restart. ---*/
 
-    solver_container[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver_container, config, config->GetRestart_Iter() - 1,
-                                                    true);
+    solver_container[MESH_0][FLOW_SOL]->LoadRestart(geometry, solver_container, config, TimeIter-1, true);
 
-    /*--- Load an additional restart file for the turbulence model ---*/
+    /*--- Load an additional restart file for the turbulence model. ---*/
     if (rans)
-      solver_container[MESH_0][TURB_SOL]->LoadRestart(geometry, solver_container, config, config->GetRestart_Iter() - 1,
-                                                      false);
+      solver_container[MESH_0][TURB_SOL]->LoadRestart(geometry, solver_container, config, TimeIter-1, false);
 
     /*--- Push back this new solution to time level N. ---*/
 
     for (unsigned short iMesh = 0; iMesh <= config->GetnMGLevels(); iMesh++) {
       solver_container[iMesh][FLOW_SOL]->GetNodes()->Set_Solution_time_n();
-      if (rans) {
-        solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
-      }
+      if (rans) solver_container[iMesh][TURB_SOL]->GetNodes()->Set_Solution_time_n();
+
+      geometry[iMesh]->nodes->SetVolume_n();
+      if (config->GetGrid_Movement()) geometry[iMesh]->nodes->SetCoord_n();
     }
   }
 }