Hybrid Parallel AD (Part 1/?)

.github/workflows/regression.yml

@@ -16,7 +16,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix: 
-        config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, BaseOMP]
+        config_set: [BaseMPI, ReverseMPI, ForwardMPI, BaseNoMPI, ReverseNoMPI, ForwardNoMPI, BaseOMP, ReverseOMP, ForwardOMP]
         include:
           - config_set: BaseMPI
             flags: '-Denable-pywrapper=true -Denable-tests=true --warnlevel=3 --werror'
@@ -32,6 +32,10 @@ jobs:
             flags: '-Denable-directdiff=true -Denable-normal=false -Dwith-mpi=disabled -Denable-tests=true --warnlevel=3 --werror'
           - config_set: BaseOMP
             flags: '-Dwith-omp=true -Denable-mixedprec=true -Denable-tecio=false --warnlevel=3 --werror'
+          - config_set: ReverseOMP
+            flags: '-Denable-autodiff=true -Denable-normal=false -Dwith-omp=true -Denable-mixedprec=true -Denable-tecio=false --warnlevel=3 --werror'
+          - config_set: ForwardOMP
+            flags: '-Denable-directdiff=true -Denable-normal=false -Dwith-omp=true -Denable-mixedprec=true -Denable-tecio=false --warnlevel=3 --werror'
     runs-on: ubuntu-latest
     steps:
       - name: Cache Object Files

.gitmodules

@@ -15,3 +15,6 @@
 [submodule "subprojects/Mutationpp"]
         path = subprojects/Mutationpp
         url = https://github.com/mutationpp/Mutationpp.git	
+[submodule "externals/opdi"]
+	path = externals/opdi
+	url = https://github.com/SciCompKL/OpDiLib

Common/include/basic_types/ad_structure.hpp

@@ -1,7 +1,7 @@
 /*!
  * \file ad_structure.hpp
  * \brief Main routines for the algorithmic differentiation (AD) structure.
- * \author T. Albring
+ * \author T. Albring, J. Blühdorn
  * \version 7.1.1 "Blackbird"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -27,7 +27,8 @@
 
 #pragma once
 
-#include "datatype_structure.hpp"
+#include "../code_config.hpp"
+#include "../parallelization/omp_structure.hpp"
 
 /*!
  * \namespace AD
@@ -278,62 +279,92 @@ namespace AD{
 
   extern int adjointVectorPosition;
 
-  /*--- Reference to the tape ---*/
-
-  extern su2double::TapeType& globalTape;
-
   extern bool Status;
 
   extern bool PreaccActive;
 
   extern bool PreaccEnabled;
 
-  extern su2double::TapeType::Position StartPosition, EndPosition;
+#ifdef HAVE_OPDI
+  using CoDiTapePosition = su2double::TapeType::Position;
+  using OpDiState = void*;
+  using TapePosition = std::pair<CoDiTapePosition, OpDiState>;
+#else
+  using TapePosition = su2double::TapeType::Position;
+#endif
+
+  extern TapePosition StartPosition, EndPosition;
 
-  extern std::vector<su2double::TapeType::Position> TapePositions;
+  extern std::vector<TapePosition> TapePositions;
 
   extern std::vector<su2double::GradientData> localInputValues;
 
   extern std::vector<su2double*> localOutputValues;
 
   extern codi::PreaccumulationHelper<su2double> PreaccHelper;
 
+  /*--- Reference to the tape. ---*/
+
+  FORCEINLINE su2double::TapeType& getGlobalTape() {
+    return su2double::getGlobalTape();
+  }
+
   FORCEINLINE void RegisterInput(su2double &data, bool push_index = true) {
-    AD::globalTape.registerInput(data);
+    AD::getGlobalTape().registerInput(data);
     if (push_index) {
       inputValues.push_back(data.getGradientData());
     }
   }
 
-  FORCEINLINE void RegisterOutput(su2double& data) {AD::globalTape.registerOutput(data);}
+  FORCEINLINE void RegisterOutput(su2double& data) {AD::getGlobalTape().registerOutput(data);}
 
   FORCEINLINE void ResetInput(su2double &data) {data.getGradientData() = su2double::GradientData();}
 
-  FORCEINLINE void StartRecording() {AD::globalTape.setActive();}
+  FORCEINLINE void StartRecording() {AD::getGlobalTape().setActive();}
 
-  FORCEINLINE void StopRecording() {AD::globalTape.setPassive();}
+  FORCEINLINE void StopRecording() {AD::getGlobalTape().setPassive();}
 
-  FORCEINLINE bool TapeActive() { return AD::globalTape.isActive(); }
+  FORCEINLINE bool TapeActive() { return AD::getGlobalTape().isActive(); }
 
-  FORCEINLINE void PrintStatistics() {AD::globalTape.printStatistics();}
+  FORCEINLINE void PrintStatistics() {AD::getGlobalTape().printStatistics();}
 
-  FORCEINLINE void ClearAdjoints() {AD::globalTape.clearAdjoints(); }
+  FORCEINLINE void ClearAdjoints() {AD::getGlobalTape().clearAdjoints(); }
 
-  FORCEINLINE void ComputeAdjoint() {AD::globalTape.evaluate(); adjointVectorPosition = 0;}
+  FORCEINLINE void ComputeAdjoint() {
+  #if defined(HAVE_OPDI)
+    opdi::logic->prepareEvaluate();
+  #endif
+    AD::getGlobalTape().evaluate();
+    adjointVectorPosition = 0;
+  }
 
   FORCEINLINE void ComputeAdjoint(unsigned short enter, unsigned short leave) {
-    AD::globalTape.evaluate(TapePositions[enter], TapePositions[leave]);
+  #if defined(HAVE_OPDI)
+    opdi::logic->recoverState(TapePositions[enter].second);
+    opdi::logic->prepareEvaluate();
+    AD::getGlobalTape().evaluate(TapePositions[enter].first, TapePositions[leave].first);
+  #else
+    AD::getGlobalTape().evaluate(TapePositions[enter], TapePositions[leave]);
+  #endif
     if (leave == 0)
       adjointVectorPosition = 0;
   }
 
   FORCEINLINE void Reset() {
-    globalTape.reset();
+    AD::getGlobalTape().reset();
+  #if defined(HAVE_OPDI)
+    opdi::logic->reset();
+  #endif
     if (inputValues.size() != 0) {
       adjointVectorPosition = 0;
       inputValues.clear();
     }
     if (TapePositions.size() != 0) {
+    #if defined(HAVE_OPDI)
+      for (TapePosition& pos : TapePositions) {
+        opdi::logic->freeState(pos.second);
+      }
+    #endif
       TapePositions.clear();
     }
   }
@@ -343,11 +374,11 @@ namespace AD{
   }
 
   FORCEINLINE void SetDerivative(int index, const double val) {
-    AD::globalTape.setGradient(index, val);
+    AD::getGlobalTape().setGradient(index, val);
   }
 
   FORCEINLINE double GetDerivative(int index) {
-    return AD::globalTape.getGradient(index);
+    return AD::getGlobalTape().getGradient(index);
   }
 
   /*--- Base case for parameter pack expansion. ---*/
@@ -361,6 +392,11 @@ namespace AD{
     SetPreaccIn(moreData...);
   }
 
+  template<class T, class... Ts, su2enable_if<std::is_same<T,su2double>::value> = 0>
+  FORCEINLINE void SetPreaccIn(T&& data, Ts&&... moreData) {
+    static_assert(!std::is_same<T,su2double>::value, "rvalues cannot be registered");
+  }
+
   template<class T>
   FORCEINLINE void SetPreaccIn(const T& data, const int size) {
     if (PreaccActive) {
@@ -384,20 +420,8 @@ namespace AD{
     }
   }
 
-  template<class T>
-  FORCEINLINE void SetPreaccIn(const T& data, const int size_x, const int size_y, const int size_z) {
-    if (!PreaccActive) return;
-    for (int i = 0; i < size_x; i++) {
-      for (int j = 0; j < size_y; j++) {
-        for (int k = 0; k < size_z; k++) {
-          if (data[i][j][k].isActive()) PreaccHelper.addInput(data[i][j][k]);
-        }
-      }
-    }
-  }
-
   FORCEINLINE void StartPreacc() {
-    if (globalTape.isActive() && PreaccEnabled) {
+    if (AD::getGlobalTape().isActive() && PreaccEnabled) {
       PreaccHelper.start();
       PreaccActive = true;
     }
@@ -438,7 +462,11 @@ namespace AD{
   }
 
   FORCEINLINE void Push_TapePosition() {
-    TapePositions.push_back(AD::globalTape.getPosition());
+  #if defined(HAVE_OPDI)
+    TapePositions.push_back({AD::getGlobalTape().getPosition(), opdi::logic->exportState()});
+  #else
+    TapePositions.push_back(AD::getGlobalTape().getPosition());
+  #endif
   }
 
   FORCEINLINE void EndPreacc(){
@@ -478,15 +506,15 @@ namespace AD{
   }
 
   FORCEINLINE void SetExtFuncOut(su2double& data) {
-    if (globalTape.isActive()) {
+    if (AD::getGlobalTape().isActive()) {
       FuncHelper->addOutput(data);
     }
   }
 
   template<class T>
   FORCEINLINE void SetExtFuncOut(T&& data, const int size) {
     for (int i = 0; i < size; i++) {
-      if (globalTape.isActive()) {
+      if (AD::getGlobalTape().isActive()) {
         FuncHelper->addOutput(data[i]);
       }
     }
@@ -496,7 +524,7 @@ namespace AD{
   FORCEINLINE void SetExtFuncOut(T&& data, const int size_x, const int size_y) {
     for (int i = 0; i < size_x; i++) {
       for (int j = 0; j < size_y; j++) {
-        if (globalTape.isActive()) {
+        if (AD::getGlobalTape().isActive()) {
           FuncHelper->addOutput(data[i][j]);
         }
       }
@@ -511,7 +539,7 @@ namespace AD{
   FORCEINLINE void EndExtFunc() { delete FuncHelper; }
 
   FORCEINLINE bool BeginPassive() {
-    if(AD::globalTape.isActive()) {
+    if(AD::getGlobalTape().isActive()) {
       StopRecording();
       return true;
     }

Common/include/basic_types/datatype_structure.hpp

@@ -30,87 +30,10 @@
 #include <iostream>
 #include <complex>
 #include <cstdio>
-#include <type_traits>
-
-#if defined(_MSC_VER)
-#define FORCEINLINE __forceinline
-#elif defined(__GNUC__) || defined(__clang__) || defined(__INTEL_COMPILER)
-#define FORCEINLINE inline __attribute__((always_inline))
-#else
-#define FORCEINLINE inline
-#endif
-
-#if defined(__GNUC__) || defined(__clang__) || defined(__INTEL_COMPILER)
-#define NEVERINLINE inline __attribute__((noinline))
-#else
-#define NEVERINLINE inline
-#endif
-
-#if defined(__INTEL_COMPILER)
-/*--- Disable warnings related to inline attributes. ---*/
-#pragma warning disable 2196
-#pragma warning disable 3415
-/*--- Disable warnings related to overloaded virtual. ---*/
-#pragma warning disable 654
-#pragma warning disable 1125
-#if defined(CODI_FORWARD_TYPE) || defined(CODI_REVERSE_TYPE)
-#pragma warning disable 1875
-#endif
-#endif
-
-/*--- Convenience SFINAE typedef to conditionally
- * enable/disable function template overloads. ---*/
-template<bool condition>
-using su2enable_if = typename std::enable_if<condition,bool>::type;
-
-/*--- Depending on the datatype defined during the configuration,
- * include the correct definition, and create the main typedef. ---*/
-
-#if defined(CODI_REVERSE_TYPE) // reverse mode AD
-#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 CODI_INDEX_TAPE
-using su2double = codi::RealReverseIndex;
-#elif CODI_PRIMAL_TAPE
-using su2double = codi::RealReversePrimal;
-#elif CODI_PRIMAL_INDEX_TAPE
-using su2double = codi::RealReversePrimalIndex;
-#else
-using su2double = codi::RealReverse;
-#endif
-
-#elif defined(CODI_FORWARD_TYPE) // forward mode AD
-#include "codi.hpp"
-using su2double = codi::RealForward;
-
-#else // primal / direct / no AD
-using su2double = double;
-#endif
 
+#include "../code_config.hpp"
 #include "ad_structure.hpp"
 
-/*--- This type can be used for (rare) compatiblity cases or for
- * computations that are intended to be (always) passive. ---*/
-using passivedouble = double;
-
-/*--- Define a type for potentially lower precision operations. ---*/
-#ifdef USE_MIXED_PRECISION
-using su2mixedfloat = float;
-#else
-using su2mixedfloat = passivedouble;
-#endif
-
 /*!
  * \namespace SU2_TYPE
  * \brief Namespace for defining the datatype wrapper routines, this acts as a base
@@ -174,11 +97,11 @@ namespace SU2_TYPE {
 
 #ifdef CODI_REVERSE_TYPE
   FORCEINLINE passivedouble GetSecondary(const su2double& data) {
-    return AD::globalTape.getGradient(AD::inputValues[AD::adjointVectorPosition++]);
+    return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]);
   }
 
   FORCEINLINE passivedouble GetDerivative(const su2double& data) {
-    return AD::globalTape.getGradient(AD::inputValues[AD::adjointVectorPosition++]);
+    return AD::getGlobalTape().getGradient(AD::inputValues[AD::adjointVectorPosition++]);
   }
 #else // forward
   FORCEINLINE passivedouble GetSecondary(const su2double& data) {return data.getGradient();}