Merge pull request #1159 from su2code/reduce_disc_adj_mem

Reduce discrete adjoint memory usage ~25%

Author: pcarruscag

Common/include/geometry/CGeometry.hpp

@@ -492,17 +492,31 @@ class CGeometry {
    * \brief Get the edge index from using the nodes of the edge.
    * \param[in] first_point - First point of the edge.
    * \param[in] second_point - Second point of the edge.
+   * \param[in] error - Throw error if edge does not exist.
    * \return Index of the edge.
    */
-  long FindEdge(unsigned long first_point, unsigned long second_point) const;
+  inline long FindEdge(unsigned long first_point, unsigned long second_point, bool error = true) const {
+    for (unsigned short iNode = 0; iNode < nodes->GetnPoint(first_point); iNode++) {
+      auto iPoint = nodes->GetPoint(first_point, iNode);
+      if (iPoint == second_point) return nodes->GetEdge(first_point, iNode);
+    }
+    if (error) {
+      char buf[100];
+      SPRINTF(buf, "Can't find the edge that connects %lu and %lu.", first_point, second_point);
+      SU2_MPI::Error(buf, CURRENT_FUNCTION);
+    }
+    return -1;
+  }
 
   /*!
    * \brief Get the edge index from using the nodes of the edge.
    * \param[in] first_point - First point of the edge.
    * \param[in] second_point - Second point of the edge.
    * \return Index of the edge.
    */
-  bool CheckEdge(unsigned long first_point, unsigned long second_point) const;
+  inline bool CheckEdge(unsigned long first_point, unsigned long second_point) const {
+    return FindEdge(first_point, second_point, false) >= 0;
+  }
 
   /*!
    * \brief Get the distance between a plane (defined by three point) and a point.
@@ -684,11 +698,6 @@ class CGeometry {
    */
   inline virtual void GatherInOutAverageValues(CConfig *config, bool allocate) {}
 
-  /*!
-   * \brief Sets CG coordinates.
-   */
-  inline virtual void SetCoord_CG(void) {}
-
   /*!
    * \brief Set max length.
    * \param[in] config - Definition of the particular problem.
@@ -705,9 +714,8 @@ class CGeometry {
   /*!
    * \brief A virtual member.
    * \param[in] config - Definition of the particular problem.
-   * \param[in] action - Allocate or not the new elements.
    */
-  inline virtual void VisualizeControlVolume(CConfig *config, unsigned short action) {}
+  inline virtual void VisualizeControlVolume(const CConfig *config) const {}
 
   /*!
    * \brief A virtual member.
@@ -732,7 +740,7 @@ class CGeometry {
    * \param[in] config - Definition of the particular problem.
    * \param[in] action - Allocate or not the new elements.
    */
-  inline virtual void SetBoundControlVolume(CConfig *config, unsigned short action) {}
+  inline virtual void SetBoundControlVolume(const CConfig *config, unsigned short action) {}
 
   /*!
    * \brief A virtual member.
@@ -1589,10 +1597,9 @@ class CGeometry {
                               const su2double *cg_elem, vector<long> &neighbours, vector<bool> &is_neighbor) const;
 
   /*!
-   * \brief Compute and store the volume of the elements.
-   * \param[in] config - Problem configuration.
+   * \brief Compute and store the volume of the primal elements.
    */
-  void SetElemVolume(CConfig *config);
+  void SetElemVolume();
 
   /*!
    * \brief Set the multigrid index for the current geometry object.
@@ -1610,7 +1617,7 @@ class CGeometry {
    * \brief A virtual member.
    * \param config - Config
    */
-  inline virtual void ComputeMeshQualityStatistics(CConfig *config) {}
+  inline virtual void ComputeMeshQualityStatistics(const CConfig *config) {}
 
   /*!
    * \brief Get the sparse pattern of "type" with given level of fill.

Common/include/geometry/CPhysicalGeometry.hpp

@@ -444,11 +444,6 @@ class CPhysicalGeometry final : public CGeometry {
    */
   void GatherInOutAverageValues(CConfig *config, bool allocate) override;
 
-  /*!
-   * \brief Set the center of gravity of the face, elements and edges.
-   */
-  void SetCoord_CG(void) override;
-
   /*!
    * \brief Set the edge structure of the control volume.
    * \param[in] config - Definition of the particular problem.
@@ -459,9 +454,8 @@ class CPhysicalGeometry final : public CGeometry {
   /*!
    * \brief Visualize the structure of the control volume(s).
    * \param[in] config - Definition of the particular problem.
-   * \param[in] action - Allocate or not the new elements.
    */
-  void VisualizeControlVolume(CConfig *config, unsigned short action) override;
+  void VisualizeControlVolume(const CConfig *config) const override;
 
   /*!
    * \brief Mach the near field boundary condition.
@@ -487,7 +481,7 @@ class CPhysicalGeometry final : public CGeometry {
    * \param[in] config - Definition of the particular problem.
    * \param[in] action - Allocate or not the new elements.
    */
-  void SetBoundControlVolume(CConfig *config, unsigned short action) override;
+  void SetBoundControlVolume(const CConfig *config, unsigned short action) override;
 
   /*!
    * \brief Set the maximum cell-center to cell-center distance for CVs.
@@ -598,7 +592,7 @@ class CPhysicalGeometry final : public CGeometry {
    * \brief Compute 3 grid quality metrics: orthogonality angle, dual cell aspect ratio, and dual cell volume ratio.
    * \param[in] config - Definition of the particular problem.
    */
-  void ComputeMeshQualityStatistics(CConfig *config) override;
+  void ComputeMeshQualityStatistics(const CConfig *config) override;
 
   /*!
    * \brief Find and store the closest neighbor to a vertex.

Common/include/geometry/dual_grid/CDualGrid.hpp

@@ -74,17 +74,16 @@ class CDualGrid{
    * \param[in] val_coord_Edge_CG - Coordinates of the centre of gravity of the edge.
    * \param[in] val_coord_FaceElem_CG - Coordinates of the centre of gravity of the face of an element.
    * \param[in] val_coord_Elem_CG - Coordinates of the centre of gravity of the element.
-   * \param[in] config - Definition of the particular problem.
    */
-  virtual void SetNodes_Coord(su2double *val_coord_Edge_CG, su2double *val_coord_FaceElem_CG, su2double *val_coord_Elem_CG) = 0;
+  virtual void SetNodes_Coord(const su2double *val_coord_Edge_CG, const su2double *val_coord_FaceElem_CG,
+                              const su2double *val_coord_Elem_CG) = 0;
 
   /*!
    * \overload
    * \param[in] val_coord_Edge_CG - Coordinates of the centre of gravity of the edge.
    * \param[in] val_coord_Elem_CG - Coordinates of the centre of gravity of the element.
-   * \param[in] config - Definition of the particular problem.
    */
-  virtual void SetNodes_Coord(su2double *val_coord_Edge_CG, su2double *val_coord_Elem_CG) = 0;
+  virtual void SetNodes_Coord(const su2double *val_coord_Edge_CG, const su2double *val_coord_Elem_CG) = 0;
 
   /*!
    * \brief A pure virtual member.

Common/include/geometry/dual_grid/CEdge.hpp

@@ -29,6 +29,8 @@
 
 #include "../../containers/C2DContainer.hpp"
 
+class CPhysicalGeometry;
+
 /*!
  * \class CEdge
  * \brief Class for defining the edges of the dual grid.
@@ -41,7 +43,8 @@ class CEdge {
   using NodeArray = C2DContainer<Index, Index, StorageType::ColumnMajor, 64, DynamicSize, 2>;
   NodeArray Nodes;           /*!< \brief Vector to store the node indices of the edge. */
   su2activematrix Normal;    /*!< \brief Normal (area) of the edge. */
-  su2activematrix Coord_CG;  /*!< \brief Center-of-gravity (mid point) of the edge. */
+
+  friend class CPhysicalGeometry;
 
 public:
   enum NodePosition : unsigned long {LEFT = 0, RIGHT = 1};
@@ -58,25 +61,6 @@ class CEdge {
    */
   CEdge() = delete;
 
-  /*!
-   * \brief Set the center of gravity of the edge.
-   * \param[in] iEdge - Edge index.
-   * \param[in] nodeCoord - Coordinates of the two nodes.
-   */
-  template<class T>
-  void SetCoord_CG(unsigned long iEdge, const T& nodeCoord) {
-    for (auto iDim = 0u; iDim < Coord_CG.cols(); ++iDim)
-      Coord_CG(iEdge,iDim) = 0.5 * (nodeCoord[0][iDim] + nodeCoord[1][iDim]);
-  }
-
-  /*!
-   * \brief Obtain the center of gravity of the edge.
-   * \param[in] iEdge - Edge index.
-   * \param[in] iDim - Dimension.
-   * \return Coordinate of the centre of gravity.
-   */
-  inline su2double GetCG(unsigned long iEdge, unsigned long iDim) const { return Coord_CG(iEdge,iDim); }
-
   /*!
    * \brief Get left/right node index defining the edge.
    * \param[in] iEdge - Edge index.

Common/include/geometry/dual_grid/CPoint.hpp

@@ -36,6 +36,7 @@
 using namespace std;
 
 class CConfig;
+class CPhysicalGeometry;
 
 /*!
  * \class CPoint
@@ -44,6 +45,8 @@ class CConfig;
  */
 class CPoint {
 private:
+  friend class CPhysicalGeometry;
+
   const unsigned long nDim = 0;
 
   su2vector<unsigned long> GlobalIndex;   /*!< \brief Global index in the parallel simulation. */

Common/include/geometry/dual_grid/CVertex.hpp

@@ -75,15 +75,16 @@ class CVertex : public CDualGrid {
    * \param[in] val_coord_Elem_CG - Coordinates of the centre of gravity of the element.
    * \return Compute the normal (dimensional) to the face that makes the vertex.
    */
-  void SetNodes_Coord(su2double *val_coord_Edge_CG, su2double *val_coord_FaceElem_CG, su2double *val_coord_Elem_CG) override;
+  void SetNodes_Coord(const su2double *val_coord_Edge_CG, const su2double *val_coord_FaceElem_CG,
+                      const su2double *val_coord_Elem_CG) override;
 
   /*!
    * \overload
    * \param[in] val_coord_Edge_CG - Coordinates of the centre of gravity of the edge.
    * \param[in] val_coord_Elem_CG - Coordinates of the centre of gravity of the element.
    * \return Compute the normal (dimensional) to the face that makes the vertex.
    */
-  void SetNodes_Coord(su2double *val_coord_Edge_CG, su2double *val_coord_Elem_CG) override;
+  void SetNodes_Coord(const su2double *val_coord_Edge_CG, const su2double *val_coord_Elem_CG) override;
 
   /*!
    * \brief Copy the the normal vector of a face.

Common/include/geometry/primal_grid/CPrimalGrid.hpp

@@ -49,9 +49,7 @@ class CPrimalGrid {
   long *Neighbor_Elements;      /*!< \brief Vector to store the elements surronding an element. */
   short *PeriodIndexNeighbors;  /*!< \brief Vector to store the periodic index of a neighbor.
                                             A -1 indicates no periodic transformation to the neighbor. */
-  su2double *Coord_CG;             /*!< \brief Coordinates of the center-of-gravity of the element. */
-  su2double **Coord_FaceElems_CG;  /*!< \brief Coordinates of the center-of-gravity of the face of the
-                                               elements. */
+  su2double Coord_CG[3] = {0.0}; /*!< \brief Coordinates of the center-of-gravity of the element. */
   static unsigned short nDim;    /*!< \brief Dimension of the element (2D or 3D) useful for triangles,
                                                quadrilateral and edges. */
   unsigned long DomainElement;     /*!< \brief Only for boundaries, in this variable the 3D elements which
@@ -170,14 +168,23 @@ class CPrimalGrid {
    * \brief Set the center of gravity of an element (including edges).
    * \param[in] val_coord - Coordinates of the element.
    */
-  void SetCoord_CG(const su2double* const* val_coord);
+  template<class T>
+  inline su2double* SetCoord_CG(const T& val_coord) {
+    for (unsigned short iDim = 0; iDim < nDim; iDim++) {
+      Coord_CG[iDim] = 0.0;
+      for (unsigned short iNode = 0; iNode < GetnNodes(); iNode++)
+        Coord_CG[iDim] += val_coord[iNode][iDim]/su2double(GetnNodes());
+    }
+    return Coord_CG;
+  }
 
   /*!
    * \brief Get the center of gravity of an element (including edges).
    * \param[in] val_dim - Coordinate of the center of gravity.
    * \return Coordinates of the center of gravity.
    */
-  inline su2double GetCG(unsigned short val_dim) { return Coord_CG[val_dim]; }
+  inline su2double GetCG(unsigned short val_dim) const { return Coord_CG[val_dim]; }
+  inline const su2double* GetCG() const { return Coord_CG; }
 
   /*!
    * \brief Set the center of gravity of an element (including edges).
@@ -192,14 +199,6 @@ class CPrimalGrid {
    */
   inline su2double GetVolume(void) const { return Volume; }
 
-  /*!
-   * \brief Get the CG of a face of an element.
-   * \param[in] val_face - Local index of the face.
-   * \param[in] val_dim - Coordinate of the center of gravity.
-   * \return Coordinates of the center of gravity.
-   */
-  inline su2double GetFaceCG(unsigned short val_face, unsigned short val_dim) { return Coord_FaceElems_CG[val_face][val_dim]; }
-
   /*!
    * \brief Get all the neighbors of an element.
    * \return List of all the neighbor of an element.
@@ -431,7 +430,7 @@ class CPrimalGrid {
    * \brief Virtual function to make available the number of donor elements for the wall function treatment.
    * \return The number of donor elements.
    */
-  inline  virtual unsigned short GetNDonorsWallFunctions(void) {return 0;}
+  inline virtual unsigned short GetNDonorsWallFunctions(void) {return 0;}
 
   /*!
    * \brief Virtual function to make available the pointer to the vector for the donor elements

Common/include/linear_algebra/vector_expressions.hpp

@@ -157,10 +157,16 @@ FORCEINLINE auto FUN(decay_t<S> u, const CVecExpr<V,S>& v)                    \
   RETURNS( EXPR<Bcast<S>,V,S>(Bcast<S>(u), v.derived())                       \
 )                                                                             \
 
-/*--- std::max/min have issues (maybe because they return by reference). ---*/
+/*--- std::max/min have issues (maybe because they return by reference).
+ * For AD codi::max/min need to be used to avoid issues in debug builds. ---*/
 
+#if defined(CODI_REVERSE_TYPE) || defined(CODI_FORWARD_TYPE)
+#define max_impl math::max
+#define min_impl math::min
+#else
 #define max_impl(a,b) a<b? Scalar(b) : Scalar(a)
 #define min_impl(a,b) b<a? Scalar(b) : Scalar(a)
+#endif
 MAKE_BINARY_FUN(max, max_, max_impl)
 MAKE_BINARY_FUN(min, min_, min_impl)
 MAKE_BINARY_FUN(pow, pow_, math::pow)

Common/src/geometry/CGeometry.cpp

@@ -1389,26 +1389,6 @@ su2double CGeometry::Point2Plane_Distance(const su2double *Coord, const su2doubl
 
 }
 
-long CGeometry::FindEdge(unsigned long first_point, unsigned long second_point) const {
-
-  for (unsigned short iNode = 0; iNode < nodes->GetnPoint(first_point); iNode++) {
-    auto iPoint = nodes->GetPoint(first_point, iNode);
-    if (iPoint == second_point) return nodes->GetEdge(first_point, iNode);
-  }
-
-  char buf[100];
-  SPRINTF(buf, "Can't find the edge that connects %lu and %lu.", first_point, second_point);
-  SU2_MPI::Error(buf, CURRENT_FUNCTION);
-  return 0;
-}
-
-bool CGeometry::CheckEdge(unsigned long first_point, unsigned long second_point) const {
-
-  for (auto iPoint : nodes->GetPoints(first_point))
-    if (iPoint == second_point) return true;
-  return false;
-}
-
 void CGeometry::SetEdges(void) {
 
   nEdge = 0;
@@ -2545,7 +2525,6 @@ void CGeometry::UpdateGeometry(CGeometry **geometry_container, CConfig *config)
     geometry_container[MESH_0]->CompleteComms(geometry_container[MESH_0], config, GRID_VELOCITY);
   }
 
-  geometry_container[MESH_0]->SetCoord_CG();
   geometry_container[MESH_0]->SetControlVolume(config, UPDATE);
   geometry_container[MESH_0]->SetBoundControlVolume(config, UPDATE);
   geometry_container[MESH_0]->SetMaxLength(config);
@@ -3499,7 +3478,7 @@ bool CGeometry::GetRadialNeighbourhood(const unsigned long iElem_global,
   return finished;
 }
 
-void CGeometry::SetElemVolume(CConfig *config)
+void CGeometry::SetElemVolume()
 {
   SU2_OMP_PARALLEL
   {