diff --git a/Common/include/linear_algebra/CPreconditioner.hpp b/Common/include/linear_algebra/CPreconditioner.hpp
index e4fc7cf159f..b73a1f8967c 100644
--- a/Common/include/linear_algebra/CPreconditioner.hpp
+++ b/Common/include/linear_algebra/CPreconditioner.hpp
@@ -114,13 +114,39 @@ class CJacobiPreconditioner final : public CPreconditioner<ScalarType> {
    * \param[out] v - CSysVector that is the result of the preconditioning
    */
   inline void operator()(const CSysVector<ScalarType>& u, CSysVector<ScalarType>& v) const override {
-    sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config);
+    if (config->GetCUDA()) {
+#ifdef HAVE_CUDA
+      sparse_matrix.GPUComputeJacobiPreconditioner(u, v, geometry, config);
+#else
+      SU2_MPI::Error(
+          "\nError in launching sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with "
+          "CUDA "
+          "options enabled in Meson to access GPU Functions",
+          CURRENT_FUNCTION);
+#endif
+    } else {
+      sparse_matrix.ComputeJacobiPreconditioner(u, v, geometry, config);
+    }
   }
 
   /*!
    * \note Request the associated matrix to build the preconditioner.
    */
-  inline void Build() override { sparse_matrix.BuildJacobiPreconditioner(); }
+  inline void Build() override {
+    if (config->GetCUDA()) {
+#ifdef HAVE_CUDA
+      sparse_matrix.GPUBuildJacobiPreconditioner();
+#else
+      SU2_MPI::Error(
+          "\nError in building sparse matrix Preconditioner Function\nENABLE_CUDA is set to YES\nPlease compile with "
+          "CUDA "
+          "options enabled in Meson to access GPU Functions",
+          CURRENT_FUNCTION);
+#endif
+    } else {
+      sparse_matrix.BuildJacobiPreconditioner();
+    }
+  }
 };
 
 /*!
diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp
index 11baf6edf8e..d11ccbdf370 100644
--- a/Common/include/linear_algebra/CSysMatrix.hpp
+++ b/Common/include/linear_algebra/CSysMatrix.hpp
@@ -176,7 +176,9 @@ class CSysMatrix {
   /*!< \brief Level structure for alternative shared memory parallelization of ILU. */
   CCompressedSparsePatternUL levels_ilu;
 
-  ScalarType* invM; /*!< \brief Inverse of (Jacobi) preconditioner. */
+  ScalarType* invM;             /*!< \brief Inverse of (Jacobi) preconditioner. */
+  ScalarType* d_invM;           /*!< \brief Inverse of (Jacobi) preconditioner on device. */
+  bool invM_is_managed = false; /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */
 
   /*--- Temporary (hence mutable) working memory used in the Linelet preconditioner, outer vector is for threads ---*/
   mutable vector<vector<const ScalarType*> >
@@ -924,6 +926,11 @@ class CSysMatrix {
    */
   void BuildJacobiPreconditioner();
 
+  /*!
+   * \brief Build the Jacobi preconditioner on GPU.
+   */
+  void GPUBuildJacobiPreconditioner();
+
   /*!
    * \brief Multiply CSysVector by the preconditioner
    * \param[in] vec - CSysVector to be multiplied by the preconditioner.
@@ -934,6 +941,16 @@ class CSysMatrix {
   void ComputeJacobiPreconditioner(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod, CGeometry* geometry,
                                    const CConfig* config) const;
 
+  /*!
+   * \brief Multiply CSysVector by the preconditioner
+   * \param[in] vec - CSysVector to be multiplied by the preconditioner.
+   * \param[out] prod - Result of the product A*vec.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   */
+  void GPUComputeJacobiPreconditioner(const CSysVector<ScalarType>& vec, CSysVector<ScalarType>& prod,
+                                      CGeometry* geometry, const CConfig* config) const;
+
   /*!
    * \brief Build the ILU preconditioner.
    */
diff --git a/Common/include/linear_algebra/CSysVector.hpp b/Common/include/linear_algebra/CSysVector.hpp
index 1498b549bbb..94aac89e113 100644
--- a/Common/include/linear_algebra/CSysVector.hpp
+++ b/Common/include/linear_algebra/CSysVector.hpp
@@ -34,6 +34,7 @@
 #include "../parallelization/mpi_structure.hpp"
 #include "../parallelization/omp_structure.hpp"
 #include "../parallelization/vectorization.hpp"
+#include "gpu_ast.hpp"
 #include "vector_expressions.hpp"
 #include "../../include/CConfig.hpp"
 
@@ -77,6 +78,14 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   unsigned long nVar = 1;       /*!< \brief Number of elements in a block. */
 
   ScalarType* d_vec_val = nullptr; /*!< \brief Device Pointer to store the vector values on the GPU. */
+  bool vec_is_managed = false;     /*!< \brief Boolean that indicates whether GPU supports Unified Memory or not */
+  bool useCuda = false;            /*!< \brief Whether CUDA is enabled. */
+  enum class GPUMemState {
+    HOST,
+    DEVICE,
+    SYNCED
+  }; /*!< \brief enumerate the possible memory synchronization states for this */
+  mutable GPUMemState memState = GPUMemState::HOST; /*!< \brief mutable GPUMemState initialized to HOST */
 
 #ifdef HAVE_OMP
   mutable std::unique_ptr<ScalarType[]>
@@ -164,7 +173,12 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \note Not defined for expressions because we do not know their sizes.
    * \param[in] u - Vector being copied.
    */
-  CSysVector(const CSysVector& u) { Initialize(u.GetNBlk(), u.GetNBlkDomain(), u.nVar, u.vec_val, true); }
+  CSysVector(const CSysVector& u) {
+#ifdef HAVE_CUDA
+    u.SyncToHost();
+#endif
+    Initialize(u.GetNBlk(), u.GetNBlkDomain(), u.nVar, u.vec_val, true);
+  }
 
   /*!
    * \brief Swap contents with another vector.
@@ -177,6 +191,7 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     std::swap(nElmDomain, other.nElmDomain);
     std::swap(nVar, other.nVar);
     std::swap(dot_scratch, other.dot_scratch);
+    std::swap(memState, other.memState);
   }
 
   /*!
@@ -210,7 +225,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   void PassiveCopy(const CSysVector<T>& other) {
     /*--- This is a method and not the overload of an operator to make sure who
      * calls it knows the consequence to the derivative information (lost) ---*/
-
+#ifdef HAVE_CUDA
+    other.SyncToHost();
+#endif
     /*--- check if self-assignment, otherwise perform deep copy ---*/
     if ((const void*)this == (const void*)&other) return;
 
@@ -220,8 +237,58 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     CSYSVEC_PARFOR
     for (auto i = 0ul; i < nElm; i++) vec_val[i] = SU2_TYPE::GetValue(other[i]);
     END_CSYSVEC_PARFOR
+#ifdef HAVE_CUDA
+    MarkHostDirty();
+#endif
   }
 
+  /*!
+   * \brief enumeratethe possible vector-scalar operators on GPU
+   */
+  enum class GPUScalarOp { SET, ADD, SUB, MUL, DIV };
+
+  /*!
+   * \brief enumerate the possible vector-vector operators on GPU
+   */
+  enum class GPUVectorOp { SET, ADD, SUB, NEG };
+
+  /*!
+   * \brief map GPU scalar Ops to the appropriate GPU Operation code in the tree
+   * \note prevents SET=assignement operations raising an error
+   */
+  inline GPUOpType ToASTCombineOp(GPUScalarOp op) {
+    switch (op) {
+      case GPUScalarOp::ADD:
+        return GPUOpType::ADD;
+      case GPUScalarOp::SUB:
+        return GPUOpType::SUB;
+      case GPUScalarOp::MUL:
+        return GPUOpType::MUL;
+      case GPUScalarOp::DIV:
+        return GPUOpType::DIV;
+      default:
+        SU2_MPI::Error("ToASTCombineOp called with GPUScalarOp::SET.", CURRENT_FUNCTION);
+        return GPUOpType::ADD;
+    }
+  }
+
+  /*!
+   * \brief appends a new vector node to the Abstract Syntax Tree for the given payload
+   */
+  int BuildAST(GPUASTPayload<ScalarType>& payload) const {
+    int idx = payload.NewNode();
+    payload.nodes[idx].op = GPUOpType::VEC;
+    payload.nodes[idx].d_ptr = this->d_vec_val;
+    return idx;
+  }
+
+  /*!
+   * \brief method to launch the generic Unary Operation Kernel and apply given functors on GPU
+   * \param[in] op - Unitary operation
+   * \param[in] val - scalar value
+   */
+  void GPUUnaryOperation(GPUScalarOp op, ScalarType val);
+
   /*!
    * \brief Performs the memory copy from host to device.
    * \param[in] trigger - boolean value that decides whether to conduct the transfer or not. True by default.
@@ -245,6 +312,16 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    */
   inline ScalarType* GetDevicePointer() const { return d_vec_val; }
 
+  /*!
+   * \brief return pointer that points to the CSysVector values in CPU memory
+   */
+  ScalarType* data() const {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
+    return vec_val;
+  }
+
   /*!
    * \brief return the number of local elements in the CSysVector
    */
@@ -275,14 +352,35 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] i - Local index to access.
    * \return Value at position i.
    */
-  inline ScalarType& operator[](unsigned long i) { return vec_val[i]; }
-  inline const ScalarType& operator[](unsigned long i) const { return vec_val[i]; }
+  inline ScalarType& operator[](unsigned long i) {
+#ifdef HAVE_CUDA
+    SyncToHost();
+    MarkHostDirty();
+#endif
+    return vec_val[i];
+  }
+  inline const ScalarType& operator[](unsigned long i) const {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
+    return vec_val[i];
+  }
 
   /*!
    * \brief Iterators for range for loops.
    */
-  inline const ScalarType* begin() const { return vec_val; }
-  inline const ScalarType* end() const { return vec_val + nElm; }
+  inline const ScalarType* begin() const {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
+    return vec_val;
+  }
+  inline const ScalarType* end() const {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
+    return vec_val + nElm;
+  }
 
   /*!
    * \brief Access operator with assignment permitted block version.
@@ -290,8 +388,17 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] iVar - Index of variable.
    * \return Value at position (i,j).
    */
-  inline ScalarType& operator()(unsigned long iPoint, unsigned long iVar) { return vec_val[iPoint * nVar + iVar]; }
+  inline ScalarType& operator()(unsigned long iPoint, unsigned long iVar) {
+#ifdef HAVE_CUDA
+    SyncToHost();
+    MarkHostDirty();
+#endif
+    return vec_val[iPoint * nVar + iVar];
+  }
   inline const ScalarType& operator()(unsigned long iPoint, unsigned long iVar) const {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
     return vec_val[iPoint * nVar + iVar];
   }
 
@@ -301,17 +408,46 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] other - Another vector.
    */
   CSysVector& operator=(const CSysVector& other) {
+#ifdef HAVE_CUDA
+    other.SyncToHost();
+#endif
     CSYSVEC_PARFOR
     for (auto i = 0ul; i < nElm; ++i) vec_val[i] = other.vec_val[i];
     END_CSYSVEC_PARFOR
+#ifdef HAVE_CUDA
+    this->MarkHostDirty();
+#endif
     return *this;
   }
 
   /*!
-   * \brief Compound assignement operations with scalars and expressions.
+   * \brief Compound assignement operations with scalars and expressions, GPU or CPU
    * \param[in] val/expr - Scalar value or expression.
    */
-#define MAKE_COMPOUND(OP)                                                 \
+#ifdef HAVE_CUDA
+#define MAKE_COMPOUND(OP, TAG)                                            \
+  CSysVector& operator OP(ScalarType val) {                               \
+    GPUUnaryOperation(GPUScalarOp::TAG, val);                             \
+    return *this;                                                         \
+  }                                                                       \
+  template <class T>                                                      \
+  CSysVector& operator OP(const VecExpr::CVecExpr<T, ScalarType>& expr) { \
+    GPUASTPayload<ScalarType> payload;                                    \
+    int rhs_idx = expr.derived().BuildAST(payload);                       \
+    int root_idx = rhs_idx;                                               \
+    if (GPUScalarOp::TAG != GPUScalarOp::SET) {                           \
+      int self_idx = this->BuildAST(payload);                             \
+      root_idx = payload.NewNode();                                       \
+      payload.nodes[root_idx].op = ToASTCombineOp(GPUScalarOp::TAG);      \
+      payload.nodes[root_idx].left = self_idx;                            \
+      payload.nodes[root_idx].right = rhs_idx;                            \
+    }                                                                     \
+    payload.root_idx = root_idx;                                          \
+    LaunchGPUAST(payload, this->nElm);                                    \
+    return *this;                                                         \
+  }
+#else
+#define MAKE_COMPOUND(OP, TAG)                                            \
   CSysVector& operator OP(ScalarType val) {                               \
     CSYSVEC_PARFOR                                                        \
     for (auto i = 0ul; i < nElm; ++i) vec_val[i] OP val;                  \
@@ -325,11 +461,13 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     END_CSYSVEC_PARFOR                                                    \
     return *this;                                                         \
   }
-  MAKE_COMPOUND(=)
-  MAKE_COMPOUND(+=)
-  MAKE_COMPOUND(-=)
-  MAKE_COMPOUND(*=)
-  MAKE_COMPOUND(/=)
+#endif
+
+  MAKE_COMPOUND(=, SET)
+  MAKE_COMPOUND(+=, ADD)
+  MAKE_COMPOUND(-=, SUB)
+  MAKE_COMPOUND(*=, MUL)
+  MAKE_COMPOUND(/=, DIV)
 #undef MAKE_COMPOUND
 
   /*!
@@ -346,9 +484,12 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   ScalarType dot(const VecExpr::CVecExpr<T, ScalarType>& expr) const {
     /*--- All threads get the same "view" of the vectors. ---*/
     SU2_OMP_BARRIER
+    ScalarType sum = 0.0;
 
+#ifdef HAVE_CUDA
+    sum = GPUDot(static_cast<const CSysVector&>(expr.derived()));  // assuming expr is a vector
+#else
     /*--- Local dot product for each thread. ---*/
-    ScalarType sum = 0.0;
 
     CSYSVEC_PARFOR
     for (auto i = 0ul; i < nElmDomain; ++i) {
@@ -357,11 +498,13 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     END_CSYSVEC_PARFOR
 
     dot_scratch[omp_get_thread_num()] = sum;
+#endif
 
     BEGIN_SU2_OMP_SAFE_GLOBAL_ACCESS {
+#ifndef HAVE_CUDA
       /*--- Reduce over all threads in an ordered way to ensure a deterministic result. ---*/
       for (int i = 1; i < omp_get_num_threads(); ++i) sum += dot_scratch[i];
-
+#endif
       /*--- Reduce across all mpi ranks, only the master thread communicates. ---*/
       const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE;
       SelectMPIWrapper<ScalarType>::W::Allreduce(&sum, &dot_scratch[0], 1, mpi_type, MPI_SUM, SU2_MPI::GetComm());
@@ -372,6 +515,20 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     return dot_scratch[0];
   }
 
+  /*!
+   * \brief Dot product between "this" and another vector.
+   * \param[in] other - the other CSysVector.
+   * \return Result of dot product
+   */
+  ScalarType GPUDot(const CSysVector& other) const;
+
+  /*!
+   * \brief launch method for the generic Abstract Syntax Tree GPU kernel
+   * \param[in] payload - payload is the abstract syntax tree
+   * \param[in] nElm - number of elements in the vector
+   */
+  void LaunchGPUAST(const GPUASTPayload<ScalarType>& payload, unsigned long nElm);
+
   /*!
    * \brief Computes the product of V^T W efficiencly, where V and W are tall matrices stored as vectors of CSysVector.
    * \param[in] V - Tall matrix.
@@ -384,6 +541,61 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
   static const su2matrix<ScalarType>& multiDot(const std::vector<CSysVector>& V, size_t i0, size_t n,
                                                const std::vector<CSysVector>& W, size_t m);
 
+  /*!
+   * \brief Computes the product of V^T W on the GPU, where V and W are tall matrices stored as vectors of CSysVector.
+   * \param[in] V - Tall matrix.
+   * \param[in] i0 - First column of V to consider.
+   * \param[in] n - Number of columns to consider from V starting at i0.
+   * \param[in] W - Tall matrix.
+   * \param[in] m - Number of columns to consider from W.
+   * \return n by m matrix with the result of the product.
+   */
+  static const su2matrix<ScalarType>& multiDotGPU(const std::vector<CSysVector<ScalarType>>& V, const size_t i0,
+                                                  const size_t n, const std::vector<CSysVector<ScalarType>>& W,
+                                                  const size_t m);
+
+  /*!
+   * \brief Computes v = vs * ws or v += vs * ws with unrolling of up to 4 iterations on the GPU
+   * \param[in] n - number of vectors to consider
+   * \param[in] ws - array of scalar weights corresponding to the device pointers to vectors
+   * \param[in] vs_ptrs - array of device pointers
+   * \param[in] v - target vector
+   * \param[in] inc - If true, adds results to target vector. If false, overwrites
+   */
+  static void LinearCombinationGPU(const unsigned long n, const std::vector<CSysVector<ScalarType>>& vs,
+                                   const ScalarType* ws, CSysVector<ScalarType>& v, bool inc = false);
+
+  /*!
+   * \brief method to mark the vector as lastly modified by device
+   */
+  void MarkDeviceDirty() const { memState = GPUMemState::DEVICE; }
+
+  /*!
+   * \brief method to mark the vector as lastly modified by host
+   */
+  void MarkHostDirty() const { memState = GPUMemState::HOST; }
+
+  /*!
+   * \brief wrapper around HtDTransfer to manage memory state
+   */
+  void SyncToDevice() const {
+    if (memState == GPUMemState::HOST) {
+      HtDTransfer();
+      memState = GPUMemState::SYNCED;
+    }
+  }
+  /*!
+   * \brief wrapper around DtHTransfer to manage memory state
+   */
+  void SyncToHost() const {
+    if (memState == GPUMemState::DEVICE) {
+      DtHTransfer();
+      memState = GPUMemState::SYNCED;
+    }
+  }
+
+  void print_memory_state() const { cout << "memory state: " << static_cast<int>(memState) << endl; }
+
   /*!
    * \brief Squared L2 norm of the vector (via dot with self).
    * \return Squared L2 norm.
@@ -401,15 +613,32 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    * \param[in] iPoint - Index of block.
    * \return Pointer to start of block.
    */
-  inline ScalarType* GetBlock(unsigned long iPoint) { return &vec_val[iPoint * nVar]; }
-  inline const ScalarType* GetBlock(unsigned long iPoint) const { return &vec_val[iPoint * nVar]; }
+  inline ScalarType* GetBlock(unsigned long iPoint) {
+#ifdef HAVE_CUDA
+    SyncToHost();
+    MarkHostDirty();
+#endif
+    return &vec_val[iPoint * nVar];
+  }
+  inline const ScalarType* GetBlock(unsigned long iPoint) const {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
+    return &vec_val[iPoint * nVar];
+  }
 
   /*!
    * \brief Set the values to zero for one block.
    * \param[in] iPoint - Index of the block being set to zero.
    */
   inline void SetBlock_Zero(unsigned long iPoint) {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
     for (auto iVar = 0ul; iVar < nVar; iVar++) vec_val[iPoint * nVar + iVar] = 0.0;
+#ifdef HAVE_CUDA
+    MarkHostDirty();
+#endif
   }
 
   /*!
@@ -421,11 +650,17 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
    */
   template <class VectorType, bool Overwrite = true>
   FORCEINLINE void SetBlock(unsigned long iPoint, const VectorType& block, ScalarType alpha = 1) {
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
     if (Overwrite) {
       for (auto i = 0ul; i < nVar; ++i) vec_val[iPoint * nVar + i] = alpha * block[i];
     } else {
       for (auto i = 0ul; i < nVar; ++i) vec_val[iPoint * nVar + i] += alpha * block[i];
     }
+#ifdef HAVE_CUDA
+    MarkHostDirty();
+#endif
   }
 
   /*!
@@ -468,13 +703,18 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     assert(nVar == this->nVar);
     ScalarType vec[N][nVar];
     UnpackBlock(vector, mask, vec);
-
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
     /*--- Update one by one skipping if mask is 0. ---*/
     for (size_t k = 0; k < N; ++k) {
       if (mask[k] == 0) continue;
       SU2_OMP_SIMD
       for (size_t i = 0; i < nVar; ++i) vec_val[iPoint[k] * nVar + i] = vec[k][i];
     }
+#ifdef HAVE_CUDA
+    MarkHostDirty();
+#endif
   }
 
   /*!
@@ -489,7 +729,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
     assert(nVar == this->nVar);
     ScalarType vec[N][nVar];
     UnpackBlock(vector, mask, vec);
-
+#ifdef HAVE_CUDA
+    SyncToHost();
+#endif
     /*--- Update one by one skipping if mask is 0. ---*/
     for (size_t k = 0; k < N; ++k) {
       if (mask[k] == 0) continue;
@@ -499,6 +741,9 @@ class CSysVector : public VecExpr::CVecExpr<CSysVector<ScalarType>, ScalarType>
         vec_val[jPoint[k] * nVar + i] -= vec[k][i];
       }
     }
+#ifdef HAVE_CUDA
+    MarkHostDirty();
+#endif
   }
 };
 
diff --git a/Common/include/linear_algebra/GPUComms.cuh b/Common/include/linear_algebra/GPUComms.cuh
index 13854381877..1f8f1202135 100644
--- a/Common/include/linear_algebra/GPUComms.cuh
+++ b/Common/include/linear_algebra/GPUComms.cuh
@@ -51,3 +51,16 @@ inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=t
 }
 
 #define gpuErrChk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
+
+/*!
+ * \brief Prefetch a CUDA Unified Memory array to a device asynchronously
+ * \param[in] ptr, pointer to the memory we want to prefetch.
+ * \param[in] size in bytes.
+ * \tparam ZeroInit, initialize memory to 0.
+ */
+template <class T>
+inline void gpu_um_prefetch(T* ptr, size_t size, int device) noexcept {
+#ifdef HAVE_CUDA
+  gpuErrChk(cudaMemPrefetchAsync((void*)ptr, size, device));
+#endif
+}
diff --git a/Common/include/linear_algebra/gpu_ast.hpp b/Common/include/linear_algebra/gpu_ast.hpp
new file mode 100644
index 00000000000..9047d5dfae5
--- /dev/null
+++ b/Common/include/linear_algebra/gpu_ast.hpp
@@ -0,0 +1,70 @@
+/*!
+ * \file gpu_ast.hpp
+ * \brief definition of the Abstract Syntax Tree for evaluating CVecExpr expression trees on GPU with a single generic
+ * interpreter kernel \author D. Di giusto \version 8.5.0 "Harrier"
+ *
+ * SU2 Project Website: https://su2code.github.io
+ *
+ * The SU2 Project is maintained by the SU2 Foundation
+ * (http://su2foundation.org)
+ *
+ * Copyright 2012-2024, SU2 Contributors (cf. AUTHORS.md)
+ *
+ * SU2 is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ * SU2 is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with SU2. If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#pragma once
+
+#include <iostream>
+#include "../parallelization/mpi_structure.hpp"
+
+/*!
+ * \brief enumeration of the possible node-operations on GPU to build abstract syntax trees
+ */
+enum class GPUOpType { VEC, SCALAR, ADD, SUB, MUL, DIV, NEG };
+
+/*!
+ * \brief GPU Abstract Syntax Tree node structure
+ */
+template <typename ScalarType>
+struct GPUASTNode {
+  GPUOpType op = GPUOpType::SCALAR;   // operation for this node
+  int left = -1;                      // index of left child node-operation, -1 for leaf
+  int right = -1;                     // index of right child node-operation, -1 for leaf
+  const ScalarType* d_ptr = nullptr;  // device pointer for vector
+  ScalarType val = ScalarType(0);     // constant value for scalar
+};
+
+/*!
+ * \brief arbitrary max limit on the number of nodes in one Abstract Syntax tree
+ */
+constexpr int MAX_AST_NODES = 16;
+
+/*!
+ * \brief struct for the abstract syntax tree payload sent from CPU to GPU during one kernel launch
+ */
+template <class ScalarType>
+struct GPUASTPayload {
+  GPUASTNode<ScalarType> nodes[MAX_AST_NODES];  // array hosting all nodes for a single expression-operator
+  int root_idx = 0;
+  int node_count = 0;
+
+  // allocation method
+  int NewNode() {
+    if (node_count >= MAX_AST_NODES) {
+      SU2_MPI::Error("GPU AST node budget exceeded, you should increase MAX_AST_NODES!", __FUNCTION__);
+    }
+    return node_count++;
+  }
+};
diff --git a/Common/include/linear_algebra/vector_expressions.hpp b/Common/include/linear_algebra/vector_expressions.hpp
index a0d0ce28901..e35379a3609 100644
--- a/Common/include/linear_algebra/vector_expressions.hpp
+++ b/Common/include/linear_algebra/vector_expressions.hpp
@@ -35,6 +35,10 @@
 #include <cmath>
 #include <cstdint>
 
+#ifdef HAVE_CUDA
+#include "gpu_ast.hpp"
+#endif
+
 namespace VecExpr {
 /// \addtogroup VecExpr
 /// @{
@@ -78,6 +82,18 @@ class Bcast : public CVecExpr<Bcast<Scalar>, Scalar> {
   static constexpr bool StoreAsRef = false;
   FORCEINLINE Bcast(const Scalar& x_) : x(x_) {}
   FORCEINLINE const Scalar& operator[](size_t) const { return x; }
+
+/*!
+ * \brief appends a new scalar node to the abstract expression tree of a given payload
+ */
+#ifdef HAVE_CUDA
+  int BuildAST(GPUASTPayload<Scalar>& payload) const {
+    int idx = payload.NewNode();
+    payload.nodes[idx].op = GPUOpType::SCALAR;
+    payload.nodes[idx].val = x;
+    return idx;
+  }
+#endif
 };
 
 /*!
@@ -117,38 +133,78 @@ namespace math = ::std;
 #define RETURNS(...) \
   ->decltype(__VA_ARGS__) { return __VA_ARGS__; }
 
+/*--- Macro to inject GPU abstract syntax trees into unary functions
+ * The corresponding NONE macro is for operators not yet existing ---*/
+
+#ifdef HAVE_CUDA
+#define GPU_UNARY_AST(OPTYPE)                          \
+  int BuildAST(GPUASTPayload<Scalar>& payload) const { \
+    int child = u.BuildAST(payload);                   \
+    int idx = payload.NewNode();                       \
+    payload.nodes[idx].op = GPUOpType::OPTYPE;         \
+    payload.nodes[idx].left = child;                   \
+    return idx;                                        \
+  }
+#define GPU_UNARY_AST_NONE()
+#else
+#define GPU_UNARY_AST(OPTYPE)
+#define GPU_UNARY_AST_NONE()
+#endif
+
 /*--- Macro to create expression classes (EXPR) and overloads (FUN) for unary
  * functions, based on their coefficient-wise implementation (IMPL). ---*/
 
-#define MAKE_UNARY_FUN(FUN, EXPR, IMPL)                             \
-  /*!--- Expression class. ---*/                                    \
-  template <class U, class Scalar>                                  \
-  class EXPR : public CVecExpr<EXPR<U, Scalar>, Scalar> {           \
-    store_t<const U> u;                                             \
-                                                                    \
-   public:                                                          \
-    static constexpr bool StoreAsRef = false;                       \
-    FORCEINLINE EXPR(const U& u_) : u(u_) {}                        \
-    FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) \
-  };                                                                \
-  /*!--- Function overload, returns an expression object. ---*/     \
-  template <class U, class S>                                       \
+#define MAKE_UNARY_FUN(FUN, EXPR, IMPL, GPU_AST)                            \
+  /*!--- Expression class. ---*/                                            \
+  template <class U, class Scalar>                                          \
+  class EXPR : public CVecExpr<EXPR<U, Scalar>, Scalar> {                   \
+    store_t<const U> u;                                                     \
+                                                                            \
+   public:                                                                  \
+    static constexpr bool StoreAsRef = false;                               \
+    FORCEINLINE EXPR(const U& u_) : u(u_) {}                                \
+    FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i])) GPU_AST \
+  };                                                                        \
+  /*!--- Function overload, returns an expression object. ---*/             \
+  template <class U, class S>                                               \
   FORCEINLINE auto FUN(const CVecExpr<U, S>& u) RETURNS(EXPR<U, S>(u.derived()))
 
 #define sign_impl(x) Scalar(1 - 2 * (x < 0))
-MAKE_UNARY_FUN(operator-, minus_, -)
-MAKE_UNARY_FUN(abs, abs_, math::abs)
-MAKE_UNARY_FUN(exp, exp_, math::exp)
-MAKE_UNARY_FUN(sqrt, sqrt_, math::sqrt)
-MAKE_UNARY_FUN(sign, sign_, sign_impl)
+MAKE_UNARY_FUN(operator-, minus_, -, GPU_UNARY_AST(NEG))
+MAKE_UNARY_FUN(abs, abs_, math::abs, GPU_UNARY_AST_NONE())
+MAKE_UNARY_FUN(exp, exp_, math::exp, GPU_UNARY_AST_NONE())
+MAKE_UNARY_FUN(sqrt, sqrt_, math::sqrt, GPU_UNARY_AST_NONE())
+MAKE_UNARY_FUN(sign, sign_, sign_impl, GPU_UNARY_AST_NONE())
 #undef sign_impl
 
 #undef MAKE_UNARY_FUN
+#undef GPU_UNARY_AST
+#undef GPU_UNARY_AST_NONE
+
+/*--- Macro to inject GPU abstract syntax trees into binary functions
+ * The corresponding NONE macro is for operators not yet existing ---*/
+
+#ifdef HAVE_CUDA
+#define GPU_BINARY_AST(OPTYPE)                         \
+  int BuildAST(GPUASTPayload<Scalar>& payload) const { \
+    int l = u.BuildAST(payload);                       \
+    int r = v.BuildAST(payload);                       \
+    int idx = payload.NewNode();                       \
+    payload.nodes[idx].op = GPUOpType::OPTYPE;         \
+    payload.nodes[idx].left = l;                       \
+    payload.nodes[idx].right = r;                      \
+    return idx;                                        \
+  }
+#define GPU_BINARY_AST_NONE()
+#else
+#define GPU_BINARY_AST(OPTYPE)
+#define GPU_BINARY_AST_NONE()
+#endif
 
 /*--- Macro to create expressions and overloads for binary functions. ---*/
 
 // clang-format off
-#define MAKE_BINARY_FUN(FUN, EXPR, IMPL)                                                                              \
+#define MAKE_BINARY_FUN(FUN, EXPR, IMPL, GPU_AST)                                                                     \
   /*!--- Expression class. ---*/                                                                                      \
   template <class U, class V, class Scalar>                                                                           \
   class EXPR : public CVecExpr<EXPR<U, V, Scalar>, Scalar> {                                                          \
@@ -159,6 +215,7 @@ MAKE_UNARY_FUN(sign, sign_, sign_impl)
     static constexpr bool StoreAsRef = false;                                                                         \
     FORCEINLINE EXPR(const U& u_, const V& v_) : u(u_), v(v_) {}                                                      \
     FORCEINLINE auto operator[](size_t i) const RETURNS(IMPL(u[i], v[i]))                                             \
+    GPU_AST                                                                                                           \
   };                                                                                                                  \
   /*!--- Vector with vector function overload. ---*/                                                                  \
   template <class U, class V, class S>                                                                                \
@@ -191,9 +248,9 @@ using std::fmax;
 using std::fmin;
 #undef MAKE_FMINMAX_OVERLOADS
 
-MAKE_BINARY_FUN(fmax, max_, fmax)
-MAKE_BINARY_FUN(fmin, min_, fmin)
-MAKE_BINARY_FUN(pow, pow_, math::pow)
+MAKE_BINARY_FUN(fmax, max_, fmax, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(fmin, min_, fmin, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(pow, pow_, math::pow, GPU_BINARY_AST_NONE())
 
 /*--- sts::plus and co. were tried, the code was horrendous (due to the forced
  * conversion between different types) and creating functions for these ops
@@ -203,10 +260,10 @@ MAKE_BINARY_FUN(pow, pow_, math::pow)
 #define sub_impl(a, b) a - b
 #define mul_impl(a, b) a* b
 #define div_impl(a, b) a / b
-MAKE_BINARY_FUN(operator+, add_, add_impl)
-MAKE_BINARY_FUN(operator-, sub_, sub_impl)
-MAKE_BINARY_FUN(operator*, mul_, mul_impl)
-MAKE_BINARY_FUN(operator/, div_, div_impl)
+MAKE_BINARY_FUN(operator+, add_, add_impl, GPU_BINARY_AST(ADD))
+MAKE_BINARY_FUN(operator-, sub_, sub_impl, GPU_BINARY_AST(SUB))
+MAKE_BINARY_FUN(operator*, mul_, mul_impl, GPU_BINARY_AST(MUL))
+MAKE_BINARY_FUN(operator/, div_, div_impl, GPU_BINARY_AST(DIV))
 #undef add_impl
 #undef sub_impl
 #undef mul_impl
@@ -224,12 +281,12 @@ MAKE_BINARY_FUN(operator/, div_, div_impl)
 #define ne_impl(a, b) TO_PASSIVE(a != b)
 #define lt_impl(a, b) TO_PASSIVE(a < b)
 #define gt_impl(a, b) TO_PASSIVE(a > b)
-MAKE_BINARY_FUN(operator<=, le_, le_impl)
-MAKE_BINARY_FUN(operator>=, ge_, ge_impl)
-MAKE_BINARY_FUN(operator==, eq_, eq_impl)
-MAKE_BINARY_FUN(operator!=, ne_, ne_impl)
-MAKE_BINARY_FUN(operator<, lt_, lt_impl)
-MAKE_BINARY_FUN(operator>, gt_, gt_impl)
+MAKE_BINARY_FUN(operator<=, le_, le_impl, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(operator>=, ge_, ge_impl, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(operator==, eq_, eq_impl, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(operator!=, ne_, ne_impl, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(operator<, lt_, lt_impl, GPU_BINARY_AST_NONE())
+MAKE_BINARY_FUN(operator>, gt_, gt_impl, GPU_BINARY_AST_NONE())
 #undef TO_PASSIVE
 #undef le_impl
 #undef ge_impl
@@ -239,6 +296,8 @@ MAKE_BINARY_FUN(operator>, gt_, gt_impl)
 #undef gt_impl
 
 #undef MAKE_BINARY_FUN
+#undef GPU_BINARY_AST
+#undef GPU_BINARY_AST_NONE
 
 /// @}
 }  // namespace VecExpr
diff --git a/Common/include/toolboxes/allocation_toolbox.hpp b/Common/include/toolboxes/allocation_toolbox.hpp
index 9c357405b75..35eeda15353 100644
--- a/Common/include/toolboxes/allocation_toolbox.hpp
+++ b/Common/include/toolboxes/allocation_toolbox.hpp
@@ -98,6 +98,51 @@ inline void aligned_free(T* ptr) noexcept {
 }  // namespace MemoryAllocation
 
 namespace GPUMemoryAllocation {
+/*!
+ * \brief Unified Memory support verification
+ * \return 1 if UM is supported, false otherwise
+ */
+inline bool UMSupported() noexcept {
+#if defined(HAVE_CUDA)
+  static int managed = -1;
+  if (managed == -1) {
+    int device = 0;  // one device per process supposed
+    gpuErrChk(cudaGetDevice(&device));
+    gpuErrChk(cudaDeviceGetAttribute(&managed, cudaDevAttrManagedMemory, device));
+  }
+  return managed == 1;
+#else
+  return false;
+#endif
+}
+inline int GetCurrentDevice() noexcept {
+#ifdef HAVE_CUDA
+  static int device = -1;
+  if (device == -1) gpuErrChk(cudaGetDevice(&device));
+  return device;
+#else
+  return -1;
+#endif
+}
+
+/*!
+ * \brief Memory allocation for variables through CUDA Unified Memory: one pointer valid on host and device
+ * \param[in] size in bytes.
+ * \tparam ZeroInit, initialize memory to 0.
+ * \return Pointer to memory, always use gpu_free to deallocate.
+ */
+template <class T, bool ZeroInit = false>
+inline T* gpu_um_alloc(size_t size) noexcept {
+  void* ptr = nullptr;
+
+#if defined(HAVE_CUDA)
+  gpuErrChk(cudaMallocManaged((void**)(&ptr), size));
+  if (ZeroInit) std::memset(ptr, 0, size);
+#else
+  return 0;
+#endif
+  return static_cast<T*>(ptr);
+}
 /*!
  * \brief Memory allocation for variables on the GPU.
  * \param[in] size in bytes.
diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp
index f02d908a059..a865859633d 100644
--- a/Common/src/linear_algebra/CSysMatrix.cpp
+++ b/Common/src/linear_algebra/CSysMatrix.cpp
@@ -104,7 +104,6 @@ CSysMatrix<ScalarType>::~CSysMatrix() {
   MemoryAllocation::aligned_free(mat.d);
   MemoryAllocation::aligned_free(mat.l);
   MemoryAllocation::aligned_free(mat.u);
-  MemoryAllocation::aligned_free(invM);
 
   if (useCuda) {
     GPUMemoryAllocation::gpu_free(gpu.d);
@@ -116,6 +115,15 @@ CSysMatrix<ScalarType>::~CSysMatrix() {
     GPUMemoryAllocation::gpu_free(gpu.col_ind_u);
   }
 
+  if (invM_is_managed) {
+    GPUMemoryAllocation::gpu_free(invM);
+  } else {
+    MemoryAllocation::aligned_free(invM);
+    if (useCuda) {
+      GPUMemoryAllocation::gpu_free(d_invM);
+    }
+  }
+
 #ifdef USE_MKL
   mkl_jit_destroy(MatrixMatrixProductJitter);
   mkl_jit_destroy(MatrixVectorProductJitterBetaZero);
@@ -241,7 +249,18 @@ void CSysMatrix<ScalarType>::Initialize(unsigned long npoint, unsigned long npoi
     allocAndInit(ilu.u, ilu.nnz_u * nVar * nEqn);
   }
 
-  if (diag_needed) allocAndInit(invM, nPointDomain * nVar * nEqn);
+  if (diag_needed) {
+    if (useCuda && GPUMemoryAllocation::UMSupported()) {
+      invM = GPUMemoryAllocation::gpu_um_alloc<ScalarType, true>(nPointDomain * nVar * nEqn * sizeof(ScalarType));
+      d_invM = invM;  // temporary alias for testing
+      invM_is_managed = true;
+    } else {
+      allocAndInit(invM, nPointDomain * nVar * nEqn);
+      if (useCuda) {
+        d_invM = GPUMemoryAllocation::gpu_alloc<ScalarType, true>(nPointDomain * nVar * nEqn * sizeof(ScalarType));
+      }
+    }
+  }
 
   /*--- Thread parallel initialization. ---*/
 
diff --git a/Common/src/linear_algebra/CSysMatrixGPU.cu b/Common/src/linear_algebra/CSysMatrixGPU.cu
index 1ebba30097c..0cd31d966c6 100644
--- a/Common/src/linear_algebra/CSysMatrixGPU.cu
+++ b/Common/src/linear_algebra/CSysMatrixGPU.cu
@@ -1,7 +1,7 @@
 /*!
  * \file CSysMatrixGPU.cu
  * \brief Implementations of Kernels and Functions for Matrix Operations on the GPU
- * \author A. Raj
+ * \author A. Raj, Jesse Li, D. Di Giusto
  * \version 8.5.0 "Harrier"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -26,7 +26,32 @@
  */
 
 #include "../../include/linear_algebra/CSysMatrix.hpp"
+#include "../../include/linear_algebra/CSysMatrix.inl"
 #include "../../include/linear_algebra/GPUComms.cuh"
+#include "../../include/linear_algebra/CSysVector.hpp"
+
+/*!
+ * \brief Matrix-vector product kernel.
+ */
+template<typename matrixType, typename vectorType>
+__global__ void GPUMatrixVectorProductKernel(matrixType *invM, vectorType* vec, vectorType* prod, unsigned long nPointDomain, unsigned long nVar)
+{
+
+  const unsigned long iPoint = blockIdx.x * blockDim.x + threadIdx.x;
+  if (iPoint >= nPointDomain) return;
+
+  const auto block = &invM[iPoint * nVar * nVar];
+  const auto rhs = &vec[iPoint * nVar];
+  auto out = &prod[iPoint * nVar];
+
+  for (auto iVar = 0; iVar < nVar; ++iVar) {
+    vectorType sum = vectorType(0);
+    for (auto jVar = 0; jVar < nVar; ++jVar) {
+      sum += block[iVar * nVar + jVar] * rhs[jVar];
+    }
+    out[iVar] = sum;
+  }
+}
 
 /*!
  * \brief Block-LDU SpMV kernel: y[iRow] = (L + D + U) * x per block-row.
@@ -82,9 +107,9 @@ void CSysMatrix<ScalarType>::GPUMatrixVectorProduct(const CSysVector<ScalarType>
     SU2_MPI::Error("CUDA CSysMatrix block-LDU SpMV requires square blocks.", CURRENT_FUNCTION);
   }
 
+  vec.SyncToDevice();
   ScalarType* d_vec = vec.GetDevicePointer();
   ScalarType* d_prod = prod.GetDevicePointer();
-  vec.HtDTransfer();
 
   dim3 blockDim(static_cast<unsigned>(nVar), 1, 1);
   dim3 gridDim(static_cast<unsigned>(nPointDomain), 1, 1);
@@ -93,7 +118,54 @@ void CSysMatrix<ScalarType>::GPUMatrixVectorProduct(const CSysVector<ScalarType>
       gpu.row_ptr_u, gpu.col_ind_u, gpu.u, d_vec, d_prod);
   gpuErrChk(cudaGetLastError());
 
-  prod.DtHTransfer();
+  prod.MarkDeviceDirty();
+}
+
+template <class ScalarType>
+void CSysMatrix<ScalarType>::GPUComputeJacobiPreconditioner(const CSysVector<ScalarType>& vec,
+                                                            CSysVector<ScalarType>& prod, CGeometry* geometry,
+                                                            const CConfig* config) const {
+  SU2_ZONE_SCOPED
+  /*--- Apply Jacobi preconditioner, y = D^{-1} * x, the inverse of the diagonal is already known and synced to device ---*/
+
+  vec.SyncToDevice();
+  ScalarType* d_vec = vec.GetDevicePointer();
+  ScalarType* d_prod = prod.GetDevicePointer();
+
+  prod.GPUSetVal(0.0);
+
+  dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1);
+  int gridx = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nPointDomain);
+  dim3 gridDim(gridx, 1, 1);
+
+  GPUMatrixVectorProductKernel<<<gridDim, blockDim>>>(d_invM, d_vec, d_prod, nPointDomain, nVar);
+  gpuErrChk( cudaPeekAtLastError() );
+
+  prod.MarkDeviceDirty(); // this triggers () operators to sync data from device when preparing mpi buffers
+
+  /*--- MPI Parallelization ---*/
+  CSysMatrixComms::Initiate(prod, geometry, config);
+  CSysMatrixComms::Complete(prod, geometry, config);
+
+  prod.SyncToDevice(); // this will trigger migration if mpi comms updated the host data
+
+}
+
+template <class ScalarType>
+void CSysMatrix<ScalarType>::GPUBuildJacobiPreconditioner() {
+  SU2_ZONE_SCOPED
+  /*--- Build Jacobi preconditioner (M = D), compute and store the inverses of the diagonal blocks. ---*/
+  SU2_OMP_FOR_DYN(omp_heavy_size)
+  for (unsigned long iPoint = 0; iPoint < nPointDomain; iPoint++)
+    InverseDiagonalBlock(iPoint, &(invM[iPoint * nVar * nVar]));
+  END_SU2_OMP_FOR
+
+  //copy to device or prefetch
+  if (invM_is_managed) {
+    gpu_um_prefetch(d_invM, nPointDomain * nVar * nVar * sizeof(ScalarType), GPUMemoryAllocation::GetCurrentDevice());
+  } else {
+    gpuErrChk(cudaMemcpy(d_invM, invM, nPointDomain * nVar * nVar * sizeof(ScalarType), cudaMemcpyHostToDevice));
+  }
 }
 
 template class CSysMatrix<su2mixedfloat>; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.
diff --git a/Common/src/linear_algebra/CSysSolve.cpp b/Common/src/linear_algebra/CSysSolve.cpp
index aae5d9ce707..14fa8d7b6de 100644
--- a/Common/src/linear_algebra/CSysSolve.cpp
+++ b/Common/src/linear_algebra/CSysSolve.cpp
@@ -33,6 +33,7 @@
 #include "../../include/linear_algebra/CSysMatrix.hpp"
 #include "../../include/linear_algebra/CMatrixVectorProduct.hpp"
 #include "../../include/linear_algebra/CPreconditioner.hpp"
+#include "../../include/linear_algebra/CSysVector.hpp"
 
 SU2_IGNORE_WARNING("-Wmaybe-uninitialized")
 #include "Eigen/Eigenvalues"
@@ -107,6 +108,15 @@ void LinearCombinationImpl(const unsigned long n, const Vectors& vs, const Weigh
 template <class ScalarType, class Weights>
 void LinearCombinationImpl(const unsigned long n, const std::vector<CSysVector<ScalarType>>& vs, const Weights& ws,
                            CSysVector<ScalarType>& v, bool inc = false) {
+#ifdef HAVE_CUDA
+  std::vector<ScalarType> ws_host(n);  // collect weights into simple host array
+  for (unsigned long i = 0; i < n; ++i) {
+    ws_host[i] = static_cast<ScalarType>(ws(i));
+  }
+  CSysVector<ScalarType>::LinearCombinationGPU(n, vs, ws_host.data(), v, inc);
+  return;
+#endif
+
   LinearCombinationImpl(
       n, [&vs](auto i) -> auto& { return vs[i]; }, ws, v, inc);
 }
diff --git a/Common/src/linear_algebra/CSysVector.cpp b/Common/src/linear_algebra/CSysVector.cpp
index f7df34633a0..cbbdf1e3e4a 100644
--- a/Common/src/linear_algebra/CSysVector.cpp
+++ b/Common/src/linear_algebra/CSysVector.cpp
@@ -52,7 +52,15 @@ void CSysVector<ScalarType>::Initialize(unsigned long numBlk, unsigned long numB
 
   if (vec_val == nullptr) vec_val = MemoryAllocation::aligned_alloc<ScalarType, true>(64, nElm * sizeof(ScalarType));
 
-  d_vec_val = GPUMemoryAllocation::gpu_alloc<ScalarType, true>(nElm * sizeof(ScalarType));
+#ifdef HAVE_CUDA
+  useCuda = true;
+#endif
+  if (useCuda) {
+    if (d_vec_val == nullptr) {
+      d_vec_val = GPUMemoryAllocation::gpu_alloc<ScalarType, true>(nElm * sizeof(ScalarType));
+    }
+    MarkHostDirty();
+  }
 
 #ifdef HAVE_OMP
   dot_scratch.reset(new ScalarType[omp_get_max_threads()]);
@@ -73,6 +81,11 @@ const su2matrix<ScalarType>& CSysVector<ScalarType>::multiDot(const std::vector<
                                                               const std::vector<CSysVector<ScalarType>>& W,
                                                               const size_t m) {
   SU2_ZONE_SCOPED
+
+#ifdef HAVE_CUDA
+  return multiDotGPU(V, i0, n, W, m);
+#endif
+
   static constexpr size_t BLOCK_SIZE = 1024;
   static su2matrix<ScalarType> shared;
 
@@ -139,7 +152,9 @@ CSysVector<ScalarType>::~CSysVector() {
   }
   MemoryAllocation::aligned_free(vec_val);
 
-  GPUMemoryAllocation::gpu_free(d_vec_val);
+  if (useCuda) {
+    GPUMemoryAllocation::gpu_free(d_vec_val);
+  }
 }
 
 /*--- Explicit instantiations ---*/
diff --git a/Common/src/linear_algebra/CSysVectorGPU.cu b/Common/src/linear_algebra/CSysVectorGPU.cu
index 94ec17bb88f..be451800599 100644
--- a/Common/src/linear_algebra/CSysVectorGPU.cu
+++ b/Common/src/linear_algebra/CSysVectorGPU.cu
@@ -1,7 +1,7 @@
 /*!
  * \file CSysVectorGPU.cu
  * \brief Implementations of Kernels and Functions for Vector Operations on the GPU
- * \author A. Raj
+ * \author A. Raj, D. Di giusto
  * \version 8.5.0 "Harrier"
  *
  * SU2 Project Website: https://su2code.github.io
@@ -27,6 +27,181 @@
 
 #include "../../include/linear_algebra/CSysVector.hpp"
 #include "../../include/linear_algebra/GPUComms.cuh"
+#include "../../include/linear_algebra/gpu_ast.hpp"
+#include <execinfo.h> // Needed for backtrace and backtrace_symbols
+#include <cstdlib>
+
+/*!
+ * \brief evaluates the abstract syntax tree iteratively.
+ * The tree is unrolled to prevent recursion and undefinite memory mapping.
+ */
+template <typename ScalarType>
+__device__ __forceinline__ ScalarType EvaluateASTIterative(const GPUASTPayload<ScalarType>& payload,
+                                                            unsigned long elem_idx) {
+  ScalarType values[MAX_AST_NODES];
+
+#pragma unroll
+  for (int idx = 0; idx < payload.node_count; ++idx) {
+    const GPUASTNode<ScalarType>& node = payload.nodes[idx];
+    switch (node.op) {
+      case GPUOpType::VEC:    values[idx] = node.d_ptr[elem_idx]; break;
+      case GPUOpType::SCALAR: values[idx] = node.val; break;
+      case GPUOpType::ADD:    values[idx] = values[node.left] + values[node.right]; break;
+      case GPUOpType::SUB:    values[idx] = values[node.left] - values[node.right]; break;
+      case GPUOpType::MUL:    values[idx] = values[node.left] * values[node.right]; break;
+      case GPUOpType::DIV:    values[idx] = values[node.left] / values[node.right]; break;
+      case GPUOpType::NEG:    values[idx] = -values[node.left]; break;
+    }
+  }
+  return values[payload.root_idx];
+}
+
+/*!
+ * \brief kernel evaluating the abstract syntax tree for each thread
+ */
+template <typename ScalarType>
+__global__ void GPUASTInterpreterKernel(ScalarType* __restrict__ d_out, GPUASTPayload<ScalarType> payload, unsigned long n) {
+  unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x;
+  if (idx < n) {
+    d_out[idx] = EvaluateASTIterative(payload, idx);
+  }
+}
+
+/*!
+ * \brief method to configure and launch the kernel executing the abstract syntax tree for a given vector operators
+ */
+template <class ScalarType>
+void CSysVector<ScalarType>::LaunchGPUAST(const GPUASTPayload<ScalarType>& payload, unsigned long nElm) {
+
+  dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1);
+  int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm);
+  dim3 gridDim(numBlocks, 1, 1);
+
+  this->SyncToDevice();
+
+  GPUASTInterpreterKernel<<<gridDim, blockDim>>>(this->d_vec_val, payload, nElm);
+  gpuErrChk(cudaPeekAtLastError());
+
+  this->MarkDeviceDirty();
+}
+
+/*!
+ * \brief namespace defining the scalar operators for GPU
+ */
+namespace {
+template <class T> struct OpSetScalar { T val; __device__ __forceinline__ T operator()(T x) const { return val; } };
+template <class T> struct OpAddScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x + val; } };
+template <class T> struct OpSubScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x - val; } };
+template <class T> struct OpMulScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x * val; } };
+template <class T> struct OpDivScalar { T val; __device__ __forceinline__ T operator()(T x) const { return x / val; } };
+}  // namespace
+
+/*!
+ * \brief generic Unary Operation Kernel to apply given functors on GPU
+ */
+template <class ScalarType, class Operator>
+__global__ void GPUUnaryOperationKernel(ScalarType* __restrict__ vec, unsigned long n, Operator op) {
+  const unsigned long idx = blockIdx.x * blockDim.x + threadIdx.x;
+  if (idx < n) vec[idx] = op(vec[idx]);
+}
+
+/*!
+ * \brief dispatching method for the GPU unary operations
+ */
+template <class ScalarType>
+void CSysVector<ScalarType>::GPUUnaryOperation(GPUScalarOp op, ScalarType val) {
+
+  dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1);
+  int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElm);
+  dim3 gridDim(numBlocks, 1, 1);
+
+  this->SyncToDevice();
+
+  switch (op) {
+    case GPUScalarOp::SET:
+      GPUUnaryOperationKernel<<<gridDim, blockDim>>>(this->d_vec_val, this->nElm, OpSetScalar<ScalarType>{val});
+      break;
+    case GPUScalarOp::ADD:
+      GPUUnaryOperationKernel<<<gridDim, blockDim>>>(this->d_vec_val, this->nElm, OpAddScalar<ScalarType>{val});
+      break;
+    case GPUScalarOp::SUB:
+      GPUUnaryOperationKernel<<<gridDim, blockDim>>>(this->d_vec_val, this->nElm, OpSubScalar<ScalarType>{val});
+      break;
+    case GPUScalarOp::MUL:
+      GPUUnaryOperationKernel<<<gridDim, blockDim>>>(this->d_vec_val, this->nElm, OpMulScalar<ScalarType>{val});
+      break;
+    case GPUScalarOp::DIV:
+      GPUUnaryOperationKernel<<<gridDim, blockDim>>>(this->d_vec_val, this->nElm, OpDivScalar<ScalarType>{val});
+      break;
+  }
+  gpuErrChk(cudaPeekAtLastError());
+
+  this->MarkDeviceDirty();
+}
+
+
+/*!
+ * \brief GPU dot prodcut kernel
+ * \brief block-level reduction of elementwise products, accumulated into a single device scalar.
+ */
+template <class ScalarType>
+__global__ void GPUDotKernel(const ScalarType* __restrict__ a, const ScalarType* __restrict__ b,
+                             unsigned long n, ScalarType* __restrict__ result) {
+  //shared memory
+  extern __shared__ unsigned char smem_raw[];
+  ScalarType* sdata = reinterpret_cast<ScalarType*>(smem_raw);
+
+  //local and global thread indexes for access and block reduction
+  const unsigned long tid = threadIdx.x;
+  unsigned long idx = blockIdx.x * blockDim.x + tid;
+  const unsigned long stride = blockDim.x * gridDim.x;
+
+  //thread reduction
+  ScalarType local = ScalarType(0);
+  for (; idx < n; idx += stride) local += a[idx] * b[idx];
+
+  sdata[tid] = local;
+  __syncthreads();
+
+  //block reduction
+  for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
+      if (tid < s) sdata[tid] += sdata[tid + s];
+      __syncthreads();
+  }
+
+  //final atomic add per block
+  if (tid == 0) atomicAdd(result, sdata[0]);
+}
+
+/*!
+ * \brief GPU dot product method between this and other CSysVector
+ * \note this is a vectors read-only kernel returning a scalar on host
+ */
+template <class ScalarType>
+ScalarType CSysVector<ScalarType>::GPUDot(const CSysVector& other) const {
+
+  dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE, 1, 1);
+  int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, this->nElmDomain);
+  dim3 gridDim(numBlocks, 1, 1);
+
+  SyncToDevice();
+  other.SyncToDevice();
+
+  // allocate and zero the result scalar
+  ScalarType* d_dot_result;
+  gpuErrChk(cudaMalloc(&d_dot_result, sizeof(ScalarType)));
+  gpuErrChk(cudaMemset(d_dot_result, 0, sizeof(ScalarType)));
+
+  const size_t sharedBytes = KernelParameters::MVP_BLOCK_SIZE * sizeof(ScalarType);
+  GPUDotKernel<<<gridDim, blockDim, sharedBytes>>>(this->d_vec_val, other.d_vec_val, this->nElmDomain, d_dot_result);
+  gpuErrChk(cudaPeekAtLastError());
+
+  ScalarType result;
+  gpuErrChk(cudaMemcpy(&result, d_dot_result, sizeof(ScalarType), cudaMemcpyDeviceToHost));
+  gpuErrChk(cudaFree(d_dot_result));
+
+  return result;
+}
 
 template<class ScalarType>
 void CSysVector<ScalarType>::HtDTransfer(bool trigger) const
@@ -46,4 +221,189 @@ void CSysVector<ScalarType>::GPUSetVal(ScalarType val, bool trigger) const
    if(trigger) gpuErrChk(cudaMemset((void*)(d_vec_val), val, (sizeof(ScalarType)*nElm)));
 }
 
+/*!
+ * \brief multi vector product CUDA kernel one line of blocks per pair V[i0+i],W[j];
+ *        Configurable multiple blocks reducing over the size of the vectors
+ */
+template <class ScalarType>
+__global__ void GPUmultiDot(const ScalarType* const* __restrict__ d_V, const size_t n,
+                            const ScalarType* const* __restrict__ d_W, const size_t m, const size_t size,
+                            ScalarType* __restrict__ d_local)
+{
+  // Map each x,y block to the specific (i,j) dot product
+  const size_t pair_idx = blockIdx.y;
+  if (pair_idx >= n * m) return;
+
+  const size_t i = pair_idx / m;
+  const size_t j = pair_idx % m;
+
+  //get the corresponding vectors
+  const ScalarType* __restrict__ vi = d_V[i];
+  const ScalarType* __restrict__ wj = d_W[j];
+
+  // grid strided loop over the vector elements
+  ScalarType local_sum = 0.0;
+  const size_t tid = blockIdx.x * blockDim.x + threadIdx.x;
+  const size_t stride = gridDim.x * blockDim.x;
+
+  for (size_t k = tid; k < size; k += stride)
+  {
+      local_sum += vi[k] * wj[k];
+  }
+
+  // shared memory reduction within the block
+  extern __shared__ char shared_mem[];
+  ScalarType* sdata = reinterpret_cast<ScalarType*>(shared_mem);
+
+  sdata[threadIdx.x] = local_sum;
+  __syncthreads();
+
+  // parallel reduction on the block
+  for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1)
+  {
+      if (threadIdx.x < s)
+      {
+        sdata[threadIdx.x] += sdata[threadIdx.x + s];
+      }
+      __syncthreads();
+  }
+
+  // atomic add of each block partial sum to the output matrix, operated by thread 0 of each block
+  if (threadIdx.x == 0)
+  {
+      atomicAdd(&d_local[i * m + j], sdata[0]);
+  }
+}
+
+/*!
+ * \brief multi vector dot produt method for GPU
+ * \note this is a vectors-read only method that returns an array of scalars
+ */
+template <class ScalarType>
+const su2matrix<ScalarType>& CSysVector<ScalarType>::multiDotGPU(const std::vector<CSysVector<ScalarType>>& V,
+                                                                 const size_t i0, const size_t n,
+                                                                 const std::vector<CSysVector<ScalarType>>& W,
+                                                                 const size_t m) {
+  SU2_ZONE_SCOPED
+
+  static su2matrix<ScalarType> shared;
+  if (n == 0 || m == 0) return shared;
+
+  const size_t size = V[0].nElmDomain;
+
+  // get all the device pointers
+  std::vector<const ScalarType*> h_V_ptrs(n), h_W_ptrs(m);
+  for (size_t i = 0; i < n; ++i){
+    V[i0 + i].SyncToDevice();
+    h_V_ptrs[i] = V[i0 + i].GetDevicePointer();
+  }
+  for (size_t j = 0; j < m; ++j){
+    W[j].SyncToDevice();
+    h_W_ptrs[j] = W[j].GetDevicePointer();
+  }
+
+  // copy the pointers to the device arrays of pointers
+  const ScalarType** d_V_ptrs;
+  const ScalarType** d_W_ptrs;
+  gpuErrChk(cudaMalloc(&d_V_ptrs, n * sizeof(ScalarType*)));
+  gpuErrChk(cudaMalloc(&d_W_ptrs, m * sizeof(ScalarType*)));
+  gpuErrChk(cudaMemcpy(d_V_ptrs, h_V_ptrs.data(), n * sizeof(ScalarType*), cudaMemcpyHostToDevice));
+  gpuErrChk(cudaMemcpy(d_W_ptrs, h_W_ptrs.data(), m * sizeof(ScalarType*), cudaMemcpyHostToDevice));
+
+  // allocate result buffer, zero it
+  ScalarType* d_local;
+  gpuErrChk(cudaMalloc(&d_local, n * m * sizeof(ScalarType)));
+  gpuErrChk(cudaMemset(d_local, 0, n * m * sizeof(ScalarType)));
+
+  dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1);
+  int numBlocksPerPair = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, size);
+  dim3 gridDim(numBlocksPerPair, n * m, 1);
+
+  GPUmultiDot<<<gridDim, blockDim, KernelParameters::MVP_BLOCK_SIZE * sizeof(ScalarType)>>>(d_V_ptrs, n, d_W_ptrs, m, size, d_local);
+  gpuErrChk(cudaGetLastError());
+
+  // copy result to host for MPI reduce
+  su2matrix<ScalarType> local(n,m);
+  gpuErrChk(cudaMemcpy(local.data(), d_local, n * m * sizeof(ScalarType), cudaMemcpyDeviceToHost));
+
+  /*--- Single AllReduce of the result, only the master thread communicates. ---*/
+  SU2_OMP_MASTER {
+    shared.resize(n, m);
+
+    const auto mpi_type = (sizeof(ScalarType) < sizeof(double)) ? MPI_FLOAT : MPI_DOUBLE;
+    SelectMPIWrapper<ScalarType>::W::Allreduce(local.data(), shared.data(), n * m, mpi_type, MPI_SUM,
+                                               SU2_MPI::GetComm());
+  }
+  END_SU2_OMP_MASTER
+
+  /*--- All threads have the same view of the result. ---*/
+  SU2_OMP_BARRIER
+
+
+  // clean allocations
+  gpuErrChk(cudaFree(d_local));
+  gpuErrChk(cudaFree(d_V_ptrs));
+  gpuErrChk(cudaFree(d_W_ptrs));
+
+  return shared;
+}
+
+template<class ScalarType, int N>
+struct WeightedVecs {
+  const ScalarType* ptrs[N];
+  ScalarType weights[N];
+};
+
+/*!
+ * \brief linear combination kernel to calculate the next vector v from weights and vectors
+ */
+template<class ScalarType, int N>
+__global__ void LinearCombinationKernel(ScalarType* __restrict__ v, WeightedVecs<ScalarType, N> wv,
+                                        int n, unsigned long nElm, bool inc)
+{
+  const unsigned long k = blockIdx.x * blockDim.x + threadIdx.x;
+  if (k >= nElm) return;
+
+  //handle overwriting or combination with existing
+  ScalarType result = inc ? v[k] : ScalarType(0);
+
+  #pragma unroll
+  for (int i = 0; i < N; ++i) // N is known at compile time (4), this unrolls to: if (i < n) result += weight[i] * vector[i][k]; i<4
+    if (i < n) result += wv.weights[i] * wv.ptrs[i][k];
+  v[k] = result;
+}
+
+/*!
+ * \brief dispatcher for the linear combination kernel on GPU
+ */
+template<class ScalarType>
+void CSysVector<ScalarType>::LinearCombinationGPU(const unsigned long n, const std::vector<CSysVector<ScalarType>>& vs, const ScalarType* ws,
+                                                  CSysVector<ScalarType>& v, bool inc)
+{
+  const unsigned long nElm = v.nElmDomain;
+  dim3 blockDim(KernelParameters::MVP_BLOCK_SIZE,1,1);
+  int numBlocks = KernelParameters::round_up_division(KernelParameters::MVP_BLOCK_SIZE, nElm);
+  dim3 gridDim(numBlocks, 1, 1);
+
+  v.SyncToDevice();
+  ScalarType* d_v = v.GetDevicePointer();
+
+  for (unsigned long i = 0; i < n; i += 4) {
+    const int rem = static_cast<int>(std::min(n - i, 4ul));
+    //prepare vectors pointers and corresponding weights, passing them by value
+    WeightedVecs<ScalarType, 4> vs_ws = {};
+    for (int j = 0; j < rem; ++j) {
+      vs[i + j].SyncToDevice();                       //ensure is on GPU, was copied in multiDot
+      vs_ws.ptrs[j] = vs[i + j].GetDevicePointer();  // get the pointer
+      vs_ws.weights[j] = ws[i + j];                  // plain array indexing, not ws(k)
+    }
+    //calculate the linear combination on GPU, handle more than 4 vectors through inc || i > 0
+    LinearCombinationKernel<ScalarType, 4><<<gridDim, blockDim>>>(d_v, vs_ws, rem, nElm, inc || i > 0);
+    gpuErrChk(cudaPeekAtLastError());
+  }
+
+  v.MarkDeviceDirty();
+
+}
+
 template class CSysVector<su2double>; //This is a temporary fix for invalid instantiations due to separating the member function from the header file the class is defined in. Will try to rectify it in coming commits.