replace AD_BEGIN_PASSIVE macro by a function

pcarruscag · pcarruscag · commit f7748994632d · 2020-06-02T14:50:54.000+01:00
diff --git a/Common/include/basic_types/ad_structure.hpp b/Common/include/basic_types/ad_structure.hpp
@@ -249,6 +249,16 @@ namespace AD{
    */
   inline void Push_TapePosition() {}
 
+  /*!
+   * \brief Start a passive region, i.e. stop recording.
+   */
+  inline void BeginPassive() {}
+
+  /*!
+   * \brief End a passive region, i.e. start recording if we were recording before.
+   */
+  inline void EndPassive() {}
+
 #else
   using CheckpointHandler = codi::DataStore;
 
@@ -488,27 +498,23 @@ namespace AD{
 
   FORCEINLINE void EndExtFunc() { delete FuncHelper; }
 
-#endif // CODI_REVERSE_TYPE
+  FORCEINLINE void BeginPassive() {
+    if(AD::globalTape.isActive()) {
+      AD::globalTape.setPassive();
+      AD::Status = true;
+    }
+  }
 
-} // namespace AD
+  FORCEINLINE void EndPassive() {
+    if(AD::Status) {
+     AD::globalTape.setActive();
+     AD::Status = false;
+    }
+  }
 
-/*--- Macro to begin and end sections with a passive tape ---*/
+#endif // CODI_REVERSE_TYPE
 
-#ifdef CODI_REVERSE_TYPE
-#define AD_BEGIN_PASSIVE         \
-  if(AD::globalTape.isActive()) {\
-     AD::globalTape.setPassive();\
-     AD::Status = true;          \
-  }
-#define AD_END_PASSIVE           \
-  if(AD::Status) {               \
-     AD::globalTape.setActive(); \
-     AD::Status = false;         \
-  }
-#else
-#define AD_BEGIN_PASSIVE
-#define AD_END_PASSIVE
-#endif
+} // namespace AD
 
 
 /*--- If we compile under OSX we have to overload some of the operators for
diff --git a/Common/include/linear_algebra/CSysSolve_b.hpp b/Common/include/linear_algebra/CSysSolve_b.hpp
@@ -27,6 +27,8 @@
 
 #pragma once
 
+#include "../basic_types/datatype_structure.hpp"
+
 #ifdef CODI_REVERSE_TYPE
 template<class ScalarType>
 struct CSysSolve_b {
diff --git a/Common/src/adt_structure.cpp b/Common/src/adt_structure.cpp
@@ -332,7 +332,7 @@ void CADTPointsOnlyClass::DetermineNearestNode_impl(vector<unsigned long>& front
                                                     unsigned long   &pointID,
                                                     int             &rankID) const {
 
-  AD_BEGIN_PASSIVE
+  AD::BeginPassive();
 
   /*--------------------------------------------------------------------------*/
   /*--- Step 1: Initialize the nearest node to the central node of the     ---*/
@@ -448,7 +448,7 @@ void CADTPointsOnlyClass::DetermineNearestNode_impl(vector<unsigned long>& front
     if(frontLeaves.size() == 0) break;
   }
 
-  AD_END_PASSIVE
+  AD::EndPassive();
 
   /* Recompute the distance to get the correct dependency if we use AD */
   coorTarget = coorPoints.data() + nDimADT*minIndex;
@@ -790,7 +790,7 @@ void CADTElemClass::DetermineNearestElement_impl(vector<CBBoxTargetClass>& BBoxT
                                                  unsigned long   &elemID,
                                                  int             &rankID) const {
 
-  AD_BEGIN_PASSIVE
+  AD::BeginPassive();
 
   /*----------------------------------------------------------------------------*/
   /*--- Step 1: Initialize the distance (squared) to the quaranteed distance ---*/
@@ -953,7 +953,7 @@ void CADTElemClass::DetermineNearestElement_impl(vector<CBBoxTargetClass>& BBoxT
     }
   }
 
-  AD_END_PASSIVE
+  AD::EndPassive();
 
   /* At the moment the square of the distance is stored in dist. Compute
      the correct value. */
diff --git a/Common/src/geometry/CGeometry.cpp b/Common/src/geometry/CGeometry.cpp
@@ -1890,7 +1890,7 @@ void CGeometry::ComputeAirfoil_Section(su2double *Plane_P0, su2double *Plane_Nor
                                        vector<su2double> &Zcoord_Airfoil, vector<su2double> &Variable_Airfoil,
                                        bool original_surface, CConfig *config) {
 
-  AD_BEGIN_PASSIVE
+  AD::BeginPassive();
 
   unsigned short iMarker, iNode, jNode, iDim, Index = 0;
   bool intersect;
@@ -2609,7 +2609,7 @@ void CGeometry::ComputeAirfoil_Section(su2double *Plane_P0, su2double *Plane_Nor
 
   }
 
-  AD_END_PASSIVE
+  AD::EndPassive();
 
 }
 
diff --git a/Common/src/geometry/CPhysicalGeometry.cpp b/Common/src/geometry/CPhysicalGeometry.cpp
@@ -7264,7 +7264,7 @@ void CPhysicalGeometry::SetMaxLength(CConfig* config) {
      * neighbor, then enable it again and recompute the distance.
      * This reduces the overhead of storing irrelevant computations. ---*/
 
-    AD_BEGIN_PASSIVE
+    AD::BeginPassive();
 
     su2double max_delta=0;
     auto max_neighbor = iPoint;
@@ -7285,7 +7285,7 @@ void CPhysicalGeometry::SetMaxLength(CConfig* config) {
       }
     }
 
-    AD_END_PASSIVE
+    AD::EndPassive();
 
     /*--- Recompute and set. ---*/
     const su2double* Coord_j = nodes->GetCoord(max_neighbor);
diff --git a/Common/src/interface_interpolation/CIsoparametric.cpp b/Common/src/interface_interpolation/CIsoparametric.cpp
@@ -384,8 +384,8 @@ int CIsoparametric::QuadrilateralIsoparameters(const su2double X[][3], const su2
 
   /*--- Finding Xi and Eta is a "third order" effect that we do not
    *    differentiate (also because we need to iterate). ---*/
-  const bool tapeActive = AD::TapeActive();
-  AD::StopRecording();
+
+  AD::BeginPassive();
 
   for (int iter = 0; iter < NITER; ++iter) {
 
@@ -431,7 +431,7 @@ int CIsoparametric::QuadrilateralIsoparameters(const su2double X[][3], const su2
     if (eps < tol) break;
   }
 
-  if (tapeActive) AD::StartRecording();
+  AD::EndPassive();
 
   int outOfBounds = 0;
 
diff --git a/Common/src/linear_algebra/CSysSolve_b.cpp b/Common/src/linear_algebra/CSysSolve_b.cpp
@@ -25,7 +25,6 @@
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
 
-#include "../../include/CConfig.hpp"
 #include "../../include/linear_algebra/CSysSolve_b.hpp"
 #include "../../include/linear_algebra/CSysSolve.hpp"
 #include "../../include/linear_algebra/CSysMatrix.hpp"
diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp
@@ -40,7 +40,12 @@
  * CSysVector and do not have the same work scheduling must use a
  * SU2_OMP_BARRIER before using the vector.
  */
+#ifndef _MSC_VER
+#define PARALLEL_FOR SU2_OMP(for simd schedule(static,omp_chunk_size) nowait)
+#else
+/*--- MSVC does not support SIMD by default. ---*/
 #define PARALLEL_FOR SU2_OMP(for schedule(static,omp_chunk_size) nowait)
+#endif
 
 template<class ScalarType>
 CSysVector<ScalarType>::CSysVector(void) {
diff --git a/SU2_CFD/include/limiters/computeLimiters_impl.hpp b/SU2_CFD/include/limiters/computeLimiters_impl.hpp
@@ -98,12 +98,9 @@ void computeLimiters_impl(CSolver* solver,
                      omp_get_max_threads(), OMP_MAX_CHUNK);
 #endif
 
-  bool tapeActive = false;
-
+  /*--- If limiters are frozen do not record the computation ---*/
   if (config.GetDiscrete_Adjoint() && config.GetFrozen_Limiter_Disc()) {
-    /*--- If limiters are frozen do not record the computation ---*/
-    tapeActive = AD::TapeActive();
-    AD::StopRecording();
+    AD::BeginPassive();
   }
 
   CLimiterDetails<LimiterKind> limiterDetails;
@@ -240,6 +237,6 @@ void computeLimiters_impl(CSolver* solver,
     solver->CompleteComms(&geometry, &config, kindMpiComm);
   }
 
-  if (tapeActive) AD::StartRecording();
+  AD::EndPassive();
 
 }
diff --git a/SU2_CFD/include/solvers/CMeshSolver.hpp b/SU2_CFD/include/solvers/CMeshSolver.hpp
@@ -39,14 +39,14 @@ class CMeshSolver final : public CFEASolver {
 
   bool stiffness_set;          /*!< \brief Element-based stiffness is set. */
 
-  su2double MinVolume_Ref,     /*!< \brief Minimum volume in reference and current (deformed) configuration. */
-            MinVolume_Curr;
+  unsigned long ElemCounter;   /*!< \brief Error (negative volume) counter. */
 
-  su2double MaxVolume_Ref,
-            MaxVolume_Curr;
-
-  su2double MinDistance;
-  su2double MaxDistance;
+  /*!
+   * \brief Minimum/Maximum distance and volume (in reference and current (deformed) coords).
+   */
+  su2double MinVolume, MinVolume_Ref, MinVolume_Curr;
+  su2double MaxVolume, MaxVolume_Ref, MaxVolume_Curr;
+  su2double MinDistance, MaxDistance;
 
   vector<CMeshElement> element; /*!< \brief Vector which stores element information for each problem. */
 
diff --git a/SU2_CFD/src/iteration_structure.cpp b/SU2_CFD/src/iteration_structure.cpp
@@ -231,12 +231,10 @@ void CIteration::SetMesh_Deformation(CGeometry **geometry,
 
   /*--- Perform the elasticity mesh movement ---*/
 
-  const bool ActiveTape = AD::TapeActive();
-
   if ((kind_recording != MESH_DEFORM) && !config->GetMultizone_Problem()) {
     /*--- In a primal run, AD::TapeActive returns a false ---*/
     /*--- In any other recordings, the tape is passive during the deformation. ---*/
-    AD::StopRecording();
+    AD::BeginPassive();
   }
 
   /*--- Set the stiffness of each element mesh into the mesh numerics ---*/
@@ -248,7 +246,7 @@ void CIteration::SetMesh_Deformation(CGeometry **geometry,
   solver[MESH_SOL]->DeformMesh(geometry, numerics[MESH_SOL], config);
 
   /*--- Continue recording. ---*/
-  if (ActiveTape) AD::StartRecording();
+  AD::EndPassive();
 
 }
 
diff --git a/SU2_CFD/src/numerics/CNumerics.cpp b/SU2_CFD/src/numerics/CNumerics.cpp
@@ -334,7 +334,7 @@ void CNumerics::GetInviscidIncProjFlux(const su2double *val_density,
 void CNumerics::GetInviscidProjJac(const su2double *val_velocity, const su2double *val_energy,
                                    const su2double *val_normal, su2double val_scale,
                                    su2double **val_Proj_Jac_Tensor) const {
-  AD_BEGIN_PASSIVE
+  AD::BeginPassive();
   unsigned short iDim, jDim;
   su2double sqvel, proj_vel, phi, a1, a2;
 
@@ -365,15 +365,15 @@ void CNumerics::GetInviscidProjJac(const su2double *val_velocity, const su2doubl
   for (iDim = 0; iDim < nDim; iDim++)
     val_Proj_Jac_Tensor[nDim+1][iDim+1] = val_scale*(val_normal[iDim]*a1-a2*val_velocity[iDim]*proj_vel);
   val_Proj_Jac_Tensor[nDim+1][nDim+1] = val_scale*Gamma*proj_vel;
-  AD_END_PASSIVE
+  AD::EndPassive();
 }
 
 
 void CNumerics::GetInviscidProjJac(const su2double *val_velocity, const su2double *val_enthalpy,
                                    const su2double *val_chi, const su2double *val_kappa,
                                    const su2double *val_normal, su2double val_scale,
                                    su2double **val_Proj_Jac_Tensor) const {
-  AD_BEGIN_PASSIVE
+  AD::BeginPassive();
   unsigned short iDim, jDim;
   su2double sqvel, proj_vel, phi, a1, a2;
 
@@ -404,15 +404,15 @@ void CNumerics::GetInviscidProjJac(const su2double *val_velocity, const su2doubl
   for (iDim = 0; iDim < nDim; iDim++)
     val_Proj_Jac_Tensor[nDim+1][iDim+1] = val_scale*(val_normal[iDim]*a1-a2*val_velocity[iDim]*proj_vel);
   val_Proj_Jac_Tensor[nDim+1][nDim+1] = val_scale*(a2+1)*proj_vel;
-  AD_END_PASSIVE
+  AD::EndPassive();
 }
 
 void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2double *val_velocity,
                                       const su2double *val_betainc2, const su2double *val_cp,
                                       const su2double *val_temperature, const su2double *val_dRhodT,
                                       const su2double *val_normal, su2double val_scale,
                                       su2double **val_Proj_Jac_Tensor) const {
-  AD_BEGIN_PASSIVE
+  AD::BeginPassive();
   unsigned short iDim;
   su2double proj_vel;
 
@@ -475,7 +475,7 @@ void CNumerics::GetInviscidIncProjJac(const su2double *val_density, const su2dou
     val_Proj_Jac_Tensor[4][4] = val_scale*((*val_cp)*((*val_temperature)*(*val_dRhodT) + (*val_density))*proj_vel);
 
   }
-  AD_END_PASSIVE
+  AD::EndPassive();
 }
 
 void CNumerics::GetPreconditioner(const su2double *val_density, const su2double *val_velocity,
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -1081,8 +1081,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, c
   const su2double simp_minstiff = config->GetSIMP_MinStiffness();
 
   /*--- Never record this method as the mass matrix is passive (but the mass residual is not). ---*/
-  const bool ActiveTape = AD::TapeActive();
-  AD::StopRecording();
+  AD::BeginPassive();
 
   /*--- Start OpenMP parallel region. ---*/
 
@@ -1157,7 +1156,7 @@ void CFEASolver::Compute_MassMatrix(CGeometry *geometry, CNumerics **numerics, c
 
   } // end SU2_OMP_PARALLEL
 
-  if (ActiveTape) AD::StartRecording();
+  AD::EndPassive();
 
 }
 
@@ -1330,8 +1329,7 @@ void CFEASolver::Compute_NodalStressRes(CGeometry *geometry, CNumerics **numeric
 void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics, const CConfig *config) {
 
   /*--- Never record this method as atm it is not differentiable. ---*/
-  const bool ActiveTape = AD::TapeActive();
-  AD::StopRecording();
+  AD::BeginPassive();
 
   const bool prestretch_fem = config->GetPrestretch();
 
@@ -1610,7 +1608,7 @@ void CFEASolver::Compute_NodalStress(CGeometry *geometry, CNumerics **numerics,
 
   }
 
-  if (ActiveTape) AD::StartRecording();
+  AD::EndPassive();
 
 }
 
diff --git a/SU2_CFD/src/solvers/CMeshSolver.cpp b/SU2_CFD/src/solvers/CMeshSolver.cpp

Original file line number	Diff line number	Diff line change
`@@ -1890,7 +1890,7 @@ void CGeometry::ComputeAirfoil_Section(su2double Plane_P0, su2double Plane_Nor`
`1890`	`1890`	`vector<su2double> &Zcoord_Airfoil, vector<su2double> &Variable_Airfoil,`
`1891`	`1891`	`bool original_surface, CConfig *config) {`
`1892`	`1892`
`1893`		`- AD_BEGIN_PASSIVE`
	`1893`	`+ AD::BeginPassive();`
`1894`	`1894`
`1895`	`1895`	`unsigned short iMarker, iNode, jNode, iDim, Index = 0;`
`1896`	`1896`	`bool intersect;`
`@@ -2609,7 +2609,7 @@ void CGeometry::ComputeAirfoil_Section(su2double Plane_P0, su2double Plane_Nor`
`2609`	`2609`
`2610`	`2610`	`}`
`2611`	`2611`
`2612`		`- AD_END_PASSIVE`
	`2612`	`+ AD::EndPassive();`
`2613`	`2613`
`2614`	`2614`	`}`
`2615`	`2615`
Original file line number	Diff line number	Diff line change
`@@ -7264,7 +7264,7 @@ void CPhysicalGeometry::SetMaxLength(CConfig* config) {`
`7264`	`7264`	`* neighbor, then enable it again and recompute the distance.`
`7265`	`7265`	`* This reduces the overhead of storing irrelevant computations. ---*/`
`7266`	`7266`
`7267`		`- AD_BEGIN_PASSIVE`
	`7267`	`+ AD::BeginPassive();`
`7268`	`7268`
`7269`	`7269`	`su2double max_delta=0;`
`7270`	`7270`	`auto max_neighbor = iPoint;`
`@@ -7285,7 +7285,7 @@ void CPhysicalGeometry::SetMaxLength(CConfig* config) {`
`7285`	`7285`	`}`
`7286`	`7286`	`}`
`7287`	`7287`
`7288`		`- AD_END_PASSIVE`
	`7288`	`+ AD::EndPassive();`
`7289`	`7289`
`7290`	`7290`	`/--- Recompute and set. ---/`
`7291`	`7291`	`const su2double* Coord_j = nodes->GetCoord(max_neighbor);`
Original file line number	Diff line number	Diff line change
`@@ -98,12 +98,9 @@ void computeLimiters_impl(CSolver* solver,`
`98`	`98`	`omp_get_max_threads(), OMP_MAX_CHUNK);`
`99`	`99`	`#endif`
`100`	`100`
`101`		`- bool tapeActive = false;`
`102`		`-`
	`101`	`+ /--- If limiters are frozen do not record the computation ---/`
`103`	`102`	`if (config.GetDiscrete_Adjoint() && config.GetFrozen_Limiter_Disc()) {`
`104`		`- /--- If limiters are frozen do not record the computation ---/`
`105`		`- tapeActive = AD::TapeActive();`
`106`		`- AD::StopRecording();`
	`103`	`+ AD::BeginPassive();`
`107`	`104`	`}`
`108`	`105`
`109`	`106`	`CLimiterDetails<LimiterKind> limiterDetails;`
`@@ -240,6 +237,6 @@ void computeLimiters_impl(CSolver* solver,`
`240`	`237`	`solver->CompleteComms(&geometry, &config, kindMpiComm);`
`241`	`238`	`}`
`242`	`239`
`243`		`- if (tapeActive) AD::StartRecording();`
	`240`	`+ AD::EndPassive();`
`244`	`241`
`245`	`242`	`}`