Merge branch 'develop' into develop

Author: kursatyurt

Common/include/CConfig.hpp

@@ -5804,6 +5804,13 @@ class CConfig {
    */
   su2double GetTranslation_Rate(unsigned short iDim) const { return Translation_Rate[iDim];}
 
+  /*!
+   * \brief Set the translational velocity of the mesh.
+   * \param[in] iDim - spatial component
+   * \return Translational velocity of the mesh.
+   */
+  void SetTranslation_Rate(unsigned short iDim, su2double val) { Translation_Rate[iDim] = val;}
+
   /*!
    * \brief Get the translational velocity of the marker.
    * \param[in] iMarkerMoving -  Index of the moving marker (as specified in Marker_Moving)

Common/include/basic_types/ad_structure.hpp

@@ -74,6 +74,11 @@ inline void RegisterInput(su2double& data, bool push_index = true) {}
  */
 inline void RegisterOutput(su2double& data) {}
 
+/*!
+ * \brief Resize the adjoint vector, for subsequent access without bounds checking.
+ */
+inline void ResizeAdjoints() {}
+
 /*!
  * \brief Sets the adjoint value at index to val
  * \param[in] index - Position in the adjoint vector.
@@ -369,11 +374,27 @@ FORCEINLINE void Reset() {
   }
 }
 
+FORCEINLINE void ResizeAdjoints() { AD::getTape().resizeAdjointVector(); }
+
 FORCEINLINE void SetIndex(int& index, const su2double& data) { index = data.getIdentifier(); }
 
-FORCEINLINE void SetDerivative(int index, const double val) { AD::getTape().setGradient(index, val); }
+// WARNING: For performance reasons, this method does not perform bounds checking.
+// When using it, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
+FORCEINLINE void SetDerivative(int index, const double val) {
+  if (index == 0)  // Allow multiple threads to "set the derivative" of passive variables without causing data races.
+    return;
 
-FORCEINLINE double GetDerivative(int index) { return AD::getTape().getGradient(index); }
+  using BoundsChecking = codi::GradientAccessTapeInterface<su2double::Gradient, su2double::Identifier>::BoundsChecking;
+  AD::getTape().setGradient(index, val, BoundsChecking::False);
+}
+
+// WARNING: For performance reasons, this method does not perform bounds checking.
+// If called after tape evaluations, the adjoints should exist.
+// Otherwise, please ensure sufficient adjoint vector size by a call to AD::ResizeAdjoints().
+FORCEINLINE double GetDerivative(int index) {
+  using BoundsChecking = codi::GradientAccessTapeInterface<su2double::Gradient, su2double::Identifier>::BoundsChecking;
+  return AD::getTape().getGradient(index, BoundsChecking::False);
+}
 
 FORCEINLINE bool IsIdentifierActive(su2double const& value) {
   return getTape().isIdentifierActive(value.getIdentifier());
@@ -523,26 +544,14 @@ FORCEINLINE void delete_handler(void* handler) {
 
 FORCEINLINE bool BeginPassive() {
   if (AD::getTape().isActive()) {
-    StopRecording();
+    AD::getTape().setPassive();
     return true;
   }
   return false;
 }
 
 FORCEINLINE void EndPassive(bool wasActive) {
-  if (wasActive) StartRecording();
-}
-
-FORCEINLINE bool PausePreaccumulation() {
-  const auto current = PreaccEnabled;
-  if (!current) return false;
-  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
-  return true;
-}
-
-FORCEINLINE void ResumePreaccumulation(bool wasActive) {
-  if (!wasActive) return;
-  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
+  if (wasActive) AD::getTape().setActive();
 }
 
 FORCEINLINE void StartNoSharedReading() {
@@ -558,6 +567,19 @@ FORCEINLINE void EndNoSharedReading() {
   opdi::logic->addReverseBarrier();
 #endif
 }
+
+FORCEINLINE bool PausePreaccumulation() {
+  const auto current = PreaccEnabled;
+  if (!current) return false;
+  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = false;)
+  return true;
+}
+
+FORCEINLINE void ResumePreaccumulation(bool wasActive) {
+  if (!wasActive) return;
+  SU2_OMP_SAFE_GLOBAL_ACCESS(PreaccEnabled = true;)
+}
+
 #endif  // CODI_REVERSE_TYPE
 
 void Initialize();

Common/include/code_config.hpp

@@ -96,7 +96,7 @@ FORCEINLINE Out su2staticcast_p(In ptr) {
 #include "codi/tools/data/externalFunctionUserData.hpp"
 
 #if defined(HAVE_OMP)
-using su2double = codi::RealReverseIndexOpenMP;
+using su2double = codi::RealReverseIndexOpenMPGen<double, double>;
 #else
 #if defined(CODI_INDEX_TAPE)
 using su2double = codi::RealReverseIndex;

Common/include/geometry/CGeometry.hpp

@@ -595,17 +595,6 @@ class CGeometry {
     return FindEdge(first_point, second_point, false) >= 0;
   }
 
-  /*!
-   * \brief Get the distance between a plane (defined by three point) and a point.
-   * \param[in] Coord - Coordinates of the point.
-   * \param[in] iCoord - Coordinates of the first point that defines the plane.
-   * \param[in] jCoord - Coordinates of the second point that defines the plane.
-   * \param[in] kCoord - Coordinates of the third point that defines the plane.
-   * \return Signed distance.
-   */
-  su2double Point2Plane_Distance(const su2double* Coord, const su2double* iCoord, const su2double* jCoord,
-                                 const su2double* kCoord);
-
   /*!
    * \brief Create a file for testing the geometry.
    */

Common/include/grid_movement/CFreeFormDefBox.hpp

@@ -642,12 +642,10 @@ class CFreeFormDefBox : public CGridMovement {
   }
 
   /*!
-   * \brief Set, at each vertex, the index of the free form FFDBox that contains the vertex.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   * \param[in] iFFDBox - Index of the FFDBox.
+   * \brief Returns true if the point is inside the FFD.
+   * \param[in] coord - Coordinate of the point to check.
    */
-  bool GetPointFFD(CGeometry* geometry, CConfig* config, unsigned long iPoint) const;
+  bool CheckPointInsideFFD(const su2double* coord) const;
 
   /*!
    * \brief Set the zone of the computational domain that is going to be deformed.

Common/include/option_structure.hpp

@@ -779,11 +779,13 @@ static const MapType<std::string, ENUM_GUST_TYPE> Gust_Type_Map = {
  */
 enum ENUM_GUST_DIR {
   X_DIR = 0,  /*!< \brief Gust direction-X. */
-  Y_DIR = 1   /*!< \brief Gust direction-Y. */
+  Y_DIR = 1,  /*!< \brief Gust direction-Y. */
+  Z_DIR = 2   /*!< \brief Gust direction-Z. */
 };
 static const MapType<std::string, ENUM_GUST_DIR> Gust_Dir_Map = {
   MakePair("X_DIR", X_DIR)
   MakePair("Y_DIR", Y_DIR)
+  MakePair("Z_DIR", Z_DIR)
 };
 
 // If you add to ENUM_CENTERED, you must also add the option to ENUM_CONVECTIVE

Common/include/parallelization/omp_structure.hpp

@@ -188,12 +188,14 @@ void omp_finalize();
  * thread, with all threads and memory views synchronized both beforehand and afterwards.
  */
 
-#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS \
-  SU2_OMP_BARRIER                        \
+#define BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS             \
+  SU2_OMP_BARRIER                                    \
+  if (omp_in_parallel()) AD::StartNoSharedReading(); \
   SU2_OMP_MASTER
 
-#define END_SU2_OMP_SAFE_GLOBAL_ACCESS \
-  END_SU2_OMP_MASTER                   \
+#define END_SU2_OMP_SAFE_GLOBAL_ACCESS             \
+  END_SU2_OMP_MASTER                               \
+  if (omp_in_parallel()) AD::EndNoSharedReading(); \
   SU2_OMP_BARRIER
 
 #define SU2_OMP_SAFE_GLOBAL_ACCESS(...) BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS{__VA_ARGS__} END_SU2_OMP_SAFE_GLOBAL_ACCESS

Common/include/toolboxes/geometry_toolbox.hpp

@@ -150,6 +150,15 @@ inline void QuadrilateralNormal(const T& coords, U* normal) {
   normal[2] *= 0.5;
 }
 
+/*! \brief Signed distance from a point to a plane defined by 3 coordinates. */
+template <class T, class U>
+inline U PointToPlaneDistance(const T& coords, const U* point) {
+  U normal[3], distance[3];
+  TriangleNormal(coords, normal);
+  Distance(3, point, coords[0], distance);
+  return DotProduct(3, distance, normal) / Norm(3, normal);
+}
+
 /*!
  * \brief Compute a 3D rotation matrix.
  * \note The implicit ordering is rotation about the x, y, and then z axis.

Common/src/geometry/CGeometry.cpp

@@ -1246,30 +1246,6 @@ void CGeometry::PostPeriodicSends(CGeometry* geometry, const CConfig* config, un
 #endif
 }
 
-su2double CGeometry::Point2Plane_Distance(const su2double* Coord, const su2double* iCoord, const su2double* jCoord,
-                                          const su2double* kCoord) {
-  su2double CrossProduct[3], iVector[3], jVector[3], distance, modulus;
-  unsigned short iDim;
-
-  for (iDim = 0; iDim < 3; iDim++) {
-    iVector[iDim] = jCoord[iDim] - iCoord[iDim];
-    jVector[iDim] = kCoord[iDim] - iCoord[iDim];
-  }
-
-  CrossProduct[0] = iVector[1] * jVector[2] - iVector[2] * jVector[1];
-  CrossProduct[1] = iVector[2] * jVector[0] - iVector[0] * jVector[2];
-  CrossProduct[2] = iVector[0] * jVector[1] - iVector[1] * jVector[0];
-
-  modulus =
-      sqrt(CrossProduct[0] * CrossProduct[0] + CrossProduct[1] * CrossProduct[1] + CrossProduct[2] * CrossProduct[2]);
-
-  distance = 0.0;
-  for (iDim = 0; iDim < 3; iDim++) distance += CrossProduct[iDim] * (Coord[iDim] - iCoord[iDim]);
-  distance /= modulus;
-
-  return distance;
-}
-
 void CGeometry::SetEdges() {
   nEdge = 0;
   for (auto iPoint = 0ul; iPoint < nPoint; iPoint++) {

Common/src/grid_movement/CFreeFormDefBox.cpp

@@ -28,6 +28,7 @@
 #include "../../include/grid_movement/CFreeFormDefBox.hpp"
 #include "../../include/grid_movement/CBezierBlending.hpp"
 #include "../../include/grid_movement/CBSplineBlending.hpp"
+#include "../../include/toolboxes/geometry_toolbox.hpp"
 
 CFreeFormDefBox::CFreeFormDefBox() : CGridMovement() {}
 
@@ -1014,13 +1015,8 @@ su2double* CFreeFormDefBox::GetParametricCoord_Iterative(unsigned long iPoint, s
   return ParamCoord;
 }
 
-bool CFreeFormDefBox::GetPointFFD(CGeometry* geometry, CConfig* config, unsigned long iPoint) const {
-  bool Inside = true;
-  bool cylindrical = (config->GetFFD_CoordSystem() == CYLINDRICAL);
-  bool spherical = (config->GetFFD_CoordSystem() == SPHERICAL);
-  bool polar = (config->GetFFD_CoordSystem() == POLAR);
-
-  /*--- indices of the FFD box. Note that the front face is labelled 0,1,2,3 and the back face is 4,5,6,7 ---*/
+bool CFreeFormDefBox::CheckPointInsideFFD(const su2double* coord) const {
+  /*--- Indices of the FFD box. Note that the front face is labelled 0,1,2,3 and the back face is 4,5,6,7 ---*/
 
   unsigned short Index[6][5] = {{0, 1, 2, 3, 0},   // front side
                                 {1, 5, 6, 2, 1},   // right side
@@ -1029,76 +1025,34 @@ bool CFreeFormDefBox::GetPointFFD(CGeometry* geometry, CConfig* config, unsigned
                                 {4, 5, 1, 0, 4},   // bottom side
                                 {4, 7, 6, 5, 4}};  // back side
 
-  /*--- The current approach is to subdivide each of the 6 faces of the hexahedral FFD box into 4 triangles
-      by defining a supporting middle point. This allows nonplanar FFD boxes.
-      Note that the definition of the FFD box is as follows: the FFD box is a 6-sided die and we are looking at the side
-     "1". The opposite side is side "6". If we are looking at side "1", we define the nodes counterclockwise. If we are
-     looking at side "6", we define the face clockwise  ---*/
-
-  unsigned short nDim = geometry->GetnDim();
-
-  su2double Coord[3] = {0.0, 0.0, 0.0};
-  for (unsigned short iDim = 0; iDim < nDim; iDim++) Coord[iDim] = geometry->nodes->GetCoord(iPoint, iDim);
-
-  su2double X_0, Y_0, Z_0, Xbar, Ybar, Zbar;
-
-  if (cylindrical) {
-    X_0 = config->GetFFD_Axis(0);
-    Y_0 = config->GetFFD_Axis(1);
-    Z_0 = config->GetFFD_Axis(2);
-
-    Xbar = Coord[0] - X_0;
-    Ybar = Coord[1] - Y_0;
-    Zbar = Coord[2] - Z_0;
-
-    Coord[0] = sqrt(Ybar * Ybar + Zbar * Zbar);
-    Coord[1] = atan2(Zbar, Ybar);
-    if (Coord[1] > PI_NUMBER / 2.0) Coord[1] -= 2.0 * PI_NUMBER;
-    Coord[2] = Xbar;
+  /*--- The current approach is to subdivide each of the 6 faces of the hexahedral FFD box into 4 triangles by defining
+   a supporting middle point. This allows nonplanar FFD boxes. Note that the definition of the FFD box is as follows:
+   the FFD box is a 6-sided die and we are looking at the side "1". The opposite side is side "6". If we are looking at
+   side "1", we define the nodes counterclockwise. If we are looking at side "6", we define the face clockwise. ---*/
 
-  }
-
-  else if (spherical || polar) {
-    X_0 = config->GetFFD_Axis(0);
-    Y_0 = config->GetFFD_Axis(1);
-    Z_0 = config->GetFFD_Axis(2);
-
-    Xbar = Coord[0] - X_0;
-    Ybar = Coord[1] - Y_0;
-    Zbar = Coord[2] - Z_0;
-
-    Coord[0] = sqrt(Xbar * Xbar + Ybar * Ybar + Zbar * Zbar);
-    Coord[1] = atan2(Zbar, Ybar);
-    if (Coord[1] > PI_NUMBER / 2.0) Coord[1] -= 2.0 * PI_NUMBER;
-    Coord[2] = acos(Xbar / Coord[0]);
-  }
+  /*--- Loop over the faces of the FFD box. ---*/
 
-  /*--- loop over the faces of the FFD box  ---*/
+  for (unsigned short iFace = 0; iFace < 6; iFace++) {
+    /*--- Every face needs an interpolated middle point for the triangles. ---*/
 
-  for (unsigned short iVar = 0; iVar < 6; iVar++) {
     su2double P[3] = {0.0, 0.0, 0.0};
-
-    /*--- every face needs an interpolated middle point for the triangles ---*/
-
     for (int p = 0; p < 4; p++) {
-      P[0] += 0.25 * Coord_Corner_Points[Index[iVar][p]][0];
-      P[1] += 0.25 * Coord_Corner_Points[Index[iVar][p]][1];
-      P[2] += 0.25 * Coord_Corner_Points[Index[iVar][p]][2];
+      P[0] += 0.25 * Coord_Corner_Points[Index[iFace][p]][0];
+      P[1] += 0.25 * Coord_Corner_Points[Index[iFace][p]][1];
+      P[2] += 0.25 * Coord_Corner_Points[Index[iFace][p]][2];
     }
 
-    /*--- loop over the 4 triangles making up the FFD box. The sign is equal for all distances ---*/
+    /*--- Loop over the 4 triangles making up the FFD box. The sign should be equal for all distances. ---*/
 
-    for (unsigned short jVar = 0; jVar < 4; jVar++) {
-      su2double Distance_Point = geometry->Point2Plane_Distance(Coord, Coord_Corner_Points[Index[iVar][jVar]],
-                                                                Coord_Corner_Points[Index[iVar][jVar + 1]], P);
-      if (Distance_Point < 0) {
-        Inside = false;
-        return Inside;
+    for (int iNode = 0; iNode < 4; iNode++) {
+      const su2double* plane[] = {P, Coord_Corner_Points[Index[iFace][iNode]],
+                                  Coord_Corner_Points[Index[iFace][iNode + 1]]};
+      if (GeometryToolbox::PointToPlaneDistance(plane, coord) < 0) {
+        return false;
       }
     }
   }
-
-  return Inside;
+  return true;
 }
 
 su2double CFreeFormDefBox::GetDerivative1(su2double* uvw, unsigned short val_diff, unsigned short* ijk,