diff --git a/src/algebra/TinyMatrix.hpp b/src/algebra/TinyMatrix.hpp
index eda3b5f19c3f7950bda6acce50a3e39b7c9f2cec..888acb06889073c1a8ae8115486a6cc38e484ab5 100644
--- a/src/algebra/TinyMatrix.hpp
+++ b/src/algebra/TinyMatrix.hpp
@@ -1,8 +1,8 @@
 #ifndef TINY_MATRIX_HPP
 #define TINY_MATRIX_HPP
 
+#include <PastisMacros.hpp>
 #include <PastisAssert.hpp>
-#include <Kokkos_Core.hpp>
 
 #include <Types.hpp>
 #include <TinyVector.hpp>
@@ -16,14 +16,14 @@ private:
   T m_values[N*N];
   static_assert((N>0),"TinyMatrix size must be strictly positive");
 
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   constexpr size_t _index(const size_t& i, const size_t& j) const noexcept
   {
     return std::move(i*N+j);
   }
 
   template <typename... Args>
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
   {
     m_values[N*N-1-sizeof...(args)] = t;
@@ -33,7 +33,7 @@ private:
   }
 
 public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix operator-() const
   {
     TinyMatrix opposed;
@@ -43,14 +43,14 @@ public:
     return std::move(opposed);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr friend TinyMatrix operator*(const T& t, const TinyMatrix& A)
   {
     TinyMatrix B = A;
     return std::move(B *= t);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator*=(const T& t)
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -59,7 +59,7 @@ public:
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix operator*(const TinyMatrix& B) const
   {
     const TinyMatrix& A = *this;
@@ -76,7 +76,7 @@ public:
     return std::move(AB);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector<N,T> operator*(const TinyVector<N,T>& x) const
   {
     const TinyMatrix& A = *this;
@@ -91,7 +91,7 @@ public:
     return std::move(Ax);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr friend std::ostream& operator<<(std::ostream& os, const TinyMatrix& A)
   {
     if constexpr(N==1) {
@@ -111,7 +111,7 @@ public:
     return os;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr bool operator==(const TinyMatrix& A) const
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -120,13 +120,13 @@ public:
     return true;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr bool operator!=(const TinyMatrix& A) const
   {
     return not this->operator==(A);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix operator+(const TinyMatrix& A) const
   {
     TinyMatrix sum;
@@ -136,7 +136,7 @@ public:
     return std::move(sum);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix operator+(TinyMatrix&& A) const
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -145,7 +145,7 @@ public:
     return std::move(A);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix operator-(const TinyMatrix& A) const
   {
     TinyMatrix difference;
@@ -155,7 +155,7 @@ public:
     return std::move(difference);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix operator-(TinyMatrix&& A) const
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -164,7 +164,7 @@ public:
     return std::move(A);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator+=(const TinyMatrix& A)
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -173,7 +173,7 @@ public:
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator-=(const TinyMatrix& A)
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -182,21 +182,21 @@ public:
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr T& operator()(const size_t& i, const size_t& j) noexcept(NO_ASSERT)
   {
     Assert((i<N) and (j<N));
     return m_values[_index(i,j)];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr const T& operator()(const size_t& i, const size_t& j) const noexcept(NO_ASSERT)
   {
     Assert((i<N) and (j<N));
     return m_values[_index(i,j)];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator=(const ZeroType& z) noexcept
   {
     static_assert(std::is_arithmetic<T>(),"Cannot assign 'zero' value for non-arithmetic types");
@@ -206,7 +206,7 @@ public:
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator=(const IdentityType& I) noexcept
   {
     static_assert(std::is_arithmetic<T>(),"Cannot assign 'identity' value for non-arithmetic types");
@@ -218,7 +218,7 @@ public:
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator=(const TinyMatrix& A) noexcept
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -227,24 +227,24 @@ public:
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix& operator=(TinyMatrix&& A) noexcept = default;
 
   template <typename... Args>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix(const T& t, Args&&... args) noexcept
   {
     static_assert(sizeof...(args)==N*N-1, "wrong number of parameters");
     this->_unpackVariadicInput(t, args...);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix() noexcept
   {
     ;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix(const ZeroType& z) noexcept
   {
     static_assert(std::is_arithmetic<T>(),"Cannot construct from 'zero' value for non-arithmetic types");
@@ -253,7 +253,7 @@ public:
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix(const IdentityType& I) noexcept
   {
     static_assert(std::is_arithmetic<T>(),"Cannot construct from 'identity' value for non-arithmetic types");
@@ -264,7 +264,7 @@ public:
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyMatrix(const TinyMatrix& A) noexcept
   {
     for (size_t i=0; i<N*N; ++i) {
@@ -272,15 +272,15 @@ public:
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   TinyMatrix(TinyMatrix&& A) noexcept = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~TinyMatrix()=default;
 };
 
 template <size_t N, typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyMatrix<N,T> tensorProduct(const TinyVector<N,T>& x,
                                         const TinyVector<N,T>& y)
 {
@@ -294,7 +294,7 @@ constexpr TinyMatrix<N,T> tensorProduct(const TinyVector<N,T>& x,
 }
 
 template <size_t N, typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr T det(const TinyMatrix<N,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
@@ -338,7 +338,7 @@ constexpr T det(const TinyMatrix<N,T>& A)
 }
 
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr T det(const TinyMatrix<1,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
@@ -346,7 +346,7 @@ constexpr T det(const TinyMatrix<1,T>& A)
 }
 
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr T det(const TinyMatrix<2,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
@@ -354,7 +354,7 @@ constexpr T det(const TinyMatrix<2,T>& A)
 }
 
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr T det(const TinyMatrix<3,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "determinent is not defined for non arithmetic types");
@@ -365,7 +365,7 @@ constexpr T det(const TinyMatrix<3,T>& A)
 }
 
 template <size_t N, typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyMatrix<N-1,T> getMinor(const TinyMatrix<N,T>& A,
                                      const size_t& I,
                                      const size_t& J)
@@ -393,11 +393,11 @@ constexpr TinyMatrix<N-1,T> getMinor(const TinyMatrix<N,T>& A,
 }
 
 template <size_t N, typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyMatrix<N,T> inverse(const TinyMatrix<N,T>& A);
 
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
@@ -408,7 +408,7 @@ constexpr TinyMatrix<1,T> inverse(const TinyMatrix<1,T>& A)
 }
 
 template <size_t N, typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr T cofactor(const TinyMatrix<N,T>& A,
                      const size_t& i,
                      const size_t& j)
@@ -420,7 +420,7 @@ constexpr T cofactor(const TinyMatrix<N,T>& A,
 }
 
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
@@ -436,7 +436,7 @@ constexpr TinyMatrix<2,T> inverse(const TinyMatrix<2,T>& A)
 
 
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyMatrix<3,T> inverse(const TinyMatrix<3,T>& A)
 {
   static_assert(std::is_arithmetic<T>::value, "inverse is not defined for non arithmetic types");
diff --git a/src/algebra/TinyVector.hpp b/src/algebra/TinyVector.hpp
index c783bb7968696c120f96ef82014a5443639871d9..988d5f0ea45e68c840199c03c626ffc61f119e90 100644
--- a/src/algebra/TinyVector.hpp
+++ b/src/algebra/TinyVector.hpp
@@ -3,10 +3,13 @@
 
 #include <iostream>
 
+#include <PastisMacros.hpp>
 #include <PastisAssert.hpp>
-#include <Kokkos_Core.hpp>
+
 #include <Types.hpp>
 
+#include <cmath>
+
 template <size_t N, typename T=double>
 class TinyVector
 {
@@ -15,7 +18,7 @@ class TinyVector
   static_assert((N>0),"TinyVector size must be strictly positive");
 
   template <typename... Args>
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   constexpr void _unpackVariadicInput(const T& t, Args&&... args) noexcept
   {
     m_values[N-1-sizeof...(args)] = t;
@@ -25,7 +28,7 @@ class TinyVector
   }
 
  public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector operator-() const
   {
     TinyVector opposed;
@@ -35,13 +38,13 @@ class TinyVector
     return std::move(opposed);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr size_t dimension() const
   {
     return N;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr bool operator==(const TinyVector& v) const
   {
     for (size_t i=0; i<N; ++i) {
@@ -50,13 +53,13 @@ class TinyVector
     return true;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr bool operator!=(const TinyVector& v) const
   {
     return not this->operator==(v);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr T operator,(const TinyVector& v) const
   {
     T t = m_values[0]*v.m_values[0];
@@ -66,7 +69,7 @@ class TinyVector
     return std::move(t);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector& operator*=(const T& t)
   {
     for (size_t i=0; i<N; ++i) {
@@ -75,21 +78,21 @@ class TinyVector
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr friend TinyVector operator*(const T& t, const TinyVector& v)
   {
     TinyVector w = v;
     return std::move(w*=t);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr friend TinyVector operator*(const T& t, TinyVector&& v)
   {
     v *= t;
     return std::move(v);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr friend std::ostream& operator<<(std::ostream& os, const TinyVector& v)
   {
     os << '(' << v.m_values[0];
@@ -100,7 +103,7 @@ class TinyVector
     return os;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector operator+(const TinyVector& v) const
   {
     TinyVector sum;
@@ -110,7 +113,7 @@ class TinyVector
     return std::move(sum);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector operator+(TinyVector&& v) const
   {
     for (size_t i=0; i<N; ++i) {
@@ -119,7 +122,7 @@ class TinyVector
     return std::move(v);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector operator-(const TinyVector& v) const
   {
     TinyVector difference;
@@ -129,7 +132,7 @@ class TinyVector
     return std::move(difference);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector operator-(TinyVector&& v) const
   {
     for (size_t i=0; i<N; ++i) {
@@ -138,7 +141,7 @@ class TinyVector
     return std::move(v);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector& operator+=(const TinyVector& v)
   {
     for (size_t i=0; i<N; ++i) {
@@ -147,7 +150,7 @@ class TinyVector
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector& operator-=(const TinyVector& v)
   {
     for (size_t i=0; i<N; ++i) {
@@ -156,21 +159,21 @@ class TinyVector
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr T& operator[](const size_t& i) noexcept(NO_ASSERT)
   {
     Assert(i<N);
     return m_values[i];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr const T& operator[](const size_t& i) const noexcept(NO_ASSERT)
   {
     Assert(i<N);
     return m_values[i];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector& operator=(const ZeroType& z) noexcept
   {
     static_assert(std::is_arithmetic<T>(),"Cannot assign 'zero' value for non-arithmetic types");
@@ -180,7 +183,7 @@ class TinyVector
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const TinyVector& operator=(const TinyVector& v) noexcept
   {
     for (size_t i=0; i<N; ++i) {
@@ -189,24 +192,24 @@ class TinyVector
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector& operator=(TinyVector&& v) noexcept = default;
 
   template <typename... Args>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector(const T& t, Args&&... args) noexcept
   {
     static_assert(sizeof...(args)==N-1, "wrong number of parameters");
     this->_unpackVariadicInput(t, args...);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector() noexcept
   {
     ;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector(const ZeroType& z) noexcept
   {
     static_assert(std::is_arithmetic<T>(),"Cannot construct from 'zero' value for non-arithmetic types");
@@ -215,7 +218,7 @@ class TinyVector
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector(const TinyVector& v) noexcept
   {
     for (size_t i=0; i<N; ++i) {
@@ -223,15 +226,15 @@ class TinyVector
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   constexpr TinyVector(TinyVector&& v) noexcept = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~TinyVector() noexcept = default;
 };
 
 template <size_t N, typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr T l2Norm(const TinyVector<N,T>& x)
 {
   static_assert(std::is_arithmetic<T>(),"Cannot compute L2 norm for non-arithmetic types");
@@ -241,7 +244,7 @@ constexpr T l2Norm(const TinyVector<N,T>& x)
 
 // Cross product is only defined for K^3 vectors
 template <typename T>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 constexpr TinyVector<3,T> crossProduct(const TinyVector<3,T>& u, const TinyVector<3,T>& v)
 {
   TinyVector<3,T> cross_product(u[1]*v[2]-u[2]*v[1],
diff --git a/src/main.cpp b/src/main.cpp
index ba6f6c504ae02f7859d722de46295e192f92fffe..69c9ec388503ef93294a153f5a87e76bdccc8cb6 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,15 +1,7 @@
 #include <iostream>
-#include <Kokkos_Core.hpp>
-#include <RevisionInfo.hpp>
-#include <rang.hpp>
-#include <FPEManager.hpp>
-#include <SignalManager.hpp>
-#include <ConsoleManager.hpp>
+#include <PastisUtils.hpp>
 
-// #include <RawKokkosAcousticSolver.hpp>
-// #include <MeshLessAcousticSolver.hpp>
-// #include <AcousticSolverClass.hpp>
-// #include <AcousticSolverTest.hpp>
+#include <rang.hpp>
 
 #include <Connectivity.hpp>
 
@@ -19,6 +11,8 @@
 
 #include <VTKWriter.hpp>
 
+#include <Timer.hpp>
+
 #include <TinyVector.hpp>
 #include <TinyMatrix.hpp>
 
@@ -28,106 +22,19 @@
 
 #include <GmshReader.hpp>
 
-#include <CLI/CLI.hpp>
 #include <limits>
 #include <map>
 
 int main(int argc, char *argv[])
 {
-  long unsigned number = 10;
-  std::string filename;
-  {
-    CLI::App app{"Pastis help"};
-
-    app.add_option("-n,--number", number, "Number of cells");//->required();
-
-    app.add_option("filename,-f,--filename", filename, "gmsh file");//->required();
-
-    int threads=-1;
-    app.add_option("--threads", threads, "Number of Kokkos threads")->check(CLI::Range(1,std::numeric_limits<decltype(threads)>::max()));
-
-    std::string colorize="auto";
-    app.add_set("--colorize", colorize, {"auto", "yes", "no"}, "Colorize console output", true);
-
-    bool disable_fpe = false;
-    app.add_flag("--no-fpe", disable_fpe, "Do not trap floating point exceptions");
-    bool disable_signals = false;
-    app.add_flag("--no-signal", disable_signals, "Do not catches signals");
-
-    std::string pause_on_error="auto";
-    app.add_set("--pause-on-error", pause_on_error, {"auto", "yes", "no"}, "Pause for debugging on unexpected error", true);
-
-    std::atexit([](){std::cout << rang::style::reset;});
-    try {
-      app.parse(argc, argv);
-    } catch (const CLI::ParseError &e) {
-      return app.exit(e);
-    }
-
-    ConsoleManager::init(colorize);
-    FPEManager::init(not disable_fpe);
-    SignalManager::setPauseForDebug(pause_on_error);
-    SignalManager::init(not disable_signals);
-  }
-
-  std::cout << "Code version: "
-            << rang::style::bold << RevisionInfo::version() << rang::style::reset << '\n';
-
-  std::cout << "-------------------- "
-            << rang::fg::green
-            << "git info"
-            << rang::fg::reset
-            <<" -------------------------"
-            << '\n';
-  std::cout << "tag:  " << rang::fg::reset
-            << rang::style::bold << RevisionInfo::gitTag() << rang::style::reset << '\n';
-  std::cout << "HEAD: " << rang::style::bold << RevisionInfo::gitHead() << rang::style::reset << '\n';
-  std::cout << "hash: " << rang::style::bold << RevisionInfo::gitHash() << rang::style::reset << "  (";
-
-  if (RevisionInfo::gitIsClean()) {
-    std::cout << rang::fgB::green << "clean" << rang::fg::reset;
-  } else {
-    std::cout << rang::fgB::red << "dirty" << rang::fg::reset;
-  }
-  std::cout << ")\n";
-  std::cout << "-------------------------------------------------------\n";
-
-  Kokkos::initialize(argc, argv);
-  Kokkos::DefaultExecutionSpace::print_configuration(std::cout);
+  std::string filename = initialize(argc, argv);
 
   std::map<std::string, double> method_cost_map;
 
-  // { // Basic function based acoustic solver
-  //   Kokkos::Timer timer;
-  //   timer.reset();
-  //   RawKokkos::AcousticSolver(number);
-  //   method_cost_map["RawKokkos"] = timer.seconds();
-  // }
-
-  // { // class for acoustic solver (mesh less)
-  //   Kokkos::Timer timer;
-  //   timer.reset();
-  //   MeshLessAcousticSolver acoustic_solver(number);
-  //   method_cost_map["MeshLessAcousticSolver"] = timer.seconds();
-  // }
-
-  // { // class for acoustic solver
-  //   Kokkos::Timer timer;
-  //   timer.reset();
-  //   AcousticSolverClass acoustic_solver(number);
-  //   method_cost_map["AcousticSolverClass"] = timer.seconds();
-  // }
-
-  // { // class for acoustic solver test
-  //   Kokkos::Timer timer;
-  //   timer.reset();
-  //   AcousticSolverTest acoustic_solver(number);
-  //   method_cost_map["AcousticSolverTest"] = timer.seconds();
-  // }
   try  {
     if (filename != "") {
       std::cout << "Reading (gmsh) " << rang::style::underline << filename << rang::style::reset << " ...\n";
-      Kokkos::Timer gmsh_timer;
+      Timer gmsh_timer;
       gmsh_timer.reset();
       GmshReader gmsh_reader(filename);
       method_cost_map["Mesh building"] = gmsh_timer.seconds();
@@ -154,7 +61,7 @@ int main(int argc, char *argv[])
 
           const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
 
-          Kokkos::Timer timer;
+          Timer timer;
           timer.reset();
           MeshDataType mesh_data(mesh);
 
@@ -263,7 +170,7 @@ int main(int argc, char *argv[])
 
           const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
 
-          Kokkos::Timer timer;
+          Timer timer;
           timer.reset();
           MeshDataType mesh_data(mesh);
 
@@ -359,7 +266,7 @@ int main(int argc, char *argv[])
 
           const MeshType& mesh = dynamic_cast<const MeshType&>(*gmsh_reader.mesh());
 
-          Kokkos::Timer timer;
+          Timer timer;
           timer.reset();
           MeshDataType mesh_data(mesh);
 
@@ -449,7 +356,7 @@ int main(int argc, char *argv[])
     std::exit(1);
   }
 
-  Kokkos::finalize();
+  finalize();
 
   std::string::size_type size=0;
   for (const auto& method_cost : method_cost_map) {
@@ -458,14 +365,14 @@ int main(int argc, char *argv[])
 
   for (const auto& method_cost : method_cost_map) {
     std::cout << "* ["
-      	      << rang::fgB::cyan
-	      << std::setw(size) << std::left
-	      << method_cost.first
-	      << rang::fg::reset
-	      << "] Execution time: "
-	      << rang::style::bold
-	      << method_cost.second
-	      << rang::style::reset << '\n';
+              << rang::fgB::cyan
+              << std::setw(size) << std::left
+              << method_cost.first
+              << rang::fg::reset
+              << "] Execution time: "
+              << rang::style::bold
+              << method_cost.second
+              << rang::style::reset << '\n';
   }
 
   return 0;
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index f84f3835abe00acb25561ada3d071a29c5adc167..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(), KOKKOS_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(), KOKKOS_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 cfcb7f67909e2c7668bb11d842eb5437e57734c2..c4cd00ab7b8a3667379571f69a6281c8538be817 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -1,10 +1,12 @@
 #ifndef CONNECTIVITY_HPP
 #define CONNECTIVITY_HPP
 
+#include <PastisMacros.hpp>
 #include <PastisAssert.hpp>
+#include <PastisUtils.hpp>
+
 #include <TinyVector.hpp>
 
-#include <Kokkos_Core.hpp>
 #include <ItemValue.hpp>
 
 #include <IConnectivity.hpp>
@@ -70,14 +72,14 @@ class ConnectivityFace<2>
     return os;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   bool operator==(const ConnectivityFace& f) const
   {
     return ((m_node0_id == f.m_node0_id) and
             (m_node1_id == f.m_node1_id));
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   bool operator<(const ConnectivityFace& f) const
   {
     return ((m_node0_id<f.m_node0_id) or
@@ -85,13 +87,13 @@ class ConnectivityFace<2>
              (m_node1_id<f.m_node1_id)));
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace& operator=(const ConnectivityFace&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace& operator=(ConnectivityFace&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
   {
     Assert(given_node_id_list.size()==2);
@@ -101,13 +103,13 @@ class ConnectivityFace<2>
     m_node1_id = max;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace(const ConnectivityFace&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace(ConnectivityFace&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~ConnectivityFace() = default;
 };
 
@@ -139,19 +141,19 @@ class ConnectivityFace<3>
     return os;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const bool& reversed() const
   {
     return m_reversed;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const std::vector<unsigned int>& nodeIdList() const
   {
     return m_node_id_list;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
   {
     const auto min_id = std::min_element(node_list.begin(), node_list.end());
@@ -172,7 +174,7 @@ class ConnectivityFace<3>
     return rotated_node_list;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
       : m_reversed(false),
         m_node_id_list(_sort(given_node_id_list))
@@ -194,7 +196,7 @@ class ConnectivityFace<3>
     return false;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   bool operator<(const ConnectivityFace& f) const
   {
     const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
@@ -205,23 +207,23 @@ class ConnectivityFace<3>
     return m_node_id_list.size() < f.m_node_id_list.size();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace& operator=(const ConnectivityFace&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace& operator=(ConnectivityFace&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace(const ConnectivityFace&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace(ConnectivityFace&&) = default;
 
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityFace() = delete;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~ConnectivityFace() = default;
 };
 
@@ -272,7 +274,7 @@ class Connectivity final
   void _computeCellFaceAndFaceNodeConnectivities();
 
   template <typename SubItemValuePerItemType>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const SubItemValuePerItemType&
   _lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const
   {
@@ -287,7 +289,7 @@ class Connectivity final
 
   friend class ConnectivityComputer;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const ConnectivityMatrix& _getMatrix(const ItemType& item_type_0,
                                        const ItemType& item_type_1) const
   {
@@ -304,7 +306,7 @@ class Connectivity final
 
  public:
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const bool& isConnectivityMatrixBuilt(const ItemType& item_type_0,
                                         const ItemType& item_type_1) const
   {
@@ -314,7 +316,7 @@ class Connectivity final
   }
   template <ItemType source_item_type,
             ItemType target_item_type>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto getItemToItemMatrix() const
   {
     return ItemToItemMatrix<source_item_type,
@@ -322,93 +324,93 @@ class Connectivity final
                                                          target_item_type));
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto cellToFaceMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto cellToEdgeMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto cellToNodeMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto faceToCellMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto faceToEdgeMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto faceToNodeMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto edgeToCellMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto edgeToFaceMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto edgeToNodeMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto nodeToCellMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto nodeToFaceMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto nodeToEdgeMatrix() const
   {
     return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>();
   }
 
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& cellFaceIsReversed() const
   {
     static_assert(dimension>1, "reversed faces makes no sense in dimension 1");
     return m_cell_face_is_reversed;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& cellLocalNumbersInTheirNodes() const
   {
     return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& cellLocalNumbersInTheirEdges() const
   {
     if constexpr (dimension>2) {
@@ -418,7 +420,7 @@ class Connectivity final
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& cellLocalNumbersInTheirFaces() const
   {
     if constexpr (dimension>1) {
@@ -428,7 +430,7 @@ class Connectivity final
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& faceLocalNumbersInTheirCells() const
   {
     if constexpr(dimension>1) {
@@ -438,21 +440,21 @@ class Connectivity final
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& faceLocalNumbersInTheirEdges() const
   {
     static_assert(dimension>2,"this function has no meaning in 1d or 2d");
     return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_edges);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& faceLocalNumbersInTheirNodes() const
   {
     static_assert(dimension>1,"this function has no meaning in 1d");
     return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& edgeLocalNumbersInTheirCells() const
   {
     if constexpr (dimension>2) {
@@ -462,34 +464,34 @@ class Connectivity final
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& edgeLocalNumbersInTheirFaces() const
   {
     static_assert(dimension>2, "this function has no meaning in 1d or 2d");
     return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_faces);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& edgeLocalNumbersInTheirNodes() const
   {
     static_assert(dimension>2, "this function has no meaning in 1d or 2d");
     return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_nodes);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& nodeLocalNumbersInTheirCells() const
   {
     return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& nodeLocalNumbersInTheirEdges() const
   {
     static_assert(dimension>2, "this function has no meaning in 1d or 2d");
     return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_edges);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto& nodeLocalNumbersInTheirFaces() const
   {
     static_assert(dimension>1,"this function has no meaning in 1d");
@@ -526,7 +528,7 @@ class Connectivity final
     return m_ref_node_list[i];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfNodes() const final
   {
     const auto& node_to_cell_matrix
@@ -534,7 +536,7 @@ class Connectivity final
     return node_to_cell_matrix.numRows();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfEdges() const final
   {
     const auto& edge_to_node_matrix
@@ -542,7 +544,7 @@ class Connectivity final
     return edge_to_node_matrix.numRows();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfFaces() const final
   {
     const auto& face_to_node_matrix
@@ -550,7 +552,7 @@ class Connectivity final
     return face_to_node_matrix.numRows();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfCells() const final
   {
     const auto& cell_to_node_matrix
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index aa6f5c56ae612133afc3040d84fd6fe3c807b507..8c087e3baf65cec60ae10c5a130e51748dfe19e5 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -5,7 +5,7 @@
 
 
 template <typename ConnectivityType>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 ConnectivityMatrix ConnectivityComputer::
 computeConnectivityMatrix(const ConnectivityType& connectivity,
                           const ItemType& item_type,
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 4387070d0bc9ffc91e88f16c160b87a65e176925..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
@@ -30,31 +30,31 @@ class ConnectivityMatrix
     return m_host_matrix.row_map;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto numEntries() const
   {
     return m_host_matrix.entries.extent(0);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto numRows() const
   {
     return m_host_matrix.numRows();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto rowConst(const size_t& j) const
   {
     return m_host_matrix.rowConst(j);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto rowMap(const size_t& j) const
   {
     return m_host_matrix.row_map[j];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer)
       : m_host_matrix{Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", initializer)},
         m_is_built{true}
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index ea003eb9b6c41d7f049c648ab62b01a01b0add26..f187da556f28f8e4d0661077f6e397d8ce11b8f7 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -1,4 +1,5 @@
 #include <GmshReader.hpp>
+#include <PastisMacros.hpp>
 
 #include <iostream>
 #include <fstream>
@@ -17,7 +18,8 @@
 #include <iomanip>
 
 template<typename T>
-inline std::string stringify(const T & t)
+PASTIS_INLINE
+std::string stringify(const T & t)
 {
   std::ostringstream oss;
   oss << t;
@@ -25,7 +27,8 @@ inline std::string stringify(const T & t)
 }
 
 template<>
-inline std::string stringify<std::string>(const std::string& t)
+PASTIS_INLINE
+std::string stringify<std::string>(const std::string& t)
 {
   return t;
 }
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index e27c2ef08bf13e7e9946c095da1cb628027e3a44..f60db961d80f939c64139f8746a674e4ce4de91e 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -33,28 +33,28 @@ class IConnectivity
 };
 
 template <>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 size_t IConnectivity::numberOf<ItemType::node>() const
 {
   return this->numberOfNodes();
 }
 
 template <>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 size_t IConnectivity::numberOf<ItemType::edge>() const
 {
   return this->numberOfEdges();
 }
 
 template <>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 size_t IConnectivity::numberOf<ItemType::face>() const
 {
   return this->numberOfFaces();
 }
 
 template <>
-KOKKOS_INLINE_FUNCTION
+PASTIS_INLINE
 size_t IConnectivity::numberOf<ItemType::cell>() const
 {
   return this->numberOfCells();
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index 510c421fc9c1a67d0565d4b5e5e3e08c91c4ccd9..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,
@@ -23,38 +23,38 @@ class ItemToItemMatrix
 
    public:
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     size_t size() const
     {
       return m_row.length;
     }
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     TargetItemId operator[](const size_t& j) const
     {
       return m_row(j);
     }
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubItemList(const RowType& row)
         : m_row{row}
     {
       ;
     }
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubItemList& operator=(const SubItemList&) = default;
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubItemList& operator=(SubItemList&&) = default;
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubItemList(const SubItemList&) = default;
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubItemList(SubItemList&&) = default;
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     ~SubItemList() = default;
   };
 
@@ -62,7 +62,7 @@ class ItemToItemMatrix
   const ConnectivityMatrix& m_connectivity_matrix;
 
  public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto operator[](const SourceItemId& source_id) const
   {
     using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
@@ -70,7 +70,7 @@ class ItemToItemMatrix
   }
 
   template <typename IndexType>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const auto operator[](const IndexType& source_id) const
   {
     static_assert(std::is_same<IndexType, SourceItemId>(),
@@ -79,26 +79,26 @@ class ItemToItemMatrix
     return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemToItemMatrix(const ConnectivityMatrix& connectivity_matrix)
       : m_connectivity_matrix{connectivity_matrix}
   {
     ;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemToItemMatrix& operator=(const ItemToItemMatrix&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemToItemMatrix& operator=(ItemToItemMatrix&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemToItemMatrix(ItemToItemMatrix&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemToItemMatrix(const ItemToItemMatrix&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~ItemToItemMatrix() = default;
 };
 
diff --git a/src/mesh/ItemType.hpp b/src/mesh/ItemType.hpp
index 965f888b1b7315f28f5132b0457b86ca1e85b210..1b5f3cb452e898dc957ff7f96fb07c4955312f16 100644
--- a/src/mesh/ItemType.hpp
+++ b/src/mesh/ItemType.hpp
@@ -13,8 +13,8 @@ enum class ItemType
   cell = 3
 };
 
-inline constexpr
-std::string_view itemName(const ItemType& item_type)
+PASTIS_INLINE
+constexpr std::string_view itemName(const ItemType& item_type)
 {
   std::string_view name;
   switch(item_type){
@@ -44,7 +44,8 @@ struct ItemTypeId {};
 template <>
 struct ItemTypeId<1>
 {
-  inline static constexpr size_t itemTId(const ItemType& item_type) {
+  PASTIS_INLINE
+  static constexpr size_t itemTId(const ItemType& item_type) {
     size_t i = std::numeric_limits<size_t>::max();
     switch(item_type) {
       case ItemType::cell: {
@@ -66,7 +67,8 @@ struct ItemTypeId<1>
 template <>
 struct ItemTypeId<2>
 {
-  inline static constexpr size_t itemTId(const ItemType& item_type) {
+  PASTIS_INLINE
+  static constexpr size_t itemTId(const ItemType& item_type) {
     size_t i = std::numeric_limits<size_t>::max();
     switch(item_type) {
       case ItemType::cell: {
@@ -91,7 +93,8 @@ struct ItemTypeId<2>
 template <>
 struct ItemTypeId<3>
 {
-  inline static constexpr size_t itemTId(const ItemType& item_type) {
+  PASTIS_INLINE
+  static constexpr size_t itemTId(const ItemType& item_type) {
     size_t i = std::numeric_limits<size_t>::max();
     switch(item_type) {
       case ItemType::cell: {
diff --git a/src/mesh/ItemValue.hpp b/src/mesh/ItemValue.hpp
index 77eb1cb37d9998fb1af6e4a27d05ff7fcdb60077..754f66ac7e450af35cf5d89e6c22b4f040c4a661 100644
--- a/src/mesh/ItemValue.hpp
+++ b/src/mesh/ItemValue.hpp
@@ -31,13 +31,13 @@ class ItemValue
                    item_type>;
 
  public:
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   const bool& isBuilt() const
   {
     return m_is_built;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t size() const
   {
     Assert(m_is_built);
@@ -46,7 +46,7 @@ class ItemValue
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   DataType& operator[](const ItemId& i) const
   {
     Assert(m_is_built);
@@ -61,7 +61,7 @@ class ItemValue
     return m_values[i];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfItems() const
   {
     Assert(m_is_built);
@@ -69,7 +69,7 @@ class ItemValue
   }
 
   template <typename DataType2>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemValue&
   operator=(const ItemValue<DataType2, item_type>& value_per_item)
   {
@@ -88,7 +88,7 @@ class ItemValue
   }
 
   template <typename DataType2>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ItemValue(const ItemValue<DataType2, item_type>& value_per_item)
   {
     this->operator=(value_per_item);
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 0188d507cb2d463b899e3e318b5a75a002b152e4..c21910e7a5bb24f08ec1106db08793912ee0bc7b 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -27,50 +27,50 @@ private:
   NodeValue<Rd> m_mutable_xr;
 
 public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const size_t meshDimension() const
   {
     return dimension;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const Connectivity& connectivity() const
   {
     return *m_connectivity;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfNodes() const
   {
     return m_connectivity->numberOfNodes();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfFaces() const
   {
     return m_connectivity->numberOfFaces();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfCells() const
   {
     return m_connectivity->numberOfCells();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const NodeValue<const Rd>& xr() const
   {
     return m_xr;
   }
 
   [[deprecated("should rework this class: quite ugly")]]
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   NodeValue<Rd> mutableXr() const
   {
     return m_mutable_xr;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Mesh(const std::shared_ptr<Connectivity>& connectivity,
        NodeValue<Rd>& xr)
       : m_connectivity{connectivity},
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 1a07f1bcc0e477a0141e0ab7d4f7065aad633336..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>
@@ -30,7 +30,7 @@ class MeshData
   CellValue<const Rd> m_xj;
   CellValue<const double> m_Vj;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void _updateCenter()
   { // Computes vertices isobarycenter
     if constexpr (dimension == 1) {
@@ -40,7 +40,7 @@ class MeshData
           = m_mesh.connectivity().cellToNodeMatrix();
 
       CellValue<Rd> xj(m_mesh.connectivity());
-      Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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(), KOKKOS_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) {
@@ -66,7 +66,7 @@ class MeshData
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void _updateVolume()
   {
     const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -74,7 +74,7 @@ class MeshData
         = m_mesh.connectivity().cellToNodeMatrix();
 
     CellValue<double> Vj(m_mesh.connectivity());
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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];
 
@@ -86,7 +86,7 @@ class MeshData
     m_Vj = Vj;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void _updateCjr() {
     if constexpr (dimension == 1) {
       // Cjr/njr/ljr are constant overtime
@@ -98,7 +98,7 @@ class MeshData
 
       {
         NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
-        Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_LAMBDA(const size_t& jr){
+        parallel_for(Cjr.numberOfValues(), PASTIS_LAMBDA(const size_t& jr){
             Cjr[jr] = zero;
           });
 
-        Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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 0a4b9ec63ff74b71e25342144871c34e3236c163..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(), KOKKOS_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(), KOKKOS_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;
@@ -96,16 +96,19 @@ class MeshFlatNodeBoundary
   const Rd m_outgoing_normal;
 
   template <typename MeshType>
-  inline Rd _getNormal(const MeshType& mesh);
+  PASTIS_INLINE
+  Rd _getNormal(const MeshType& mesh);
 
   template <typename MeshType>
-  inline void _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
-                                   const TinyVector<2,double>& xmin,
-                                   const TinyVector<2,double>& xmax,
-                                   const MeshType& mesh) const;
+  PASTIS_INLINE
+  void _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
+                            const TinyVector<2,double>& xmin,
+                            const TinyVector<2,double>& xmax,
+                            const MeshType& mesh) const;
 
   template <typename MeshType>
-  inline Rd _getOutgoingNormal(const MeshType& mesh);
+  PASTIS_INLINE
+  Rd _getOutgoingNormal(const MeshType& mesh);
  public:
   const Rd& outgoingNormal() const
   {
@@ -155,7 +158,7 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
 
   const NodeValue<const R2>& xr = mesh.xr();
 
-  Kokkos::parallel_for(m_node_list.size(), KOKKOS_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";
@@ -166,7 +169,8 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
 
 template <>
 template <typename MeshType>
-inline TinyVector<1,double>
+PASTIS_INLINE
+TinyVector<1,double>
 MeshFlatNodeBoundary<1>::
 _getNormal(const MeshType& mesh)
 {
@@ -183,7 +187,8 @@ _getNormal(const MeshType& mesh)
 
 template <>
 template <typename MeshType>
-inline TinyVector<2,double>
+PASTIS_INLINE
+TinyVector<2,double>
 MeshFlatNodeBoundary<2>::
 _getNormal(const MeshType& mesh)
 {
@@ -225,7 +230,8 @@ _getNormal(const MeshType& mesh)
 
 template <>
 template <typename MeshType>
-inline TinyVector<3,double>
+PASTIS_INLINE
+TinyVector<3,double>
 MeshFlatNodeBoundary<3>::
 _getNormal(const MeshType& mesh)
 {
@@ -309,7 +315,8 @@ _getNormal(const MeshType& mesh)
 
 template <>
 template <typename MeshType>
-inline TinyVector<1,double>
+PASTIS_INLINE
+TinyVector<1,double>
 MeshFlatNodeBoundary<1>::
 _getOutgoingNormal(const MeshType& mesh)
 {
@@ -344,7 +351,8 @@ _getOutgoingNormal(const MeshType& mesh)
 
 template <>
 template <typename MeshType>
-inline TinyVector<2,double>
+PASTIS_INLINE
+TinyVector<2,double>
 MeshFlatNodeBoundary<2>::
 _getOutgoingNormal(const MeshType& mesh)
 {
@@ -379,7 +387,8 @@ _getOutgoingNormal(const MeshType& mesh)
 
 template <>
 template <typename MeshType>
-inline TinyVector<3,double>
+PASTIS_INLINE
+TinyVector<3,double>
 MeshFlatNodeBoundary<3>::
 _getOutgoingNormal(const MeshType& mesh)
 {
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 3853038a1822ccd5bc6e4eb6d5c2d2696c8876ba..afe1ee4fa99277a199cafeee9aceac935dee4400 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -53,25 +53,25 @@ class SubItemValuePerItem<DataType,
     using data_type = DataType;
 
    private:
-    KOKKOS_RESTRICT DataType* const m_sub_values;
+    PASTIS_RESTRICT DataType* const m_sub_values;
     const size_t m_size;
 
    public:
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     const DataType& operator[](const size_t& i) const
     {
       Assert(i<m_size);
       return m_sub_values[i];
     }
 
-    KOKKOS_FORCEINLINE_FUNCTION
+    PASTIS_FORCEINLINE
     DataType& operator[](const size_t& i)
     {
       Assert(i<m_size);
       return m_sub_values[i];
     }
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     const size_t& size() const
     {
       return m_size;
@@ -79,10 +79,10 @@ class SubItemValuePerItem<DataType,
 
     SubView(const SubView&) = delete;
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubView(SubView&&) = default;
 
-    KOKKOS_INLINE_FUNCTION
+    PASTIS_INLINE
     SubView(const Array<DataType>& values,
             const size_t& begin,
             const size_t& end)
@@ -94,7 +94,7 @@ class SubItemValuePerItem<DataType,
     }
   };
 
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   const bool& isBuilt() const
   {
     return m_is_built;
@@ -102,7 +102,7 @@ class SubItemValuePerItem<DataType,
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   DataType& operator()(const ItemId& j, const size_t& r) const
   {
     Assert(m_is_built);
@@ -110,7 +110,7 @@ class SubItemValuePerItem<DataType,
   }
 
   template <typename IndexType>
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   DataType& operator()(const IndexType& j, const size_t& r) const
   {
     static_assert(std::is_same<IndexType, size_t>(),
@@ -118,7 +118,7 @@ class SubItemValuePerItem<DataType,
     return m_values[m_host_row_map(size_t{j})+r];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfValues() const
   {
     Assert(m_is_built);
@@ -127,7 +127,7 @@ class SubItemValuePerItem<DataType,
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
-  KOKKOS_FORCEINLINE_FUNCTION
+  PASTIS_FORCEINLINE
   DataType& operator[](const size_t& i) const
   {
     Assert(m_is_built);
@@ -142,7 +142,7 @@ class SubItemValuePerItem<DataType,
     return m_values[i];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfItems() const
   {
     Assert(m_is_built);
@@ -150,14 +150,14 @@ class SubItemValuePerItem<DataType,
     return m_host_row_map.extent(0);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t numberOfSubValues(const size_t& i_cell) const
   {
     Assert(m_is_built);
     return m_host_row_map(i_cell+1)-m_host_row_map(i_cell);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   SubView itemValues(const size_t& i_cell)
   {
     Assert(m_is_built);
@@ -168,7 +168,7 @@ class SubItemValuePerItem<DataType,
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   SubView itemValues(const size_t& i_cell) const
   {
     Assert(m_is_built);
@@ -178,7 +178,7 @@ class SubItemValuePerItem<DataType,
   }
 
   template <typename DataType2>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   SubItemValuePerItem&
   operator=(const SubItemValuePerItem<DataType2, sub_item_type, item_type>& sub_item_value_per_item)
   {
@@ -198,7 +198,7 @@ class SubItemValuePerItem<DataType,
   }
 
   template <typename DataType2>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   SubItemValuePerItem(const SubItemValuePerItem<DataType2, sub_item_type, item_type>& sub_item_value_per_item)
   {
     this->operator=(sub_item_value_per_item);
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 375828f41a86eef02f1d7fe87e48bb6ff4388223..3afc434969159d4d8a547040be195669e2360f52 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -34,24 +34,24 @@ class AcousticSolver
   using Rdd = TinyMatrix<dimension>;
 
  private:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const CellValue<const double>
   computeRhoCj(const CellValue<const double>& rhoj,
                const CellValue<const double>& cj)
   {
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
         m_rhocj[j] = rhoj[j]*cj[j];
       });
     return m_rhocj;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void computeAjr(const CellValue<const double>& rhocj,
                   const NodeValuePerCell<const Rd>& Cjr,
                   const NodeValuePerCell<const double>& ljr,
                   const NodeValuePerCell<const Rd>& njr)
   {
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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) {
@@ -60,7 +60,7 @@ class AcousticSolver
       });
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const NodeValue<const Rdd>
   computeAr(const NodeValuePerCell<const Rdd>& Ajr) {
     const auto& node_to_cell_matrix
@@ -68,7 +68,7 @@ class AcousticSolver
     const auto& node_local_numbers_in_their_cells
         = m_connectivity.nodeLocalNumbersInTheirCells();
 
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_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
@@ -85,7 +85,7 @@ class AcousticSolver
     return m_Ar;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   const NodeValue<const Rd>
   computeBr(const NodeValuePerCell<Rdd>& Ajr,
             const NodeValuePerCell<const Rd>& Cjr,
@@ -96,7 +96,7 @@ class AcousticSolver
     const auto& node_local_numbers_in_their_cells
         = m_connectivity.nodeLocalNumbersInTheirCells();
 
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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,12 +191,12 @@ class AcousticSolver
   void inverse(const NodeValue<const Rdd>& A,
                NodeValue<Rdd>& inv_A) const
   {
-    Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
+    parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
         inv_A[r] = ::inverse(A[r]);
       });
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void computeExplicitFluxes(const NodeValue<const Rd>& xr,
                              const CellValue<const Rd>& xj,
                              const CellValue<const double>& rhoj,
@@ -250,7 +250,7 @@ class AcousticSolver
     ;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   double acoustic_dt(const CellValue<const double>& Vj,
                      const CellValue<const double>& cj) const
   {
@@ -258,7 +258,7 @@ class AcousticSolver
     const auto& cell_to_node_matrix
         = m_mesh.connectivity().cellToNodeMatrix();
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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 704787756bc97ac60b48699c8215c9444cb87bf2..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, KOKKOS_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, KOKKOS_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, KOKKOS_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 8747cb0c5cff48ca787e4186638cb03902966ca7..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, KOKKOS_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;
@@ -125,16 +125,16 @@ public:
     return *m_boundary_condition;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   BoundaryConditionHandler& operator=(BoundaryConditionHandler&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   BoundaryConditionHandler(const BoundaryConditionHandler&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   BoundaryConditionHandler(BoundaryConditionHandler&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   BoundaryConditionHandler(std::shared_ptr<BoundaryCondition> boundary_condition)
     : m_boundary_condition(boundary_condition)
   {
diff --git a/src/scheme/FiniteVolumesEulerUnknowns.hpp b/src/scheme/FiniteVolumesEulerUnknowns.hpp
index 2eaaf474dbbe520da30cc2c5b0ad5b4737e12538..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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_LAMBDA(const CellId& j){
+    parallel_for(m_mesh.numberOfCells(), PASTIS_LAMBDA(const CellId& j){
         m_uj[j] = zero;
       });
 
-    Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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(), KOKKOS_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 e72efb042c5c9006960a15dfd20c87ed29c29d82..e7273ecee82142014b378a55980d66c0305f016a 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -1,8 +1,10 @@
 #ifndef ARRAY_HPP
 #define ARRAY_HPP
 
+#include <PastisMacros.hpp>
+#include <PastisUtils.hpp>
+
 #include <PastisAssert.hpp>
-#include <Kokkos_Core.hpp>
 
 template <typename DataType>
 class Array
@@ -18,20 +20,20 @@ class Array
   friend Array<std::add_const_t<DataType>>;
 
  public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   size_t size() const
   {
     return m_values.extent(0);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   DataType& operator[](const index_type& i) const
   {
     Assert(i<m_values.extent(0));
     return m_values[i];
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array(const size_t& size)
       : m_values("anonymous", size)
   {
@@ -40,7 +42,7 @@ class Array
   }
 
   template <typename DataType2>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array& operator=(const Array<DataType2>& array)
   {
     // ensures that DataType is the same as source DataType2
@@ -54,29 +56,29 @@ class Array
     return *this;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array& operator=(const Array&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array& operator=(Array&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array() = default;
 
   template <typename DataType2>
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array(const Array<DataType2>& array)
   {
     this->operator=(array);
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array(Array&&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   Array(const Array&) = default;
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~Array() = default;
 };
 
diff --git a/src/utils/ArrayUtils.hpp b/src/utils/ArrayUtils.hpp
index 9a55abfebde94b4ab53adff473c52907ec5385f8..b64452d992b9728dfc2e3a117f790d8893ec155c 100644
--- a/src/utils/ArrayUtils.hpp
+++ b/src/utils/ArrayUtils.hpp
@@ -1,7 +1,8 @@
 #ifndef ARRAY_UTILS_HPP
 #define ARRAY_UTILS_HPP
 
-#include <Kokkos_Core.hpp>
+#include <PastisMacros.hpp>
+#include <PastisUtils.hpp>
 
 template <typename ArrayType>
 class ReduceMin
@@ -12,15 +13,15 @@ class ReduceMin
   using index_type = typename ArrayType::index_type;
 
  public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   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;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void operator()(const index_type& i, data_type& data) const
   {
     if (m_array[i] < data) {
@@ -28,7 +29,7 @@ class ReduceMin
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void join(volatile data_type& dst,
             const volatile data_type& src) const
   {
@@ -37,20 +38,20 @@ class ReduceMin
     }
   }
 
-  KOKKOS_INLINE_FUNCTION void
-  init(data_type& value) const
+  PASTIS_INLINE
+  void init(data_type& value) const
   {
     value = std::numeric_limits<data_type>::max();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ReduceMin(const ArrayType& array)
       : m_array(array)
   {
     ;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~ReduceMin() = default;
 };
 
@@ -64,15 +65,15 @@ class ReduceMax
   using index_type = typename ArrayType::index_type;
 
  public:
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   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;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void operator()(const index_type& i, data_type& data) const
   {
     if (m_array[i] > data) {
@@ -80,7 +81,7 @@ class ReduceMax
     }
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   void join(volatile data_type& dst,
             const volatile data_type& src) const
   {
@@ -89,20 +90,20 @@ class ReduceMax
     }
   }
 
-  KOKKOS_INLINE_FUNCTION void
-  init(data_type& value) const
+  PASTIS_INLINE
+  void init(data_type& value) const
   {
     value = -std::numeric_limits<data_type>::max();
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ReduceMax(const ArrayType& array)
       : m_array(array)
   {
     ;
   }
 
-  KOKKOS_INLINE_FUNCTION
+  PASTIS_INLINE
   ~ReduceMax() = default;
 };
 
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index ad047d0ea29d6f6bff9001ca880d1419f2595d20..3599d7ae878000f1c7f5251a050f8ca54bd986d7 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -8,6 +8,7 @@ add_library(
   BacktraceManager.cpp
   ConsoleManager.cpp
   FPEManager.cpp
+  PastisUtils.cpp
   RevisionInfo.cpp
   SignalManager.cpp)
 
@@ -35,7 +36,7 @@ add_custom_command(TARGET PastisGitRevison
   DEPENDS ${CMAKE_CURRENT_BINARY_DIR}/pastis_git_revision.hpp
   COMMENT ""
   )
-  
+
 add_custom_command(OUTPUT ${CMAKE_CURRENT_BINARY_DIR}/pastis_git_revision.hpp
   PRE_BUILD
   COMMAND ${CMAKE_COMMAND} -E copy_if_different
diff --git a/src/utils/FPEManager.cpp b/src/utils/FPEManager.cpp
index 253ef99c37c34493296a71dfd96d47cf4e2d1c76..d4a4247682afa9f72c771933dcf091e00010d4bc 100644
--- a/src/utils/FPEManager.cpp
+++ b/src/utils/FPEManager.cpp
@@ -1,4 +1,5 @@
 #include <FPEManager.hpp>
+#include <PastisMacros.hpp>
 #include <pastis_config.hpp>
 #include <rang.hpp>
 
@@ -12,7 +13,8 @@
 
 // Public domain polyfill for feenableexcept on OS X
 // http://www-personal.umich.edu/~williams/archive/computation/fe-handling-example.c
-inline int feenableexcept(unsigned int excepts)
+PASTIS_INLINE
+int feenableexcept(unsigned int excepts)
 {
   static fenv_t fenv;
   unsigned int new_excepts = excepts & FE_ALL_EXCEPT;
@@ -31,7 +33,8 @@ inline int feenableexcept(unsigned int excepts)
   return fesetenv(&fenv) ? -1 : old_excepts;
 }
 
-inline int fedisableexcept(unsigned int excepts)
+PASTIS_INLINE
+int fedisableexcept(unsigned int excepts)
 {
   static fenv_t fenv;
   unsigned int new_excepts = excepts & FE_ALL_EXCEPT;
diff --git a/src/utils/PastisMacros.hpp b/src/utils/PastisMacros.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..59d2a579512a67a637f8c47f5745db89efea0a24
--- /dev/null
+++ b/src/utils/PastisMacros.hpp
@@ -0,0 +1,13 @@
+#ifndef PASTIS_MACROS_HPP
+#define PASTIS_MACROS_HPP
+
+#include <Kokkos_Macros.hpp>
+
+#define PASTIS_RESTRICT KOKKOS_RESTRICT
+
+#define PASTIS_INLINE KOKKOS_INLINE_FUNCTION
+#define PASTIS_FORCEINLINE KOKKOS_FORCEINLINE_FUNCTION
+
+#define PASTIS_LAMBDA KOKKOS_LAMBDA
+
+#endif // PASTIS_MACROS_HPP
diff --git a/src/utils/PastisUtils.cpp b/src/utils/PastisUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4149fabe36f0d0ee86c18600e866d30529ad6fbd
--- /dev/null
+++ b/src/utils/PastisUtils.cpp
@@ -0,0 +1,92 @@
+#include <PastisUtils.hpp>
+#include <Kokkos_Core.hpp>
+
+#include <RevisionInfo.hpp>
+#include <rang.hpp>
+
+#include <FPEManager.hpp>
+#include <SignalManager.hpp>
+#include <ConsoleManager.hpp>
+
+#include <CLI/CLI.hpp>
+
+std::string initialize(int& argc, char* argv[])
+{
+  long unsigned number = 10;
+  std::string filename;
+
+  std::cout << "Pastis version: "
+            << rang::style::bold << RevisionInfo::version() << rang::style::reset << '\n';
+
+  std::cout << "-------------------- "
+            << rang::fg::green
+            << "git info"
+            << rang::fg::reset
+            <<" -------------------------"
+            << '\n';
+  std::cout << "tag:  " << rang::fg::reset
+            << rang::style::bold << RevisionInfo::gitTag() << rang::style::reset << '\n';
+  std::cout << "HEAD: " << rang::style::bold << RevisionInfo::gitHead() << rang::style::reset << '\n';
+  std::cout << "hash: " << rang::style::bold << RevisionInfo::gitHash() << rang::style::reset << "  (";
+
+  if (RevisionInfo::gitIsClean()) {
+    std::cout << rang::fgB::green << "clean" << rang::fg::reset;
+  } else {
+    std::cout << rang::fgB::red << "dirty" << rang::fg::reset;
+  }
+  std::cout << ")\n";
+  std::cout << "-------------------------------------------------------\n";
+  {
+    CLI::App app{"Pastis help"};
+
+    app.add_option("-n,--number", number, "Number of cells");//->required();
+
+    app.add_option("filename,-f,--filename", filename, "gmsh file");//->required();
+
+    int threads=-1;
+    app.add_option("--threads", threads, "Number of Kokkos threads")->check(CLI::Range(1,std::numeric_limits<decltype(threads)>::max()));
+
+    std::string colorize="auto";
+    app.add_set("--colorize", colorize, {"auto", "yes", "no"}, "Colorize console output", true);
+
+    bool disable_fpe = false;
+    app.add_flag("--no-fpe", disable_fpe, "Do not trap floating point exceptions");
+    bool disable_signals = false;
+    app.add_flag("--no-signal", disable_signals, "Do not catches signals");
+
+    std::string pause_on_error="auto";
+    app.add_set("--pause-on-error", pause_on_error, {"auto", "yes", "no"}, "Pause for debugging on unexpected error", true);
+
+    std::atexit([](){std::cout << rang::style::reset;});
+    try {
+      app.parse(argc, argv);
+    } catch (const CLI::ParseError &e) {
+      std::exit(app.exit(e));
+    }
+
+    ConsoleManager::init(colorize);
+    FPEManager::init(not disable_fpe);
+    SignalManager::setPauseForDebug(pause_on_error);
+    SignalManager::init(not disable_signals);
+  }
+
+  Kokkos::initialize(argc,argv);
+  std::cout << "-------------------- "
+            << rang::fg::green
+            << "exec info"
+            << rang::fg::reset
+            <<" ------------------------"
+            << '\n';
+
+  std::cout << rang::style::bold;
+  Kokkos::DefaultExecutionSpace::print_configuration(std::cout);
+  std::cout << rang::style::reset;
+  std::cout << "-------------------------------------------------------\n";
+
+  return filename;
+}
+
+void finalize()
+{
+  Kokkos::finalize();
+}
diff --git a/src/utils/PastisUtils.hpp b/src/utils/PastisUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d0a68cfdb91af7acbb5393cc00fdc0c2611a4848
--- /dev/null
+++ b/src/utils/PastisUtils.hpp
@@ -0,0 +1,33 @@
+#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);
+}
+
+std::string initialize(int& argc, char* argv[]);
+
+void finalize();
+
+#endif // PASTIS_UTILS_HPP
diff --git a/src/utils/Timer.hpp b/src/utils/Timer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..44d4a0c51868a9c7d877c6e3be74082d118918fb
--- /dev/null
+++ b/src/utils/Timer.hpp
@@ -0,0 +1,8 @@
+#ifndef TIMER_HPP
+#define TIMER_HPP
+
+#include <Kokkos_Timer.hpp>
+
+using Timer = Kokkos::Timer;
+
+#endif // TIMER_HPP
diff --git a/tests/test_ItemType.cpp b/tests/test_ItemType.cpp
index 096dbdaf5b8a85099c8ec892a15c7f6ceb1000ca..600b0882ed26b1f506484a913e1168b26fa173fd 100644
--- a/tests/test_ItemType.cpp
+++ b/tests/test_ItemType.cpp
@@ -1,4 +1,5 @@
 #include <catch.hpp>
+#include <PastisMacros.hpp>
 
 #include <ItemType.hpp>