Merge branch 'develop' into NEMO_opt_init

.github/workflows/regression.yml

@@ -60,7 +60,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix: 
-        testscript: ['tutorials.py', 'parallel_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py']
+        testscript: ['tutorials.py', 'parallel_regression.py', 'parallel_regression_AD.py', 'serial_regression.py', 'serial_regression_AD.py', 'hybrid_regression.py', 'hybrid_regression_AD.py']
         include:
           - testscript: 'tutorials.py'
             tag: MPI
@@ -74,6 +74,8 @@ jobs:
             tag: NoMPI
           - testscript: 'hybrid_regression.py'
             tag: OMP
+          - testscript: 'hybrid_regression_AD.py'
+            tag: OMP
     steps:
       - name: Download All artifact
         uses: actions/download-artifact@v2

Common/include/basic_types/ad_structure.hpp

@@ -252,7 +252,7 @@ namespace AD{
 
   /*!
    * \brief Start a passive region, i.e. stop recording.
-   * \return True is tape was active.
+   * \return True if tape was active.
    */
   inline bool BeginPassive() { return false; }
 
@@ -262,6 +262,28 @@ namespace AD{
    */
   inline void EndPassive(bool wasActive) {}
 
+  /*!
+   * \brief Pause the use of preaccumulation.
+   * \return True if preaccumulation was active.
+   */
+  inline bool PausePreaccumulation() { return false; }
+
+  /*!
+   * \brief Resume the use of preaccumulation.
+   * \param[in] wasActive - Whether preaccumulation was active before pausing.
+   */
+  inline void ResumePreaccumulation(bool wasActive) {}
+
+  /*!
+   * \brief Begin a hybrid parallel adjoint evaluation mode that assumes an inherently safe reverse path.
+   */
+  inline void StartNoSharedReading() {}
+
+  /*!
+   * \brief End the "no shared reading" adjoint evaluation mode.
+   */
+  inline void EndNoSharedReading() {}
+
 #else
   using CheckpointHandler = codi::DataStore;
 
@@ -271,9 +293,10 @@ namespace AD{
 
   extern ExtFuncHelper* FuncHelper;
 
-  extern bool Status;
-
   extern bool PreaccActive;
+#ifdef HAVE_OPDI
+  SU2_OMP(threadprivate(PreaccActive))
+#endif
 
   extern bool PreaccEnabled;
 
@@ -290,6 +313,9 @@ namespace AD{
   extern std::vector<TapePosition> TapePositions;
 
   extern codi::PreaccumulationHelper<su2double> PreaccHelper;
+#ifdef HAVE_OPDI
+  SU2_OMP(threadprivate(PreaccHelper))
+#endif
 
   /*--- Reference to the tape. ---*/
 
@@ -446,6 +472,7 @@ namespace AD{
   FORCEINLINE void EndPreacc(){
     if (PreaccActive) {
       PreaccHelper.finish(false);
+      PreaccActive = false;
     }
   }
 
@@ -522,6 +549,39 @@ namespace AD{
 
   FORCEINLINE void EndPassive(bool wasActive) { if(wasActive) StartRecording(); }
 
+  FORCEINLINE bool PausePreaccumulation() {
+    const auto current = PreaccEnabled;
+    if (!current) return false;
+    SU2_OMP_BARRIER
+    SU2_OMP_MASTER
+    PreaccEnabled = false;
+    END_SU2_OMP_MASTER
+    SU2_OMP_BARRIER
+    return true;
+  }
+
+  FORCEINLINE void ResumePreaccumulation(bool wasActive) {
+    if (!wasActive) return;
+    SU2_OMP_BARRIER
+    SU2_OMP_MASTER
+    PreaccEnabled = true;
+    END_SU2_OMP_MASTER
+    SU2_OMP_BARRIER
+  }
+
+  FORCEINLINE void StartNoSharedReading() {
+#ifdef HAVE_OPDI
+    opdi::logic->setAdjointAccessMode(opdi::LogicInterface::AdjointAccessMode::Classical);
+    opdi::logic->addReverseBarrier();
+#endif
+  }
+
+  FORCEINLINE void EndNoSharedReading() {
+#ifdef HAVE_OPDI
+    opdi::logic->setAdjointAccessMode(opdi::LogicInterface::AdjointAccessMode::Atomic);
+    opdi::logic->addReverseBarrier();
+#endif
+  }
 #endif // CODI_REVERSE_TYPE
 
 } // namespace AD

Common/include/code_config.hpp

@@ -79,25 +79,15 @@ using su2conditional_t = typename su2conditional<B,T,F>::type;
 #include "codi.hpp"
 #include "codi/tools/dataStore.hpp"
 
-#ifndef CODI_INDEX_TAPE
-#define CODI_INDEX_TAPE 0
-#endif
-#ifndef CODI_PRIMAL_TAPE
-#define CODI_PRIMAL_TAPE 0
-#endif
-#ifndef CODI_PRIMAL_INDEX_TAPE
-#define CODI_PRIMAL_INDEX_TAPE 0
-#endif
-
 #if defined(HAVE_OMP)
 using su2double = codi::RealReverseIndexParallel;
 #else
-#if CODI_INDEX_TAPE
+#if defined(CODI_INDEX_TAPE)
 using su2double = codi::RealReverseIndex;
-#elif CODI_PRIMAL_TAPE
-using su2double = codi::RealReversePrimal;
-#elif CODI_PRIMAL_INDEX_TAPE
-using su2double = codi::RealReversePrimalIndex;
+//#elif defined(CODI_PRIMAL_TAPE)
+//using su2double = codi::RealReversePrimal;
+//#elif defined(CODI_PRIMAL_INDEX_TAPE)
+//using su2double = codi::RealReversePrimalIndex;
 #else
 using su2double = codi::RealReverse;
 #endif

Common/include/linear_algebra/CSysSolve.hpp

@@ -256,6 +256,7 @@ class CSysSolve {
   void HandleTemporariesOut(CSysVector<OtherType>& LinSysSol) {
 
     /*--- Reset the pointers. ---*/
+    SU2_OMP_BARRIER
     SU2_OMP_MASTER {
       LinSysRes_ptr = nullptr;
       LinSysSol_ptr = nullptr;
@@ -276,6 +277,7 @@ class CSysSolve {
     LinSysSol.PassiveCopy(LinSysSol_tmp);
 
     /*--- Reset the pointers. ---*/
+    SU2_OMP_BARRIER
     SU2_OMP_MASTER {
       LinSysRes_ptr = nullptr;
       LinSysSol_ptr = nullptr;

Common/include/toolboxes/graph_toolbox.hpp

@@ -527,7 +527,7 @@ T createNaturalColoring(Index_t numInnerIndexes)
  * \param[out] indexColor - Optional, vector with colors given to the outer indices.
  * \return Coloring in the same type of the input pattern.
  */
-template<class T, typename Color_t = char, size_t MaxColors = 32, size_t MaxMB = 128>
+template<class T, typename Color_t = char, size_t MaxColors = 64, size_t MaxMB = 128>
 T colorSparsePattern(const T& pattern, size_t groupSize = 1, bool balanceColors = false,
                      std::vector<Color_t>* indexColor = nullptr)
 {

Common/include/toolboxes/printing_toolbox.hpp

@@ -71,7 +71,6 @@ namespace PrintingToolbox {
 class CTablePrinter{
 public:
   CTablePrinter(std::ostream * output, const std::string & separator = "|");
-  ~CTablePrinter();
 
   enum alignment {
     CENTER,

Common/src/CConfig.cpp

@@ -4477,11 +4477,7 @@ void CConfig::SetPostprocessing(SU2_COMPONENT val_software, unsigned short val_i
 #if defined CODI_REVERSE_TYPE
   AD_Mode = YES;
 
-#if defined HAVE_OMP
-  AD::PreaccEnabled = false;
-#else
   AD::PreaccEnabled = AD_Preaccumulation;
-#endif
 
 #else
   if (AD_Mode == YES) {

Common/src/basic_types/ad_structure.cpp

@@ -35,9 +35,16 @@ namespace AD {
   std::vector<TapePosition> TapePositions;
 
   bool PreaccActive = false;
+#ifdef HAVE_OPDI
+  SU2_OMP(threadprivate(PreaccActive))
+#endif
+
   bool PreaccEnabled = true;
 
   codi::PreaccumulationHelper<su2double> PreaccHelper;
+#ifdef HAVE_OPDI
+  SU2_OMP(threadprivate(PreaccHelper))
+#endif
 
   ExtFuncHelper* FuncHelper;
 

Common/src/geometry/CPhysicalGeometry.cpp

@@ -7701,7 +7701,8 @@ void CPhysicalGeometry::SetBoundControlVolume(const CConfig *config, unsigned sh
 
       const auto nNodes = bound[iMarker][iElem]->GetnNodes();
 
-      AD::StartPreacc();
+      /*--- Cannot preaccumulate if hybrid parallel due to shared reading. ---*/
+      if (omp_get_num_threads() == 1) AD::StartPreacc();
 
       /*--- Get pointers to the coordinates of all the element nodes ---*/
       array<const su2double*, N_POINTS_MAXIMUM> Coord;

Common/src/linear_algebra/CSysMatrix.cpp

@@ -185,10 +185,17 @@ void CSysMatrix<ScalarType>::Initialize(unsigned long npoint, unsigned long npoi
   /*--- This is akin to the row_ptr. ---*/
   omp_partitions = new unsigned long [omp_num_parts+1];
 
-  /// TODO: Use a work estimate to produce more balanced partitions.
-  auto pts_per_part = roundUpDiv(nPointDomain, omp_num_parts);
-  for(auto part = 0ul; part < omp_num_parts; ++part)
-    omp_partitions[part] = part * pts_per_part;
+  /*--- Work estimate based on non-zeros to produce balanced partitions. ---*/
+
+  const auto row_ptr_prec = ilu_needed? row_ptr_ilu : row_ptr;
+  const auto nnz_prec = row_ptr_prec[nPointDomain];
+
+  const auto nnz_per_part = roundUpDiv(nnz_prec, omp_num_parts);
+
+  for (auto iPoint = 0ul, part = 0ul; iPoint < nPointDomain; ++iPoint) {
+    if (row_ptr_prec[iPoint] >= part*nnz_per_part)
+      omp_partitions[part++] = iPoint;
+  }
   omp_partitions[omp_num_parts] = nPointDomain;
 
   /*--- Generate MKL Kernels ---*/

Common/src/toolboxes/printing_toolbox.cpp

@@ -43,10 +43,6 @@ PrintingToolbox::CTablePrinter::CTablePrinter(std::ostream * output, const std::
   precision_ = 6;
 }
 
-PrintingToolbox::CTablePrinter::~CTablePrinter(){
-
-}
-
 int PrintingToolbox::CTablePrinter::GetNumColumns() const {
   return (int)column_headers_.size();
 }
@@ -102,8 +98,6 @@ void PrintingToolbox::CTablePrinter::PrintHorizontalLine() {
 
 void PrintingToolbox::CTablePrinter::PrintHeader(){
 
-
-
   if (print_header_top_line_) PrintHorizontalLine();
   *out_stream_ << separator_;
   int indent = 0;
@@ -140,4 +134,3 @@ void PrintingToolbox::CTablePrinter::PrintHeader(){
 void PrintingToolbox::CTablePrinter::PrintFooter(){
   PrintHorizontalLine();
 }
-

SU2_CFD/include/gradients/computeGradientsGreenGauss.hpp

@@ -76,7 +76,8 @@ void computeGradientsGreenGauss(CSolver* solver,
   {
     auto nodes = geometry.nodes;
 
-    AD::StartPreacc();
+    /*--- Cannot preaccumulate if hybrid parallel due to shared reading. ---*/
+    if (omp_get_num_threads() == 1) AD::StartPreacc();
     AD::SetPreaccIn(nodes->GetVolume(iPoint));
     AD::SetPreaccIn(nodes->GetPeriodicVolume(iPoint));
 

SU2_CFD/include/gradients/computeGradientsLeastSquares.hpp

@@ -203,7 +203,8 @@ void computeGradientsLeastSquares(CSolver* solver,
     auto nodes = geometry.nodes;
     const auto coord_i = nodes->GetCoord(iPoint);
 
-    AD::StartPreacc();
+    /*--- Cannot preaccumulate if hybrid parallel due to shared reading. ---*/
+    if (omp_get_num_threads() == 1) AD::StartPreacc();
     AD::SetPreaccIn(coord_i, nDim);
 
     for (size_t iVar = varBegin; iVar < varEnd; ++iVar)

SU2_CFD/include/limiters/computeLimiters_impl.hpp

@@ -132,7 +132,8 @@ void computeLimiters_impl(CSolver* solver,
     auto nodes = geometry.nodes;
     const auto coord_i = nodes->GetCoord(iPoint);
 
-    AD::StartPreacc();
+    /*--- Cannot preaccumulate if hybrid parallel due to shared reading. ---*/
+    if (omp_get_num_threads() == 1) AD::StartPreacc();
     AD::SetPreaccIn(coord_i, nDim);
 
     for (size_t iVar = varBegin; iVar < varEnd; ++iVar)

SU2_CFD/include/output/CFVMOutput.hpp

@@ -34,12 +34,7 @@ class CFVMOutput : public COutput{
   /*!
    * \brief Constructor of the class
    */
-  CFVMOutput(CConfig *config, unsigned short nDim, bool femOutput);
-
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CFVMOutput(void) = default;
+  CFVMOutput(const CConfig *config, unsigned short nDim, bool femOutput);
 
   /*!
    * \brief Add Coordinates to output.

SU2_CFD/include/output/CFlowCompOutput.hpp

@@ -47,7 +47,7 @@ class CFlowCompOutput final: public CFlowOutput {
    * \brief Constructor of the class
    * \param[in] config - Definition of the particular problem.
    */
-  CFlowCompOutput(CConfig *config, unsigned short nDim);
+  CFlowCompOutput(const CConfig *config, unsigned short nDim);
 
   /*!
    * \brief Load the history output field values

SU2_CFD/include/output/CFlowOutput.hpp

@@ -38,22 +38,22 @@ class CFlowOutput : public CFVMOutput{
    * \brief Constructor of the class
    * \param[in] config - Definition of the particular problem.
    */
-  CFlowOutput(CConfig *config, unsigned short nDim, bool femOutput);
+  CFlowOutput(const CConfig *config, unsigned short nDim, bool femOutput);
 
   /*!
    * \brief Add flow surface output fields
    * \param[in] config - Definition of the particular problem.
    */
-  void AddAnalyzeSurfaceOutput(CConfig *config);
+  void AddAnalyzeSurfaceOutput(const CConfig *config);
 
   /*!
    * \brief Set flow surface output field values
    * \param[in] solver - The container holding all solution data.
    * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
+   * \param[in,out] config - Definition of the particular problem.
    * \param[in] output - Boolean indicating whether information should be written to screen
    */
-  void SetAnalyzeSurface(CSolver *solver, CGeometry *geometry, CConfig *config, bool output);
+  void SetAnalyzeSurface(const CSolver *solver, const CGeometry *geometry, CConfig *config, bool output);
 
   /*!
    * \brief Add aerodynamic coefficients as output fields
@@ -152,10 +152,9 @@ class CFlowOutput : public CFVMOutput{
   /*!
    * \brief Write the forces breakdown file
    * \param[in] config - Definition of the particular problem per zone.
-   * \param[in] geometry - Geometrical definition of the problem.
-   * \param[in] solver_container - The container holding all solution data.
+   * \param[in] flow_solver - The container holding all solution data.
    */
-  void WriteForcesBreakdown(CConfig *config, CGeometry *geometry, CSolver **solver_container);
+  void WriteForcesBreakdown(const CConfig *config, const CSolver *flow_solver) const;
 
   /*!
    * \brief Set the time averaged output fields.

SU2_CFD/include/output/CMultizoneOutput.hpp

@@ -61,12 +61,7 @@ class CMultizoneOutput final: public COutput {
   /*!
    * \brief Constructor of the class.
    */
-  CMultizoneOutput(CConfig *driver_config, CConfig** config, unsigned short nDim);
-
-  /*!
-   * \brief Destructor of the class.
-   */
-  ~CMultizoneOutput(void) override;
+  CMultizoneOutput(const CConfig *driver_config, const CConfig* const* config, unsigned short nDim);
 
   /*!
    * \brief Load the multizone history output field values