diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index 8ef511660f1fbc3760ab82b63f3603c517634524..888acb06889073c1a8ae8115486a6cc38e484ab5 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -3,7 +3,6 @@
 
 #include <PastisMacros.hpp>
 #include <PastisAssert.hpp>
-#include <Kokkos_Core.hpp>
 
 #include <Types.hpp>
 #include <TinyVector.hpp>
diff --git a/src/main.cpp b/src/main.cpp
index ba6f6c504ae02f7859d722de46295e192f92fffe..064a382270a6834a2e27bf6e53b29e1cbc059c9f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -6,11 +6,6 @@
 #include <SignalManager.hpp>
 #include <ConsoleManager.hpp>
 
-// #include <RawKokkosAcousticSolver.hpp>
-// #include <MeshLessAcousticSolver.hpp>
-// #include <AcousticSolverClass.hpp>
-// #include <AcousticSolverTest.hpp>
-
 #include <Connectivity.hpp>
 
 #include <Mesh.hpp>
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 9d95f2896a073345f6d163534795143ca4e47487..90024775defe4ee51cbc6691328e4752ac1429b2 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -228,14 +228,14 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
 
   {
     CellValue<CellType> cell_type(*this);
-    Kokkos::parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         cell_type[j] = cell_type_vector[j];
       });
     m_cell_type = cell_type;
   }
   {
     CellValue<double> inv_cell_nb_nodes(*this);
-    Kokkos::parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
         inv_cell_nb_nodes[j] = 1./cell_nodes.length;
       });
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 4c83d4adc02cc241d328eb8b22f06c6c66a1ef1f..c4cd00ab7b8a3667379571f69a6281c8538be817 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -3,9 +3,10 @@
 
 #include <PastisMacros.hpp>
 #include <PastisAssert.hpp>
+#include <PastisUtils.hpp>
+
 #include <TinyVector.hpp>
 
-#include <Kokkos_Core.hpp>
 #include <ItemValue.hpp>
 
 #include <IConnectivity.hpp>
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index b8fa1e009a6f5f32a00c2e4defdc742b9872b899..512c2cfdfda322f1b410543a75eea2af7dc19061 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -1,7 +1,7 @@
 #ifndef CONNECTIVITY_MATRIX_HPP
 #define CONNECTIVITY_MATRIX_HPP
 
-#include <Kokkos_Core.hpp>
+#include <PastisUtils.hpp>
 #include <Kokkos_StaticCrsGraph.hpp>
 
 class ConnectivityMatrix
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index 70240b8efa1a4b96739b1b6814778d2fdb65e84d..38f8234c7ab4a1a1801644655a0d933a568cfced 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -4,7 +4,7 @@
 #include <ItemType.hpp>
 #include <ItemId.hpp>
 
-#include <Kokkos_Core.hpp>
+#include <PastisUtils.hpp>
 #include <ConnectivityMatrix.hpp>
 
 template <ItemType SourceItemType,
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 0f016f0e95ec577fbf5e51561a97c45a64254731..6bd7e892576d329a82de8e10519f5b29fabb6e57 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -1,7 +1,7 @@
 #ifndef MESH_DATA_HPP
 #define MESH_DATA_HPP
 
-#include <Kokkos_Core.hpp>
+#include <PastisUtils.hpp>
 #include <TinyVector.hpp>
 
 #include <ItemValue.hpp>
@@ -40,7 +40,7 @@ class MeshData
           = m_mesh.connectivity().cellToNodeMatrix();
 
       CellValue<Rd> xj(m_mesh.connectivity());
-      Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+      parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
           const auto& cell_nodes = cell_to_node_matrix[j];
           xj[j] = 0.5*(xr[cell_nodes[0]]+xr[cell_nodes[1]]);
         });
@@ -54,7 +54,7 @@ class MeshData
       const auto& cell_to_node_matrix
           = m_mesh.connectivity().cellToNodeMatrix();
       CellValue<Rd> xj(m_mesh.connectivity());
-      Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+      parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
           Rd X = zero;
           const auto& cell_nodes = cell_to_node_matrix[j];
           for (size_t R=0; R<cell_nodes.size(); ++R) {
@@ -74,7 +74,7 @@ class MeshData
         = m_mesh.connectivity().cellToNodeMatrix();
 
     CellValue<double> Vj(m_mesh.connectivity());
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         double sum_cjr_xr = 0;
         const auto& cell_nodes = cell_to_node_matrix[j];
 
@@ -98,7 +98,7 @@ class MeshData
 
       {
         NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+        parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
             const auto& cell_nodes = cell_to_node_matrix[j];
             for (size_t R=0; R<cell_nodes.size(); ++R) {
               int Rp1 = (R+1)%cell_nodes.size();
@@ -112,7 +112,7 @@ class MeshData
 
       {
         NodeValuePerCell<double> ljr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
+        parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             ljr[jr] = l2Norm(m_Cjr[jr]);
           });
         m_ljr = ljr;
@@ -120,7 +120,7 @@ class MeshData
 
       {
         NodeValuePerCell<Rd> njr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
+        parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
           });
         m_njr = njr;
@@ -132,7 +132,7 @@ class MeshData
       const auto& face_to_node_matrix
           = m_mesh.connectivity().faceToNodeMatrix();
 
-      Kokkos::parallel_for(m_mesh.numberOfFaces(), PASTIS_LAMBDA(const FaceId& l) {
+      parallel_for(m_mesh.numberOfFaces(), PASTIS_LAMBDA(const FaceId& l) {
           const auto& face_nodes = face_to_node_matrix[l];
           const size_t nb_nodes = face_nodes.size();
           std::vector<Rd> dxr(nb_nodes);
@@ -164,11 +164,11 @@ class MeshData
 
       {
         NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
-        Kokkos::parallel_for(Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
+        parallel_for(Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             Cjr[jr] = zero;
           });
 
-        Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
             const auto& cell_nodes = cell_to_node_matrix[j];
 
             const auto& cell_faces = cell_to_face_matrix[j];
@@ -209,7 +209,7 @@ class MeshData
 
       {
         NodeValuePerCell<double> ljr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
+        parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             ljr[jr] = l2Norm(m_Cjr[jr]);
           });
         m_ljr = ljr;
@@ -217,7 +217,7 @@ class MeshData
 
       {
         NodeValuePerCell<Rd> njr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
+        parallel_for(m_Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             njr[jr] = (1./m_ljr[jr])*m_Cjr[jr];
           });
         m_njr = njr;
@@ -271,7 +271,7 @@ class MeshData
       // in 1d Cjr are computed once for all
       {
         NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+        parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
             Cjr(j,0)=-1;
             Cjr(j,1)= 1;
           });
@@ -281,7 +281,7 @@ class MeshData
       m_njr = m_Cjr;
       {
         NodeValuePerCell<double> ljr(m_mesh.connectivity());
-        Kokkos::parallel_for(ljr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
+        parallel_for(ljr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             ljr[jr] = 1;
         });
         m_ljr = ljr;
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index a8026cc8609a178c891ca2a67e3739bb897c5067..9f945e07b2c8b50889e014732c370cc74d1cb45a 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -37,7 +37,7 @@ class MeshNodeBoundary
         = mesh.connectivity().faceToCellMatrix();
 
     const Array<const FaceId>& face_list = ref_face_list.faceList();
-    Kokkos::parallel_for(face_list.size(), PASTIS_LAMBDA(const int& l){
+    parallel_for(face_list.size(), PASTIS_LAMBDA(const int& l){
         const auto& face_cells = face_to_cell_matrix[face_list[l]];
         if (face_cells.size()>1) {
           std::cerr << "internal faces cannot be used to define mesh boundaries\n";
@@ -64,7 +64,7 @@ class MeshNodeBoundary
     node_ids.resize(std::distance(node_ids.begin(), last));
 
     Array<NodeId> node_list(node_ids.size());
-    Kokkos::parallel_for(node_ids.size(), PASTIS_LAMBDA(const int& r){
+    parallel_for(node_ids.size(), PASTIS_LAMBDA(const int& r){
         node_list[r] = node_ids[r];
       });
     m_node_list = node_list;
@@ -158,7 +158,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
 
   const NodeValue<const R2>& xr = mesh.xr();
 
-  Kokkos::parallel_for(m_node_list.size(), PASTIS_LAMBDA(const size_t& r) {
+  parallel_for(m_node_list.size(), PASTIS_LAMBDA(const size_t& r) {
       const R2& x = xr[m_node_list[r]];
       if ((x-origin,normal)>1E-13*length) {
         std::cerr << "this FlatBoundary is not flat!\n";
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 389e916e0c3f9c0a5b1364fc6e1bf06b7cfb9a12..3afc434969159d4d8a547040be195669e2360f52 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -39,7 +39,7 @@ class AcousticSolver
   computeRhoCj(const CellValue<const double>& rhoj,
                const CellValue<const double>& cj)
   {
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
         m_rhocj[j] = rhoj[j]*cj[j];
       });
     return m_rhocj;
@@ -51,7 +51,7 @@ class AcousticSolver
                   const NodeValuePerCell<const double>& ljr,
                   const NodeValuePerCell<const Rd>& njr)
   {
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
         const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
         const double& rho_c = rhocj[j];
         for (size_t r=0; r<nb_nodes; ++r) {
@@ -68,7 +68,7 @@ class AcousticSolver
     const auto& node_local_numbers_in_their_cells
         = m_connectivity.nodeLocalNumbersInTheirCells();
 
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
         Rdd sum = zero;
         const auto& node_to_cell = node_to_cell_matrix[r];
         const auto& node_local_number_in_its_cells
@@ -96,7 +96,7 @@ class AcousticSolver
     const auto& node_local_numbers_in_their_cells
         = m_connectivity.nodeLocalNumbersInTheirCells();
 
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
         Rd& br = m_br[r];
         br = zero;
         const auto& node_to_cell = node_to_cell_matrix[r];
@@ -144,7 +144,7 @@ class AcousticSolver
 
           const Array<const NodeId>& node_list
               = symmetry_bc.nodeList();
-          Kokkos::parallel_for(symmetry_bc.numberOfNodes(), PASTIS_LAMBDA(const int& r_number) {
+          parallel_for(symmetry_bc.numberOfNodes(), PASTIS_LAMBDA(const int& r_number) {
               const NodeId r = node_list[r_number];
 
               m_Ar[r] = P*m_Ar[r]*P + nxn;
@@ -162,7 +162,7 @@ class AcousticSolver
   {
     inverse(Ar, m_inv_Ar);
     const NodeValue<const Rdd> invAr = m_inv_Ar;
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
         m_ur[r]=invAr[r]*br[r];
       });
 
@@ -179,7 +179,7 @@ class AcousticSolver
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().cellToNodeMatrix();
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         for (size_t r=0; r<cell_nodes.size(); ++r) {
@@ -191,7 +191,7 @@ class AcousticSolver
   void inverse(const NodeValue<const Rdd>& A,
                NodeValue<Rdd>& inv_A) const
   {
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
         inv_A[r] = ::inverse(A[r]);
       });
   }
@@ -258,7 +258,7 @@ class AcousticSolver
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().cellToNodeMatrix();
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         double S = 0;
@@ -297,7 +297,7 @@ class AcousticSolver
         = m_mesh.connectivity().cellToNodeMatrix();
 
     const CellValue<const double> inv_mj = unknowns.invMj();
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         Rd momentum_fluxes = zero;
@@ -311,18 +311,18 @@ class AcousticSolver
         Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
       });
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
         ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
       });
 
     NodeValue<Rd> mutable_xr = m_mesh.mutableXr();
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r){
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r){
         mutable_xr[r] += dt*ur[r];
       });
     m_mesh_data.updateAllData();
 
     const CellValue<const double> mj = unknowns.mj();
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         rhoj[j] = mj[j]/Vj[j];
       });
   }
diff --git a/src/scheme/BlockPerfectGas.hpp b/src/scheme/BlockPerfectGas.hpp
index 13dbfd81660ddd40d1e63229f6da78eaa5a05a05..f3ea76520c656d7d38ea0f104a4a01cca9c95a85 100644
--- a/src/scheme/BlockPerfectGas.hpp
+++ b/src/scheme/BlockPerfectGas.hpp
@@ -34,7 +34,7 @@ public:
     const CellValue<const double>& e = m_ej;
     const CellValue<const double>& gamma = m_gammaj;
 
-    Kokkos::parallel_for(nj, PASTIS_LAMBDA(const CellId& j){
+    parallel_for(nj, PASTIS_LAMBDA(const CellId& j){
         const double gamma_minus_one = gamma[j]-1;
         m_pj[j] = gamma_minus_one*rho[j]*e[j];
         m_cj[j] = std::sqrt(gamma[j]*gamma_minus_one*e[j]);
@@ -48,12 +48,12 @@ public:
     const CellValue<const double>& p = m_pj;
     const CellValue<const double>& gamma = m_gammaj;
 
-    Kokkos::parallel_for(nj, PASTIS_LAMBDA(const CellId& j){
+    parallel_for(nj, PASTIS_LAMBDA(const CellId& j){
         m_ej[j] = p[j]/(rho[j]*(gamma[j]-1));
       });
 
     const CellValue<const double>& e = m_ej;
-    Kokkos::parallel_for(nj, PASTIS_LAMBDA(const CellId& j){
+    parallel_for(nj, PASTIS_LAMBDA(const CellId& j){
         m_cj[j] = std::sqrt(gamma[j]*(gamma[j]-1)*e[j]);
       });
   }
diff --git a/src/scheme/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
index ed73e14bcbd4569227dd790b9c1425ac76ced28a..0380f97c021b637a623146cd1be7f18524913e2d 100644
--- a/src/scheme/BoundaryCondition.hpp
+++ b/src/scheme/BoundaryCondition.hpp
@@ -69,7 +69,7 @@ public:
       m_number_of_faces(faces.size())
   {
     Array<unsigned int> face_list(faces.size());
-    Kokkos::parallel_for(m_number_of_faces, PASTIS_LAMBDA(const int& f){
+    parallel_for(m_number_of_faces, PASTIS_LAMBDA(const int& f){
         face_list[f]=faces[f];
       });
     m_face_list = face_list;
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index fc1178b77ed6b6b9f8f04860cd875af67ae90b5f..109e93d5a13acd2ce5c23b2b83576af21665c3c7 100644
--- a/src/scheme/FiniteVolumesEulerUnknowns.hpp
+++ b/src/scheme/FiniteVolumesEulerUnknowns.hpp
@@ -124,7 +124,7 @@ public:
   {
     const CellValue<const Rd>& xj = m_mesh_data.xj();
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         if (xj[j][0]<0.5) {
           m_rhoj[j]=1;
         } else {
@@ -132,7 +132,7 @@ public:
         }
       });
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         if (xj[j][0]<0.5) {
           m_pj[j]=1;
         } else {
@@ -140,27 +140,27 @@ public:
         }
       });
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         m_uj[j] = zero;
       });
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         m_gammaj[j] = 1.4;
       });
 
     BlockPerfectGas block_eos(m_rhoj, m_ej, m_pj, m_gammaj, m_cj);
     block_eos.updateEandCFromRhoP();
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         m_Ej[j] = m_ej[j]+0.5*(m_uj[j],m_uj[j]);
       });
 
     const CellValue<const double>& Vj = m_mesh_data.Vj();
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         m_mj[j] = m_rhoj[j] * Vj[j];
       });
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         m_inv_mj[j] = 1./m_mj[j];
       });
   }
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index bfcb6231dca981bf2a10898fecce79343d5a0167..e7273ecee82142014b378a55980d66c0305f016a 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -2,7 +2,7 @@
 #define ARRAY_HPP
 
 #include <PastisMacros.hpp>
-#include <Kokkos_Core.hpp>
+#include <PastisUtils.hpp>
 
 #include <PastisAssert.hpp>
 
diff --git a/src/utils/ArrayUtils.hpp b/src/utils/ArrayUtils.hpp
index 26f4b623307497b6bae7b7f0416b5eef322284d7..b64452d992b9728dfc2e3a117f790d8893ec155c 100644
--- a/src/utils/ArrayUtils.hpp
+++ b/src/utils/ArrayUtils.hpp
@@ -2,7 +2,7 @@
 #define ARRAY_UTILS_HPP
 
 #include <PastisMacros.hpp>
-#include <Kokkos_Core.hpp>
+#include <PastisUtils.hpp>
 
 template <typename ArrayType>
 class ReduceMin
@@ -17,7 +17,7 @@ class ReduceMin
   operator data_type()
   {
     data_type reduced_value;
-    Kokkos::parallel_reduce(m_array.size(), *this, reduced_value);
+    parallel_reduce(m_array.size(), *this, reduced_value);
     return reduced_value;
   }
 
@@ -69,7 +69,7 @@ class ReduceMax
   operator data_type()
   {
     data_type reduced_value;
-    Kokkos::parallel_reduce(m_array.size(), *this, reduced_value);
+    parallel_reduce(m_array.size(), *this, reduced_value);
     return reduced_value;
   }
 
diff --git a/src/utils/PastisUtils.hpp b/src/utils/PastisUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..adf12176fa47847b439aa334a4f7f9d2ada1468f
--- /dev/null
+++ b/src/utils/PastisUtils.hpp
@@ -0,0 +1,29 @@
+#ifndef PASTIS_UTILS_HPP
+#define PASTIS_UTILS_HPP
+
+#include <Kokkos_Core.hpp>
+#include <PastisMacros.hpp>
+
+#include <functional>
+#include <string>
+
+template <typename FunctionType>
+PASTIS_FORCEINLINE
+void parallel_for(const size_t& size,
+                  const FunctionType& lambda,
+                  const std::string& label = "")
+{
+  Kokkos::parallel_for(size, lambda, label);
+}
+
+template <typename ArrayType, typename ReturnType>
+PASTIS_FORCEINLINE
+void parallel_reduce(const size_t& size,
+                     const ArrayType& array,
+                     ReturnType& value,
+                     const std::string& label = "")
+{
+  Kokkos::parallel_reduce(label, size, array, value);
+}
+
+#endif // PASTIS_UTILS_HPP