diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp index e2e9488a36f4..be6d247ad5b2 100644 --- a/Common/include/geometry/CGeometry.hpp +++ b/Common/include/geometry/CGeometry.hpp @@ -492,9 +492,21 @@ 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. @@ -502,7 +514,9 @@ class CGeometry { * \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 &neighbours, vector &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. diff --git a/Common/include/geometry/CPhysicalGeometry.hpp b/Common/include/geometry/CPhysicalGeometry.hpp index 84e1acdd89c8..c833c843a6de 100644 --- a/Common/include/geometry/CPhysicalGeometry.hpp +++ b/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. diff --git a/Common/include/geometry/dual_grid/CDualGrid.hpp b/Common/include/geometry/dual_grid/CDualGrid.hpp index 3044f8556cd1..8d4d6120c556 100644 --- a/Common/include/geometry/dual_grid/CDualGrid.hpp +++ b/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. diff --git a/Common/include/geometry/dual_grid/CEdge.hpp b/Common/include/geometry/dual_grid/CEdge.hpp index b9d0de79d428..845bfda95c7e 100644 --- a/Common/include/geometry/dual_grid/CEdge.hpp +++ b/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; 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 - 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. diff --git a/Common/include/geometry/dual_grid/CPoint.hpp b/Common/include/geometry/dual_grid/CPoint.hpp index 499c9829fa95..6e504444b190 100644 --- a/Common/include/geometry/dual_grid/CPoint.hpp +++ b/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 GlobalIndex; /*!< \brief Global index in the parallel simulation. */ diff --git a/Common/include/geometry/dual_grid/CVertex.hpp b/Common/include/geometry/dual_grid/CVertex.hpp index a0b200115d7a..39522889e54c 100644 --- a/Common/include/geometry/dual_grid/CVertex.hpp +++ b/Common/include/geometry/dual_grid/CVertex.hpp @@ -75,7 +75,8 @@ 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 @@ -83,7 +84,7 @@ 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_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. diff --git a/Common/include/geometry/primal_grid/CPrimalGrid.hpp b/Common/include/geometry/primal_grid/CPrimalGrid.hpp index 82b49169a708..39cea851d5f4 100644 --- a/Common/include/geometry/primal_grid/CPrimalGrid.hpp +++ b/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 + 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 diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp index 55d819eceb30..a44534dedc4b 100644 --- a/Common/include/linear_algebra/vector_expressions.hpp +++ b/Common/include/linear_algebra/vector_expressions.hpp @@ -157,10 +157,16 @@ FORCEINLINE auto FUN(decay_t u, const CVecExpr& v) \ RETURNS( EXPR,V,S>(Bcast(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) aGetnPoint(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 { diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp index 6f2e5526c08a..e6e7e729115b 100644 --- a/Common/src/geometry/CPhysicalGeometry.cpp +++ b/Common/src/geometry/CPhysicalGeometry.cpp @@ -6677,148 +6677,6 @@ void CPhysicalGeometry::GatherInOutAverageValues(CConfig *config, bool allocate) } - -void CPhysicalGeometry::SetCoord_CG(void) { - - unsigned short iMarker, iNode; - unsigned long elem_poin, edge_poin, iElem, iEdge; - - /*--- Buffer of pointers to node coordinates ---*/ - array Coord; - - /*--- Compute the center of gravity for elements ---*/ - - SU2_OMP_FOR_STAT(roundUpDiv(nElem,2*omp_get_max_threads())) - for (iElem = 0; iElemGetnNodes() <= N_POINTS_MAXIMUM && "Insufficient N_POINTS_MAXIMUM"); - - /*--- Store the coordinates for all the element nodes ---*/ - for (iNode = 0; iNode < elem[iElem]->GetnNodes(); iNode++) { - elem_poin = elem[iElem]->GetNode(iNode); - Coord[iNode] = nodes->GetCoord(elem_poin); - } - - /*--- Compute the element CG coordinates ---*/ - elem[iElem]->SetCoord_CG(Coord.data()); - } - - /*--- Center of gravity for face elements ---*/ - - SU2_OMP_FOR_DYN(1) - for (iMarker = 0; iMarker < nMarker; iMarker++) { - for (iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) { - - /*--- Store the coordinates for all the element nodes ---*/ - for (iNode = 0; iNode < bound[iMarker][iElem]->GetnNodes(); iNode++) { - elem_poin = bound[iMarker][iElem]->GetNode(iNode); - Coord[iNode] = nodes->GetCoord(elem_poin); - } - - /*--- Compute the element CG coordinates ---*/ - bound[iMarker][iElem]->SetCoord_CG(Coord.data()); - } - } - - /*--- Center of gravity for edges ---*/ - - SU2_OMP_FOR_STAT(roundUpDiv(nEdge,2*omp_get_max_threads())) - for (iEdge = 0; iEdge < nEdge; iEdge++) { - - /*--- Store the coordinates for all the element nodes ---*/ - for (iNode = 0; iNode < edges->GetnNodes(); iNode++) { - edge_poin=edges->GetNode(iEdge,iNode); - Coord[iNode] = nodes->GetCoord(edge_poin); - } - - /*--- Compute the edge CG coordinates ---*/ - edges->SetCoord_CG(iEdge, Coord.data()); - } - -} - -void CPhysicalGeometry::SetBoundControlVolume(CConfig *config, unsigned short action) { - - SU2_OMP_MASTER { - - unsigned short Neighbor_Node, iMarker, iNode, iNeighbor_Nodes, iDim; - unsigned long Neighbor_Point, iVertex, iPoint, iElem; - long iEdge; - su2double Area, *NormalFace = nullptr; - - /*--- Update values of faces of the edge ---*/ - - if (action != ALLOCATE) - for (iMarker = 0; iMarker < nMarker; iMarker++) - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) - vertex[iMarker][iVertex]->SetZeroValues(); - - su2double *Coord_Edge_CG = new su2double [nDim]; - su2double *Coord_Elem_CG = new su2double [nDim]; - su2double *Coord_Vertex = new su2double [nDim]; - - /*--- Loop over all the markers ---*/ - - for (iMarker = 0; iMarker < nMarker; iMarker++) - - /*--- Loop over all the boundary elements ---*/ - - for (iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) - - /*--- Loop over all the nodes of the boundary ---*/ - - for (iNode = 0; iNode < bound[iMarker][iElem]->GetnNodes(); iNode++) { - iPoint = bound[iMarker][iElem]->GetNode(iNode); - iVertex = nodes->GetVertex(iPoint, iMarker); - - /*--- Loop over the neighbor nodes, there is a face for each one ---*/ - - for (iNeighbor_Nodes = 0; iNeighbor_Nodes < bound[iMarker][iElem]->GetnNeighbor_Nodes(iNode); iNeighbor_Nodes++) { - Neighbor_Node = bound[iMarker][iElem]->GetNeighbor_Nodes(iNode, iNeighbor_Nodes); - Neighbor_Point = bound[iMarker][iElem]->GetNode(Neighbor_Node); - - /*--- Shared edge by the Neighbor Point and the point ---*/ - - iEdge = FindEdge(iPoint, Neighbor_Point); - for (iDim = 0; iDim < nDim; iDim++) { - Coord_Edge_CG[iDim] = edges->GetCG(iEdge,iDim); - Coord_Elem_CG[iDim] = bound[iMarker][iElem]->GetCG(iDim); - Coord_Vertex[iDim] = nodes->GetCoord(iPoint, iDim); - } - switch (nDim) { - case 2: - - /*--- Store the 2D face ---*/ - - if (iNode == 0) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Elem_CG, Coord_Vertex); - if (iNode == 1) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Vertex, Coord_Elem_CG); - break; - case 3: - - /*--- Store the 3D face ---*/ - - if (iNeighbor_Nodes == 0) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Elem_CG, Coord_Edge_CG, Coord_Vertex); - if (iNeighbor_Nodes == 1) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Edge_CG, Coord_Elem_CG, Coord_Vertex); - break; - } - } - } - - delete[] Coord_Edge_CG; - delete[] Coord_Elem_CG; - delete[] Coord_Vertex; - - /*--- Check if there is a normal with null area ---*/ - - for (iMarker = 0; iMarker < nMarker; iMarker ++) - for (iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { - NormalFace = vertex[iMarker][iVertex]->GetNormal(); - Area = GeometryToolbox::Norm(nDim, NormalFace); - if (Area == 0.0) for (iDim = 0; iDim < nDim; iDim++) NormalFace[iDim] = EPS*EPS; - } - - } SU2_OMP_BARRIER -} - void CPhysicalGeometry::SetMaxLength(CConfig* config) { SU2_OMP_FOR_STAT(roundUpDiv(nPointDomain,omp_get_max_threads())) @@ -7637,59 +7495,98 @@ void CPhysicalGeometry::MatchPeriodic(CConfig *config, void CPhysicalGeometry::SetControlVolume(CConfig *config, unsigned short action) { - SU2_OMP_MASTER { - - unsigned long face_iPoint = 0, face_jPoint = 0, iPoint, iElem; - unsigned short nEdgesFace = 1, iFace, iEdgesFace, iDim; - su2double Coord_Edge_CG[3] = {0.0}, Coord_FaceElem_CG[3] = {0.0}, Coord_Elem_CG[3] = {0.0}; - /*--- Update values of faces of the edge ---*/ if (action != ALLOCATE) { - edges->SetZeroValues(); - for (iPoint = 0; iPoint < nPoint; iPoint++) + su2double ZeroArea[MAXNDIM] = {0.0}; + + SU2_OMP_FOR_STAT(1024) + for (auto iEdge = 0ul; iEdge < nEdge; iEdge++) + edges->SetNormal(iEdge, ZeroArea); + + SU2_OMP_FOR_STAT(1024) + for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) nodes->SetVolume(iPoint, 0.0); } + SU2_OMP_MASTER { /*--- The following is difficult to parallelize with threads. ---*/ + su2double my_DomainVolume = 0.0; - for (iElem = 0; iElem < nElem; iElem++) { - for (iFace = 0; iFace < elem[iElem]->GetnFaces(); iFace++) { + for (auto iElem = 0ul; iElem < nElem; iElem++) { + + const auto nNodes = elem[iElem]->GetnNodes(); + + /*--- To make preaccumulation more effective, use as few inputs + as possible, recomputing intermediate quantities as needed. ---*/ + AD::StartPreacc(); + + /*--- Get pointers to the coordinates of all the element nodes ---*/ + array Coord; + + for (unsigned short iNode = 0; iNode < nNodes; iNode++) { + auto iPoint = elem[iElem]->GetNode(iNode); + Coord[iNode] = nodes->GetCoord(iPoint); +#ifdef CODI_REVERSE_TYPE + /*--- The same points and edges will be referenced multiple times as they are common + to many of the element's faces, therefore they are "registered" here only once. ---*/ + AD::SetPreaccIn(nodes->Volume(iPoint)); + for (unsigned short jNode = iNode+1; jNode < nNodes; jNode++) { + auto jPoint = elem[iElem]->GetNode(jNode); + auto iEdge = FindEdge(iPoint, jPoint, false); + if (iEdge >= 0) AD::SetPreaccIn(edges->Normal[iEdge], nDim); + } +#endif + } + AD::SetPreaccIn(Coord, nNodes, nDim); + + /*--- Compute the element median CG coordinates ---*/ + auto Coord_Elem_CG = elem[iElem]->SetCoord_CG(Coord); + AD::SetPreaccOut(Coord_Elem_CG, nDim); + + for (unsigned short iFace = 0; iFace < elem[iElem]->GetnFaces(); iFace++) { /*--- In 2D all the faces have only one edge ---*/ - if (nDim == 2) nEdgesFace = 1; - /*--- In 3D the number of edges per face is the same as the number of point per face ---*/ - if (nDim == 3) nEdgesFace = elem[iElem]->GetnNodesFace(iFace); + unsigned short nEdgesFace = 1; + + /*--- In 3D the number of edges per face is the same as the number of point + per face and the median CG of the face is needed. ---*/ + su2double Coord_FaceElem_CG[MAXNDIM] = {0.0}; + if (nDim == 3) { + nEdgesFace = elem[iElem]->GetnNodesFace(iFace); + + for (unsigned short iNode = 0; iNode < nEdgesFace; iNode++) { + auto NodeFace = elem[iElem]->GetFaces(iFace, iNode); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Coord_FaceElem_CG[iDim] += Coord[NodeFace][iDim]/nEdgesFace; + } + } /*-- Loop over the edges of a face ---*/ - for (iEdgesFace = 0; iEdgesFace < nEdgesFace; iEdgesFace++) { + for (unsigned short iEdgesFace = 0; iEdgesFace < nEdgesFace; iEdgesFace++) { + + const auto face_iNode = elem[iElem]->GetFaces(iFace,iEdgesFace); + unsigned short face_jNode; - /*--- In 2D only one edge (two points) per edge ---*/ if (nDim == 2) { - face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,0)); - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,1)); + /*--- In 2D only one edge (two points) per edge ---*/ + face_jNode = elem[iElem]->GetFaces(iFace,1); } - - /*--- In 3D there are several edges in each face ---*/ - if (nDim == 3) { - face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, iEdgesFace)); - if (iEdgesFace != nEdgesFace-1) - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, iEdgesFace+1)); - else - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,0)); + else { + /*--- In 3D we "circle around" the face ---*/ + face_jNode = elem[iElem]->GetFaces(iFace, (iEdgesFace+1)%nEdgesFace); } + const auto face_iPoint = elem[iElem]->GetNode(face_iNode); + const auto face_jPoint = elem[iElem]->GetNode(face_jNode); + /*--- We define a direction (from the smalest index to the greatest) --*/ const bool change_face_orientation = (face_iPoint > face_jPoint); const auto iEdge = FindEdge(face_iPoint, face_jPoint); - for (iDim = 0; iDim < nDim; iDim++) { - Coord_Edge_CG[iDim] = edges->GetCG(iEdge,iDim); - Coord_Elem_CG[iDim] = elem[iElem]->GetCG(iDim); - Coord_FaceElem_CG[iDim] = elem[iElem]->GetFaceCG(iFace, iDim); + su2double Coord_Edge_CG[MAXNDIM] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Coord_Edge_CG[iDim] = 0.5 * (Coord[face_iNode][iDim] + Coord[face_jNode][iDim]); } - const su2double* Coord_FaceiPoint = nodes->GetCoord(face_iPoint); - const su2double* Coord_FacejPoint = nodes->GetCoord(face_jPoint); - su2double Volume_i, Volume_j; if (nDim == 2) { @@ -7699,8 +7596,8 @@ void CPhysicalGeometry::SetControlVolume(CConfig *config, unsigned short action) else edges->SetNodes_Coord(iEdge, Coord_Edge_CG, Coord_Elem_CG); - Volume_i = CEdge::GetVolume(Coord_FaceiPoint, Coord_Edge_CG, Coord_Elem_CG); - Volume_j = CEdge::GetVolume(Coord_FacejPoint, Coord_Edge_CG, Coord_Elem_CG); + Volume_i = CEdge::GetVolume(Coord[face_iNode], Coord_Edge_CG, Coord_Elem_CG); + Volume_j = CEdge::GetVolume(Coord[face_jNode], Coord_Edge_CG, Coord_Elem_CG); } else { /*--- Three dimensional problem ---*/ @@ -7709,8 +7606,8 @@ void CPhysicalGeometry::SetControlVolume(CConfig *config, unsigned short action) else edges->SetNodes_Coord(iEdge, Coord_Edge_CG, Coord_FaceElem_CG, Coord_Elem_CG); - Volume_i = CEdge::GetVolume(Coord_FaceiPoint, Coord_Edge_CG, Coord_FaceElem_CG, Coord_Elem_CG); - Volume_j = CEdge::GetVolume(Coord_FacejPoint, Coord_Edge_CG, Coord_FaceElem_CG, Coord_Elem_CG); + Volume_i = CEdge::GetVolume(Coord[face_iNode], Coord_Edge_CG, Coord_FaceElem_CG, Coord_Elem_CG); + Volume_j = CEdge::GetVolume(Coord[face_jNode], Coord_Edge_CG, Coord_FaceElem_CG, Coord_Elem_CG); } nodes->AddVolume(face_iPoint, Volume_i); @@ -7719,15 +7616,19 @@ void CPhysicalGeometry::SetControlVolume(CConfig *config, unsigned short action) my_DomainVolume += Volume_i+Volume_j; } } - } - /*--- Check if there is a normal with null area ---*/ - for (auto iEdge = 0ul; iEdge < nEdge; iEdge++) { - const auto Area2 = GeometryToolbox::SquaredNorm(nDim, edges->GetNormal(iEdge)); - if (Area2 == 0.0) { - su2double DefaultArea[3] = {EPS*EPS}; - edges->SetNormal(iEdge, DefaultArea); +#ifdef CODI_REVERSE_TYPE + for (unsigned short iNode = 0; iNode < nNodes; iNode++) { + auto iPoint = elem[iElem]->GetNode(iNode); + AD::SetPreaccOut(nodes->Volume(iPoint)); + for (unsigned short jNode = iNode+1; jNode < nNodes; jNode++) { + auto jPoint = elem[iElem]->GetNode(jNode); + auto iEdge = FindEdge(iPoint, jPoint, false); + if (iEdge >= 0) AD::SetPreaccOut(edges->Normal[iEdge], nDim); + } } +#endif + AD::EndPreacc(); } su2double DomainVolume; @@ -7740,77 +7641,157 @@ void CPhysicalGeometry::SetControlVolume(CConfig *config, unsigned short action) } } SU2_OMP_BARRIER + + /*--- Check if there is a normal with null area ---*/ + SU2_OMP_FOR_STAT(1024) + for (auto iEdge = 0ul; iEdge < nEdge; iEdge++) { + const auto Area2 = GeometryToolbox::SquaredNorm(nDim, edges->GetNormal(iEdge)); + su2double DefaultArea[MAXNDIM] = {EPS*EPS}; + if (Area2 == 0.0) edges->SetNormal(iEdge, DefaultArea); + } } -void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short action) { +void CPhysicalGeometry::SetBoundControlVolume(const CConfig *config, unsigned short action) { - /*--- This routine is only meant for visualization in serial currently ---*/ -#ifndef HAVE_MPI + /*--- Clear normals ---*/ - unsigned long face_iPoint = 0, face_jPoint = 0, iElem, iPoint_Viz; - long iEdge; - unsigned short nEdgesFace = 1, iFace, iEdgesFace, iDim; - su2double *Coord_Edge_CG, *Coord_FaceElem_CG, *Coord_Elem_CG, *Coord_FaceiPoint, - *Coord_FacejPoint; - int counter = 0; - char cstr[MAX_STRING_SIZE], buffer[50]; - ofstream Tecplot_File; - string mesh_filename; - vector X, Y, Z, X_n, Y_n, Z_n; - su2double r1[3], r2[3], CrossProduct[3]; + if (action != ALLOCATE) { + SU2_OMP_FOR_DYN(1) + for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) + for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) + vertex[iMarker][iVertex]->SetZeroValues(); + } - /*--- Access the point number for control volume we want to vizualize ---*/ + /*--- Loop over all the boundary elements ---*/ + + SU2_OMP_FOR_DYN(1) + for (unsigned short iMarker = 0; iMarker < nMarker; iMarker++) { + for (unsigned long iElem = 0; iElem < nElem_Bound[iMarker]; iElem++) { + + const auto nNodes = bound[iMarker][iElem]->GetnNodes(); + + AD::StartPreacc(); + + /*--- Get pointers to the coordinates of all the element nodes ---*/ + array Coord; + + for (unsigned short iNode = 0; iNode < nNodes; iNode++) { + const auto iPoint = bound[iMarker][iElem]->GetNode(iNode); + const auto iVertex = nodes->GetVertex(iPoint, iMarker); + Coord[iNode] = nodes->GetCoord(iPoint); + AD::SetPreaccIn(vertex[iMarker][iVertex]->GetNormal(), nDim); + } + AD::SetPreaccIn(Coord, nNodes, nDim); + + /*--- Compute the element CG coordinates ---*/ + auto Coord_Elem_CG = bound[iMarker][iElem]->SetCoord_CG(Coord); + AD::SetPreaccOut(Coord_Elem_CG, nDim); + + /*--- Loop over all the nodes of the boundary element ---*/ - iPoint_Viz = config->GetVisualize_CV(); + for (unsigned short iNode = 0; iNode < nNodes; iNode++) { + const auto iPoint = bound[iMarker][iElem]->GetNode(iNode); + const auto iVertex = nodes->GetVertex(iPoint, iMarker); + auto Coord_Vertex = Coord[iNode]; - /*--- Allocate some structures for building the dual CVs ---*/ + /*--- Loop over the neighbor nodes, there is a face for each one ---*/ + + for (unsigned short iNeighbor = 0; iNeighbor < bound[iMarker][iElem]->GetnNeighbor_Nodes(iNode); iNeighbor++) { + if (nDim == 2) { + /*--- Store the 2D face ---*/ + if (iNode == 0) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Elem_CG, Coord_Vertex); + if (iNode == 1) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Vertex, Coord_Elem_CG); + } + else { + const auto Neighbor_Node = bound[iMarker][iElem]->GetNeighbor_Nodes(iNode, iNeighbor); + auto Neighbor_Coord = Coord[Neighbor_Node]; + + /*--- Store the 3D face ---*/ + su2double Coord_Edge_CG[MAXNDIM] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Coord_Edge_CG[iDim] = 0.5 * (Coord_Vertex[iDim] + Neighbor_Coord[iDim]); + + if (iNeighbor == 0) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Elem_CG, Coord_Edge_CG, Coord_Vertex); + if (iNeighbor == 1) vertex[iMarker][iVertex]->SetNodes_Coord(Coord_Edge_CG, Coord_Elem_CG, Coord_Vertex); + } + } + } + + for (unsigned short iNode = 0; iNode < nNodes; iNode++) { + const auto iPoint = bound[iMarker][iElem]->GetNode(iNode); + const auto iVertex = nodes->GetVertex(iPoint, iMarker); + AD::SetPreaccOut(vertex[iMarker][iVertex]->GetNormal(), nDim); + } + AD::EndPreacc(); + } + } + + /*--- Check if there is a normal with null area ---*/ + + SU2_OMP_FOR_DYN(1) + for (unsigned short iMarker = 0; iMarker < nMarker; iMarker ++) { + for (unsigned long iVertex = 0; iVertex < nVertex[iMarker]; iVertex++) { + auto Area2 = GeometryToolbox::SquaredNorm(nDim, vertex[iMarker][iVertex]->GetNormal()); + su2double DefaultArea[MAXNDIM] = {EPS*EPS}; + if (Area2 == 0.0) vertex[iMarker][iVertex]->SetNormal(DefaultArea); + } + } +} + +void CPhysicalGeometry::VisualizeControlVolume(const CConfig *config) const { - Coord_Edge_CG = new su2double [nDim]; - Coord_FaceElem_CG = new su2double [nDim]; - Coord_Elem_CG = new su2double [nDim]; - Coord_FaceiPoint = new su2double [nDim]; - Coord_FacejPoint = new su2double [nDim]; + /*--- Access the point number for control volume we want to vizualize ---*/ + + auto iPoint = GetGlobal_to_Local_Point(config->GetVisualize_CV()); + if (iPoint < 0) return; + const unsigned long iPoint_Viz = iPoint; + + int counter = 0; + vector X, Y, Z; /*--- Loop over each face of each element ---*/ - CrossProduct[0] = 0.0; CrossProduct[1] = 0.0; CrossProduct[2] = 0.0; + for (auto iElem = 0ul; iElem < nElem; iElem++) { - for (iElem = 0; iElem < nElem; iElem++) { + /*--- Get pointers to the coordinates of all the element nodes ---*/ + array Coord; - for (iFace = 0; iFace < elem[iElem]->GetnFaces(); iFace++) { + for (unsigned short iNode = 0; iNode < elem[iElem]->GetnNodes(); iNode++) { + auto elem_poin = elem[iElem]->GetNode(iNode); + Coord[iNode] = nodes->GetCoord(elem_poin); + } + + for (unsigned short iFace = 0; iFace < elem[iElem]->GetnFaces(); iFace++) { /*--- In 2D all the faces have only one edge ---*/ - if (nDim == 2) nEdgesFace = 1; - /*--- In 3D the number of edges per face is the same as the number of point per face ---*/ - if (nDim == 3) nEdgesFace = elem[iElem]->GetnNodesFace(iFace); + unsigned short nEdgesFace = 1; - /*-- Loop over the edges of a face ---*/ - for (iEdgesFace = 0; iEdgesFace < nEdgesFace; iEdgesFace++) { + /*--- In 3D the number of edges per face is the same as the number of point + per face and the median CG of the face is needed. ---*/ + su2double Coord_FaceElem_CG[MAXNDIM] = {0.0}; + if (nDim == 3) { + nEdgesFace = elem[iElem]->GetnNodesFace(iFace); - /*--- In 2D only one edge (two points) per edge ---*/ - if (nDim == 2) { - face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,0)); - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,1)); + for (unsigned short iNode = 0; iNode < nEdgesFace; iNode++) { + auto NodeFace = elem[iElem]->GetFaces(iFace, iNode); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Coord_FaceElem_CG[iDim] += Coord[NodeFace][iDim]/nEdgesFace; } + } - /*--- In 3D there are several edges in each face ---*/ - if (nDim == 3) { - face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, iEdgesFace)); - if (iEdgesFace != nEdgesFace-1) - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, iEdgesFace+1)); - else - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,0)); - } + /*-- Loop over the edges of a face ---*/ + for (unsigned short iEdgesFace = 0; iEdgesFace < nEdgesFace; iEdgesFace++) { - /*--- We define a direction (from the smallest index to the greatest) --*/ - iEdge = FindEdge(face_iPoint, face_jPoint); + const auto face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,iEdgesFace)); + unsigned long face_jPoint; - for (iDim = 0; iDim < nDim; iDim++) { - Coord_Edge_CG[iDim] = edges->GetCG(iEdge,iDim); - Coord_Elem_CG[iDim] = elem[iElem]->GetCG(iDim); - Coord_FaceElem_CG[iDim] = elem[iElem]->GetFaceCG(iFace, iDim); - Coord_FaceiPoint[iDim] = nodes->GetCoord(face_iPoint, iDim); - Coord_FacejPoint[iDim] = nodes->GetCoord(face_jPoint, iDim); + if (nDim == 2) { + /*--- In 2D only one edge (two points) per edge ---*/ + face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,1)); + } + else { + /*--- In 3D we "circle around" the face ---*/ + face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, (iEdgesFace+1)%nEdgesFace)); } /*--- Print out the coordinates for a set of triangles making @@ -7818,21 +7799,22 @@ void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short a if (face_iPoint == iPoint_Viz || face_jPoint == iPoint_Viz) { + const su2double* Coord_FaceiPoint = nodes->GetCoord(face_iPoint); + const su2double* Coord_FacejPoint = nodes->GetCoord(face_jPoint); + const su2double* Coord_Elem_CG = elem[iElem]->GetCG(); + + su2double Coord_Edge_CG[MAXNDIM] = {0.0}; + for (unsigned short iDim = 0; iDim < nDim; iDim++) { + Coord_Edge_CG[iDim] = 0.5 * (Coord_FaceiPoint[iDim] + Coord_FacejPoint[iDim]); + } + if (nDim == 2) { X.push_back(Coord_Elem_CG[0]); X.push_back(Coord_Edge_CG[0]); Y.push_back(Coord_Elem_CG[1]); Y.push_back(Coord_Edge_CG[1]); - } else if (nDim == 3) { + } else { X.push_back(Coord_FaceElem_CG[0]); X.push_back(Coord_Edge_CG[0]); X.push_back(Coord_Elem_CG[0]); Y.push_back(Coord_FaceElem_CG[1]); Y.push_back(Coord_Edge_CG[1]); Y.push_back(Coord_Elem_CG[1]); Z.push_back(Coord_FaceElem_CG[2]); Z.push_back(Coord_Edge_CG[2]); Z.push_back(Coord_Elem_CG[2]); - - for (iDim = 0; iDim < nDim; iDim++) { - r1[iDim] = Coord_FaceElem_CG[iDim]-Coord_Elem_CG[iDim]; - r2[iDim] = Coord_Edge_CG[iDim]-Coord_Elem_CG[iDim]; - } - CrossProduct[0] += 0.5*(r1[1]*r2[2] - r1[2]*r2[1]); - CrossProduct[1] += 0.5*(r1[2]*r2[0] - r1[0]*r2[2]); - CrossProduct[2] += 0.5*(r1[0]*r2[1] - r1[1]*r2[0]); } counter++; } @@ -7842,10 +7824,9 @@ void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short a /*--- Write a Tecplot file to visualize the CV ---*/ - strcpy(cstr,"dual_cv"); - SPRINTF (buffer, "_%d.dat", SU2_TYPE::Int(iPoint_Viz)); - strcat(cstr, buffer); - + ofstream Tecplot_File; + char cstr[MAX_STRING_SIZE]; + SPRINTF(cstr, "dual_cv_%d.dat", config->GetVisualize_CV()); Tecplot_File.open(cstr); Tecplot_File << "TITLE= \"Visualization of the control volume\"" << endl; @@ -7853,7 +7834,7 @@ void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short a Tecplot_File << "VARIABLES = \"x\",\"y\" " << endl; Tecplot_File << "ZONE NODES= "<< counter*2 <<", ELEMENTS= "; Tecplot_File << counter <<", DATAPACKING=POINT, ZONETYPE=FEQUADRILATERAL"<< endl; - } if (nDim == 3) { + } else { Tecplot_File << "VARIABLES = \"x\",\"y\",\"z\" " << endl; Tecplot_File << "ZONE NODES= "<< counter*3 <<", ELEMENTS= "; Tecplot_File << counter <<", DATAPACKING=POINT, ZONETYPE=FEBRICK"<< endl; @@ -7862,7 +7843,7 @@ void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short a /*--- Write coordinates for the nodes in the order that they were found for each of the edges/triangles making up a dual control volume. ---*/ - for (vector::size_type i = 0; i != X.size(); i++) { + for (size_t i = 0; i != X.size(); i++) { Tecplot_File << X[i] << "\t" << Y[i]; if (nDim == 3) Tecplot_File << "\t" << Z[i]; Tecplot_File << "\n"; @@ -7870,31 +7851,16 @@ void CPhysicalGeometry::VisualizeControlVolume(CConfig *config, unsigned short a /*--- Create a new connectivity table in the order the faces were found ---*/ - int j; - for (int i= 0; i < counter; i++) { + for (int i = 0, j; i < counter; i++) { if (nDim == 2) { j = i*2; Tecplot_File << j+1 <<"\t"< Coord; - vector Coord_Edge_CG(nDim); - vector Coord_FaceElem_CG(nDim); - vector Coord_Elem_CG(nDim); - vector Coord_FaceiPoint(nDim); - vector Coord_FacejPoint(nDim); + for (unsigned short iNode = 0; iNode < elem[iElem]->GetnNodes(); iNode++) { + auto elem_poin = elem[iElem]->GetNode(iNode); + Coord[iNode] = nodes->GetCoord(elem_poin); + } + + const su2double* Coord_Elem_CG = elem[iElem]->GetCG(); - for (unsigned long iElem = 0; iElem < nElem; iElem++) for (unsigned short iFace = 0; iFace < elem[iElem]->GetnFaces(); iFace++) { /*--- In 2D all the faces have only one edge ---*/ + unsigned short nEdgesFace = 1; - if (nDim == 2) nEdgesFace = 1; - - /*--- In 3D the number of edges per face is the same - as the number of points per face. ---*/ + /*--- In 3D the number of edges per face is the same as the number of point + per face and the median CG of the face is needed. ---*/ + su2double Coord_FaceElem_CG[MAXNDIM] = {0.0}; + if (nDim == 3) { + nEdgesFace = elem[iElem]->GetnNodesFace(iFace); - if (nDim == 3) nEdgesFace = elem[iElem]->GetnNodesFace(iFace); + for (unsigned short iNode = 0; iNode < nEdgesFace; iNode++) { + auto NodeFace = elem[iElem]->GetFaces(iFace, iNode); + for (unsigned short iDim = 0; iDim < nDim; iDim++) + Coord_FaceElem_CG[iDim] += Coord[NodeFace][iDim]/nEdgesFace; + } + } /*-- Loop over the edges of a face ---*/ - for (unsigned short iEdgesFace = 0; iEdgesFace < nEdgesFace; iEdgesFace++) { - /*--- In 2D only one edge (two points) per edge ---*/ + const auto face_iNode = elem[iElem]->GetFaces(iFace,iEdgesFace); + unsigned short face_jNode; if (nDim == 2) { - face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,0)); - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,1)); + /*--- In 2D only one edge (two points) per edge ---*/ + face_jNode = elem[iElem]->GetFaces(iFace,1); } - - /*--- In 3D there are several edges in each face ---*/ - - if (nDim == 3) { - face_iPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, iEdgesFace)); - if (iEdgesFace != nEdgesFace-1) - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace, iEdgesFace+1)); - else - face_jPoint = elem[iElem]->GetNode(elem[iElem]->GetFaces(iFace,0)); + else { + /*--- In 3D we "circle around" the face ---*/ + face_jNode = elem[iElem]->GetFaces(iFace, (iEdgesFace+1)%nEdgesFace); } - /*--- Locate the edge for the two selected points. ---*/ - - const unsigned long iEdge = FindEdge(face_iPoint, face_jPoint); - - /*--- Collect the CG and coordinates for this sub-element face. ---*/ + const auto face_iPoint = elem[iElem]->GetNode(face_iNode); + const auto face_jPoint = elem[iElem]->GetNode(face_jNode); + su2double Coord_Edge_CG[MAXNDIM] = {0.0}; for (unsigned short iDim = 0; iDim < nDim; iDim++) { - Coord_Edge_CG[iDim] = edges->GetCG(iEdge,iDim); - Coord_Elem_CG[iDim] = elem[iElem]->GetCG(iDim); - Coord_FaceElem_CG[iDim] = elem[iElem]->GetFaceCG(iFace, iDim); - Coord_FaceiPoint[iDim] = nodes->GetCoord(face_iPoint, iDim); - Coord_FacejPoint[iDim] = nodes->GetCoord(face_jPoint, iDim); + Coord_Edge_CG[iDim] = 0.5 * (Coord[face_iNode][iDim] + Coord[face_jNode][iDim]); } /*--- Access the sub-volume of the element separately in 2D or 3D. ---*/ su2double Volume_i, Volume_j; - switch (nDim) { - case 2: - Volume_i = CEdge::GetVolume(Coord_FaceiPoint.data(), - Coord_Edge_CG.data(), - Coord_Elem_CG.data()); - - Volume_j = CEdge::GetVolume(Coord_FacejPoint.data(), - Coord_Edge_CG.data(), - Coord_Elem_CG.data()); - break; - case 3: - Volume_i = CEdge::GetVolume(Coord_FaceiPoint.data(), - Coord_Edge_CG.data(), - Coord_FaceElem_CG.data(), - Coord_Elem_CG.data()); - - Volume_j = CEdge::GetVolume(Coord_FacejPoint.data(), - Coord_Edge_CG.data(), - Coord_FaceElem_CG.data(), - Coord_Elem_CG.data()); - break; + if (nDim == 2) { + Volume_i = CEdge::GetVolume(Coord[face_iNode], Coord_Edge_CG, Coord_Elem_CG); + Volume_j = CEdge::GetVolume(Coord[face_jNode], Coord_Edge_CG, Coord_Elem_CG); + } + else { + Volume_i = CEdge::GetVolume(Coord[face_iNode], Coord_Edge_CG, + Coord_FaceElem_CG, Coord_Elem_CG); + Volume_j = CEdge::GetVolume(Coord[face_jNode], Coord_Edge_CG, + Coord_FaceElem_CG, Coord_Elem_CG); } /*--- Check if sub-elem volume is the min or max for iPoint. ---*/ @@ -8558,6 +8509,7 @@ void CPhysicalGeometry::ComputeMeshQualityStatistics(CConfig *config) { } } + } /*--- Compute the metrics with a final loop over the vertices. Also compute the local min and max values here for reporting. ---*/ diff --git a/Common/src/geometry/dual_grid/CEdge.cpp b/Common/src/geometry/dual_grid/CEdge.cpp index f2995644e025..c207ec075d64 100644 --- a/Common/src/geometry/dual_grid/CEdge.cpp +++ b/Common/src/geometry/dual_grid/CEdge.cpp @@ -36,7 +36,6 @@ CEdge::CEdge(unsigned long nEdge, unsigned long nDim) { const auto nEdgeSIMD = nextMultiple(nEdge, simd::preferredLen()); Nodes.resize(nEdgeSIMD,2) = 0; Normal.resize(nEdgeSIMD,nDim) = su2double(0.0); - Coord_CG.resize(nEdgeSIMD,nDim) = su2double(0.0); } void CEdge::SetZeroValues(void) { @@ -52,24 +51,13 @@ su2double CEdge::GetVolume(const su2double *coord_Edge_CG, su2double vec_a[nDim] = {0.0}, vec_b[nDim] = {0.0}, vec_c[nDim] = {0.0}, vec_d[nDim] = {0.0}; - AD::StartPreacc(); - AD::SetPreaccIn(coord_Edge_CG, nDim); - AD::SetPreaccIn(coord_Elem_CG, nDim); - AD::SetPreaccIn(coord_FaceElem_CG, nDim); - AD::SetPreaccIn(coord_Point, nDim); - Distance(nDim, coord_Edge_CG, coord_Point, vec_a); Distance(nDim, coord_FaceElem_CG, coord_Point, vec_b); Distance(nDim, coord_Elem_CG, coord_Point, vec_c); CrossProduct(vec_a, vec_b, vec_d); - su2double Local_Volume = fabs(DotProduct(nDim, vec_c, vec_d)) / 6.0; - - AD::SetPreaccOut(Local_Volume); - AD::EndPreacc(); - - return Local_Volume; + return fabs(DotProduct(nDim, vec_c, vec_d)) / 6.0; } su2double CEdge::GetVolume(const su2double *coord_Edge_CG, @@ -80,20 +68,10 @@ su2double CEdge::GetVolume(const su2double *coord_Edge_CG, su2double vec_a[nDim] = {0.0}, vec_b[nDim] = {0.0}; - AD::StartPreacc(); - AD::SetPreaccIn(coord_Edge_CG, nDim); - AD::SetPreaccIn(coord_Elem_CG, nDim); - AD::SetPreaccIn(coord_Point, nDim); - Distance(nDim, coord_Elem_CG, coord_Point, vec_a); Distance(nDim, coord_Edge_CG, coord_Point, vec_b); - su2double Local_Volume = 0.5 * fabs(vec_a[0]*vec_b[1] - vec_a[1]*vec_b[0]); - - AD::SetPreaccOut(Local_Volume); - AD::EndPreacc(); - - return Local_Volume; + return 0.5 * fabs(vec_a[0]*vec_b[1] - vec_a[1]*vec_b[0]); } void CEdge::SetNodes_Coord(unsigned long iEdge, @@ -105,12 +83,6 @@ void CEdge::SetNodes_Coord(unsigned long iEdge, su2double vec_a[nDim] = {0.0}, vec_b[nDim] = {0.0}, Dim_Normal[nDim]; - AD::StartPreacc(); - AD::SetPreaccIn(coord_Edge_CG, nDim); - AD::SetPreaccIn(coord_Elem_CG, nDim); - AD::SetPreaccIn(coord_FaceElem_CG, nDim); - AD::SetPreaccIn(Normal[iEdge], nDim); - Distance(nDim, coord_Elem_CG, coord_Edge_CG, vec_a); Distance(nDim, coord_FaceElem_CG, coord_Edge_CG, vec_b); @@ -118,26 +90,12 @@ void CEdge::SetNodes_Coord(unsigned long iEdge, for (auto iDim = 0ul; iDim < nDim; ++iDim) Normal(iEdge,iDim) += 0.5 * Dim_Normal[iDim]; - - AD::SetPreaccOut(Normal[iEdge], nDim); - AD::EndPreacc(); } void CEdge::SetNodes_Coord(unsigned long iEdge, const su2double *coord_Edge_CG, const su2double *coord_Elem_CG) { - constexpr unsigned long nDim = 2; - - AD::StartPreacc(); - AD::SetPreaccIn(coord_Elem_CG, nDim); - AD::SetPreaccIn(coord_Edge_CG, nDim); - AD::SetPreaccIn(Normal[iEdge], nDim); - Normal(iEdge,0) += coord_Elem_CG[1] - coord_Edge_CG[1]; Normal(iEdge,1) -= coord_Elem_CG[0] - coord_Edge_CG[0]; - - AD::SetPreaccOut(Normal[iEdge], nDim); - AD::EndPreacc(); - } diff --git a/Common/src/geometry/dual_grid/CVertex.cpp b/Common/src/geometry/dual_grid/CVertex.cpp index 3bce07b55ae9..06fb710f3767 100644 --- a/Common/src/geometry/dual_grid/CVertex.cpp +++ b/Common/src/geometry/dual_grid/CVertex.cpp @@ -26,50 +26,34 @@ */ #include "../../../include/geometry/dual_grid/CVertex.hpp" +#include "../../../include/toolboxes/geometry_toolbox.hpp" + +using namespace GeometryToolbox; CVertex::CVertex(unsigned long val_point, unsigned short val_nDim) : CDualGrid(val_nDim) { Nodes[0] = val_point; } -void CVertex::SetNodes_Coord(su2double *val_coord_Edge_CG, su2double *val_coord_FaceElem_CG, su2double *val_coord_Elem_CG) { - - su2double vec_a[3] = {0.0,0.0,0.0}, vec_b[3] = {0.0,0.0,0.0}; - unsigned short iDim; - - assert(nDim == 3); - - AD::StartPreacc(); - AD::SetPreaccIn(val_coord_Edge_CG, nDim); - AD::SetPreaccIn(val_coord_Elem_CG, nDim); - AD::SetPreaccIn(val_coord_FaceElem_CG, nDim); - AD::SetPreaccIn(Normal, nDim); +void CVertex::SetNodes_Coord(const su2double *coord_Edge_CG, + const su2double *coord_FaceElem_CG, + const su2double *coord_Elem_CG) { - for (iDim = 0; iDim < nDim; iDim++) { - vec_a[iDim] = val_coord_Elem_CG[iDim]-val_coord_Edge_CG[iDim]; - vec_b[iDim] = val_coord_FaceElem_CG[iDim]-val_coord_Edge_CG[iDim]; - } + constexpr unsigned long nDim = 3; + su2double vec_a[nDim] = {0.0}, vec_b[nDim] = {0.0}, Dim_Normal[nDim]; - Normal[0] += 0.5 * ( vec_a[1] * vec_b[2] - vec_a[2] * vec_b[1]); - Normal[1] -= 0.5 * ( vec_a[0] * vec_b[2] - vec_a[2] * vec_b[0]); - Normal[2] += 0.5 * ( vec_a[0] * vec_b[1] - vec_a[1] * vec_b[0]); + Distance(nDim, coord_Elem_CG, coord_Edge_CG, vec_a); + Distance(nDim, coord_FaceElem_CG, coord_Edge_CG, vec_b); - AD::SetPreaccOut(Normal, nDim); - AD::EndPreacc(); + CrossProduct(vec_a, vec_b, Dim_Normal); + for (auto iDim = 0ul; iDim < nDim; ++iDim) + Normal[iDim] += 0.5 * Dim_Normal[iDim]; } -void CVertex::SetNodes_Coord(su2double *val_coord_Edge_CG, su2double *val_coord_Elem_CG) { - - AD::StartPreacc(); - AD::SetPreaccIn(val_coord_Elem_CG, nDim); - AD::SetPreaccIn(val_coord_Edge_CG, nDim); - AD::SetPreaccIn(Normal, nDim); +void CVertex::SetNodes_Coord(const su2double *val_coord_Edge_CG, + const su2double *val_coord_Elem_CG) { Normal[0] += val_coord_Elem_CG[1]-val_coord_Edge_CG[1]; - Normal[1] -= (val_coord_Elem_CG[0]-val_coord_Edge_CG[0]); - - AD::SetPreaccOut(Normal, nDim); - AD::EndPreacc(); - + Normal[1] -= val_coord_Elem_CG[0]-val_coord_Edge_CG[0]; } diff --git a/Common/src/geometry/primal_grid/CHexahedron.cpp b/Common/src/geometry/primal_grid/CHexahedron.cpp index 3bb2a8a5c379..dd82981d9a87 100644 --- a/Common/src/geometry/primal_grid/CHexahedron.cpp +++ b/Common/src/geometry/primal_grid/CHexahedron.cpp @@ -49,19 +49,9 @@ CHexahedron::CHexahedron(unsigned long val_point_0, unsigned long val_point_1, unsigned long val_point_2, unsigned long val_point_3, unsigned long val_point_4, unsigned long val_point_5, unsigned long val_point_6, unsigned long val_point_7) : CPrimalGrid() { - unsigned short iDim, iFace, iNeighbor_Elements; + unsigned short iNeighbor_Elements; - /*--- Allocate center-of-gravity coordinates ---*/ nDim = 3; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_CG[iDim] = 0.0; - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; @@ -79,14 +69,7 @@ CHexahedron::CHexahedron(unsigned long val_point_0, unsigned long val_point_1, } -CHexahedron::~CHexahedron() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CHexahedron::~CHexahedron() {} void CHexahedron::Change_Orientation(void) { swap(Nodes[1], Nodes[3]); diff --git a/Common/src/geometry/primal_grid/CLine.cpp b/Common/src/geometry/primal_grid/CLine.cpp index 27d57798beed..3f903dae9462 100644 --- a/Common/src/geometry/primal_grid/CLine.cpp +++ b/Common/src/geometry/primal_grid/CLine.cpp @@ -48,20 +48,7 @@ unsigned short CLine::maxNodesFace = 2; CLine::CLine(unsigned long val_point_0, unsigned long val_point_1, unsigned short val_nDim) : CPrimalGrid() { - unsigned short iDim, iFace; - - /*--- Allocate CG coordinates ---*/ - nDim = val_nDim; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_CG[iDim] = 0.0; - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } /*--- Allocate and define face structure of the element ---*/ @@ -71,11 +58,4 @@ CLine::CLine(unsigned long val_point_0, unsigned long val_point_1, } -CLine::~CLine() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CLine::~CLine() {} diff --git a/Common/src/geometry/primal_grid/CPrimalGrid.cpp b/Common/src/geometry/primal_grid/CPrimalGrid.cpp index 5c08bdf2b989..7481ff27b47b 100644 --- a/Common/src/geometry/primal_grid/CPrimalGrid.cpp +++ b/Common/src/geometry/primal_grid/CPrimalGrid.cpp @@ -36,8 +36,6 @@ CPrimalGrid::CPrimalGrid(void) { Neighbor_Elements = nullptr; ElementOwnsFace = nullptr; PeriodIndexNeighbors = nullptr; - Coord_CG = nullptr; - Coord_FaceElems_CG = nullptr; JacobianFaceIsConstant = nullptr; GlobalIndex = 0; @@ -46,40 +44,12 @@ CPrimalGrid::CPrimalGrid(void) { CPrimalGrid::~CPrimalGrid() { delete[] Nodes; - delete[] Coord_CG; delete[] Neighbor_Elements; delete[] ElementOwnsFace; delete[] PeriodIndexNeighbors; delete[] JacobianFaceIsConstant; } -void CPrimalGrid::SetCoord_CG(const su2double* const* val_coord) { - unsigned short iDim, iNode, NodeFace, iFace; - - AD::StartPreacc(); - AD::SetPreaccIn(val_coord, GetnNodes(), nDim); - - for (iDim = 0; iDim < nDim; iDim++) { - Coord_CG[iDim] = 0.0; - for (iNode = 0; iNode < GetnNodes(); iNode++) - Coord_CG[iDim] += val_coord[iNode][iDim]/su2double(GetnNodes()); - } - - for (iFace = 0; iFace < GetnFaces(); iFace++) - for (iDim = 0; iDim < nDim; iDim++) { - Coord_FaceElems_CG[iFace][iDim] = 0.0; - for (iNode = 0; iNode < GetnNodesFace(iFace); iNode++) { - NodeFace = GetFaces(iFace, iNode); - Coord_FaceElems_CG[iFace][iDim] += val_coord[NodeFace][iDim]/su2double(GetnNodesFace(iFace)); - } - } - - AD::SetPreaccOut(Coord_CG, nDim); - AD::SetPreaccOut(Coord_FaceElems_CG, GetnFaces(), nDim); - AD::EndPreacc(); - -} - void CPrimalGrid::GetAllNeighbor_Elements() { cout << "( "; for (unsigned short iFace = 0; iFace < GetnFaces(); iFace++) diff --git a/Common/src/geometry/primal_grid/CPrism.cpp b/Common/src/geometry/primal_grid/CPrism.cpp index ca7a2e665a0b..6c858ef30669 100644 --- a/Common/src/geometry/primal_grid/CPrism.cpp +++ b/Common/src/geometry/primal_grid/CPrism.cpp @@ -48,19 +48,9 @@ unsigned short CPrism::maxNodesFace = 4; CPrism::CPrism(unsigned long val_point_0, unsigned long val_point_1, unsigned long val_point_2, unsigned long val_point_3, unsigned long val_point_4, unsigned long val_point_5) : CPrimalGrid() { - unsigned short iDim, iFace, iNeighbor_Elements; + unsigned short iNeighbor_Elements; - /*--- Allocate CG coordinates ---*/ nDim = 3; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) Coord_CG[iDim] = 0.0; - - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; @@ -80,14 +70,7 @@ CPrism::CPrism(unsigned long val_point_0, unsigned long val_point_1, } -CPrism::~CPrism() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CPrism::~CPrism() {} void CPrism::Change_Orientation(void) { swap(Nodes[0], Nodes[1]); diff --git a/Common/src/geometry/primal_grid/CPyramid.cpp b/Common/src/geometry/primal_grid/CPyramid.cpp index 02ead89dd266..eed37cfdf3d9 100644 --- a/Common/src/geometry/primal_grid/CPyramid.cpp +++ b/Common/src/geometry/primal_grid/CPyramid.cpp @@ -48,19 +48,9 @@ unsigned short CPyramid::maxNodesFace = 4; CPyramid::CPyramid(unsigned long val_point_0, unsigned long val_point_1, unsigned long val_point_2, unsigned long val_point_3, unsigned long val_point_4) : CPrimalGrid() { - unsigned short iDim, iFace, iNeighbor_Elements; + unsigned short iNeighbor_Elements; - /*--- Allocate CG coordinates ---*/ nDim = 3; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_CG[iDim] = 0.0; - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; @@ -79,14 +69,7 @@ CPyramid::CPyramid(unsigned long val_point_0, unsigned long val_point_1, } -CPyramid::~CPyramid() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CPyramid::~CPyramid() {} void CPyramid::Change_Orientation(void) { swap(Nodes[1],Nodes[3]); diff --git a/Common/src/geometry/primal_grid/CQuadrilateral.cpp b/Common/src/geometry/primal_grid/CQuadrilateral.cpp index 7fd629a63325..f732c9594363 100644 --- a/Common/src/geometry/primal_grid/CQuadrilateral.cpp +++ b/Common/src/geometry/primal_grid/CQuadrilateral.cpp @@ -48,19 +48,9 @@ unsigned short CQuadrilateral::maxNodesFace = 2; CQuadrilateral::CQuadrilateral(unsigned long val_point_0, unsigned long val_point_1, unsigned long val_point_2, unsigned long val_point_3, unsigned short val_nDim) : CPrimalGrid() { - unsigned short iDim, iFace, iNeighbor_Elements; + unsigned short iNeighbor_Elements; - /*--- Allocate CG coordinates ---*/ nDim = val_nDim; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_CG[iDim] = 0.0; - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; @@ -78,14 +68,7 @@ CQuadrilateral::CQuadrilateral(unsigned long val_point_0, unsigned long val_poin } -CQuadrilateral::~CQuadrilateral() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CQuadrilateral::~CQuadrilateral() {} void CQuadrilateral::Change_Orientation(void) { swap(Nodes[1], Nodes[3]); diff --git a/Common/src/geometry/primal_grid/CTetrahedron.cpp b/Common/src/geometry/primal_grid/CTetrahedron.cpp index bb99d5d0d47c..a747ac9d2ff0 100644 --- a/Common/src/geometry/primal_grid/CTetrahedron.cpp +++ b/Common/src/geometry/primal_grid/CTetrahedron.cpp @@ -47,19 +47,9 @@ unsigned short CTetrahedron::maxNodesFace = 3; CTetrahedron::CTetrahedron(unsigned long val_point_0, unsigned long val_point_1, unsigned long val_point_2, unsigned long val_point_3) : CPrimalGrid() { - unsigned short iDim, iFace, iNeighbor_Elements; + unsigned short iNeighbor_Elements; - /*--- Allocate CG coordinates ---*/ nDim = 3; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_CG[iDim] = 0.0; - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; @@ -77,14 +67,7 @@ CTetrahedron::CTetrahedron(unsigned long val_point_0, unsigned long val_point_1, } -CTetrahedron::~CTetrahedron() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CTetrahedron::~CTetrahedron() {} void CTetrahedron::Change_Orientation(void) { swap(Nodes[0],Nodes[1]); diff --git a/Common/src/geometry/primal_grid/CTriangle.cpp b/Common/src/geometry/primal_grid/CTriangle.cpp index 4a251d543169..7ddea102c19b 100644 --- a/Common/src/geometry/primal_grid/CTriangle.cpp +++ b/Common/src/geometry/primal_grid/CTriangle.cpp @@ -47,19 +47,10 @@ unsigned short CTriangle::maxNodesFace = 2; CTriangle::CTriangle(unsigned long val_point_0, unsigned long val_point_1, unsigned long val_point_2, unsigned short val_nDim) : CPrimalGrid() { - unsigned short iDim, iFace, iNeighbor_Elements; + unsigned short iNeighbor_Elements; - /*--- Allocate CG coordinates ---*/ nDim = val_nDim; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_CG[iDim] = 0.0; - Coord_FaceElems_CG = new su2double* [nFaces]; - for (iFace = 0; iFace < nFaces; iFace++) { - Coord_FaceElems_CG[iFace] = new su2double [nDim]; - for (iDim = 0; iDim < nDim; iDim++) - Coord_FaceElems_CG[iFace][iDim] = 0.0; - } + /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; Nodes[0] = val_point_0; @@ -75,11 +66,4 @@ CTriangle::CTriangle(unsigned long val_point_0, unsigned long val_point_1, } -CTriangle::~CTriangle() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CTriangle::~CTriangle() {} diff --git a/Common/src/geometry/primal_grid/CVertexMPI.cpp b/Common/src/geometry/primal_grid/CVertexMPI.cpp index a97331295de9..6401d60dcb19 100644 --- a/Common/src/geometry/primal_grid/CVertexMPI.cpp +++ b/Common/src/geometry/primal_grid/CVertexMPI.cpp @@ -38,12 +38,8 @@ unsigned short CVertexMPI::VTK_Type = 1; unsigned short CVertexMPI::maxNodesFace = 0; CVertexMPI::CVertexMPI(unsigned long val_point, unsigned short val_nDim) : CPrimalGrid() { - unsigned short iDim; - /*--- Allocate CG coordinates ---*/ nDim = val_nDim; - Coord_CG = new su2double[nDim]; - for (iDim = 0; iDim < nDim; iDim++) Coord_CG[iDim] = 0.0; /*--- Allocate and define face structure of the element ---*/ Nodes = new unsigned long[nNodes]; @@ -54,13 +50,6 @@ CVertexMPI::CVertexMPI(unsigned long val_point, unsigned short val_nDim) : CPrim } -CVertexMPI::~CVertexMPI() { - unsigned short iFaces; - - for (iFaces = 0; iFaces < nFaces; iFaces++) - if (Coord_FaceElems_CG[iFaces] != nullptr) delete[] Coord_FaceElems_CG[iFaces]; - delete[] Coord_FaceElems_CG; - -} +CVertexMPI::~CVertexMPI() {} void CVertexMPI::Change_Orientation(void) { } diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp index 0aaff36590bd..68b35eb613a1 100644 --- a/Common/src/grid_movement/CVolumetricMovement.cpp +++ b/Common/src/grid_movement/CVolumetricMovement.cpp @@ -90,7 +90,6 @@ void CVolumetricMovement::UpdateDualGrid(CGeometry *geometry, CConfig *config) { /*--- After moving all nodes, update the dual mesh. Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/ - geometry->SetCoord_CG(); geometry->SetControlVolume(config, UPDATE); geometry->SetBoundControlVolume(config, UPDATE); geometry->SetMaxLength(config); diff --git a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl index 84e11f019e9a..26b4bb7850b7 100644 --- a/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl +++ b/SU2_CFD/include/solvers/CFVMFlowSolverBase.inl @@ -495,6 +495,7 @@ CFVMFlowSolverBase::~CFVMFlowSolverBase() { } delete nodes; + delete edgeNumerics; } template diff --git a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp index 06172348f96d..61516cce31c2 100644 --- a/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjMultizoneDriver.cpp @@ -564,11 +564,20 @@ void CDiscAdjMultizoneDriver::SetRecording(unsigned short kind_recording, Kind_T } } - if (rank == MASTER_NODE) { - if(kind_recording != NONE && config_container[record_zone]->GetWrt_AD_Statistics()) { - AD::PrintStatistics(); + if (kind_recording != NONE && driver_config->GetWrt_AD_Statistics()) { + if (rank == MASTER_NODE) AD::PrintStatistics(); +#ifdef CODI_REVERSE_TYPE + if (size > SINGLE_NODE) { + su2double myMem = AD::globalTape.getTapeValues().getUsedMemorySize(), totMem = 0.0; + SU2_MPI::Allreduce(&myMem, &totMem, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + if (rank == MASTER_NODE) { + cout << "MPI\n"; + cout << "-------------------------------------\n"; + cout << " Total memory used : " << totMem << " MB\n"; + cout << "-------------------------------------\n" << endl; + } } - cout << "-------------------------------------------------------------------------\n" << endl; +#endif } AD::StopRecording(); diff --git a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp index 495876b2c6af..b1561b048f94 100644 --- a/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp +++ b/SU2_CFD/src/drivers/CDiscAdjSinglezoneDriver.cpp @@ -291,9 +291,20 @@ void CDiscAdjSinglezoneDriver::SetRecording(unsigned short kind_recording){ SetObjFunction(); - if (rank == MASTER_NODE && kind_recording != NONE && config_container[ZONE_0]->GetWrt_AD_Statistics()) { - AD::PrintStatistics(); - cout << "-------------------------------------------------------------------------\n" << endl; + if (kind_recording != NONE && config_container[ZONE_0]->GetWrt_AD_Statistics()) { + if (rank == MASTER_NODE) AD::PrintStatistics(); +#ifdef CODI_REVERSE_TYPE + if (size > SINGLE_NODE) { + su2double myMem = AD::globalTape.getTapeValues().getUsedMemorySize(), totMem = 0.0; + SU2_MPI::Allreduce(&myMem, &totMem, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD); + if (rank == MASTER_NODE) { + cout << "MPI\n"; + cout << "-------------------------------------\n"; + cout << " Total memory used : " << totMem << " MB\n"; + cout << "-------------------------------------\n" << endl; + } + } +#endif } AD::StopRecording(); diff --git a/SU2_CFD/src/drivers/CDriver.cpp b/SU2_CFD/src/drivers/CDriver.cpp index 2a3448f1b256..841a0d6cb822 100644 --- a/SU2_CFD/src/drivers/CDriver.cpp +++ b/SU2_CFD/src/drivers/CDriver.cpp @@ -804,24 +804,19 @@ void CDriver::Geometrical_Preprocessing_FVM(CConfig *config, CGeometry **&geomet geometry[MESH_0]->SetEdges(); geometry[MESH_0]->SetVertex(config); - /*--- Compute cell center of gravity ---*/ - - if ((rank == MASTER_NODE) && (!fea)) cout << "Computing centers of gravity." << endl; - SU2_OMP_PARALLEL { - geometry[MESH_0]->SetCoord_CG(); - } - /*--- Create the control volume structures ---*/ if ((rank == MASTER_NODE) && (!fea)) cout << "Setting the control volume structure." << endl; - geometry[MESH_0]->SetControlVolume(config, ALLOCATE); - geometry[MESH_0]->SetBoundControlVolume(config, ALLOCATE); + SU2_OMP_PARALLEL { + geometry[MESH_0]->SetControlVolume(config, ALLOCATE); + geometry[MESH_0]->SetBoundControlVolume(config, ALLOCATE); + } /*--- Visualize a dual control volume if requested ---*/ if ((config->GetVisualize_CV() >= 0) && - (config->GetVisualize_CV() < (long)geometry[MESH_0]->GetnPointDomain())) - geometry[MESH_0]->VisualizeControlVolume(config, UPDATE); + (config->GetVisualize_CV() < (long)geometry[MESH_0]->GetGlobal_nPointDomain())) + geometry[MESH_0]->VisualizeControlVolume(config); /*--- Identify closest normal neighbor ---*/ diff --git a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp index 5b49e9f64eec..d889f036b4ce 100644 --- a/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp +++ b/SU2_CFD/src/iteration/CDiscAdjFEAIteration.cpp @@ -385,7 +385,7 @@ void CDiscAdjFEAIteration::SetDependencies(CSolver***** solver, CGeometry**** ge /*--- We only differentiate wrt to this variable in the adjoint secondary recording. ---*/ if (config[iZone]->GetTopology_Optimization() && (kind_recording == MESH_COORDS)) { /*--- The filter may require the volumes of the elements. ---*/ - structural_geometry->SetElemVolume(config[iZone]); + structural_geometry->SetElemVolume(); /// TODO: Ideally there would be a way to capture this dependency without the `static_cast`, but /// making it a virtual method of CSolver does not feel "right" as its purpose could be confused. static_cast(dir_solver)->FilterElementDensities(structural_geometry, config[iZone]); diff --git a/SU2_CFD/src/solvers/CEulerSolver.cpp b/SU2_CFD/src/solvers/CEulerSolver.cpp index 792b80b96eff..d05bbfc0bc3a 100644 --- a/SU2_CFD/src/solvers/CEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CEulerSolver.cpp @@ -2171,8 +2171,7 @@ void CEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_con bool cont_adjoint = config->GetContinuous_Adjoint(); bool disc_adjoint = config->GetDiscrete_Adjoint(); bool implicit = (config->GetKind_TimeIntScheme() == EULER_IMPLICIT); - bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || - (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); + bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); bool center_jst = (config->GetKind_Centered_Flow() == JST) && (iMesh == MESH_0); bool center_jst_ke = (config->GetKind_Centered_Flow() == JST_KE) && (iMesh == MESH_0); bool center_jst_mat = (config->GetKind_Centered_Flow() == JST_MAT) && (iMesh == MESH_0); @@ -2275,12 +2274,11 @@ void CEulerSolver::CommonPreprocessing(CGeometry *geometry, CSolver **solver_con void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - unsigned long InnerIter = config->GetInnerIter(); - bool cont_adjoint = config->GetContinuous_Adjoint(); - bool muscl = (config->GetMUSCL_Flow() || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == ROE)); - bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); - bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED) || (cont_adjoint && config->GetKind_ConvNumScheme_AdjFlow() == SPACE_CENTERED); - bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; + const auto InnerIter = config->GetInnerIter(); + const bool muscl = config->GetMUSCL_Flow() && (iMesh == MESH_0); + const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); + const bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); + const bool van_albada = (config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE); /*--- Common preprocessing steps. ---*/ @@ -2288,7 +2286,7 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container /*--- Upwind second order reconstruction ---*/ - if ((muscl && !center) && (iMesh == MESH_0) && !Output) { + if (!Output && muscl && !center) { /*--- Gradient computation for MUSCL reconstruction. ---*/ @@ -2303,10 +2301,8 @@ void CEulerSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container /*--- Limiter computation ---*/ - if (limiter && (iMesh == MESH_0) && !Output && !van_albada) - SetPrimitive_Limiter(geometry, config); + if (limiter && !van_albada) SetPrimitive_Limiter(geometry, config); } - } unsigned long CEulerSolver::SetPrimitive_Variables(CSolver **solver_container, CConfig *config, bool Output) { @@ -10112,7 +10108,6 @@ void CEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConfig /*--- Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/ - geometry[MESH_0]->SetCoord_CG(); geometry[MESH_0]->SetControlVolume(config, UPDATE); geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry[MESH_0]->SetMaxLength(config); diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp index cdc7f9ff6a1d..fe10390bf57c 100644 --- a/SU2_CFD/src/solvers/CFEASolver.cpp +++ b/SU2_CFD/src/solvers/CFEASolver.cpp @@ -782,7 +782,7 @@ void CFEASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, * This only needs to be done for the undeformed (initial) shape. */ if (topology_mode && !topol_filter_applied) { - geometry->SetElemVolume(config); + geometry->SetElemVolume(); FilterElementDensities(geometry, config); topol_filter_applied = true; } diff --git a/SU2_CFD/src/solvers/CIncEulerSolver.cpp b/SU2_CFD/src/solvers/CIncEulerSolver.cpp index 7891a1c35c6e..50c62b0635b6 100644 --- a/SU2_CFD/src/solvers/CIncEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CIncEulerSolver.cpp @@ -3661,7 +3661,6 @@ void CIncEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CConf /*--- Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/ - geometry[MESH_0]->SetCoord_CG(); geometry[MESH_0]->SetControlVolume(config, UPDATE); geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry[MESH_0]->SetMaxLength(config); diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp index 01649dfffaae..2c2bd1ad4f50 100644 --- a/SU2_CFD/src/solvers/CMeshSolver.cpp +++ b/SU2_CFD/src/solvers/CMeshSolver.cpp @@ -558,7 +558,6 @@ void CMeshSolver::UpdateDualGrid(CGeometry *geometry, CConfig *config){ /*--- After moving all nodes, update the dual mesh. Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/ - geometry->SetCoord_CG(); geometry->SetControlVolume(config, UPDATE); geometry->SetBoundControlVolume(config, UPDATE); geometry->SetMaxLength(config); diff --git a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp index 9345cfeab89f..40e4b6dfae41 100644 --- a/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp +++ b/SU2_CFD/src/solvers/CNEMOEulerSolver.cpp @@ -3147,7 +3147,6 @@ void CNEMOEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CCon /*--- Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/ - geometry[MESH_0]->SetCoord_CG(); geometry[MESH_0]->SetControlVolume(config, UPDATE); geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry[MESH_0]->SetMaxLength(config); @@ -3173,7 +3172,6 @@ void CNEMOEulerSolver::LoadRestart(CGeometry **geometry, CSolver ***solver, CCon /*--- Recompute the edges and dual mesh control volumes in the domain and on the boundaries. ---*/ - geometry[MESH_0]->SetCoord_CG(); geometry[MESH_0]->SetControlVolume(config, UPDATE); geometry[MESH_0]->SetBoundControlVolume(config, UPDATE); geometry[MESH_0]->SetMaxLength(config); diff --git a/SU2_CFD/src/solvers/CNSSolver.cpp b/SU2_CFD/src/solvers/CNSSolver.cpp index f24e21fd4b47..c82ff6ebf6ff 100644 --- a/SU2_CFD/src/solvers/CNSSolver.cpp +++ b/SU2_CFD/src/solvers/CNSSolver.cpp @@ -98,21 +98,29 @@ CNSSolver::~CNSSolver(void) { void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - unsigned long InnerIter = config->GetInnerIter(); - bool cont_adjoint = config->GetContinuous_Adjoint(); - bool limiter_flow = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); - bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); - bool limiter_adjflow = (cont_adjoint && (config->GetKind_SlopeLimit_AdjFlow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter())); - bool van_albada = config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE; - bool wall_functions = config->GetWall_Functions(); + const auto InnerIter = config->GetInnerIter(); + const bool muscl = config->GetMUSCL_Flow() && (iMesh == MESH_0); + const bool center = (config->GetKind_ConvNumScheme_Flow() == SPACE_CENTERED); + const bool limiter = (config->GetKind_SlopeLimit_Flow() != NO_LIMITER) && (InnerIter <= config->GetLimiterIter()); + const bool van_albada = (config->GetKind_SlopeLimit_Flow() == VAN_ALBADA_EDGE); + const bool wall_functions = config->GetWall_Functions(); /*--- Common preprocessing steps (implemented by CEulerSolver) ---*/ CommonPreprocessing(geometry, solver_container, config, iMesh, iRKStep, RunTime_EqSystem, Output); - /*--- Compute gradient for MUSCL reconstruction. ---*/ + /*--- Compute gradient for MUSCL reconstruction, for output (i.e. the + turbulence solver) only density and velocity are needed ---*/ - if (config->GetReconstructionGradientRequired() && (iMesh == MESH_0)) { + const auto nPrimVarGrad_bak = nPrimVarGrad; + if (Output) { + SU2_OMP_BARRIER + SU2_OMP_MASTER + nPrimVarGrad = 1+nDim; + SU2_OMP_BARRIER + } + + if (config->GetReconstructionGradientRequired() && muscl && !center) { switch (config->GetKind_Gradient_Method_Recon()) { case GREEN_GAUSS: SetPrimitive_Gradient_GG(geometry, config, true); break; @@ -132,10 +140,15 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C SetPrimitive_Gradient_LS(geometry, config); } - /*--- Compute the limiter in case we need it in the turbulence model or to limit the - * viscous terms (check this logic with JST and 2nd order turbulence model) ---*/ + if (Output) { + SU2_OMP_MASTER + nPrimVarGrad = nPrimVarGrad_bak; + SU2_OMP_BARRIER + } - if ((iMesh == MESH_0) && (limiter_flow || limiter_turb || limiter_adjflow) && !Output && !van_albada) { + /*--- Compute the limiters ---*/ + + if (muscl && !center && limiter && !van_albada && !Output) { SetPrimitive_Limiter(geometry, config); } @@ -150,21 +163,17 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C nodes->SetVorticity_StrainMag(); + /*--- Min and Max are not really differentiable ---*/ + const bool wasActive = AD::BeginPassive(); + su2double strainMax = 0.0, omegaMax = 0.0; SU2_OMP(for schedule(static,omp_chunk_size) nowait) for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++) { - - su2double StrainMag = nodes->GetStrainMag(iPoint); - const su2double* Vorticity = nodes->GetVorticity(iPoint); - su2double Omega = sqrt(Vorticity[0]*Vorticity[0]+ Vorticity[1]*Vorticity[1]+ Vorticity[2]*Vorticity[2]); - - strainMax = max(strainMax, StrainMag); - omegaMax = max(omegaMax, Omega); - + strainMax = max(strainMax, nodes->GetStrainMag(iPoint)); + omegaMax = max(omegaMax, GeometryToolbox::Norm(3, nodes->GetVorticity(iPoint))); } - SU2_OMP_CRITICAL - { + SU2_OMP_CRITICAL { StrainMag_Max = max(StrainMag_Max, strainMax); Omega_Max = max(Omega_Max, omegaMax); } @@ -182,6 +191,8 @@ void CNSSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, C SU2_OMP_BARRIER } + AD::EndPassive(wasActive); + /*--- Compute the TauWall from the wall functions ---*/ if (wall_functions) { diff --git a/SU2_CFD/src/solvers/CTurbSASolver.cpp b/SU2_CFD/src/solvers/CTurbSASolver.cpp index 59bb5ddc4e7c..4209d0a11865 100644 --- a/SU2_CFD/src/solvers/CTurbSASolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSASolver.cpp @@ -243,12 +243,10 @@ CTurbSASolver::~CTurbSASolver(void) { void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && - (config->GetInnerIter() <= config->GetLimiterIter()); - unsigned short kind_hybridRANSLES = config->GetKind_HybridRANSLES(); - const su2double* const* PrimGrad_Flow = nullptr; - const su2double* Vorticity = nullptr; - su2double Laminar_Viscosity = 0.0; + const bool muscl = config->GetMUSCL_Turb(); + const bool limiter = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && + (config->GetInnerIter() <= config->GetLimiterIter()); + const auto kind_hybridRANSLES = config->GetKind_HybridRANSLES(); /*--- Clear residual and system matrix, not needed for * reducer strategy as we write over the entire matrix. ---*/ @@ -259,7 +257,7 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe /*--- Upwind second order reconstruction and gradients ---*/ - if (config->GetReconstructionGradientRequired()) { + if (config->GetReconstructionGradientRequired() && muscl) { if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) SetSolution_Gradient_GG(geometry, config, true); if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) @@ -274,7 +272,7 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) SetSolution_Gradient_LS(geometry, config); - if (limiter_turb) SetSolution_Limiter(geometry, config); + if (limiter && muscl) SetSolution_Limiter(geometry, config); if (kind_hybridRANSLES != NO_HYBRIDRANSLES) { @@ -283,9 +281,9 @@ void CTurbSASolver::Preprocessing(CGeometry *geometry, CSolver **solver_containe if (kind_hybridRANSLES == SA_EDDES){ SU2_OMP_FOR_STAT(omp_chunk_size) for (unsigned long iPoint = 0; iPoint < nPoint; iPoint++){ - Vorticity = solver_container[FLOW_SOL]->GetNodes()->GetVorticity(iPoint); - PrimGrad_Flow = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); - Laminar_Viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); + auto Vorticity = solver_container[FLOW_SOL]->GetNodes()->GetVorticity(iPoint); + auto PrimGrad_Flow = solver_container[FLOW_SOL]->GetNodes()->GetGradient_Primitive(iPoint); + auto Laminar_Viscosity = solver_container[FLOW_SOL]->GetNodes()->GetLaminarViscosity(iPoint); nodes->SetVortex_Tilting(iPoint, PrimGrad_Flow, Vorticity, Laminar_Viscosity); } } diff --git a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp index e9bdb40207d4..335f7539b7a3 100644 --- a/SU2_CFD/src/solvers/CTurbSSTSolver.cpp +++ b/SU2_CFD/src/solvers/CTurbSSTSolver.cpp @@ -242,8 +242,9 @@ CTurbSSTSolver::~CTurbSSTSolver(void) { void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_container, CConfig *config, unsigned short iMesh, unsigned short iRKStep, unsigned short RunTime_EqSystem, bool Output) { - const bool limiter_turb = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && - (config->GetInnerIter() <= config->GetLimiterIter()); + const bool muscl = config->GetMUSCL_Turb(); + const bool limiter = (config->GetKind_SlopeLimit_Turb() != NO_LIMITER) && + (config->GetInnerIter() <= config->GetLimiterIter()); /*--- Clear residual and system matrix, not needed for * reducer strategy as we write over the entire matrix. ---*/ @@ -254,7 +255,7 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain /*--- Upwind second order reconstruction and gradients ---*/ - if (config->GetReconstructionGradientRequired()) { + if (config->GetReconstructionGradientRequired() && muscl) { if (config->GetKind_Gradient_Method_Recon() == GREEN_GAUSS) SetSolution_Gradient_GG(geometry, config, true); if (config->GetKind_Gradient_Method_Recon() == LEAST_SQUARES) @@ -269,7 +270,7 @@ void CTurbSSTSolver::Preprocessing(CGeometry *geometry, CSolver **solver_contain if (config->GetKind_Gradient_Method() == WEIGHTED_LEAST_SQUARES) SetSolution_Gradient_LS(geometry, config); - if (limiter_turb) SetSolution_Limiter(geometry, config); + if (limiter && muscl) SetSolution_Limiter(geometry, config); } diff --git a/SU2_DEF/src/SU2_DEF.cpp b/SU2_DEF/src/SU2_DEF.cpp index c300e7bb10b9..225a74d42221 100644 --- a/SU2_DEF/src/SU2_DEF.cpp +++ b/SU2_DEF/src/SU2_DEF.cpp @@ -195,11 +195,6 @@ int main(int argc, char *argv[]) { if (config_container[iZone]->GetDesign_Variable(0) != NO_DEFORMATION) { - /*--- Compute center of gravity ---*/ - - if (rank == MASTER_NODE) cout << "Computing centers of gravity." << endl; - geometry_container[iZone]->SetCoord_CG(); - /*--- Create the dual control volume structures ---*/ if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; @@ -394,7 +389,6 @@ int main(int argc, char *argv[]) { geometry_container[iZone]->SetBoundVolume(); geometry_container[iZone]->SetEdges(); geometry_container[iZone]->SetVertex(config_container[iZone]); - geometry_container[iZone]->SetCoord_CG(); geometry_container[iZone]->SetControlVolume(config_container[iZone], ALLOCATE); geometry_container[iZone]->SetBoundControlVolume(config_container[iZone], ALLOCATE); diff --git a/SU2_DOT/src/SU2_DOT.cpp b/SU2_DOT/src/SU2_DOT.cpp index 934471b5bfa3..4960c2a40716 100644 --- a/SU2_DOT/src/SU2_DOT.cpp +++ b/SU2_DOT/src/SU2_DOT.cpp @@ -241,11 +241,6 @@ int main(int argc, char *argv[]) { if (rank == MASTER_NODE) cout << "Identify edges and vertices." << endl; geometry_container[iZone][INST_0]->SetEdges(); geometry_container[iZone][INST_0]->SetVertex(config_container[iZone]); - /*--- Compute center of gravity ---*/ - - if (rank == MASTER_NODE) cout << "Computing centers of gravity." << endl; - geometry_container[iZone][INST_0]->SetCoord_CG(); - /*--- Create the dual control volume structures ---*/ if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; diff --git a/SU2_GEO/src/SU2_GEO.cpp b/SU2_GEO/src/SU2_GEO.cpp index 3b9ea6c6f7f8..8ac58865cb9e 100644 --- a/SU2_GEO/src/SU2_GEO.cpp +++ b/SU2_GEO/src/SU2_GEO.cpp @@ -6,7 +6,7 @@ * * SU2 Project Website: https://su2code.github.io * - * The SU2 Project is maintained by the SU2 Foundation + * The SU2 Project is maintained by the SU2 Foundation * (http://su2foundation.org) * * Copyright 2012-2020, SU2 Contributors (cf. AUTHORS.md) @@ -30,21 +30,21 @@ using namespace std; int main(int argc, char *argv[]) { - + unsigned short iZone, nZone = SINGLE_ZONE; su2double StartTime = 0.0, StopTime = 0.0, UsedTime = 0.0; unsigned short iDV, iFFDBox, iPlane, nPlane, iVar; su2double *ObjectiveFunc, *ObjectiveFunc_New, *Gradient, delta_eps, **Plane_P0, **Plane_Normal, - + Fuselage_Volume = 0.0, Fuselage_WettedArea = 0.0, Fuselage_MinWidth = 0.0, Fuselage_MaxWidth = 0.0, Fuselage_MinWaterLineWidth = 0.0, Fuselage_MaxWaterLineWidth = 0.0, Fuselage_MinHeight = 0.0, Fuselage_MaxHeight = 0.0, Fuselage_MaxCurvature = 0.0, Fuselage_Volume_New = 0.0, Fuselage_WettedArea_New = 0.0, Fuselage_MinWidth_New = 0.0, Fuselage_MaxWidth_New = 0.0, Fuselage_MinWaterLineWidth_New = 0.0, Fuselage_MaxWaterLineWidth_New = 0.0, Fuselage_MinHeight_New = 0.0, Fuselage_MaxHeight_New = 0.0, Fuselage_MaxCurvature_New = 0.0, Fuselage_Volume_Grad = 0.0, Fuselage_WettedArea_Grad = 0.0, Fuselage_MinWidth_Grad = 0.0, Fuselage_MaxWidth_Grad = 0.0, Fuselage_MinWaterLineWidth_Grad = 0.0, Fuselage_MaxWaterLineWidth_Grad = 0.0, Fuselage_MinHeight_Grad = 0.0, Fuselage_MaxHeight_Grad = 0.0, Fuselage_MaxCurvature_Grad = 0.0, - + Wing_Volume = 0.0, Wing_MinThickness = 0.0, Wing_MaxThickness = 0.0, Wing_MinChord = 0.0, Wing_MaxChord = 0.0, Wing_MinLERadius = 0.0, Wing_MaxLERadius = 0.0, Wing_MinToC = 0.0, Wing_MaxToC = 0.0, Wing_ObjFun_MinToC = 0.0, Wing_MaxTwist = 0.0, Wing_MaxCurvature = 0.0, Wing_MaxDihedral = 0.0, Wing_Volume_New = 0.0, Wing_MinThickness_New = 0.0, Wing_MaxThickness_New = 0.0, Wing_MinChord_New = 0.0, Wing_MaxChord_New = 0.0, Wing_MinLERadius_New = 0.0, Wing_MaxLERadius_New = 0.0, Wing_MinToC_New = 0.0, Wing_MaxToC_New = 0.0, Wing_ObjFun_MinToC_New = 0.0, Wing_MaxTwist_New = 0.0, Wing_MaxCurvature_New = 0.0, Wing_MaxDihedral_New = 0.0, Wing_Volume_Grad = 0.0, Wing_MinThickness_Grad = 0.0, Wing_MaxThickness_Grad = 0.0, Wing_MinChord_Grad = 0.0, Wing_MaxChord_Grad = 0.0, Wing_MinLERadius_Grad = 0.0, Wing_MaxLERadius_Grad = 0.0, Wing_MinToC_Grad = 0.0, Wing_MaxToC_Grad = 0.0, Wing_ObjFun_MinToC_Grad = 0.0, Wing_MaxTwist_Grad = 0.0, Wing_MaxCurvature_Grad = 0.0, Wing_MaxDihedral_Grad = 0.0, - + Nacelle_Volume = 0.0, Nacelle_MinThickness = 0.0, Nacelle_MaxThickness = 0.0, Nacelle_MinChord = 0.0, Nacelle_MaxChord = 0.0, Nacelle_MinLERadius = 0.0, Nacelle_MaxLERadius = 0.0, Nacelle_MinToC = 0.0, Nacelle_MaxToC = 0.0, Nacelle_ObjFun_MinToC = 0.0, Nacelle_MaxTwist = 0.0, Nacelle_Volume_New = 0.0, Nacelle_MinThickness_New = 0.0, Nacelle_MaxThickness_New = 0.0, Nacelle_MinChord_New = 0.0, Nacelle_MaxChord_New = 0.0, Nacelle_MinLERadius_New = 0.0, Nacelle_MaxLERadius_New = 0.0, Nacelle_MinToC_New = 0.0, Nacelle_MaxToC_New = 0.0, Nacelle_ObjFun_MinToC_New = 0.0, Nacelle_MaxTwist_New = 0.0, Nacelle_Volume_Grad = 0.0, Nacelle_MinThickness_Grad = 0.0, Nacelle_MaxThickness_Grad = 0.0, Nacelle_MinChord_Grad = 0.0, Nacelle_MaxChord_Grad = 0.0, Nacelle_MinLERadius_Grad = 0.0, Nacelle_MaxLERadius_Grad = 0.0, Nacelle_MinToC_Grad = 0.0, Nacelle_MaxToC_Grad = 0.0, Nacelle_ObjFun_MinToC_Grad = 0.0, Nacelle_MaxTwist_Grad = 0.0; @@ -55,9 +55,9 @@ int main(int argc, char *argv[]) { bool Local_MoveSurface, MoveSurface = false; ofstream Gradient_file, ObjFunc_file; int rank, size; - + /*--- MPI initialization ---*/ - + #ifdef HAVE_MPI SU2_MPI::Init(&argc,&argv); SU2_MPI::Comm MPICommunicator(MPI_COMM_WORLD); @@ -67,17 +67,17 @@ int main(int argc, char *argv[]) { rank = SU2_MPI::GetRank(); size = SU2_MPI::GetSize(); - + /*--- Pointer to different structures that will be used throughout the entire code ---*/ - + CConfig **config_container = nullptr; CGeometry **geometry_container = nullptr; CSurfaceMovement *surface_movement = nullptr; CFreeFormDefBox** FFDBox = nullptr; - + /*--- Load in the number of zones and spatial dimensions in the mesh file (if no config file is specified, default.cfg is used) ---*/ - + if (argc == 2) { strcpy(config_file_name,argv[1]); } else { strcpy(config_file_name, "default.cfg"); } @@ -91,141 +91,136 @@ int main(int argc, char *argv[]) { nZone = config->GetnZone(); /*--- Definition of the containers per zones ---*/ - + config_container = new CConfig*[nZone]; geometry_container = new CGeometry*[nZone]; - + for (iZone = 0; iZone < nZone; iZone++) { config_container[iZone] = nullptr; geometry_container[iZone] = nullptr; } - + /*--- Loop over all zones to initialize the various classes. In most cases, nZone is equal to one. This represents the solution of a partial differential equation on a single block, unstructured mesh. ---*/ - + for (iZone = 0; iZone < nZone; iZone++) { - + /*--- Definition of the configuration option class for all zones. In this constructor, the input configuration file is parsed and all options are read and stored. ---*/ - + config_container[iZone] = new CConfig(config_file_name, SU2_GEO, true); config_container[iZone]->SetMPICommunicator(MPICommunicator); - + /*--- Definition of the geometry class to store the primal grid in the partitioning process. ---*/ - + CGeometry *geometry_aux = nullptr; - + /*--- All ranks process the grid and call ParMETIS for partitioning ---*/ - + geometry_aux = new CPhysicalGeometry(config_container[iZone], iZone, nZone); - + /*--- Color the initial grid and set the send-receive domains (ParMETIS) ---*/ - + geometry_aux->SetColorGrid_Parallel(config_container[iZone]); - + /*--- Build the grid data structures using the ParMETIS coloring. ---*/ geometry_container[iZone] = new CPhysicalGeometry(geometry_aux, config_container[iZone]); - + /*--- Deallocate the memory of geometry_aux ---*/ - + delete geometry_aux; - + /*--- Add the Send/Receive boundaries ---*/ - + geometry_container[iZone]->SetSendReceive(config_container[iZone]); - + /*--- Add the Send/Receive boundaries ---*/ - + geometry_container[iZone]->SetBoundaries(config_container[iZone]); - + } - + bool tabTecplot = config_container[ZONE_0]->GetTabular_FileFormat() == TAB_TECPLOT; - + /*--- Set up a timer for performance benchmarking (preprocessing time is included) ---*/ StartTime = SU2_MPI::Wtime(); /*--- Evaluation of the objective function ---*/ - + if (rank == MASTER_NODE) cout << endl <<"----------------------- Preprocessing computations ----------------------" << endl; - + /*--- Set the number of sections, and allocate the memory ---*/ - + if (geometry_container[ZONE_0]->GetnDim() == 2) nPlane = 1; else nPlane = config_container[ZONE_0]->GetnLocationStations(); - + Xcoord_Airfoil = new vector[nPlane]; Ycoord_Airfoil = new vector[nPlane]; Zcoord_Airfoil = new vector[nPlane]; Variable_Airfoil = new vector[nPlane]; - + Plane_P0 = new su2double*[nPlane]; Plane_Normal = new su2double*[nPlane]; for(iPlane = 0; iPlane < nPlane; iPlane++ ) { Plane_P0[iPlane] = new su2double[3]; Plane_Normal[iPlane] = new su2double[3]; } - + ObjectiveFunc = new su2double[nPlane*20]; ObjectiveFunc_New = new su2double[nPlane*20]; Gradient = new su2double[nPlane*20]; - + for (iVar = 0; iVar < nPlane*20; iVar++) { ObjectiveFunc[iVar] = 0.0; ObjectiveFunc_New[iVar] = 0.0; Gradient[iVar] = 0.0; } - - + + /*--- Compute elements surrounding points, points surrounding points ---*/ - + if (rank == MASTER_NODE) cout << "Setting local point connectivity." <SetPoint_Connectivity(); - + /*--- Check the orientation before computing geometrical quantities ---*/ - + if (config_container[ZONE_0]->GetReorientElements()) { if (rank == MASTER_NODE) cout << "Checking the numerical grid orientation of the interior elements." <Check_IntElem_Orientation(config_container[ZONE_0]); } - + /*--- Create the edge structure ---*/ - + if (rank == MASTER_NODE) cout << "Identify edges and vertices." <SetEdges(); geometry_container[ZONE_0]->SetVertex(config_container[ZONE_0]); - - /*--- Compute center of gravity ---*/ - - if (rank == MASTER_NODE) cout << "Computing centers of gravity." << endl; - geometry_container[ZONE_0]->SetCoord_CG(); - + /*--- Create the dual control volume structures ---*/ - + if (rank == MASTER_NODE) cout << "Setting the bound control volume structure." << endl; geometry_container[ZONE_0]->SetBoundControlVolume(config_container[ZONE_0], ALLOCATE); - + /*--- Compute the surface curvature ---*/ - + if (rank == MASTER_NODE) cout << "Compute the surface curvature." << endl; geometry_container[ZONE_0]->ComputeSurf_Curvature(config_container[ZONE_0]); - + /*--- Computation of positive surface area in the z-plane ---*/ - + if (rank == MASTER_NODE) cout << "Setting reference area and span." << endl; geometry_container[ZONE_0]->SetPositive_ZArea(config_container[ZONE_0]); - + /*--- Create the point-to-point MPI communication structures. ---*/ geometry_container[ZONE_0]->PreprocessP2PComms(geometry_container[ZONE_0], config_container[ZONE_0]); - + /*--- Create plane structure ---*/ - + if (rank == MASTER_NODE) cout << "Set plane structure." << endl; - + if (geometry_container[ZONE_0]->GetnDim() == 2) { Plane_Normal[0][0] = 0.0; Plane_P0[0][0] = 0.0; Plane_Normal[0][1] = 1.0; Plane_P0[0][1] = 0.0; @@ -245,29 +240,29 @@ int main(int argc, char *argv[]) { Plane_Normal[iPlane][0] = 0.0; Plane_Normal[iPlane][1] = -sin(config_container[ZONE_0]->GetLocationStations(iPlane)*PI_NUMBER/180.0); Plane_Normal[iPlane][2] = cos(config_container[ZONE_0]->GetLocationStations(iPlane)*PI_NUMBER/180.0); - + /*--- Apply tilt angle to the plane ---*/ - + su2double Tilt_Angle = config_container[ZONE_0]->GetNacelleLocation(3)*PI_NUMBER/180; su2double Plane_NormalX_Tilt = Plane_Normal[iPlane][0]*cos(Tilt_Angle) + Plane_Normal[iPlane][2]*sin(Tilt_Angle); su2double Plane_NormalY_Tilt = Plane_Normal[iPlane][1]; su2double Plane_NormalZ_Tilt = Plane_Normal[iPlane][2]*cos(Tilt_Angle) - Plane_Normal[iPlane][0]*sin(Tilt_Angle); - + /*--- Apply toe angle to the plane ---*/ - + su2double Toe_Angle = config_container[ZONE_0]->GetNacelleLocation(4)*PI_NUMBER/180; su2double Plane_NormalX_Tilt_Toe = Plane_NormalX_Tilt*cos(Toe_Angle) - Plane_NormalY_Tilt*sin(Toe_Angle); su2double Plane_NormalY_Tilt_Toe = Plane_NormalX_Tilt*sin(Toe_Angle) + Plane_NormalY_Tilt*cos(Toe_Angle); su2double Plane_NormalZ_Tilt_Toe = Plane_NormalZ_Tilt; - + /*--- Update normal vector ---*/ - + Plane_Normal[iPlane][0] = Plane_NormalX_Tilt_Toe; Plane_Normal[iPlane][1] = Plane_NormalY_Tilt_Toe; Plane_Normal[iPlane][2] = Plane_NormalZ_Tilt_Toe; - + /*--- Point in the plane ---*/ - + Plane_P0[iPlane][0] = config_container[ZONE_0]->GetNacelleLocation(0); Plane_P0[iPlane][1] = config_container[ZONE_0]->GetNacelleLocation(1); Plane_P0[iPlane][2] = config_container[ZONE_0]->GetNacelleLocation(2); @@ -276,28 +271,28 @@ int main(int argc, char *argv[]) { Plane_Normal[iPlane][1] = 1.0; Plane_P0[iPlane][1] = config_container[ZONE_0]->GetLocationStations(iPlane); } - + } } - + /*--- Compute the wing and fan description (only 3D). ---*/ - + if (geometry_container[ZONE_0]->GetnDim() == 3) { - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { - + if (rank == MASTER_NODE) { cout << "Computing the fuselage continuous description." << endl << endl; } - + geometry_container[ZONE_0]->Compute_Fuselage(config_container[ZONE_0], true, Fuselage_Volume, Fuselage_WettedArea, Fuselage_MinWidth, Fuselage_MaxWidth, Fuselage_MinWaterLineWidth, Fuselage_MaxWaterLineWidth, Fuselage_MinHeight, Fuselage_MaxHeight, Fuselage_MaxCurvature); - + /*--- Screen output for the wing definition ---*/ - + if (rank == MASTER_NODE) { if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Fuselage volume: " << Fuselage_Volume << " in^3. "; else cout << "Fuselage volume: " << Fuselage_Volume << " m^3. "; @@ -318,22 +313,22 @@ int main(int argc, char *argv[]) { if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Fuselage max. curvature: " << Fuselage_MaxCurvature << " 1/in. " << endl; else cout << "Fuselage max. curvature: " << Fuselage_MaxCurvature << " 1/m. " << endl; } - + } - + else if (config_container[ZONE_0]->GetGeo_Description() == NACELLE) { - + if (rank == MASTER_NODE) { cout << "Computing the nacelle continuous description." << endl << endl; } - + geometry_container[ZONE_0]->Compute_Nacelle(config_container[ZONE_0], true, Nacelle_Volume, Nacelle_MinThickness, Nacelle_MaxThickness, Nacelle_MinChord, Nacelle_MaxChord, Nacelle_MinLERadius, Nacelle_MaxLERadius, Nacelle_MinToC, Nacelle_MaxToC, Nacelle_ObjFun_MinToC, Nacelle_MaxTwist); /*--- Screen output for the wing definition ---*/ - + if (rank == MASTER_NODE) { if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Nacelle volume: " << Nacelle_Volume << " in^3. "; else cout << "Nacelle volume: " << Nacelle_Volume << " m^3. "; @@ -354,22 +349,22 @@ int main(int argc, char *argv[]) { cout << "Nacelle delta ToC: " << Nacelle_ObjFun_MinToC << ". "; cout << "Nacelle max. twist: " << Nacelle_MaxTwist << " deg. "<< endl; } - + } - + else { - + if (rank == MASTER_NODE) { cout << "Computing the wing continuous description." << endl << endl; } - + geometry_container[ZONE_0]->Compute_Wing(config_container[ZONE_0], true, Wing_Volume, Wing_MinThickness, Wing_MaxThickness, Wing_MinChord, Wing_MaxChord, Wing_MinLERadius, Wing_MaxLERadius, Wing_MinToC, Wing_MaxToC, Wing_ObjFun_MinToC, Wing_MaxTwist, Wing_MaxCurvature, Wing_MaxDihedral); - + /*--- Screen output for the wing definition ---*/ - + if (rank == MASTER_NODE) { if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Wing volume: " << Wing_Volume << " in^3. "; else cout << "Wing volume: " << Wing_Volume << " m^3. "; @@ -393,29 +388,29 @@ int main(int argc, char *argv[]) { else cout << "Wing max. curvature: " << Wing_MaxCurvature << " 1/m. "; cout << "Wing max. dihedral: " << Wing_MaxDihedral << " deg." << endl; } - + } - + } - + for (iPlane = 0; iPlane < nPlane; iPlane++) { - + geometry_container[ZONE_0]->ComputeAirfoil_Section(Plane_P0[iPlane], Plane_Normal[iPlane], -1E6, 1E6, -1E6, 1E6, -1E6, 1E6, nullptr, Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane], Variable_Airfoil[iPlane], true, config_container[ZONE_0]); } - + if (rank == MASTER_NODE) cout << endl <<"-------------------- Objective function evaluation ----------------------" << endl; - + if (rank == MASTER_NODE) { - + /*--- Evaluate objective function ---*/ - + for (iPlane = 0; iPlane < nPlane; iPlane++) { - + if (Xcoord_Airfoil[iPlane].size() > 1) { - + cout << "\nStation " << (iPlane+1); if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << ". XCoord: " << Plane_P0[iPlane][0] << " in, "; @@ -437,7 +432,7 @@ int main(int argc, char *argv[]) { ObjectiveFunc[2*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_Width(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); ObjectiveFunc[3*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_WaterLineWidth(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); ObjectiveFunc[4*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_Height(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); - + if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Area: " << ObjectiveFunc[0*nPlane+iPlane] << " in^2, "; else cout << "Area: " << ObjectiveFunc[0*nPlane+iPlane] << " m^2, "; if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Length: " << ObjectiveFunc[1*nPlane+iPlane] << " in, "; @@ -456,7 +451,7 @@ int main(int argc, char *argv[]) { ObjectiveFunc[3*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_LERadius(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); ObjectiveFunc[4*nPlane+iPlane] = ObjectiveFunc[1*nPlane+iPlane]/ObjectiveFunc[2*nPlane+iPlane]; ObjectiveFunc[5*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_Twist(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); - + if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Area: " << ObjectiveFunc[0*nPlane+iPlane] << " in^2, "; else cout << "Area: " << ObjectiveFunc[0*nPlane+iPlane] << " m^2, "; if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Thickness: " << ObjectiveFunc[1*nPlane+iPlane] << " in, " << endl; @@ -476,7 +471,7 @@ int main(int argc, char *argv[]) { ObjectiveFunc[3*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_LERadius(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); ObjectiveFunc[4*nPlane+iPlane] = ObjectiveFunc[1*nPlane+iPlane]/ObjectiveFunc[2*nPlane+iPlane]; ObjectiveFunc[5*nPlane+iPlane] = geometry_container[ZONE_0]->Compute_Twist(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); - + if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Area: " << ObjectiveFunc[0*nPlane+iPlane] << " in^2, "; else cout << "Area: " << ObjectiveFunc[0*nPlane+iPlane] << " m^2, "; if (config_container[ZONE_0]->GetSystemMeasurements() == US) cout << "Thickness: " << ObjectiveFunc[1*nPlane+iPlane] << " in, " << endl; @@ -489,11 +484,11 @@ int main(int argc, char *argv[]) { if (geometry_container[ZONE_0]->GetnDim() == 2) cout << "Alpha: " << ObjectiveFunc[5*nPlane+iPlane] <<" deg."; else if (geometry_container[ZONE_0]->GetnDim() == 3) cout << "Twist angle: " << ObjectiveFunc[5*nPlane+iPlane] <<" deg."; } - + } - + } - + /*--- Write the objective function in a external file ---*/ string filename = config_container[ZONE_0]->GetObjFunc_Value_FileName(); unsigned short lastindex = filename.find_last_of("."); @@ -502,15 +497,15 @@ int main(int argc, char *argv[]) { else filename += ".csv"; ObjFunc_file.open(filename.c_str(), ios::out); if (tabTecplot) ObjFunc_file << "TITLE = \"SU2_GEO Evaluation\"" << endl; - + if (geometry_container[ZONE_0]->GetnDim() == 2) { if (tabTecplot) ObjFunc_file << "VARIABLES =//" << endl; ObjFunc_file << "\"AIRFOIL_AREA\",\"AIRFOIL_THICKNESS\",\"AIRFOIL_CHORD\",\"AIRFOIL_LE_RADIUS\",\"AIRFOIL_TOC\",\"AIRFOIL_ALPHA\""; } else if (geometry_container[ZONE_0]->GetnDim() == 3) { - + if (tabTecplot) ObjFunc_file << "VARIABLES = //" << endl; - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { ObjFunc_file << "\"FUSELAGE_VOLUME\",\"FUSELAGE_WETTED_AREA\",\"FUSELAGE_MIN_WIDTH\",\"FUSELAGE_MAX_WIDTH\",\"FUSELAGE_MIN_WATERLINE_WIDTH\",\"FUSELAGE_MAX_WATERLINE_WIDTH\",\"FUSELAGE_MIN_HEIGHT\",\"FUSELAGE_MAX_HEIGHT\",\"FUSELAGE_MAX_CURVATURE\","; for (iPlane = 0; iPlane < nPlane; iPlane++) ObjFunc_file << "\"STATION"<< (iPlane+1) << "_AREA\","; @@ -546,12 +541,12 @@ int main(int argc, char *argv[]) { if (iPlane != nPlane-1) ObjFunc_file << ","; } } - + } - + if (tabTecplot) ObjFunc_file << "\nZONE T= \"Geometrical variables (value)\"" << endl; else ObjFunc_file << endl; - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { if (geometry_container[ZONE_0]->GetnDim() == 3) { ObjFunc_file << Fuselage_Volume <<", "<< Fuselage_WettedArea <<", "<< Fuselage_MinWidth <<", "<< Fuselage_MaxWidth <<", "<< Fuselage_MinWaterLineWidth <<", "<< Fuselage_MaxWaterLineWidth<<", "<< Fuselage_MinHeight <<", "<< Fuselage_MaxHeight <<", "<< Fuselage_MaxCurvature <<", "; @@ -579,26 +574,26 @@ int main(int argc, char *argv[]) { if (iPlane != (nPlane*6)-1) ObjFunc_file <<", "; } } - + ObjFunc_file.close(); - + } - + if (config_container[ZONE_0]->GetGeometryMode() == GRADIENT) { - + /*--- Definition of the Class for surface deformation ---*/ surface_movement = new CSurfaceMovement(); - + /*--- Copy coordinates to the surface structure ---*/ surface_movement->CopyBoundary(geometry_container[ZONE_0], config_container[ZONE_0]); - + /*--- Definition of the FFD deformation class ---*/ FFDBox = new CFreeFormDefBox*[MAX_NUMBER_FFD]; for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) FFDBox[iFFDBox] = nullptr; - + if (rank == MASTER_NODE) cout << endl << endl << "------------- Gradient evaluation using finite differences --------------" << endl; - + /*--- Write the gradient in a external file ---*/ if (rank == MASTER_NODE) { string filename = config_container[ZONE_0]->GetObjFunc_Grad_FileName(); @@ -608,11 +603,11 @@ int main(int argc, char *argv[]) { else filename += ".csv"; Gradient_file.open(filename.c_str(), ios::out); } - + for (iDV = 0; iDV < config_container[ZONE_0]->GetnDV(); iDV++) { /*--- Free Form deformation based ---*/ - + if ((config_container[ZONE_0]->GetDesign_Variable(iDV) == FFD_CONTROL_POINT_2D) || (config_container[ZONE_0]->GetDesign_Variable(iDV) == FFD_CAMBER_2D) || (config_container[ZONE_0]->GetDesign_Variable(iDV) == FFD_THICKNESS_2D) || @@ -624,19 +619,19 @@ int main(int argc, char *argv[]) { (config_container[ZONE_0]->GetDesign_Variable(iDV) == FFD_ROTATION) || (config_container[ZONE_0]->GetDesign_Variable(iDV) == FFD_CAMBER) || (config_container[ZONE_0]->GetDesign_Variable(iDV) == FFD_THICKNESS) ) { - + /*--- Read the FFD information in the first iteration ---*/ - + if (iDV == 0) { - + if (rank == MASTER_NODE) cout << "Read the FFD information from mesh file." << endl; - + /*--- Read the FFD information from the grid file ---*/ - + surface_movement->ReadFFDInfo(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox, config_container[ZONE_0]->GetMesh_FileName()); - + /*--- Modify the control points for polar based computations ---*/ - + if (config_container[ZONE_0]->GetFFD_CoordSystem() == CYLINDRICAL) { for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { FFDBox[iFFDBox]->SetCart2Cyl_ControlPoints(config_container[ZONE_0]); @@ -652,40 +647,40 @@ int main(int argc, char *argv[]) { FFDBox[iFFDBox]->SetCart2Sphe_ControlPoints(config_container[ZONE_0]); } } - + /*--- If the FFDBox was not defined in the input file ---*/ - + if (!surface_movement->GetFFDBoxDefinition()) { SU2_MPI::Error("The input grid doesn't have the entire FFD information!", CURRENT_FUNCTION); } - + for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { - + if (rank == MASTER_NODE) cout << "Checking FFD box dimension." << endl; surface_movement->CheckFFDDimension(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], iFFDBox); - - + + if (rank == MASTER_NODE) cout << "Check the FFD box intersections with the solid surfaces." << endl; surface_movement->CheckFFDIntersections(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], iFFDBox); - + } - + if (rank == MASTER_NODE) cout <<"-------------------------------------------------------------------------" << endl; - + } - + if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; cout << "Perform 3D deformation of the surface." << endl; } - + /*--- Apply the control point change ---*/ - + MoveSurface = false; - + for (iFFDBox = 0; iFFDBox < surface_movement->GetnFFDBox(); iFFDBox++) { - + switch ( config_container[ZONE_0]->GetDesign_Variable(iDV) ) { case FFD_CONTROL_POINT_2D : Local_MoveSurface = surface_movement->SetFFDCPChange_2D(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], FFDBox, iDV, true); break; case FFD_CAMBER_2D : Local_MoveSurface = surface_movement->SetFFDCamber_2D(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], FFDBox, iDV, true); break; @@ -700,20 +695,20 @@ int main(int argc, char *argv[]) { case FFD_THICKNESS : Local_MoveSurface = surface_movement->SetFFDThickness(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], FFDBox, iDV, true); break; case FFD_CONTROL_SURFACE : Local_MoveSurface = surface_movement->SetFFDControl_Surface(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], FFDBox, iDV, true); break; } - + /*--- Recompute cartesian coordinates using the new control points position ---*/ - + if (Local_MoveSurface) { MoveSurface = true; surface_movement->SetCartesianCoord(geometry_container[ZONE_0], config_container[ZONE_0], FFDBox[iFFDBox], iFFDBox, true); } - + } - + } - + /*--- Hicks Henne design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == HICKS_HENNE) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -722,9 +717,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetHicksHenne(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true); } - + /*--- Surface bump design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == SURFACE_BUMP) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -733,9 +728,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetSurface_Bump(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true); } - + /*--- CST design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == CST) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -744,9 +739,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetCST(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true); } - + /*--- Translation design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == TRANSLATION) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -755,9 +750,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetTranslation(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true); } - + /*--- Scale design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == SCALE) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -766,9 +761,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetScale(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true); } - + /*--- Rotation design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == ROTATION) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -777,9 +772,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetRotation(geometry_container[ZONE_0], config_container[ZONE_0], iDV, true); } - + /*--- NACA_4Digits design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == NACA_4DIGITS) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -788,9 +783,9 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetNACA_4Digits(geometry_container[ZONE_0], config_container[ZONE_0]); } - + /*--- Parabolic design variable ---*/ - + else if (config_container[ZONE_0]->GetDesign_Variable(iDV) == PARABOLIC) { if (rank == MASTER_NODE) { cout << endl << "Design variable number "<< iDV <<"." << endl; @@ -799,21 +794,21 @@ int main(int argc, char *argv[]) { MoveSurface = true; surface_movement->SetParabolic(geometry_container[ZONE_0], config_container[ZONE_0]); } - + /*--- Design variable not implement ---*/ - + else { if (rank == MASTER_NODE) cout << "Design Variable not implemented yet" << endl; } - + if (MoveSurface) { - + /*--- Compute the gradient for the volume. In 2D this is just the gradient of the area. ---*/ - + if (geometry_container[ZONE_0]->GetnDim() == 3) { - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { geometry_container[ZONE_0]->Compute_Fuselage(config_container[ZONE_0], false, Fuselage_Volume_New, Fuselage_WettedArea_New, Fuselage_MinWidth_New, Fuselage_MaxWidth_New, @@ -833,31 +828,31 @@ int main(int argc, char *argv[]) { Wing_MaxChord_New, Wing_MinLERadius_New, Wing_MaxLERadius_New, Wing_MinToC_New, Wing_MaxToC_New, Wing_ObjFun_MinToC_New, Wing_MaxTwist_New, Wing_MaxCurvature_New, Wing_MaxDihedral_New); } - + } - + /*--- Create airfoil structure ---*/ - + for (iPlane = 0; iPlane < nPlane; iPlane++) { geometry_container[ZONE_0]->ComputeAirfoil_Section(Plane_P0[iPlane], Plane_Normal[iPlane], -1E6, 1E6, -1E6, 1E6, -1E6, 1E6, nullptr, Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane], Variable_Airfoil[iPlane], false, config_container[ZONE_0]); } - + } - + /*--- Compute gradient ---*/ - + if (rank == MASTER_NODE) { - + delta_eps = config_container[ZONE_0]->GetDV_Value(iDV); - + if (delta_eps == 0) { SU2_MPI::Error("The finite difference steps is zero!!", CURRENT_FUNCTION); } - + if (MoveSurface) { - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { Fuselage_Volume_Grad = (Fuselage_Volume_New - Fuselage_Volume) / delta_eps; Fuselage_WettedArea_Grad = (Fuselage_WettedArea_New - Fuselage_WettedArea) / delta_eps; @@ -868,7 +863,7 @@ int main(int argc, char *argv[]) { Fuselage_MinHeight_Grad = (Fuselage_MinHeight_New - Fuselage_MinHeight) / delta_eps; Fuselage_MaxHeight_Grad = (Fuselage_MaxHeight_New - Fuselage_MaxHeight) / delta_eps; Fuselage_MaxCurvature_Grad = (Fuselage_MaxCurvature_New - Fuselage_MaxCurvature) / delta_eps; - + } else if (config_container[ZONE_0]->GetGeo_Description() == NACELLE) { Nacelle_Volume_Grad = (Nacelle_Volume_New - Nacelle_Volume) / delta_eps; @@ -898,80 +893,80 @@ int main(int argc, char *argv[]) { Wing_MaxCurvature_Grad = (Wing_MaxCurvature_New - Wing_MaxCurvature) / delta_eps; Wing_MaxDihedral_Grad = (Wing_MaxDihedral_New - Wing_MaxDihedral) / delta_eps; } - + for (iPlane = 0; iPlane < nPlane; iPlane++) { if (Xcoord_Airfoil[iPlane].size() > 1) { - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { - + ObjectiveFunc_New[0*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Area(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[0*nPlane + iPlane] = (ObjectiveFunc_New[0*nPlane + iPlane] - ObjectiveFunc[0*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[1*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Length(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[1*nPlane + iPlane] = (ObjectiveFunc_New[1*nPlane + iPlane] - ObjectiveFunc[1*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[2*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Width(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[2*nPlane + iPlane] = (ObjectiveFunc_New[2*nPlane + iPlane] - ObjectiveFunc[2*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[3*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_WaterLineWidth(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[3*nPlane + iPlane] = (ObjectiveFunc_New[3*nPlane + iPlane] - ObjectiveFunc[3*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[4*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Height(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[4*nPlane + iPlane] = (ObjectiveFunc_New[4*nPlane + iPlane] - ObjectiveFunc[4*nPlane + iPlane]) / delta_eps; - + } - + else if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { - + ObjectiveFunc_New[0*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Area(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[0*nPlane + iPlane] = (ObjectiveFunc_New[0*nPlane + iPlane] - ObjectiveFunc[0*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[1*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_MaxThickness(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[1*nPlane + iPlane] = (ObjectiveFunc_New[1*nPlane + iPlane] - ObjectiveFunc[1*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[2*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Chord(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[2*nPlane + iPlane] = (ObjectiveFunc_New[2*nPlane + iPlane] - ObjectiveFunc[2*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[3*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_LERadius(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[3*nPlane + iPlane] = (ObjectiveFunc_New[3*nPlane + iPlane] - ObjectiveFunc[3*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[4*nPlane + iPlane] = ObjectiveFunc_New[1*nPlane + iPlane] / ObjectiveFunc_New[2*nPlane + iPlane]; Gradient[4*nPlane + iPlane] = (ObjectiveFunc_New[4*nPlane + iPlane] - ObjectiveFunc[4*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[5*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Twist(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[5*nPlane + iPlane] = (ObjectiveFunc_New[5*nPlane + iPlane] - ObjectiveFunc[5*nPlane + iPlane]) / delta_eps; } - + else { - + ObjectiveFunc_New[0*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Area(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[0*nPlane + iPlane] = (ObjectiveFunc_New[0*nPlane + iPlane] - ObjectiveFunc[0*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[1*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_MaxThickness(Plane_P0[iPlane], Plane_Normal[iPlane], config_container[ZONE_0], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[1*nPlane + iPlane] = (ObjectiveFunc_New[1*nPlane + iPlane] - ObjectiveFunc[1*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[2*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Chord(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[2*nPlane + iPlane] = (ObjectiveFunc_New[2*nPlane + iPlane] - ObjectiveFunc[2*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[3*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_LERadius(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[3*nPlane + iPlane] = (ObjectiveFunc_New[3*nPlane + iPlane] - ObjectiveFunc[3*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[4*nPlane + iPlane] = ObjectiveFunc_New[1*nPlane + iPlane] / ObjectiveFunc_New[2*nPlane + iPlane]; Gradient[4*nPlane + iPlane] = (ObjectiveFunc_New[4*nPlane + iPlane] - ObjectiveFunc[4*nPlane + iPlane]) / delta_eps; - + ObjectiveFunc_New[5*nPlane + iPlane] = geometry_container[ZONE_0]->Compute_Twist(Plane_P0[iPlane], Plane_Normal[iPlane], Xcoord_Airfoil[iPlane], Ycoord_Airfoil[iPlane], Zcoord_Airfoil[iPlane]); Gradient[5*nPlane + iPlane] = (ObjectiveFunc_New[5*nPlane + iPlane] - ObjectiveFunc[5*nPlane + iPlane]) / delta_eps; - + } - + } } - + } - + else { - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { Fuselage_Volume_Grad = 0.0; Fuselage_WettedArea_Grad = 0.0; @@ -982,7 +977,7 @@ int main(int argc, char *argv[]) { Fuselage_MinHeight_Grad = 0.0; Fuselage_MaxHeight_Grad = 0.0; Fuselage_MaxCurvature_Grad = 0.0; - + for (iPlane = 0; iPlane < nPlane; iPlane++) { Gradient[0*nPlane + iPlane] = 0.0; Gradient[1*nPlane + iPlane] = 0.0; @@ -990,7 +985,7 @@ int main(int argc, char *argv[]) { Gradient[3*nPlane + iPlane] = 0.0; Gradient[4*nPlane + iPlane] = 0.0; } - + } else if (config_container[ZONE_0]->GetGeo_Description() == NACELLE) { Nacelle_Volume_Grad = 0.0; @@ -1004,7 +999,7 @@ int main(int argc, char *argv[]) { Nacelle_MaxToC_Grad = 0.0; Nacelle_ObjFun_MinToC_Grad = 0.0; Nacelle_MaxTwist_Grad = 0.0; - + for (iPlane = 0; iPlane < nPlane; iPlane++) { Gradient[0*nPlane + iPlane] = 0.0; Gradient[1*nPlane + iPlane] = 0.0; @@ -1028,7 +1023,7 @@ int main(int argc, char *argv[]) { Wing_MaxTwist_Grad = 0.0; Wing_MaxCurvature_Grad = 0.0; Wing_MaxDihedral_Grad = 0.0; - + for (iPlane = 0; iPlane < nPlane; iPlane++) { Gradient[0*nPlane + iPlane] = 0.0; Gradient[1*nPlane + iPlane] = 0.0; @@ -1038,11 +1033,11 @@ int main(int argc, char *argv[]) { Gradient[5*nPlane + iPlane] = 0.0; } } - + } - + /*--- Screen output ---*/ - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { if (geometry_container[ZONE_0]->GetnDim() == 3) { cout << "\nFuselage volume grad.: " << Fuselage_Volume_Grad << ". "; @@ -1055,7 +1050,7 @@ int main(int argc, char *argv[]) { cout << "Fuselage max. height grad.: " << Fuselage_MaxHeight_Grad << ". "; cout << "Fuselage max. curv. grad.: " << Fuselage_MaxCurvature_Grad << "."; } - + for (iPlane = 0; iPlane < nPlane; iPlane++) { if (Xcoord_Airfoil[iPlane].size() > 1) { cout << "\nStation " << (iPlane+1) << ". XCoord: " << Plane_P0[iPlane][0] << ". "; @@ -1081,7 +1076,7 @@ int main(int argc, char *argv[]) { cout << "Nacelle delta ToC grad.: " << Nacelle_ObjFun_MinToC_Grad << "." << endl; cout << "Nacelle max. twist grad.: " << Nacelle_MaxTwist_Grad << ". "; } - + for (iPlane = 0; iPlane < nPlane; iPlane++) { if (Xcoord_Airfoil[iPlane].size() > 1) { cout << "\nStation " << (iPlane+1) << ". YCoord: " << Plane_P0[iPlane][1] << ". "; @@ -1110,7 +1105,7 @@ int main(int argc, char *argv[]) { cout << "Wing max. curv. grad.: " << Wing_MaxCurvature_Grad << ". "; cout << "Wing max. dihedral grad.: " << Wing_MaxDihedral_Grad << "." << endl; } - + for (iPlane = 0; iPlane < nPlane; iPlane++) { if (Xcoord_Airfoil[iPlane].size() > 1) { cout << "\nStation " << (iPlane+1) << ". YCoord: " << Plane_P0[iPlane][1] << ". "; @@ -1123,19 +1118,19 @@ int main(int argc, char *argv[]) { } } } - + cout << endl; - - + + if (iDV == 0) { if (tabTecplot) Gradient_file << "TITLE = \"SU2_GEO Gradient\"" << endl; if (tabTecplot) Gradient_file << "VARIABLES = //" << endl; - + if (geometry_container[ZONE_0]->GetnDim() == 2) { Gradient_file << "\"DESIGN_VARIABLE\",\"AIRFOIL_AREA\",\"AIRFOIL_THICKNESS\",\"AIRFOIL_CHORD\",\"AIRFOIL_LE_RADIUS\",\"AIRFOIL_TOC\",\"AIRFOIL_ALPHA\""; } else if (geometry_container[ZONE_0]->GetnDim() == 3) { - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { Gradient_file << "\"DESIGN_VARIABLE\","; Gradient_file << "\"FUSELAGE_VOLUME\",\"FUSELAGE_WETTED_AREA\",\"FUSELAGE_MIN_WIDTH\",\"FUSELAGE_MAX_WIDTH\",\"FUSELAGE_MIN_WATERLINE_WIDTH\",\"FUSELAGE_MAX_WATERLINE_WIDTH\",\"FUSELAGE_MIN_HEIGHT\",\"FUSELAGE_MAX_HEIGHT\",\"FUSELAGE_MAX_CURVATURE\","; @@ -1174,15 +1169,15 @@ int main(int argc, char *argv[]) { if (iPlane != nPlane-1) Gradient_file << ","; } } - + } - + if (tabTecplot) Gradient_file << "\nZONE T= \"Geometrical variables (gradient)\"" << endl; else Gradient_file << endl; } - + Gradient_file << (iDV) <<","; - + if (config_container[ZONE_0]->GetGeo_Description() == FUSELAGE) { if (geometry_container[ZONE_0]->GetnDim() == 3) { Gradient_file << Fuselage_Volume_Grad <<","<< Fuselage_WettedArea_Grad <<","<< Fuselage_MinWidth_Grad <<","<< Fuselage_MaxWidth_Grad <<","<< Fuselage_MinWaterLineWidth_Grad <<","<< Fuselage_MaxWaterLineWidth_Grad <<","<< Fuselage_MinHeight_Grad <<","<< Fuselage_MaxHeight_Grad <<","<< Fuselage_MaxCurvature_Grad <<","; @@ -1210,40 +1205,40 @@ int main(int argc, char *argv[]) { if (iPlane != (nPlane*6)-1) Gradient_file <<","; } } - + Gradient_file << endl; - + if (iDV != (config_container[ZONE_0]->GetnDV()-1)) cout <<"-------------------------------------------------------------------------" << endl; - + } - + } - + if (rank == MASTER_NODE) Gradient_file.close(); - + } if (rank == MASTER_NODE) cout << endl <<"------------------------- Solver Postprocessing -------------------------" << endl; - + delete [] Xcoord_Airfoil; delete [] Ycoord_Airfoil; delete [] Zcoord_Airfoil; - + delete [] ObjectiveFunc; delete [] ObjectiveFunc_New; delete [] Gradient; - + for(iPlane = 0; iPlane < nPlane; iPlane++ ) { delete Plane_P0[iPlane]; delete Plane_Normal[iPlane]; } delete [] Plane_P0; delete [] Plane_Normal; - + delete config; config = nullptr; if (rank == MASTER_NODE) cout << "Deleted main variables." << endl; - - + + if (geometry_container != nullptr) { for (iZone = 0; iZone < nZone; iZone++) { if (geometry_container[iZone] != nullptr) { @@ -1256,7 +1251,7 @@ int main(int argc, char *argv[]) { delete surface_movement; if (rank == MASTER_NODE) cout << "Deleted CSurfaceMovement class." << endl; - + if (FFDBox != nullptr) { for (iFFDBox = 0; iFFDBox < MAX_NUMBER_FFD; iFFDBox++) { if (FFDBox[iFFDBox] != nullptr) { @@ -1266,7 +1261,7 @@ int main(int argc, char *argv[]) { delete [] FFDBox; } if (rank == MASTER_NODE) cout << "Deleted CFreeFormDefBox class." << endl; - + if (config_container != nullptr) { for (iZone = 0; iZone < nZone; iZone++) { if (config_container[iZone] != nullptr) { @@ -1276,32 +1271,32 @@ int main(int argc, char *argv[]) { delete [] config_container; } if (rank == MASTER_NODE) cout << "Deleted CConfig container." << endl; - + /*--- Synchronization point after a single solver iteration. Compute the wall clock time required. ---*/ StopTime = SU2_MPI::Wtime(); /*--- Compute/print the total time for performance benchmarking. ---*/ - + UsedTime = StopTime-StartTime; if (rank == MASTER_NODE) { cout << "\n\nCompleted in " << fixed << UsedTime << " seconds on "<< size; if (size == 1) cout << " core." << endl; else cout << " cores." << endl; } - + /*--- Exit the solver cleanly ---*/ - + if (rank == MASTER_NODE) cout << endl <<"------------------------- Exit Success (SU2_GEO) ------------------------" << endl << endl; - - + + /*--- Finalize MPI parallelization ---*/ - + #ifdef HAVE_MPI SU2_MPI::Finalize(); #endif - + return EXIT_SUCCESS; - + } diff --git a/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref b/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref index 5e90b70e61f9..da2e3f10472c 100644 --- a/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref +++ b/TestCases/multiple_ffd/naca0012/of_grad_directdiff.dat.ref @@ -1,3 +1,3 @@ VARIABLES="VARIABLE" , "DRAG" , "EFFICIENCY" , "FORCE_X" , "FORCE_Y" , "FORCE_Z" , "LIFT" , "MOMENT_X" , "MOMENT_Y" , "MOMENT_Z" , "SIDEFORCE" - 0 , 0.05682885157 , -0.4029984608 , 0.05415019059 , 0.123381226 , 0.0 , 0.1221705845 , 0.0 , 0.0 , 0.009249712749 , 0.0 - 1 , -0.08020054725 , 8.95834958 , -0.08614128407 , 0.2713852273 , 0.0 , 0.2731998072 , 0.0 , 0.0 , 0.00847998417 , 0.0 + 0 , 0.05685349197 , -0.4081397694 , 0.05417489043 , 0.1233787706 , 0.0 , 0.1221675908 , 0.0 , 0.0 , 0.009249264651 , 0.0 + 1 , -0.08020555475 , 8.976256797 , -0.08614647972 , 0.2713937953 , 0.0 , 0.2732084865 , 0.0 , 0.0 , 0.00847918122 , 0.0 diff --git a/UnitTests/Common/geometry/CGeometry_test.cpp b/UnitTests/Common/geometry/CGeometry_test.cpp index ef4408945a68..1fe19682dea9 100644 --- a/UnitTests/Common/geometry/CGeometry_test.cpp +++ b/UnitTests/Common/geometry/CGeometry_test.cpp @@ -125,27 +125,14 @@ TEST_CASE("Set vertex", "[Geometry]"){ } -TEST_CASE("Set center of gravity", "[Geometry]"){ +TEST_CASE("Set control volume", "[Geometry]"){ - TestCase->geometry->SetCoord_CG(); + TestCase->geometry->SetControlVolume(TestCase->config.get(), ALLOCATE); CHECK(TestCase->geometry->elem[42]->GetCG(0) == 0.625); CHECK(TestCase->geometry->elem[3]->GetCG(1) == 0.125); CHECK(TestCase->geometry->elem[25]->GetCG(2) == 0.375); - CHECK(TestCase->geometry->bound[1][4]->GetCG(0) == 1.0); - CHECK(TestCase->geometry->bound[3][2]->GetCG(1) == 1.0); - CHECK(TestCase->geometry->bound[4][3]->GetCG(2) == 0.0); - - CHECK(TestCase->geometry->edges->GetCG(10, 0) == 0.75); - CHECK(TestCase->geometry->edges->GetCG(3, 1) == 0.125); - CHECK(TestCase->geometry->edges->GetCG(22, 2) == 0.0); -} - -TEST_CASE("Set control volume", "[Geometry]"){ - - TestCase->geometry->SetControlVolume(TestCase->config.get(), ALLOCATE); - CHECK(TestCase->geometry->nodes->GetVolume(42) == Approx(0.015625)); CHECK(TestCase->geometry->edges->GetNormal(32)[0] == 0.03125); @@ -160,6 +147,10 @@ TEST_CASE("Set bound control volume", "[Geometry]"){ TestCase->geometry->SetBoundControlVolume(TestCase->config.get(), ALLOCATE); + CHECK(TestCase->geometry->bound[1][4]->GetCG(0) == 1.0); + CHECK(TestCase->geometry->bound[3][2]->GetCG(1) == 1.0); + CHECK(TestCase->geometry->bound[4][3]->GetCG(2) == 0.0); + CHECK(TestCase->geometry->vertex[0][4]->GetNormal()[0] == -0.0625); CHECK(TestCase->geometry->vertex[3][2]->GetNormal()[1] == -0.0625); CHECK(TestCase->geometry->vertex[5][3]->GetNormal()[2] == 0.03125); diff --git a/UnitTests/SU2_CFD/gradients.cpp b/UnitTests/SU2_CFD/gradients.cpp index b5d30b09fef9..b78e7595a119 100644 --- a/UnitTests/SU2_CFD/gradients.cpp +++ b/UnitTests/SU2_CFD/gradients.cpp @@ -87,7 +87,6 @@ struct GradientTestBase { geometry->Check_BoundElem_Orientation(config.get()); geometry->SetEdges(); geometry->SetVertex(config.get()); - geometry->SetCoord_CG(); geometry->SetControlVolume(config.get(), ALLOCATE); geometry->SetBoundControlVolume(config.get(), ALLOCATE); geometry->FindNormal_Neighbor(config.get()); diff --git a/UnitTests/UnitQuadTestCase.hpp b/UnitTests/UnitQuadTestCase.hpp index 07cdf9211b0f..8fad8d7892f1 100644 --- a/UnitTests/UnitQuadTestCase.hpp +++ b/UnitTests/UnitQuadTestCase.hpp @@ -99,7 +99,6 @@ struct UnitQuadTestCase { geometry->Check_BoundElem_Orientation(config.get()); geometry->SetEdges(); geometry->SetVertex(config.get()); - geometry->SetCoord_CG(); geometry->SetControlVolume(config.get(), ALLOCATE); geometry->SetBoundControlVolume(config.get(), ALLOCATE); geometry->FindNormal_Neighbor(config.get());