diff --git a/Common/include/ad_structure.hpp b/Common/include/ad_structure.hpp
index c97647652ccb..512d81d21447 100644
--- a/Common/include/ad_structure.hpp
+++ b/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();
 
 }
 
diff --git a/Common/include/ad_structure.inl b/Common/include/ad_structure.inl
index 81be57c2aac8..0f741bb85913 100644
--- a/Common/include/ad_structure.inl
+++ b/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,7 +130,7 @@ namespace AD{
 
   inline void StartPreacc() {
     if (globalTape.isActive() && PreaccEnabled) {
-      StartPosition = globalTape.getPosition();
+      PreaccHelper.start();
       PreaccActive = true;
     }
   }
@@ -132,7 +138,7 @@ namespace AD{
   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
 }
 
diff --git a/Common/include/datatypes/codi_reverse_structure.hpp b/Common/include/datatypes/codi_reverse_structure.hpp
index eab5214738e1..8221b805cfff 100644
--- a/Common/include/datatypes/codi_reverse_structure.hpp
+++ b/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
diff --git a/Common/include/grid_movement_structure.hpp b/Common/include/grid_movement_structure.hpp
index b5d7176f3eaa..7d4fa1ffe508 100644
--- a/Common/include/grid_movement_structure.hpp
+++ b/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;
diff --git a/Common/include/linear_solvers_structure.hpp b/Common/include/linear_solvers_structure.hpp
index fdfe4ceaf221..8f6336ca979e 100644
--- a/Common/include/linear_solvers_structure.hpp
+++ b/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
diff --git a/Common/include/linear_solvers_structure_b.hpp b/Common/include/linear_solvers_structure_b.hpp
index 6bcfd45264b2..877b1f758645 100644
--- a/Common/include/linear_solvers_structure_b.hpp
+++ b/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
diff --git a/Common/include/mpi_structure.hpp b/Common/include/mpi_structure.hpp
index 1528cf9e3c7e..af8a68cd5978 100644
--- a/Common/include/mpi_structure.hpp
+++ b/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)
diff --git a/Common/include/mpi_structure.inl b/Common/include/mpi_structure.inl
index 1fb314bd07ee..b1a65e3fe788 100644
--- a/Common/include/mpi_structure.inl
+++ b/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
diff --git a/Common/src/ad_structure.cpp b/Common/src/ad_structure.cpp
index 8f84885c6717..8ed974bd02ef 100644
--- a/Common/src/ad_structure.cpp
+++ b/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
 }
diff --git a/Common/src/linear_solvers_structure.cpp b/Common/src/linear_solvers_structure.cpp
index a57cddb162f0..6d0d0eee3c9d 100644
--- a/Common/src/linear_solvers_structure.cpp
+++ b/Common/src/linear_solvers_structure.cpp
@@ -643,11 +643,12 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS
 
   if (config->GetDiscrete_Adjoint()) {
 #ifdef CODI_REVERSE_TYPE
-
-   /*--- Check whether the tape is active, i.e. if it is recording and store the status ---*/
-
+    
     TapeActive = AD::globalTape.isActive();
 
+    AD::StartExtFunc(false, false);
+
+    AD::SetExtFuncIn(&LinSysRes[0], LinSysRes.GetLocSize());
 
     /*--- Stop the recording for the linear solver ---*/
 
@@ -660,17 +661,17 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS
   if (KindSolver == BCGSTAB || KindSolver == CONJUGATE_GRADIENT ||
       KindSolver == FGMRES  || KindSolver == RESTARTED_FGMRES ) {
     
-    mat_vec = new CSysMatrixVectorProduct(Jacobian, geometry, config);
+    mat_vec   = new CSysMatrixVectorProduct(Jacobian, geometry, config);    
     CPreconditioner* precond = NULL;
     
     switch (KindPrecond) {
       case JACOBI:
         Jacobian.BuildJacobiPreconditioner();
-        precond = new CJacobiPreconditioner(Jacobian, geometry, config);
+        precond   = new CJacobiPreconditioner(Jacobian, geometry, config);
         break;
       case ILU:
         Jacobian.BuildILUPreconditioner();
-        precond = new CILUPreconditioner(Jacobian, geometry, config);
+        precond   = new CILUPreconditioner(Jacobian, geometry, config);
         break;
       case LU_SGS:
         precond = new CLU_SGSPreconditioner(Jacobian, geometry, config);
@@ -745,80 +746,26 @@ unsigned long CSysSolve::Solve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CS
 
 
   if(TapeActive) {
+
+    bool RequiresTranspose = !mesh_deform; // jacobian is symmetric
+    if (!mesh_deform) KindPrecond = config->GetKind_DiscAdj_Linear_Prec();
+    else              KindPrecond = config->GetKind_Deform_Linear_Solver_Prec();
+	
     /*--- Start recording if it was stopped for the linear solver ---*/
 
     AD::StartRecording();
-
-    /*--- Prepare the externally differentiated linear solver ---*/
-
-    SetExternalSolve(Jacobian, LinSysRes, LinSysSol, geometry, config);
-
-  }
-
-  return IterLinSol;
-  
-}
-
-void CSysSolve::SetExternalSolve(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config) {
-
+    
+    AD::SetExtFuncOut(&LinSysSol[0], LinSysSol.GetLocSize());
+   
 #ifdef CODI_REVERSE_TYPE
-  
-  unsigned long size = LinSysRes.GetLocSize();
-  unsigned long i, nBlk = LinSysRes.GetNBlk(),
-                nVar = LinSysRes.GetNVar(),
-                nBlkDomain = LinSysRes.GetNBlkDomain();
-
-  bool RequiresTranspose = !mesh_deform; // jacobian is symmetric
-  
-  unsigned short KindPrecond;
-  
-  if (!mesh_deform) KindPrecond = config->GetKind_DiscAdj_Linear_Prec();
-  else              KindPrecond = config->GetKind_Deform_Linear_Solver_Prec();
-
-  /*--- Arrays to store the indices of the input/output of the linear solver.
-     * Note: They will be deleted in the CSysSolve_b::Delete_b routine. ---*/
-
-  su2double::GradientData *LinSysRes_Indices = new su2double::GradientData[size];
-  su2double::GradientData *LinSysSol_Indices = new su2double::GradientData[size];
-#if CODI_PRIMAL_INDEX_TAPE
-  su2double::Real *oldValues = new su2double::Real[size];
-#endif
-
-  for (i = 0; i < size; i++) {
-
-    /*--- Register the solution of the linear system (could already be registered when using multigrid) ---*/
-
-    if (!LinSysSol[i].isActive()) {
-#if CODI_PRIMAL_INDEX_TAPE
-      oldValues[i] = AD::globalTape.registerExtFunctionOutput(LinSysSol[i]);
-#else
-      AD::globalTape.registerInput(LinSysSol[i]);
-#endif
-    }
-
-    /*--- Store the indices ---*/
-
-    LinSysRes_Indices[i] = LinSysRes[i].getGradientData();
-    LinSysSol_Indices[i] = LinSysSol[i].getGradientData();
-  }
-
-  /*--- Push the data to the checkpoint handler for access in the reverse sweep ---*/
-
-  AD::CheckpointHandler* dataHandler = new AD::CheckpointHandler;
-
-  dataHandler->addData(LinSysRes_Indices);
-  dataHandler->addData(LinSysSol_Indices);
-#if CODI_PRIMAL_INDEX_TAPE
-  dataHandler->addData(oldValues);
-#endif
-  dataHandler->addData(size);
-  dataHandler->addData(nBlk);
-  dataHandler->addData(nVar);
-  dataHandler->addData(nBlkDomain);
-  dataHandler->addData(&Jacobian);
-  dataHandler->addData(geometry);
-  dataHandler->addData(config);
-  dataHandler->addData(this);
+    AD::FuncHelper->addUserData(&LinSysRes);
+    AD::FuncHelper->addUserData(&LinSysSol);
+    AD::FuncHelper->addUserData(&Jacobian);
+    AD::FuncHelper->addUserData(geometry);
+    AD::FuncHelper->addUserData(config);
+    AD::FuncHelper->addUserData(this);
+    AD::FuncHelper->addToTape(CSysSolve_b::Solve_b);
+#endif   
 
   /*--- Build preconditioner for the transposed Jacobian ---*/
 
@@ -833,14 +780,16 @@ void CSysSolve::SetExternalSolve(CSysMatrix & Jacobian, CSysVector & LinSysRes,
       SU2_MPI::Error("The specified preconditioner is not yet implemented for the discrete adjoint method.", CURRENT_FUNCTION);
       break;
   }
+    
+    AD::EndExtFunc();
 
-  /*--- Push the external function to the AD tape ---*/
-
-  AD::globalTape.pushExternalFunction(&CSysSolve_b::Solve_b, dataHandler, &CSysSolve_b::Delete_b);
+  }
 
-#endif
+  return IterLinSol;
+  
 }
 
+
 unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes, CSysVector & LinSysSol, CGeometry *geometry, CConfig *config) {
 #ifdef CODI_REVERSE_TYPE
 
@@ -923,3 +872,4 @@ unsigned long CSysSolve::Solve_b(CSysMatrix & Jacobian, CSysVector & LinSysRes,
   return 0;
 #endif
 }
+
diff --git a/Common/src/linear_solvers_structure_b.cpp b/Common/src/linear_solvers_structure_b.cpp
index 60848e99cd65..033dbc71f415 100644
--- a/Common/src/linear_solvers_structure_b.cpp
+++ b/Common/src/linear_solvers_structure_b.cpp
@@ -41,109 +41,40 @@
 #include "../include/matrix_structure.hpp"
 
 #ifdef CODI_REVERSE_TYPE
-void CSysSolve_b::Solve_b(AD::Tape* tape, AD::CheckpointHandler* data) {
+void CSysSolve_b::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){
+  
+  CSysVector* LinSysRes_b = NULL;
+  d->getData(LinSysRes_b);
+  
+  CSysVector* LinSysSol_b = NULL;
+  d->getData(LinSysSol_b);
   
-  /*--- Extract data from the checkpoint handler ---*/
-
-  su2double::GradientData *LinSysRes_Indices;
-  su2double::GradientData *LinSysSol_Indices;
-#if CODI_PRIMAL_INDEX_TAPE
-  su2double::Real *oldValues;
-#endif
-
-  data->getData(LinSysRes_Indices);
-  data->getData(LinSysSol_Indices);
-#if CODI_PRIMAL_INDEX_TAPE
-  data->getData(oldValues);
-#endif
-
-  unsigned long nBlk = 0, nVar = 0, nBlkDomain = 0, size = 0, i = 0;
-
-  data->getData(size);
-  data->getData(nBlk);
-  data->getData(nVar);
-  data->getData(nBlkDomain);
-
   CSysMatrix* Jacobian = NULL;
-  data->getData(Jacobian);
+  d->getData(Jacobian);
 
   CGeometry* geometry  = NULL;
-  data->getData(geometry);
+  d->getData(geometry);
 
   CConfig* config      = NULL;
-  data->getData(config);
-
+  d->getData(config);
+  
   CSysSolve* solver;
-  data->getData(solver);
+  d->getData(solver);
+  
 
   /*--- Initialize the right-hand side with the gradient of the solution of the primal linear system ---*/
 
-  CSysVector LinSysRes_b(nBlk, nBlkDomain, nVar, 0.0);
-  CSysVector LinSysSol_b(nBlk, nBlkDomain, nVar, 0.0);
-
-  for (i = 0; i < size; i ++) {
-    su2double::GradientData& index = LinSysSol_Indices[i];
-    LinSysRes_b[i] = AD::globalTape.getGradient(index);
-    LinSysSol_b[i] = 0.0;
-    AD::globalTape.gradient(index) = 0.0;
-  }
-
-  /*--- Solve the system ---*/
-
-  solver->Solve_b(*Jacobian, LinSysRes_b, LinSysSol_b, geometry, config);
-
-  /*--- Update the gradients of the right-hand side of the primal linear system ---*/
-
-  for (i = 0; i < size; i ++) {
-    su2double::GradientData& index = LinSysRes_Indices[i];
-    AD::globalTape.gradient(index) += SU2_TYPE::GetValue(LinSysSol_b[i]);
+  for (unsigned long i = 0; i < n; i ++) {
+    (*LinSysRes_b)[i] = y_b[i];
+    (*LinSysSol_b)[i] = 0.0;
   }
-
-#if CODI_PRIMAL_INDEX_TAPE
-  /*--- Set the old values that have been overwritten ---*/
-  for (i = 0; i < size; i ++) {
-    AD::globalTape.setExternalValueChange(LinSysSol_Indices[i], oldValues[i]);
+  
+  solver->Solve_b(*Jacobian, *LinSysRes_b, *LinSysSol_b, geometry, config);
+    
+  for (unsigned long i = 0; i < n; i ++) {
+    x_b[i] = SU2_TYPE::GetValue(LinSysSol_b->operator [](i));
   }
-#endif
+  
 }
 
-void CSysSolve_b::Delete_b(AD::Tape* tape, AD::CheckpointHandler* data) {
-
-  su2double::GradientData *LinSysRes_Indices = NULL;
-  su2double::GradientData *LinSysSol_Indices = NULL;
-#if CODI_PRIMAL_INDEX_TAPE
-  su2double::Real *oldValues;
-#endif
-
-  data->getData(LinSysRes_Indices);
-  data->getData(LinSysSol_Indices);
-#if CODI_PRIMAL_INDEX_TAPE
-  data->getData(oldValues);
-#endif
-
-  delete [] LinSysRes_Indices;
-  delete [] LinSysSol_Indices;
-#if CODI_PRIMAL_INDEX_TAPE
-  delete [] oldValues;
-#endif
-
-  unsigned long nBlk, nVar, nBlkDomain, size;
-
-  data->getData(size);
-  data->getData(nBlk);
-  data->getData(nVar);
-  data->getData(nBlkDomain);
-
-  CSysMatrix* Jacobian;
-  data->getData(Jacobian);
-
-  CGeometry* geometry;
-  data->getData(geometry);
-
-  CConfig* config;
-  data->getData(config);
-
-  CSysSolve* solver;
-  data->getData(solver);
-}
 #endif
diff --git a/externals/codi b/externals/codi
index bd4a639c2fe6..501dcf0305df 160000
--- a/externals/codi
+++ b/externals/codi
@@ -1 +1 @@
-Subproject commit bd4a639c2fe625a80946c8365bd2976a2868cf46
+Subproject commit 501dcf0305df147481630f20ce37c2e624fb351f
diff --git a/externals/medi b/externals/medi
index 46a97e1d6e8f..a95a23ce7585 160000
--- a/externals/medi
+++ b/externals/medi
@@ -1 +1 @@
-Subproject commit 46a97e1d6e8fdd3cb42b06534cff6acad2a49693
+Subproject commit a95a23ce7585905c3a731b28c1bb512028fc02bb