Update of CoDiPack and MeDiPack versions

Common/include/ad_structure.hpp

@@ -148,6 +148,60 @@ namespace AD{
    * for each output variable to the AD tape.
    */
   void EndPreacc();
+
+  /*!
+   * \brief Initializes an externally differentiated function. Input and output variables are set with SetExtFuncIn/SetExtFuncOut
+   * \param[in] storePrimalInput - Specifies whether the primal input values are stored for the reverse call of the external function.
+   * \param[in] storePrimalOutput - Specifies whether the primal output values are stored for the reverse call of the external function.
+   */
+  void StartExtFunc(bool storePrimalInput, bool storePrimalOutput);
+
+  /*!
+   * \brief Sets the scalar input of a externally differentiated function.
+   * \param[in] data - the scalar input variable.
+   */
+  void SetExtFuncIn(su2double &data);
+
+  /*!
+   * \brief Sets the input variables of a externally differentiated function using a 1D array.
+   * \param[in] data - the input 1D array.
+   * \param[in] size - number of rows.
+   */
+  void SetExtFuncIn(const su2double* data, const int size);
+
+  /*!
+  * \brief  Sets the input variables of a externally differentiated function using a 2D array.
+  * \param[in] data - the input 2D array.
+  * \param[in] size_x - number of rows.
+  * \param[in] size_y - number of columns.
+  */
+  void SetExtFuncIn(const su2double* const *data, const int size_x, const int size_y);
+
+  /*!
+   * \brief Sets the scalar output of a externally differentiated function.
+   * \param[in] data - the scalar output variable.
+   */
+  void SetExtFuncOut(su2double &data);
+
+  /*!
+   * \brief Sets the output variables of a externally differentiated function using a 1D array.
+   * \param[in] data - the output 1D array.
+   * \param[in] size - number of rows.
+   */
+  void SetExtFuncOut(su2double* data, const int size);
+
+  /*!
+  * \brief  Sets the output variables of a externally differentiated function using a 2D array.
+  * \param[in] data - the output 2D array.
+  * \param[in] size_x - number of rows.
+  * \param[in] size_y - number of columns.
+  */
+  void SetExtFuncOut(su2double** data, const int size_x, const int size_y);
+
+  /*!
+   * \brief Ends an external function section by deleting the structures.
+   */
+  void EndExtFunc();
 
 }
 

Common/include/ad_structure.inl

@@ -44,6 +44,10 @@ namespace AD{
 
   typedef su2double::TapeType Tape;
 
+  typedef codi::ExternalFunctionHelper<su2double> ExtFuncHelper;
+
+  extern ExtFuncHelper* FuncHelper;
+
   /*--- Stores the indices of the input variables (they might be overwritten) ---*/
 
   extern std::vector<su2double::GradientData> inputValues;
@@ -68,6 +72,8 @@ namespace AD{
 
   extern std::vector<su2double*> localOutputValues;
 
+  extern codi::PreaccumulationHelper<su2double> PreaccHelper;  
+
   inline void RegisterInput(su2double &data) {AD::globalTape.registerInput(data);
                                              inputValues.push_back(data.getGradientData());}
 
@@ -95,7 +101,7 @@ namespace AD{
   inline void SetPreaccIn(const su2double &data) {
     if (PreaccActive) {
       if (data.isActive()) {
-        localInputValues.push_back(data.getGradientData());
+        PreaccHelper.addInput(data);       
       }
     }
   }
@@ -104,7 +110,7 @@ namespace AD{
     if (PreaccActive) {
       for (unsigned short i = 0; i < size; i++) {
         if (data[i].isActive()) {
-          localInputValues.push_back(data[i].getGradientData());
+          PreaccHelper.addInput(data[i]);
         }
       }
     }
@@ -115,7 +121,7 @@ namespace AD{
       for (unsigned short i = 0; i < size_x; i++) {
         for (unsigned short j = 0; j < size_y; j++) {
           if (data[i][j].isActive()) {
-            localInputValues.push_back(data[i][j].getGradientData());
+            PreaccHelper.addInput(data[i][j]);
           }
         }
       }
@@ -124,15 +130,15 @@ namespace AD{
 
   inline void StartPreacc() {
     if (globalTape.isActive() && PreaccEnabled) {
-      StartPosition = globalTape.getPosition();
+      PreaccHelper.start();
       PreaccActive = true;
     }
   }
 
   inline void SetPreaccOut(su2double& data) {
     if (PreaccActive) {
       if (data.isActive()) {
-        localOutputValues.push_back(&data);
+        PreaccHelper.addOutput(data);
       }
     }
   }
@@ -141,7 +147,7 @@ namespace AD{
     if (PreaccActive) {
       for (unsigned short i = 0; i < size; i++) {
         if (data[i].isActive()) {
-          localOutputValues.push_back(&data[i]);
+          PreaccHelper.addOutput(data[i]);
         }
       }
     }
@@ -152,18 +158,79 @@ namespace AD{
       for (unsigned short i = 0; i < size_x; i++) {
         for (unsigned short j = 0; j < size_y; j++) {
           if (data[i][j].isActive()) {
-            localOutputValues.push_back(&data[i][j]);
+            PreaccHelper.addOutput(data[i][j]);
           }
         }
       }
     }
   }
 
+  inline void EndPreacc(){
+    if (PreaccActive) {
+      PreaccHelper.finish(false);
+    }
+  }
+
+  inline void StartExtFunc(bool storePrimalInput, bool storePrimalOutput){
+    FuncHelper = new ExtFuncHelper(true);
+    if (!storePrimalInput){
+      FuncHelper->disableInputPrimalStore();
+    }
+    if (!storePrimalOutput){
+      FuncHelper->disableOutputPrimalStore();
+    }
+  }
+
+  inline void SetExtFuncIn(const su2double &data) {
+    FuncHelper->addInput(data);       
+  }
+
+  inline void SetExtFuncIn(const su2double* data, const int size) {
+    for (int i = 0; i < size; i++) {
+      FuncHelper->addInput(data[i]);
+    }
+
+  }
+
+  inline void SetExtFuncIn(const su2double* const *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++) {
+        FuncHelper->addInput(data[i][j]);
+      }
+    }
+  }
+
+  inline void SetExtFuncOut(su2double& data) {
+    if (globalTape.isActive()) {
+      FuncHelper->addOutput(data);
+    }
+  }
+
+  inline void SetExtFuncOut(su2double* data, const int size) {
+    for (int i = 0; i < size; i++) {
+      if (globalTape.isActive()) {
+        FuncHelper->addOutput(data[i]);
+      }
+    }
+  }
+
+  inline void SetExtFuncOut(su2double** 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()) {
+          FuncHelper->addOutput(data[i][j]);
+        }
+      }
+    }
+  }
 
   inline void delete_handler(void *handler) {
     CheckpointHandler *checkpoint = static_cast<CheckpointHandler*>(handler);
     checkpoint->clear();
   }
+
+  inline void EndExtFunc(){delete FuncHelper;}
+
 #else
 
   /*--- Default implementation if reverse mode is disabled ---*/
@@ -199,6 +266,22 @@ namespace AD{
   inline void StartPreacc() {}
 
   inline void EndPreacc() {}
+
+  inline void StartExtFunc(bool storePrimalInput, bool storePrimalOutput){}
+
+  inline void SetExtFuncIn(const su2double &data) {}
+
+  inline void SetExtFuncIn(const su2double* data, const int size) {}
+
+  inline void SetExtFuncIn(const su2double* const *data, const int size_x, const int size_y) {}
+
+  inline void SetExtFuncOut(su2double& data) {}
+
+  inline void SetExtFuncOut(su2double* data, const int size) {}
+
+  inline void SetExtFuncOut(su2double** data, const int size_x, const int size_y) {}
+
+  inline void EndExtFunc(){}
 #endif
 }
 

Common/include/datatypes/codi_reverse_structure.hpp

@@ -37,7 +37,7 @@
 #pragma once
 
 #include "codi.hpp"
-#include "tools/dataStore.hpp"
+#include "codi/tools/dataStore.hpp"
 
 #ifndef CODI_INDEX_TAPE
 #  define CODI_INDEX_TAPE 0

Common/include/grid_movement_structure.hpp

@@ -53,6 +53,7 @@
 #include "matrix_structure.hpp"
 #include "vector_structure.hpp"
 #include "linear_solvers_structure.hpp"
+#include "linear_solvers_structure_b.hpp"
 #include "element_structure.hpp"
 
 using namespace std;

Common/include/linear_solvers_structure.hpp

@@ -245,17 +245,7 @@ class CSysSolve {
    * \param[in] config - Definition of the particular problem.
    */
   unsigned long Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config);
-
-  /*!
-   * \brief Prepare the linear solve during the reverse interpretation of the AD tape.
-   * \param[in] Jacobian - Jacobian Matrix for the linear system
-   * \param[in] LinSysRes - Linear system residual
-   * \param[in] LinSysSol - Linear system solution
-   * \param[in] geometry -  Geometrical definition of the problem.
-   * \param[in] config - Definition of the particular problem.
-   */
-  void SetExternalSolve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config);
-
+
   /*!
    * \brief Get the final residual.
    * \return The residual at the end of Solve

Common/include/linear_solvers_structure_b.hpp

@@ -43,7 +43,6 @@
 class CSysSolve_b{
 
 public:
-  static void Solve_b(AD::Tape* tape, AD::CheckpointHandler *data);
-  static void Delete_b(AD::Tape* tape, AD::CheckpointHandler *data);
+  static void Solve_b(const codi::RealReverse::Real* x, codi::RealReverse::Real* x_b, size_t m, const codi::RealReverse::Real* y, const codi::RealReverse::Real* y_b, size_t n, codi::DataStore* d);
 };
 #endif

Common/include/mpi_structure.hpp

@@ -60,14 +60,14 @@ class CMediMPIWrapper;
 typedef CMediMPIWrapper SU2_MPI;
 
 #if defined CODI_REVERSE_TYPE
-#include <medi/codiMediPackTypes.hpp>
+#include <codi/externals/codiMediPackTypes.hpp>
 #if CODI_PRIMAL_INDEX_TAPE
 typedef CoDiPackToolPrimalRestore<su2double> MediTool;
 #else
 typedef CoDiPackTool<su2double> MediTool;
 #endif // defined CODI_REVERSE_TYPE
 #elif defined CODI_FORWARD_TYPE
-#include <medi/codiForwardMediPackTypes.hpp>
+#include <codi/externals/codiForwardMediPackTypes.hpp>
 typedef CoDiPackForwardTool<su2double> MediTool;
 #endif // defined CODI_FORWARD_TYPE
 #define AMPI_ADOUBLE ((medi::MpiTypeInterface*)MediTool::MPI_TYPE)

Common/include/mpi_structure.inl

@@ -34,7 +34,7 @@
  * You should have received a copy of the GNU Lesser General Public
  * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
  */
-
+#include "mpi_structure.hpp"
 #pragma once
 
 #ifdef HAVE_MPI

Common/src/ad_structure.cpp

@@ -54,78 +54,9 @@ namespace AD {
   bool PreaccActive = false;
   bool PreaccEnabled = true;
 
-  void EndPreacc() {
+  codi::PreaccumulationHelper<su2double> PreaccHelper;
 
-    if(PreaccActive) {
-      unsigned short iVarOut, iVarIn;
-      unsigned short nVarOut, nVarIn;
-      su2double::GradientData index_out, index_in;
+  ExtFuncHelper* FuncHelper;
 
-      nVarOut = localOutputValues.size();
-      nVarIn  = localInputValues.size();
-
-      /*--- Store the current position of the tape ---*/
-
-      EndPosition = globalTape.getPosition();
-
-      /*--- Allocate local memory ---*/
-
-      passivedouble* local_jacobi     = new passivedouble[nVarOut*nVarIn];
-      unsigned short* nNonzero        = new unsigned short[nVarOut];
-
-      /*--- Compute the local Jacobi matrix of the code between the start and end position
-       * using the inputs and outputs declared with StartPreacc(...)/EndPreacc(...) ---*/
-
-      for (iVarOut = 0; iVarOut < nVarOut; iVarOut++) {
-        nNonzero[iVarOut] = 0;
-        index_out = localOutputValues[iVarOut]->getGradientData();
-
-        globalTape.setGradient(index_out, 1.0);
-        globalTape.evaluate(EndPosition, StartPosition);
-
-        for (iVarIn= 0; iVarIn < nVarIn; iVarIn++) {
-          index_in =  localInputValues[iVarIn];
-          local_jacobi[iVarOut*nVarIn+iVarIn] = globalTape.getGradient(index_in);
-          if (local_jacobi[iVarOut*nVarIn+iVarIn] != 0.0) {
-            nNonzero[iVarOut]++;
-          }
-          globalTape.setGradient(index_in, 0.0);
-        }
-        globalTape.setGradient(index_out, 0.0);
-        globalTape.clearAdjoints(EndPosition, StartPosition);
-      }
-
-      /*--- Reset the tape to the starting position (to reuse the part of the tape) ---*/
-
-      if (nVarOut > 0) {
-        globalTape.reset(StartPosition);
-      }
-
-      /*--- For each output create a statement on the tape and push the corresponding Jacobi entries.
-       * Note that the output variables need a new index since we did a reset of the tape section. ---*/
-
-      for (iVarOut = 0; iVarOut < nVarOut; iVarOut++) {
-        if (nNonzero[iVarOut] != 0){
-          globalTape.store(localOutputValues[iVarOut]->getValue(), localOutputValues[iVarOut]->getGradientData(), nNonzero[iVarOut]);
-          for (iVarIn = 0; iVarIn < nVarIn; iVarIn++) {
-            index_in =  localInputValues[iVarIn];
-           globalTape.pushJacobi(local_jacobi[iVarOut*nVarIn+iVarIn],
-               local_jacobi[iVarOut*nVarIn+iVarIn], local_jacobi[iVarOut*nVarIn+iVarIn], index_in);
-          }
-        }
-      }
-
-      /*--- Clear local vectors and reset indicator ---*/
-
-
-      localInputValues.clear();
-      localOutputValues.clear();
-
-      delete [] local_jacobi;
-      delete [] nNonzero;
-
-      PreaccActive = false;
-    }
-  }
 #endif
 }