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,5 +266,21 @@ 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

@@ -212,28 +212,6 @@ class CSysSolve {
    */
   unsigned long Solve(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 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_Mesh(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config);
-
-
 };
 
 #include "linear_solvers_structure.inl"

Common/include/linear_solvers_structure_b.hpp

@@ -43,8 +43,6 @@
 class CSysSolve_b{
 
 public:
-  static void Solve_b(AD::Tape* tape, AD::CheckpointHandler *data);
-  static void Solve_g(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