diff --git a/benchmarks/deepcopy.cpp b/benchmarks/deepcopy.cpp
index 2edafe00b..cb6c67668 100644
--- a/benchmarks/deepcopy.cpp
+++ b/benchmarks/deepcopy.cpp
@@ -133,4 +133,16 @@ BENCHMARK(memcpy_2d)->Args({large_dim1_2D, large_dim2_2D});
 BENCHMARK(deepcopy_2d)->Args({large_dim1_2D, large_dim2_2D});
 BENCHMARK(deepcopy_subchunk_2d)->Args({large_dim1_2D, large_dim2_2D});
 
-BENCHMARK_MAIN();
+int main(int argc, char** argv)
+{
+    ::benchmark::Initialize(&argc, argv);
+    if (::benchmark::ReportUnrecognizedArguments(argc, argv)) {
+        return 1;
+    }
+    {
+        ddc::ScopeGuard const guard;
+        ::benchmark::RunSpecifiedBenchmarks();
+    }
+    ::benchmark::Shutdown();
+    return 0;
+}
diff --git a/include/ddc/chunk_common.hpp b/include/ddc/chunk_common.hpp
index afedaf428..6b91ddf22 100644
--- a/include/ddc/chunk_common.hpp
+++ b/include/ddc/chunk_common.hpp
@@ -132,43 +132,43 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     mdomain_type m_domain;
 
 public:
-    static constexpr int rank() noexcept
+    static KOKKOS_FUNCTION constexpr int rank() noexcept
     {
         return extents_type::rank();
     }
 
-    static constexpr int rank_dynamic() noexcept
+    static KOKKOS_FUNCTION constexpr int rank_dynamic() noexcept
     {
         return extents_type::rank_dynamic();
     }
 
-    static constexpr size_type static_extent(std::size_t r) noexcept
+    static KOKKOS_FUNCTION constexpr size_type static_extent(std::size_t r) noexcept
     {
         return extents_type::static_extent(r);
     }
 
-    static constexpr bool is_always_unique() noexcept
+    static KOKKOS_FUNCTION constexpr bool is_always_unique() noexcept
     {
         return mapping_type::is_always_unique();
     }
 
-    static constexpr bool is_always_exhaustive() noexcept
+    static KOKKOS_FUNCTION constexpr bool is_always_exhaustive() noexcept
     {
         return mapping_type::is_always_exhaustive();
     }
 
-    static constexpr bool is_always_strided() noexcept
+    static KOKKOS_FUNCTION constexpr bool is_always_strided() noexcept
     {
         return mapping_type::is_always_strided();
     }
 
 public:
-    constexpr accessor_type accessor() const
+    KOKKOS_FUNCTION constexpr accessor_type accessor() const
     {
         return m_internal_mdspan.accessor();
     }
 
-    constexpr DiscreteVector<DDims...> extents() const noexcept
+    KOKKOS_FUNCTION constexpr DiscreteVector<DDims...> extents() const noexcept
     {
         return DiscreteVector<DDims...>(
                 (m_internal_mdspan.extent(type_seq_rank_v<DDims, detail::TypeSeq<DDims...>>)
@@ -176,39 +176,39 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     }
 
     template <class QueryDDim>
-    constexpr size_type extent() const noexcept
+    KOKKOS_FUNCTION constexpr size_type extent() const noexcept
     {
         return m_internal_mdspan.extent(type_seq_rank_v<QueryDDim, detail::TypeSeq<DDims...>>)
                - front<QueryDDim>(m_domain).uid();
     }
 
-    constexpr size_type size() const noexcept
+    KOKKOS_FUNCTION constexpr size_type size() const noexcept
     {
         return allocation_mdspan().size();
     }
 
-    constexpr mapping_type mapping() const noexcept
+    KOKKOS_FUNCTION constexpr mapping_type mapping() const noexcept
     {
         return allocation_mdspan().mapping();
     }
 
-    constexpr bool is_unique() const noexcept
+    KOKKOS_FUNCTION constexpr bool is_unique() const noexcept
     {
         return allocation_mdspan().is_unique();
     }
 
-    constexpr bool is_exhaustive() const noexcept
+    KOKKOS_FUNCTION constexpr bool is_exhaustive() const noexcept
     {
         return allocation_mdspan().is_exhaustive();
     }
 
-    constexpr bool is_strided() const noexcept
+    KOKKOS_FUNCTION constexpr bool is_strided() const noexcept
     {
         return allocation_mdspan().is_strided();
     }
 
     template <class QueryDDim>
-    constexpr size_type stride() const
+    KOKKOS_FUNCTION constexpr size_type stride() const
     {
         return m_internal_mdspan.stride(type_seq_rank_v<QueryDDim, detail::TypeSeq<DDims...>>);
     }
@@ -216,7 +216,7 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     /** Provide access to the domain on which this chunk is defined
      * @return the domain on which this chunk is defined
      */
-    constexpr mdomain_type domain() const noexcept
+    KOKKOS_FUNCTION constexpr mdomain_type domain() const noexcept
     {
         return m_domain;
     }
@@ -225,20 +225,22 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
      * @return the domain on which this chunk is defined
      */
     template <class... QueryDDims>
-    constexpr DiscreteDomain<QueryDDims...> domain() const noexcept
+    KOKKOS_FUNCTION constexpr DiscreteDomain<QueryDDims...> domain() const noexcept
     {
         return select<QueryDDims...>(domain());
     }
 
 protected:
     /// Empty ChunkCommon
-    constexpr ChunkCommon() = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkCommon() = default;
 
     /** Constructs a new ChunkCommon from scratch
      * @param internal_mdspan
      * @param domain
      */
-    constexpr ChunkCommon(internal_mdspan_type internal_mdspan, mdomain_type const& domain) noexcept
+    KOKKOS_FUNCTION constexpr ChunkCommon(
+            internal_mdspan_type internal_mdspan,
+            mdomain_type const& domain) noexcept
         : m_internal_mdspan(std::move(internal_mdspan))
         , m_domain(domain)
     {
@@ -251,7 +253,7 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     template <
             class Mapping = mapping_type,
             std::enable_if_t<std::is_constructible_v<Mapping, extents_type>, int> = 0>
-    constexpr ChunkCommon(ElementType* ptr, mdomain_type const& domain)
+    KOKKOS_FUNCTION constexpr ChunkCommon(ElementType* ptr, mdomain_type const& domain)
         : ChunkCommon {ptr, domain, make_mapping_for(domain)}
     {
     }
@@ -259,29 +261,31 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     /** Constructs a new ChunkCommon by copy, yields a new view to the same data
      * @param other the ChunkCommon to copy
      */
-    constexpr ChunkCommon(ChunkCommon const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkCommon(ChunkCommon const& other) = default;
 
     /** Constructs a new ChunkCommon by move
      * @param other the ChunkCommon to move
      */
-    constexpr ChunkCommon(ChunkCommon&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkCommon(ChunkCommon&& other) = default;
+
+    KOKKOS_DEFAULTED_FUNCTION ~ChunkCommon() = default;
 
     /** Copy-assigns a new value to this ChunkCommon, yields a new view to the same data
      * @param other the ChunkCommon to copy
      * @return *this
      */
-    constexpr ChunkCommon& operator=(ChunkCommon const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkCommon& operator=(ChunkCommon const& other) = default;
 
     /** Move-assigns a new value to this ChunkCommon
      * @param other the ChunkCommon to move
      * @return *this
      */
-    constexpr ChunkCommon& operator=(ChunkCommon&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkCommon& operator=(ChunkCommon&& other) = default;
 
     /** Access to the underlying allocation pointer
      * @return allocation pointer
      */
-    constexpr ElementType* data_handle() const
+    KOKKOS_FUNCTION constexpr ElementType* data_handle() const
     {
         return &m_internal_mdspan(front<DDims>(m_domain).uid()...);
     }
@@ -289,7 +293,7 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     /** Provide a modifiable view of the data
      * @return a modifiable view of the data
      */
-    constexpr internal_mdspan_type internal_mdspan() const
+    KOKKOS_FUNCTION constexpr internal_mdspan_type internal_mdspan() const
     {
         return m_internal_mdspan;
     }
@@ -297,7 +301,7 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
     /** Provide a modifiable view of the data
      * @return a modifiable view of the data
      */
-    constexpr allocation_mdspan_type allocation_mdspan() const
+    KOKKOS_FUNCTION constexpr allocation_mdspan_type allocation_mdspan() const
     {
         DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
         extents_type extents_s(::ddc::extents<DDims>(m_domain).value()...);
@@ -316,7 +320,7 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
      * @param domain the domain that sustains the view
      */
     template <class Mapping = mapping_type>
-    constexpr std::enable_if_t<
+    KOKKOS_FUNCTION constexpr std::enable_if_t<
             std::is_constructible_v<Mapping, extents_type>,
             std::experimental::layout_stride::mapping<extents_type>>
     make_mapping_for(mdomain_type const& domain)
@@ -334,7 +338,7 @@ class ChunkCommon<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy>
      * @param ptr the allocation pointer to the data_handle()
      * @param domain the domain that sustains the view
      */
-    constexpr ChunkCommon(
+    KOKKOS_FUNCTION constexpr ChunkCommon(
             ElementType* ptr,
             mdomain_type const& domain,
             std::experimental::layout_stride::mapping<extents_type>&& s_domain)
diff --git a/include/ddc/chunk_span.hpp b/include/ddc/chunk_span.hpp
index 6c4188798..6f10c1794 100644
--- a/include/ddc/chunk_span.hpp
+++ b/include/ddc/chunk_span.hpp
@@ -89,7 +89,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
 
 protected:
     template <class QueryDDim, class... ODDims>
-    constexpr auto get_slicer_for(DiscreteElement<ODDims...> const& c) const
+    KOKKOS_FUNCTION constexpr auto get_slicer_for(DiscreteElement<ODDims...> const& c) const
     {
         DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
         if constexpr (in_tags_v<QueryDDim, detail::TypeSeq<ODDims...>>) {
@@ -101,7 +101,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
     }
 
     template <class QueryDDim, class... ODDims>
-    constexpr auto get_slicer_for(DiscreteDomain<ODDims...> const& c) const
+    KOKKOS_FUNCTION constexpr auto get_slicer_for(DiscreteDomain<ODDims...> const& c) const
     {
         DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
         if constexpr (in_tags_v<QueryDDim, detail::TypeSeq<ODDims...>>) {
@@ -116,23 +116,24 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
 
 public:
     /// Empty ChunkSpan
-    constexpr ChunkSpan() = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkSpan() = default;
 
     /** Constructs a new ChunkSpan by copy, yields a new view to the same data
      * @param other the ChunkSpan to copy
      */
-    constexpr ChunkSpan(ChunkSpan const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkSpan(ChunkSpan const& other) = default;
 
     /** Constructs a new ChunkSpan by move
      * @param other the ChunkSpan to move
      */
-    constexpr ChunkSpan(ChunkSpan&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkSpan(ChunkSpan&& other) = default;
 
     /** Constructs a new ChunkSpan from a Chunk, yields a new view to the same data
      * @param other the Chunk to view
      */
     template <class OElementType, class Allocator>
-    constexpr ChunkSpan(Chunk<OElementType, mdomain_type, Allocator>& other) noexcept
+    KOKKOS_FUNCTION constexpr ChunkSpan(
+            Chunk<OElementType, mdomain_type, Allocator>& other) noexcept
         : base_type(other.m_internal_mdspan, other.m_domain)
     {
     }
@@ -146,7 +147,8 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
             class SFINAEElementType = ElementType,
             class = std::enable_if_t<std::is_const_v<SFINAEElementType>>,
             class Allocator>
-    constexpr ChunkSpan(Chunk<OElementType, mdomain_type, Allocator> const& other) noexcept
+    KOKKOS_FUNCTION constexpr ChunkSpan(
+            Chunk<OElementType, mdomain_type, Allocator> const& other) noexcept
         : base_type(other.m_internal_mdspan, other.m_domain)
     {
     }
@@ -155,7 +157,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
      * @param other the ChunkSpan to move
      */
     template <class OElementType>
-    constexpr ChunkSpan(
+    KOKKOS_FUNCTION constexpr ChunkSpan(
             ChunkSpan<OElementType, mdomain_type, layout_type, MemorySpace> const& other) noexcept
         : base_type(other.m_internal_mdspan, other.m_domain)
     {
@@ -168,7 +170,8 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
     template <
             class Mapping = mapping_type,
             std::enable_if_t<std::is_constructible_v<Mapping, extents_type>, int> = 0>
-    constexpr ChunkSpan(ElementType* const ptr, mdomain_type const& domain) : base_type(ptr, domain)
+    KOKKOS_FUNCTION constexpr ChunkSpan(ElementType* const ptr, mdomain_type const& domain)
+        : base_type(ptr, domain)
     {
     }
 
@@ -196,29 +199,32 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
      * @param domain the domain that sustains the view
      */
     template <class KokkosView, class = std::enable_if_t<Kokkos::is_view<KokkosView>::value>>
-    constexpr ChunkSpan(KokkosView const& view, mdomain_type const& domain) noexcept
+    KOKKOS_FUNCTION constexpr ChunkSpan(KokkosView const& view, mdomain_type const& domain) noexcept
         : ChunkSpan(
                 detail::build_mdspan(view, std::make_index_sequence<sizeof...(DDims)> {}),
                 domain)
     {
     }
 
+    KOKKOS_DEFAULTED_FUNCTION ~ChunkSpan() = default;
+
     /** Copy-assigns a new value to this ChunkSpan, yields a new view to the same data
      * @param other the ChunkSpan to copy
      * @return *this
      */
-    constexpr ChunkSpan& operator=(ChunkSpan const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkSpan& operator=(ChunkSpan const& other) = default;
 
     /** Move-assigns a new value to this ChunkSpan
      * @param other the ChunkSpan to move
      * @return *this
      */
-    constexpr ChunkSpan& operator=(ChunkSpan&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr ChunkSpan& operator=(ChunkSpan&& other) = default;
 
     /** Slice out some dimensions
      */
     template <class... QueryDDims>
-    constexpr auto operator[](DiscreteElement<QueryDDims...> const& slice_spec) const
+    KOKKOS_FUNCTION constexpr auto operator[](
+            DiscreteElement<QueryDDims...> const& slice_spec) const
     {
         auto subview = std::experimental::
                 submdspan(allocation_mdspan(), get_slicer_for<DDims>(slice_spec)...);
@@ -234,7 +240,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
     /** Restrict to a subdomain
      */
     template <class... QueryDDims>
-    constexpr auto operator[](DiscreteDomain<QueryDDims...> const& odomain) const
+    KOKKOS_FUNCTION constexpr auto operator[](DiscreteDomain<QueryDDims...> const& odomain) const
     {
         auto subview = std::experimental::
                 submdspan(allocation_mdspan(), get_slicer_for<DDims>(odomain)...);
@@ -259,7 +265,8 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
        * @return const-reference to this element
        */
     template <class... ODDims>
-    constexpr reference operator()(DiscreteElement<ODDims...> const& delems) const noexcept
+    KOKKOS_FUNCTION constexpr reference operator()(
+            DiscreteElement<ODDims...> const& delems) const noexcept
     {
         static_assert(sizeof...(ODDims) == sizeof...(DDims), "Invalid number of dimensions");
         assert(((select<ODDims>(delems) >= front<ODDims>(this->m_domain)) && ...));
@@ -283,7 +290,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
     /** Access to the underlying allocation pointer
      * @return allocation pointer
      */
-    constexpr ElementType* data_handle() const
+    KOKKOS_FUNCTION constexpr ElementType* data_handle() const
     {
         return base_type::data_handle();
     }
@@ -291,7 +298,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
     /** Provide a mdspan on the memory allocation
      * @return allocation mdspan
      */
-    constexpr allocation_mdspan_type allocation_mdspan() const
+    KOKKOS_FUNCTION constexpr allocation_mdspan_type allocation_mdspan() const
     {
         return base_type::allocation_mdspan();
     }
@@ -299,7 +306,7 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
     /** Provide a mdspan on the memory allocation
      * @return allocation mdspan
      */
-    constexpr auto allocation_kokkos_view() const
+    KOKKOS_FUNCTION constexpr auto allocation_kokkos_view() const
     {
         auto s = this->allocation_mdspan();
         auto kokkos_layout = detail::build_kokkos_layout(
@@ -312,12 +319,12 @@ class ChunkSpan<ElementType, DiscreteDomain<DDims...>, LayoutStridedPolicy, Memo
                 MemorySpace>(s.data_handle(), kokkos_layout);
     }
 
-    constexpr view_type span_cview() const
+    KOKKOS_FUNCTION constexpr view_type span_cview() const
     {
         return view_type(*this);
     }
 
-    constexpr span_type span_view() const
+    KOKKOS_FUNCTION constexpr span_type span_view() const
     {
         return *this;
     }
diff --git a/include/ddc/coordinate_md.hpp b/include/ddc/coordinate_md.hpp
index 05bbdf672..658910965 100644
--- a/include/ddc/coordinate_md.hpp
+++ b/include/ddc/coordinate_md.hpp
@@ -9,7 +9,7 @@
 namespace ddc {
 
 template <class... DDim, std::enable_if_t<(sizeof...(DDim) > 1), int> = 0>
-DDC_INLINE_FUNCTION Coordinate<typename DDim::continuous_dimension_type...> coordinate(
+KOKKOS_FUNCTION Coordinate<typename DDim::continuous_dimension_type...> coordinate(
         DiscreteElement<DDim...> const& c)
 {
     return Coordinate<typename DDim::continuous_dimension_type...>(coordinate(select<DDim>(c))...);
diff --git a/include/ddc/detail/kokkos.hpp b/include/ddc/detail/kokkos.hpp
index 57fdbbdb2..6de7d2771 100644
--- a/include/ddc/detail/kokkos.hpp
+++ b/include/ddc/detail/kokkos.hpp
@@ -105,7 +105,7 @@ using kokkos_to_mdspan_element_t = typename kokkos_to_mdspan_element<ET>::type;
 
 
 template <std::size_t... Is>
-Kokkos::LayoutStride make_layout_stride(
+KOKKOS_FUNCTION Kokkos::LayoutStride make_layout_stride(
         std::array<std::size_t, sizeof...(Is)> const& interleaved_extents_strides,
         std::index_sequence<Is...>)
 {
@@ -113,7 +113,7 @@ Kokkos::LayoutStride make_layout_stride(
 }
 
 template <class EP, class MP, std::size_t... Is>
-mdspan_to_kokkos_layout_t<typename MP::layout_type> build_kokkos_layout(
+KOKKOS_FUNCTION mdspan_to_kokkos_layout_t<typename MP::layout_type> build_kokkos_layout(
         EP const& ep,
         MP const& mapping,
         std::index_sequence<Is...>)
@@ -138,7 +138,9 @@ mdspan_to_kokkos_layout_t<typename MP::layout_type> build_kokkos_layout(
 }
 
 template <class DataType, class... Properties, std::size_t... Is>
-auto build_mdspan(Kokkos::View<DataType, Properties...> const view, std::index_sequence<Is...>)
+KOKKOS_FUNCTION auto build_mdspan(
+        Kokkos::View<DataType, Properties...> const view,
+        std::index_sequence<Is...>)
 {
     DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
     using element_type = kokkos_to_mdspan_element_t<DataType>;
diff --git a/include/ddc/detail/tagged_vector.hpp b/include/ddc/detail/tagged_vector.hpp
index dace14d5b..705309015 100644
--- a/include/ddc/detail/tagged_vector.hpp
+++ b/include/ddc/detail/tagged_vector.hpp
@@ -35,20 +35,21 @@ inline constexpr bool is_tagged_vector_v = IsTaggedVector<T>::value;
 
 
 template <class QueryTag, class ElementType, class... Tags>
-inline constexpr ElementType const& get(
+KOKKOS_FUNCTION constexpr ElementType const& get(
         detail::TaggedVector<ElementType, Tags...> const& tuple) noexcept
 {
     return tuple.template get<QueryTag>();
 }
 
 template <class QueryTag, class ElementType, class... Tags>
-inline constexpr ElementType& get(detail::TaggedVector<ElementType, Tags...>& tuple) noexcept
+KOKKOS_FUNCTION constexpr ElementType& get(
+        detail::TaggedVector<ElementType, Tags...>& tuple) noexcept
 {
     return tuple.template get<QueryTag>();
 }
 
 template <class QueryTag, class ElementType, class... Tags>
-inline constexpr ElementType const& get_or(
+KOKKOS_FUNCTION constexpr ElementType const& get_or(
         detail::TaggedVector<ElementType, Tags...> const& tuple,
         ElementType const& default_value) noexcept
 {
@@ -60,14 +61,14 @@ namespace detail {
 /// Unary operators: +, -
 
 template <class ElementType, class... Tags>
-constexpr inline detail::TaggedVector<ElementType, Tags...> operator+(
+KOKKOS_FUNCTION constexpr detail::TaggedVector<ElementType, Tags...> operator+(
         detail::TaggedVector<ElementType, Tags...> const& x)
 {
     return x;
 }
 
 template <class ElementType, class... Tags>
-constexpr inline detail::TaggedVector<ElementType, Tags...> operator-(
+KOKKOS_FUNCTION constexpr detail::TaggedVector<ElementType, Tags...> operator-(
         detail::TaggedVector<ElementType, Tags...> const& x)
 {
     return detail::TaggedVector<ElementType, Tags...>((-get<Tags>(x))...);
@@ -76,7 +77,7 @@ constexpr inline detail::TaggedVector<ElementType, Tags...> operator-(
 /// Internal binary operators: +, -
 
 template <class ElementType, class... Tags, class OElementType, class... OTags>
-constexpr inline auto operator+(
+KOKKOS_FUNCTION constexpr auto operator+(
         detail::TaggedVector<ElementType, Tags...> const& lhs,
         detail::TaggedVector<OElementType, OTags...> const& rhs)
 {
@@ -91,7 +92,7 @@ template <
         class OElementType,
         class = std::enable_if_t<!detail::is_tagged_vector_v<OElementType>>,
         class = std::enable_if_t<std::is_convertible_v<OElementType, ElementType>>>
-constexpr inline auto operator+(
+KOKKOS_FUNCTION constexpr auto operator+(
         detail::TaggedVector<ElementType, Tag> const& lhs,
         OElementType const& rhs)
 {
@@ -105,7 +106,7 @@ template <
         class OElementType,
         class = std::enable_if_t<!detail::is_tagged_vector_v<OElementType>>,
         class = std::enable_if_t<std::is_convertible_v<ElementType, OElementType>>>
-constexpr inline auto operator+(
+KOKKOS_FUNCTION constexpr auto operator+(
         OElementType const& lhs,
         detail::TaggedVector<ElementType, Tag> const& rhs)
 {
@@ -114,7 +115,7 @@ constexpr inline auto operator+(
 }
 
 template <class ElementType, class... Tags, class OElementType, class... OTags>
-constexpr inline auto operator-(
+KOKKOS_FUNCTION constexpr auto operator-(
         detail::TaggedVector<ElementType, Tags...> const& lhs,
         detail::TaggedVector<OElementType, OTags...> const& rhs)
 {
@@ -129,7 +130,7 @@ template <
         class OElementType,
         class = std::enable_if_t<!detail::is_tagged_vector_v<OElementType>>,
         class = std::enable_if_t<std::is_convertible_v<OElementType, ElementType>>>
-constexpr inline auto operator-(
+KOKKOS_FUNCTION constexpr auto operator-(
         detail::TaggedVector<ElementType, Tag> const& lhs,
         OElementType const& rhs)
 {
@@ -143,7 +144,7 @@ template <
         class OElementType,
         class = std::enable_if_t<!detail::is_tagged_vector_v<OElementType>>,
         class = std::enable_if_t<std::is_convertible_v<ElementType, OElementType>>>
-constexpr inline auto operator-(
+KOKKOS_FUNCTION constexpr auto operator-(
         OElementType const& lhs,
         detail::TaggedVector<ElementType, Tag> const& rhs)
 {
@@ -159,7 +160,7 @@ template <
         class... Tags,
         class = std::enable_if_t<!detail::is_tagged_vector_v<OElementType>>,
         class = std::enable_if_t<std::is_convertible_v<ElementType, OElementType>>>
-constexpr inline auto operator*(
+KOKKOS_FUNCTION constexpr auto operator*(
         ElementType const& lhs,
         detail::TaggedVector<OElementType, Tags...> const& rhs)
 {
@@ -170,14 +171,14 @@ constexpr inline auto operator*(
 } // namespace detail
 
 template <class... QueryTags, class ElementType, class... Tags>
-inline constexpr detail::TaggedVector<ElementType, QueryTags...> select(
+KOKKOS_FUNCTION constexpr detail::TaggedVector<ElementType, QueryTags...> select(
         detail::TaggedVector<ElementType, Tags...> const& arr) noexcept
 {
     return detail::TaggedVector<ElementType, QueryTags...>(arr);
 }
 
 template <class... QueryTags, class ElementType, class... Tags>
-inline constexpr detail::TaggedVector<ElementType, QueryTags...> select(
+KOKKOS_FUNCTION constexpr detail::TaggedVector<ElementType, QueryTags...> select(
         detail::TaggedVector<ElementType, Tags...>&& arr) noexcept
 {
     return detail::TaggedVector<ElementType, QueryTags...>(std::move(arr));
@@ -186,7 +187,7 @@ inline constexpr detail::TaggedVector<ElementType, QueryTags...> select(
 namespace detail {
 
 template <class QueryTag, class ElementType, class HeadTag, class... TailTags>
-constexpr detail::TaggedVector<ElementType, QueryTag> const& take_impl(
+KOKKOS_FUNCTION constexpr detail::TaggedVector<ElementType, QueryTag> const& take_impl(
         detail::TaggedVector<ElementType, HeadTag> const& head,
         detail::TaggedVector<ElementType, TailTags> const&... tags)
 {
@@ -204,7 +205,7 @@ constexpr detail::TaggedVector<ElementType, QueryTag> const& take_impl(
 }
 
 template <class QueryTag, class ElementType, class... Tags>
-constexpr detail::TaggedVector<ElementType, QueryTag> const& take(
+KOKKOS_FUNCTION constexpr detail::TaggedVector<ElementType, QueryTag> const& take(
         detail::TaggedVector<ElementType, Tags> const&... tags)
 {
     return take_impl<QueryTag>(tags...);
@@ -220,12 +221,12 @@ template <class ElementType, class Tag>
 class ConversionOperators<TaggedVector<ElementType, Tag>>
 {
 public:
-    constexpr inline operator ElementType const &() const noexcept
+    KOKKOS_FUNCTION constexpr operator ElementType const &() const noexcept
     {
         return static_cast<TaggedVector<ElementType, Tag> const*>(this)->m_values[0];
     }
 
-    constexpr inline operator ElementType&() noexcept
+    KOKKOS_FUNCTION constexpr operator ElementType&() noexcept
     {
         return static_cast<TaggedVector<ElementType, Tag>*>(this)->m_values[0];
     }
@@ -242,26 +243,29 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     std::array<ElementType, sizeof...(Tags)> m_values;
 
 public:
-    static constexpr std::size_t size() noexcept
+    using value_type = ElementType;
+
+    static KOKKOS_FUNCTION constexpr std::size_t size() noexcept
     {
         return sizeof...(Tags);
     }
 
 public:
-    inline constexpr TaggedVector() = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr TaggedVector() = default;
 
-    inline constexpr TaggedVector(TaggedVector const&) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr TaggedVector(TaggedVector const&) = default;
 
-    inline constexpr TaggedVector(TaggedVector&&) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr TaggedVector(TaggedVector&&) = default;
 
     template <class... OTags>
-    inline constexpr TaggedVector(TaggedVector<ElementType, OTags> const&... other) noexcept
+    KOKKOS_FUNCTION constexpr TaggedVector(
+            TaggedVector<ElementType, OTags> const&... other) noexcept
         : m_values {take<Tags>(other...).value()...}
     {
     }
 
     template <class OElementType, class... OTags>
-    explicit inline constexpr TaggedVector(
+    explicit KOKKOS_FUNCTION constexpr TaggedVector(
             TaggedVector<OElementType, OTags...> const& other) noexcept
         : m_values {(static_cast<ElementType>(other.template get<Tags>()))...}
     {
@@ -272,17 +276,19 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
             class = std::enable_if_t<(std::is_convertible_v<Params, ElementType> && ...)>,
             class = std::enable_if_t<(!is_tagged_vector_v<Params> && ...)>,
             class = std::enable_if_t<sizeof...(Params) == sizeof...(Tags)>>
-    explicit inline constexpr TaggedVector(Params const&... params) noexcept
+    explicit KOKKOS_FUNCTION constexpr TaggedVector(Params const&... params) noexcept
         : m_values {static_cast<ElementType>(params)...}
     {
     }
 
-    constexpr inline TaggedVector& operator=(TaggedVector const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION ~TaggedVector() = default;
+
+    KOKKOS_DEFAULTED_FUNCTION TaggedVector& operator=(TaggedVector const& other) = default;
 
-    constexpr inline TaggedVector& operator=(TaggedVector&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION TaggedVector& operator=(TaggedVector&& other) = default;
 
     template <class... OTags>
-    constexpr inline TaggedVector& operator=(
+    KOKKOS_FUNCTION constexpr TaggedVector& operator=(
             TaggedVector<ElementType, OTags...> const& other) noexcept
     {
         m_values = other.m_values;
@@ -290,48 +296,51 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     }
 
     template <class... OTags>
-    constexpr inline TaggedVector& operator=(TaggedVector<ElementType, OTags...>&& other) noexcept
+    KOKKOS_FUNCTION constexpr TaggedVector& operator=(
+            TaggedVector<ElementType, OTags...>&& other) noexcept
     {
         m_values = std::move(other.m_values);
         return *this;
     }
 
     /// Returns a reference to the underlying `std::array`
-    constexpr inline std::array<ElementType, sizeof...(Tags)>& array() noexcept
+    KOKKOS_FUNCTION constexpr std::array<ElementType, sizeof...(Tags)>& array() noexcept
     {
         return m_values;
     }
 
     /// Returns a const reference to the underlying `std::array`
-    constexpr inline std::array<ElementType, sizeof...(Tags)> const& array() const noexcept
+    KOKKOS_FUNCTION constexpr std::array<ElementType, sizeof...(Tags)> const& array() const noexcept
     {
         return m_values;
     }
 
-    constexpr inline ElementType& operator[](size_t pos)
+    KOKKOS_FUNCTION constexpr ElementType& operator[](size_t pos)
     {
         return m_values[pos];
     }
 
-    constexpr inline ElementType const& operator[](size_t pos) const
+    KOKKOS_FUNCTION constexpr ElementType const& operator[](size_t pos) const
     {
         return m_values[pos];
     }
 
     template <class OElementType, class... OTags>
-    constexpr inline bool operator==(TaggedVector<OElementType, OTags...> const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr bool operator==(
+            TaggedVector<OElementType, OTags...> const& rhs) const noexcept
     {
         return ((m_values[type_seq_rank_v<Tags, tags_seq>] == rhs.template get<Tags>()) && ...);
     }
 
     template <class OElementType, class... OTags>
-    constexpr inline bool operator!=(TaggedVector<OElementType, OTags...> const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr bool operator!=(
+            TaggedVector<OElementType, OTags...> const& rhs) const noexcept
     {
         return !(*this == rhs);
     }
 
     template <class QueryTag>
-    inline constexpr ElementType& get() noexcept
+    KOKKOS_FUNCTION constexpr ElementType& get() noexcept
     {
         using namespace detail;
         static_assert(in_tags_v<QueryTag, tags_seq>, "requested Tag absent from TaggedVector");
@@ -339,7 +348,7 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     }
 
     template <class QueryTag>
-    inline constexpr ElementType const& get() const noexcept
+    KOKKOS_FUNCTION constexpr ElementType const& get() const noexcept
     {
         using namespace detail;
         static_assert(in_tags_v<QueryTag, tags_seq>, "requested Tag absent from TaggedVector");
@@ -347,7 +356,7 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     }
 
     template <class QueryTag>
-    ElementType const& get_or(ElementType const& default_value) const&
+    KOKKOS_FUNCTION constexpr ElementType const& get_or(ElementType const& default_value) const&
     {
         DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
         if constexpr (in_tags_v<QueryTag, tags_seq>) {
@@ -359,13 +368,14 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     }
 
     template <std::size_t N = sizeof...(Tags)>
-    constexpr inline std::enable_if_t<N == 1, ElementType const&> value() const noexcept
+    KOKKOS_FUNCTION constexpr std::enable_if_t<N == 1, ElementType const&> value() const noexcept
     {
         return m_values[0];
     }
 
     template <class OElementType, class... OTags>
-    constexpr inline TaggedVector& operator+=(TaggedVector<OElementType, OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr TaggedVector& operator+=(
+            TaggedVector<OElementType, OTags...> const& rhs)
     {
         static_assert(type_seq_same_v<tags_seq, detail::TypeSeq<OTags...>>);
         ((m_values[type_seq_rank_v<Tags, tags_seq>] += rhs.template get<Tags>()), ...);
@@ -373,14 +383,15 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     }
 
     template <class OElementType>
-    constexpr inline TaggedVector& operator+=(OElementType const& rhs)
+    KOKKOS_FUNCTION constexpr TaggedVector& operator+=(OElementType const& rhs)
     {
         ((m_values[type_seq_rank_v<Tags, tags_seq>] += rhs), ...);
         return *this;
     }
 
     template <class OElementType, class... OTags>
-    constexpr inline TaggedVector& operator-=(TaggedVector<OElementType, OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr TaggedVector& operator-=(
+            TaggedVector<OElementType, OTags...> const& rhs)
     {
         static_assert(type_seq_same_v<tags_seq, detail::TypeSeq<OTags...>>);
         ((m_values[type_seq_rank_v<Tags, tags_seq>] -= rhs.template get<Tags>()), ...);
@@ -388,14 +399,15 @@ class TaggedVector : public ConversionOperators<TaggedVector<ElementType, Tags..
     }
 
     template <class OElementType>
-    constexpr inline TaggedVector& operator-=(OElementType const& rhs)
+    KOKKOS_FUNCTION constexpr TaggedVector& operator-=(OElementType const& rhs)
     {
         ((m_values[type_seq_rank_v<Tags, tags_seq>] -= rhs), ...);
         return *this;
     }
 
     template <class OElementType, class... OTags>
-    constexpr inline TaggedVector& operator*=(TaggedVector<OElementType, OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr TaggedVector& operator*=(
+            TaggedVector<OElementType, OTags...> const& rhs)
     {
         static_assert(type_seq_same_v<tags_seq, detail::TypeSeq<OTags...>>);
         ((m_values[type_seq_rank_v<Tags, tags_seq>] *= rhs.template get<Tags>()), ...);
diff --git a/include/ddc/discrete_domain.hpp b/include/ddc/discrete_domain.hpp
index 7d2f2e2b9..af96066b6 100644
--- a/include/ddc/discrete_domain.hpp
+++ b/include/ddc/discrete_domain.hpp
@@ -36,16 +36,16 @@ class DiscreteDomain
     DiscreteElement<DDims...> m_element_end;
 
 public:
-    static constexpr std::size_t rank()
+    static KOKKOS_FUNCTION constexpr std::size_t rank()
     {
         return sizeof...(DDims);
     }
 
-    DiscreteDomain() = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain() = default;
 
     /// Construct a DiscreteDomain by copies and merge of domains
     template <class... DD>
-    explicit constexpr DiscreteDomain(DD const&... domains)
+    explicit KOKKOS_FUNCTION constexpr DiscreteDomain(DD const&... domains)
         : m_element_begin(domains.front()...)
         , m_element_end((domains.front() + domains.extents())...)
     {
@@ -55,24 +55,26 @@ class DiscreteDomain
      * @param element_begin the lower bound in each direction
      * @param size the number of points in each direction
      */
-    constexpr DiscreteDomain(discrete_element_type const& element_begin, mlength_type const& size)
+    KOKKOS_FUNCTION constexpr DiscreteDomain(
+            discrete_element_type const& element_begin,
+            mlength_type const& size)
         : m_element_begin(element_begin)
         , m_element_end(element_begin + size)
     {
     }
 
-    DiscreteDomain(DiscreteDomain const& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain(DiscreteDomain const& x) = default;
 
-    DiscreteDomain(DiscreteDomain&& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain(DiscreteDomain&& x) = default;
 
-    ~DiscreteDomain() = default;
+    KOKKOS_DEFAULTED_FUNCTION ~DiscreteDomain() = default;
 
-    DiscreteDomain& operator=(DiscreteDomain const& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain& operator=(DiscreteDomain const& x) = default;
 
-    DiscreteDomain& operator=(DiscreteDomain&& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain& operator=(DiscreteDomain&& x) = default;
 
     template <class... ODims>
-    constexpr bool operator==(DiscreteDomain<ODims...> const& other) const
+    KOKKOS_FUNCTION constexpr bool operator==(DiscreteDomain<ODims...> const& other) const
     {
         return m_element_begin == other.m_element_begin && m_element_end == other.m_element_end;
     }
@@ -80,66 +82,66 @@ class DiscreteDomain
 #if __cplusplus <= 201703L
     // Shall not be necessary anymore in C++20
     // `a!=b` shall be translated by the compiler to `!(a==b)`
-    constexpr bool operator!=(DiscreteDomain const& other) const
+    KOKKOS_FUNCTION constexpr bool operator!=(DiscreteDomain const& other) const
     {
         return !(*this == other);
     }
 #endif
 
-    constexpr std::size_t size() const
+    KOKKOS_FUNCTION constexpr std::size_t size() const
     {
         return (1ul * ... * (uid<DDims>(m_element_end) - uid<DDims>(m_element_begin)));
     }
 
-    constexpr mlength_type extents() const noexcept
+    KOKKOS_FUNCTION constexpr mlength_type extents() const noexcept
     {
         return mlength_type((uid<DDims>(m_element_end) - uid<DDims>(m_element_begin))...);
     }
 
     template <class QueryDDim>
-    inline constexpr DiscreteVector<QueryDDim> extent() const noexcept
+    KOKKOS_FUNCTION constexpr DiscreteVector<QueryDDim> extent() const noexcept
     {
         return DiscreteVector<QueryDDim>(
                 uid<QueryDDim>(m_element_end) - uid<QueryDDim>(m_element_begin));
     }
 
-    constexpr discrete_element_type front() const noexcept
+    KOKKOS_FUNCTION constexpr discrete_element_type front() const noexcept
     {
         return m_element_begin;
     }
 
-    constexpr discrete_element_type back() const noexcept
+    KOKKOS_FUNCTION constexpr discrete_element_type back() const noexcept
     {
         return discrete_element_type((uid<DDims>(m_element_end) - 1)...);
     }
 
-    constexpr DiscreteDomain take_first(mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain take_first(mlength_type n) const
     {
         return DiscreteDomain(front(), n);
     }
 
-    constexpr DiscreteDomain take_last(mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain take_last(mlength_type n) const
     {
         return DiscreteDomain(front() + (extents() - n), n);
     }
 
-    constexpr DiscreteDomain remove_first(mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain remove_first(mlength_type n) const
     {
         return DiscreteDomain(front() + n, extents() - n);
     }
 
-    constexpr DiscreteDomain remove_last(mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain remove_last(mlength_type n) const
     {
         return DiscreteDomain(front(), extents() - n);
     }
 
-    constexpr DiscreteDomain remove(mlength_type n1, mlength_type n2) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain remove(mlength_type n1, mlength_type n2) const
     {
         return DiscreteDomain(front() + n1, extents() - n1 - n2);
     }
 
     template <class... ODDims>
-    constexpr auto restrict(DiscreteDomain<ODDims...> const& odomain) const
+    KOKKOS_FUNCTION constexpr auto restrict(DiscreteDomain<ODDims...> const& odomain) const
     {
         assert(((uid<ODDims>(m_element_begin) <= uid<ODDims>(odomain.m_element_begin)) && ...));
         assert(((uid<ODDims>(m_element_end) >= uid<ODDims>(odomain.m_element_end)) && ...));
@@ -151,12 +153,12 @@ class DiscreteDomain
                 DiscreteVector<DDims...>((get_or<DDims>(oextents, get<DDims>(myextents)))...));
     }
 
-    constexpr bool empty() const noexcept
+    KOKKOS_FUNCTION constexpr bool empty() const noexcept
     {
         return size() == 0;
     }
 
-    constexpr explicit operator bool()
+    KOKKOS_FUNCTION constexpr explicit operator bool()
     {
         return !empty();
     }
@@ -164,7 +166,7 @@ class DiscreteDomain
     template <
             std::size_t N = sizeof...(DDims),
             class DDim0 = std::enable_if_t<N == 1, std::tuple_element_t<0, std::tuple<DDims...>>>>
-    auto begin() const
+    KOKKOS_FUNCTION auto begin() const
     {
         return DiscreteDomainIterator<DDim0>(front());
     }
@@ -172,7 +174,7 @@ class DiscreteDomain
     template <
             std::size_t N = sizeof...(DDims),
             class DDim0 = std::enable_if_t<N == 1, std::tuple_element_t<0, std::tuple<DDims...>>>>
-    auto end() const
+    KOKKOS_FUNCTION auto end() const
     {
         return DiscreteDomainIterator<DDim0>(m_element_end);
     }
@@ -180,7 +182,7 @@ class DiscreteDomain
     template <
             std::size_t N = sizeof...(DDims),
             class DDim0 = std::enable_if_t<N == 1, std::tuple_element_t<0, std::tuple<DDims...>>>>
-    auto cbegin() const
+    KOKKOS_FUNCTION auto cbegin() const
     {
         return DiscreteDomainIterator<DDim0>(front());
     }
@@ -188,7 +190,7 @@ class DiscreteDomain
     template <
             std::size_t N = sizeof...(DDims),
             class DDim0 = std::enable_if_t<N == 1, std::tuple_element_t<0, std::tuple<DDims...>>>>
-    auto cend() const
+    KOKKOS_FUNCTION auto cend() const
     {
         return DiscreteDomainIterator<DDim0>(m_element_end);
     }
@@ -196,7 +198,7 @@ class DiscreteDomain
     template <
             std::size_t N = sizeof...(DDims),
             class = std::enable_if_t<N == 1, std::tuple_element_t<0, std::tuple<DDims...>>>>
-    constexpr decltype(auto) operator[](std::size_t n)
+    KOKKOS_FUNCTION constexpr decltype(auto) operator[](std::size_t n)
     {
         return begin()[n];
     }
@@ -204,7 +206,7 @@ class DiscreteDomain
     template <
             std::size_t N = sizeof...(DDims),
             class = std::enable_if_t<N == 1, std::tuple_element_t<0, std::tuple<DDims...>>>>
-    constexpr decltype(auto) operator[](std::size_t n) const
+    KOKKOS_FUNCTION constexpr decltype(auto) operator[](std::size_t n) const
     {
         return begin()[n];
     }
@@ -221,7 +223,7 @@ class DiscreteDomain<>
 
     using mlength_type = DiscreteVector<>;
 
-    static constexpr std::size_t rank()
+    static KOKKOS_FUNCTION constexpr std::size_t rank()
     {
         return 0;
     }
@@ -230,7 +232,8 @@ class DiscreteDomain<>
 
     // Construct a DiscreteDomain from a reordered copy of `domain`
     template <class... ODDims>
-    explicit constexpr DiscreteDomain([[maybe_unused]] DiscreteDomain<ODDims...> const& domain)
+    explicit KOKKOS_FUNCTION constexpr DiscreteDomain(
+            [[maybe_unused]] DiscreteDomain<ODDims...> const& domain)
     {
     }
 
@@ -238,23 +241,23 @@ class DiscreteDomain<>
      * @param element_begin the lower bound in each direction
      * @param size the number of points in each direction
      */
-    constexpr DiscreteDomain(
+    KOKKOS_FUNCTION constexpr DiscreteDomain(
             [[maybe_unused]] discrete_element_type const& element_begin,
             [[maybe_unused]] mlength_type const& size)
     {
     }
 
-    constexpr DiscreteDomain(DiscreteDomain const& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain(DiscreteDomain const& x) = default;
 
-    constexpr DiscreteDomain(DiscreteDomain&& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain(DiscreteDomain&& x) = default;
 
-    ~DiscreteDomain() = default;
+    KOKKOS_DEFAULTED_FUNCTION ~DiscreteDomain() = default;
 
-    DiscreteDomain& operator=(DiscreteDomain const& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain& operator=(DiscreteDomain const& x) = default;
 
-    DiscreteDomain& operator=(DiscreteDomain&& x) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomain& operator=(DiscreteDomain&& x) = default;
 
-    constexpr bool operator==([[maybe_unused]] DiscreteDomain const& other) const
+    KOKKOS_FUNCTION constexpr bool operator==([[maybe_unused]] DiscreteDomain const& other) const
     {
         return true;
     }
@@ -262,53 +265,53 @@ class DiscreteDomain<>
 #if __cplusplus <= 201703L
     // Shall not be necessary anymore in C++20
     // `a!=b` shall be translated by the compiler to `!(a==b)`
-    constexpr bool operator!=(DiscreteDomain const& other) const
+    KOKKOS_FUNCTION constexpr bool operator!=(DiscreteDomain const& other) const
     {
         return !(*this == other);
     }
 #endif
 
-    constexpr std::size_t size() const
+    KOKKOS_FUNCTION constexpr std::size_t size() const
     {
         return 1;
     }
 
-    constexpr mlength_type extents() const noexcept
+    KOKKOS_FUNCTION constexpr mlength_type extents() const noexcept
     {
         return {};
     }
 
-    constexpr discrete_element_type front() const noexcept
+    KOKKOS_FUNCTION constexpr discrete_element_type front() const noexcept
     {
         return {};
     }
 
-    constexpr discrete_element_type back() const noexcept
+    KOKKOS_FUNCTION constexpr discrete_element_type back() const noexcept
     {
         return {};
     }
 
-    constexpr DiscreteDomain take_first([[maybe_unused]] mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain take_first([[maybe_unused]] mlength_type n) const
     {
         return *this;
     }
 
-    constexpr DiscreteDomain take_last([[maybe_unused]] mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain take_last([[maybe_unused]] mlength_type n) const
     {
         return *this;
     }
 
-    constexpr DiscreteDomain remove_first([[maybe_unused]] mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain remove_first([[maybe_unused]] mlength_type n) const
     {
         return *this;
     }
 
-    constexpr DiscreteDomain remove_last([[maybe_unused]] mlength_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain remove_last([[maybe_unused]] mlength_type n) const
     {
         return *this;
     }
 
-    constexpr DiscreteDomain remove(
+    KOKKOS_FUNCTION constexpr DiscreteDomain remove(
             [[maybe_unused]] mlength_type n1,
             [[maybe_unused]] mlength_type n2) const
     {
@@ -316,24 +319,25 @@ class DiscreteDomain<>
     }
 
     template <class... ODims>
-    constexpr DiscreteDomain restrict(DiscreteDomain<ODims...> const&) const
+    KOKKOS_FUNCTION constexpr DiscreteDomain restrict(DiscreteDomain<ODims...> const&) const
     {
         return *this;
     }
 
-    constexpr bool empty() const noexcept
+    KOKKOS_FUNCTION constexpr bool empty() const noexcept
     {
         return false;
     }
 
-    constexpr explicit operator bool()
+    KOKKOS_FUNCTION constexpr explicit operator bool()
     {
         return true;
     }
 };
 
 template <class... QueryDDims, class... DDims>
-constexpr DiscreteDomain<QueryDDims...> select(DiscreteDomain<DDims...> const& domain)
+KOKKOS_FUNCTION constexpr DiscreteDomain<QueryDDims...> select(
+        DiscreteDomain<DDims...> const& domain)
 {
     return DiscreteDomain<QueryDDims...>(
             select<QueryDDims...>(domain.front()),
@@ -385,7 +389,7 @@ using cartesian_prod_t = typename cartesian_prod<DDom...>::type;
 
 // Computes the substraction DDom_a - DDom_b in the sense of linear spaces(retained dimensions are those in DDom_a which are not in DDom_b)
 template <class... DDimsA, class... DDimsB>
-constexpr auto remove_dims_of(
+KOKKOS_FUNCTION constexpr auto remove_dims_of(
         DiscreteDomain<DDimsA...> const& DDom_a,
         [[maybe_unused]] DiscreteDomain<DDimsB...> const& DDom_b) noexcept
 {
@@ -399,7 +403,7 @@ constexpr auto remove_dims_of(
 
 // Checks if dimension of DDom_a is DDim1. If not, returns restriction to DDim2 of DDom_b. May not be usefull in its own, it helps for replace_dim_of
 template <typename DDim1, typename DDim2, typename DDimA, typename... DDimsB>
-constexpr std::conditional_t<
+KOKKOS_FUNCTION constexpr std::conditional_t<
         std::is_same_v<DDimA, DDim1>,
         ddc::DiscreteDomain<DDim2>,
         ddc::DiscreteDomain<DDimA>>
@@ -416,7 +420,7 @@ replace_dim_of_1d(
 
 // Replace in DDom_a the dimension Dim1 by the dimension Dim2 of DDom_b
 template <typename DDim1, typename DDim2, typename... DDimsA, typename... DDimsB>
-constexpr auto replace_dim_of(
+KOKKOS_FUNCTION constexpr auto replace_dim_of(
         DiscreteDomain<DDimsA...> const& DDom_a,
         [[maybe_unused]] DiscreteDomain<DDimsB...> const& DDom_b) noexcept
 {
@@ -435,19 +439,22 @@ constexpr auto replace_dim_of(
 }
 
 template <class... QueryDDims, class... DDims>
-constexpr DiscreteVector<QueryDDims...> extents(DiscreteDomain<DDims...> const& domain) noexcept
+KOKKOS_FUNCTION constexpr DiscreteVector<QueryDDims...> extents(
+        DiscreteDomain<DDims...> const& domain) noexcept
 {
     return DiscreteVector<QueryDDims...>(select<QueryDDims>(domain).size()...);
 }
 
 template <class... QueryDDims, class... DDims>
-constexpr DiscreteElement<QueryDDims...> front(DiscreteDomain<DDims...> const& domain) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElement<QueryDDims...> front(
+        DiscreteDomain<DDims...> const& domain) noexcept
 {
     return DiscreteElement<QueryDDims...>(select<QueryDDims>(domain).front()...);
 }
 
 template <class... QueryDDims, class... DDims>
-constexpr DiscreteElement<QueryDDims...> back(DiscreteDomain<DDims...> const& domain) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElement<QueryDDims...> back(
+        DiscreteDomain<DDims...> const& domain) noexcept
 {
     return DiscreteElement<QueryDDims...>(select<QueryDDims>(domain).back()...);
 }
@@ -461,7 +468,7 @@ template <class... QueryDDims>
 struct Selection<detail::TypeSeq<QueryDDims...>>
 {
     template <class Domain>
-    static constexpr auto select(Domain const& domain)
+    static KOKKOS_FUNCTION constexpr auto select(Domain const& domain)
     {
         return ddc::select<QueryDDims...>(domain);
     }
@@ -470,7 +477,7 @@ struct Selection<detail::TypeSeq<QueryDDims...>>
 } // namespace detail
 
 template <class QueryDDimSeq, class... DDims>
-constexpr auto select_by_type_seq(DiscreteDomain<DDims...> const& domain)
+KOKKOS_FUNCTION constexpr auto select_by_type_seq(DiscreteDomain<DDims...> const& domain)
 {
     return detail::Selection<QueryDDimSeq>::select(domain);
 }
@@ -488,42 +495,45 @@ struct DiscreteDomainIterator
 
     using difference_type = std::ptrdiff_t;
 
-    DiscreteDomainIterator() = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteDomainIterator() = default;
 
-    constexpr explicit DiscreteDomainIterator(DiscreteElement<DDim> value) : m_value(value) {}
+    KOKKOS_FUNCTION constexpr explicit DiscreteDomainIterator(DiscreteElement<DDim> value)
+        : m_value(value)
+    {
+    }
 
-    constexpr DiscreteElement<DDim> operator*() const noexcept
+    KOKKOS_FUNCTION constexpr DiscreteElement<DDim> operator*() const noexcept
     {
         return m_value;
     }
 
-    constexpr DiscreteDomainIterator& operator++()
+    KOKKOS_FUNCTION constexpr DiscreteDomainIterator& operator++()
     {
         ++m_value.uid();
         return *this;
     }
 
-    constexpr DiscreteDomainIterator operator++(int)
+    KOKKOS_FUNCTION constexpr DiscreteDomainIterator operator++(int)
     {
         auto tmp = *this;
         ++*this;
         return tmp;
     }
 
-    constexpr DiscreteDomainIterator& operator--()
+    KOKKOS_FUNCTION constexpr DiscreteDomainIterator& operator--()
     {
         --m_value.uid();
         return *this;
     }
 
-    constexpr DiscreteDomainIterator operator--(int)
+    KOKKOS_FUNCTION constexpr DiscreteDomainIterator operator--(int)
     {
         auto tmp = *this;
         --*this;
         return tmp;
     }
 
-    constexpr DiscreteDomainIterator& operator+=(difference_type n)
+    KOKKOS_FUNCTION constexpr DiscreteDomainIterator& operator+=(difference_type n)
     {
         if (n >= difference_type(0))
             m_value.uid() += static_cast<DiscreteElementType>(n);
@@ -532,7 +542,7 @@ struct DiscreteDomainIterator
         return *this;
     }
 
-    constexpr DiscreteDomainIterator& operator-=(difference_type n)
+    KOKKOS_FUNCTION constexpr DiscreteDomainIterator& operator-=(difference_type n)
     {
         if (n >= difference_type(0))
             m_value.uid() -= static_cast<DiscreteElementType>(n);
@@ -541,69 +551,75 @@ struct DiscreteDomainIterator
         return *this;
     }
 
-    constexpr DiscreteElement<DDim> operator[](difference_type n) const
+    KOKKOS_FUNCTION constexpr DiscreteElement<DDim> operator[](difference_type n) const
     {
         return m_value + n;
     }
 
-    friend constexpr bool operator==(
+    friend KOKKOS_FUNCTION constexpr bool operator==(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
         return xx.m_value == yy.m_value;
     }
 
-    friend constexpr bool operator!=(
+    friend KOKKOS_FUNCTION constexpr bool operator!=(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
         return xx.m_value != yy.m_value;
     }
 
-    friend constexpr bool operator<(
+    friend KOKKOS_FUNCTION constexpr bool operator<(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
         return xx.m_value < yy.m_value;
     }
 
-    friend constexpr bool operator>(
+    friend KOKKOS_FUNCTION constexpr bool operator>(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
         return yy < xx;
     }
 
-    friend constexpr bool operator<=(
+    friend KOKKOS_FUNCTION constexpr bool operator<=(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
         return !(yy < xx);
     }
 
-    friend constexpr bool operator>=(
+    friend KOKKOS_FUNCTION constexpr bool operator>=(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
         return !(xx < yy);
     }
 
-    friend constexpr DiscreteDomainIterator operator+(DiscreteDomainIterator i, difference_type n)
+    friend KOKKOS_FUNCTION constexpr DiscreteDomainIterator operator+(
+            DiscreteDomainIterator i,
+            difference_type n)
     {
         return i += n;
     }
 
-    friend constexpr DiscreteDomainIterator operator+(difference_type n, DiscreteDomainIterator i)
+    friend KOKKOS_FUNCTION constexpr DiscreteDomainIterator operator+(
+            difference_type n,
+            DiscreteDomainIterator i)
     {
         return i += n;
     }
 
-    friend constexpr DiscreteDomainIterator operator-(DiscreteDomainIterator i, difference_type n)
+    friend KOKKOS_FUNCTION constexpr DiscreteDomainIterator operator-(
+            DiscreteDomainIterator i,
+            difference_type n)
     {
         return i -= n;
     }
 
-    friend constexpr difference_type operator-(
+    friend KOKKOS_FUNCTION constexpr difference_type operator-(
             DiscreteDomainIterator const& xx,
             DiscreteDomainIterator const& yy)
     {
diff --git a/include/ddc/discrete_element.hpp b/include/ddc/discrete_element.hpp
index 0d0e438e8..3907f1b40 100644
--- a/include/ddc/discrete_element.hpp
+++ b/include/ddc/discrete_element.hpp
@@ -35,31 +35,32 @@ inline constexpr bool is_discrete_element_v = IsDiscreteElement<T>::value;
 using DiscreteElementType = std::size_t;
 
 template <class Tag>
-inline constexpr DiscreteElementType const& uid(DiscreteElement<Tag> const& tuple) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElementType const& uid(DiscreteElement<Tag> const& tuple) noexcept
 {
     return tuple.uid();
 }
 
 template <class Tag>
-inline constexpr DiscreteElementType& uid(DiscreteElement<Tag>& tuple) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElementType& uid(DiscreteElement<Tag>& tuple) noexcept
 {
     return tuple.uid();
 }
 
 template <class QueryTag, class... Tags>
-inline constexpr DiscreteElementType const& uid(DiscreteElement<Tags...> const& tuple) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElementType const& uid(
+        DiscreteElement<Tags...> const& tuple) noexcept
 {
     return tuple.template uid<QueryTag>();
 }
 
 template <class QueryTag, class... Tags>
-inline constexpr DiscreteElementType& uid(DiscreteElement<Tags...>& tuple) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElementType& uid(DiscreteElement<Tags...>& tuple) noexcept
 {
     return tuple.template uid<QueryTag>();
 }
 
 template <class QueryTag, class... Tags>
-inline constexpr DiscreteElementType const& uid_or(
+KOKKOS_FUNCTION constexpr DiscreteElementType const& uid_or(
         DiscreteElement<Tags...> const& tuple,
         DiscreteElementType const& default_value) noexcept
 {
@@ -67,19 +68,21 @@ inline constexpr DiscreteElementType const& uid_or(
 }
 
 template <class... QueryTags, class... Tags>
-inline constexpr DiscreteElement<QueryTags...> select(DiscreteElement<Tags...> const& arr) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElement<QueryTags...> select(
+        DiscreteElement<Tags...> const& arr) noexcept
 {
     return DiscreteElement<QueryTags...>(arr);
 }
 
 template <class... QueryTags, class... Tags>
-inline constexpr DiscreteElement<QueryTags...> select(DiscreteElement<Tags...>&& arr) noexcept
+KOKKOS_FUNCTION constexpr DiscreteElement<QueryTags...> select(
+        DiscreteElement<Tags...>&& arr) noexcept
 {
     return DiscreteElement<QueryTags...>(std::move(arr));
 }
 
 template <class QueryTag, class... HeadTag, class... DETail>
-constexpr DiscreteElement<QueryTag> take(
+KOKKOS_FUNCTION constexpr DiscreteElement<QueryTag> take(
         DiscreteElement<HeadTag...> const& head,
         DETail const&... tail)
 {
@@ -101,7 +104,7 @@ namespace detail {
 
 /// Returns a reference to the underlying `std::array`
 template <class... Tags>
-constexpr inline std::array<DiscreteElementType, sizeof...(Tags)>& array(
+KOKKOS_FUNCTION constexpr std::array<DiscreteElementType, sizeof...(Tags)>& array(
         DiscreteElement<Tags...>& v) noexcept
 {
     return v.m_values;
@@ -109,7 +112,7 @@ constexpr inline std::array<DiscreteElementType, sizeof...(Tags)>& array(
 
 /// Returns a reference to the underlying `std::array`
 template <class... Tags>
-constexpr inline std::array<DiscreteElementType, sizeof...(Tags)> const& array(
+KOKKOS_FUNCTION constexpr std::array<DiscreteElementType, sizeof...(Tags)> const& array(
         DiscreteElement<Tags...> const& v) noexcept
 {
     return v.m_values;
@@ -126,11 +129,11 @@ class DiscreteElement
 {
     using tags_seq = detail::TypeSeq<Tags...>;
 
-    friend constexpr std::array<DiscreteElementType, sizeof...(Tags)>& detail::array<Tags...>(
-            DiscreteElement<Tags...>& v) noexcept;
+    friend KOKKOS_FUNCTION constexpr std::array<DiscreteElementType, sizeof...(Tags)>& detail::
+            array<Tags...>(DiscreteElement<Tags...>& v) noexcept;
 
-    friend constexpr std::array<DiscreteElementType, sizeof...(Tags)> const& detail::array<Tags...>(
-            DiscreteElement<Tags...> const& v) noexcept;
+    friend KOKKOS_FUNCTION constexpr std::array<DiscreteElementType, sizeof...(Tags)> const&
+    detail::array<Tags...>(DiscreteElement<Tags...> const& v) noexcept;
 
 private:
     std::array<DiscreteElementType, sizeof...(Tags)> m_values;
@@ -138,26 +141,27 @@ class DiscreteElement
 public:
     using value_type = DiscreteElementType;
 
-    static constexpr std::size_t size() noexcept
+    static KOKKOS_FUNCTION constexpr std::size_t size() noexcept
     {
         return sizeof...(Tags);
     }
 
 public:
-    inline constexpr DiscreteElement() = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteElement() = default;
 
-    inline constexpr DiscreteElement(DiscreteElement const&) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteElement(DiscreteElement const&) = default;
 
-    inline constexpr DiscreteElement(DiscreteElement&&) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteElement(DiscreteElement&&) = default;
 
     template <class... OTags>
-    explicit inline constexpr DiscreteElement(DiscreteElement<OTags...> const& other) noexcept
+    explicit KOKKOS_FUNCTION constexpr DiscreteElement(
+            DiscreteElement<OTags...> const& other) noexcept
         : m_values {other.template uid<Tags>()...}
     {
     }
 
     template <class... DE, class = std::enable_if_t<(is_discrete_element_v<DE> && ...)>>
-    explicit inline constexpr DiscreteElement(DE const&... other) noexcept
+    explicit KOKKOS_FUNCTION constexpr DiscreteElement(DE const&... other) noexcept
         : m_values {take<Tags>(other...).uid()...}
     {
     }
@@ -167,31 +171,34 @@ class DiscreteElement
             class = std::enable_if_t<(std::is_integral_v<Params> && ...)>,
             class = std::enable_if_t<(!is_discrete_element_v<Params> && ...)>,
             class = std::enable_if_t<sizeof...(Params) == sizeof...(Tags)>>
-    explicit inline constexpr DiscreteElement(Params const&... params) noexcept
+    explicit KOKKOS_FUNCTION constexpr DiscreteElement(Params const&... params) noexcept
         : m_values {static_cast<value_type>(params)...}
     {
     }
 
-    constexpr inline DiscreteElement& operator=(DiscreteElement const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION ~DiscreteElement() = default;
 
-    constexpr inline DiscreteElement& operator=(DiscreteElement&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION DiscreteElement& operator=(DiscreteElement const& other) = default;
+
+    KOKKOS_DEFAULTED_FUNCTION DiscreteElement& operator=(DiscreteElement&& other) = default;
 
     template <class... OTags>
-    constexpr inline DiscreteElement& operator=(DiscreteElement<OTags...> const& other) noexcept
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator=(
+            DiscreteElement<OTags...> const& other) noexcept
     {
         m_values = other.m_values;
         return *this;
     }
 
     template <class... OTags>
-    constexpr inline DiscreteElement& operator=(DiscreteElement<OTags...>&& other) noexcept
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator=(DiscreteElement<OTags...>&& other) noexcept
     {
         m_values = std::move(other.m_values);
         return *this;
     }
 
     template <class QueryTag>
-    value_type const& uid_or(value_type const& default_value) const&
+    KOKKOS_FUNCTION constexpr value_type const& uid_or(value_type const& default_value) const&
     {
         DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
         if constexpr (in_tags_v<QueryTag, tags_seq>) {
@@ -203,7 +210,7 @@ class DiscreteElement
     }
 
     template <class QueryTag>
-    inline constexpr value_type& uid() noexcept
+    KOKKOS_FUNCTION constexpr value_type& uid() noexcept
     {
         using namespace detail;
         static_assert(in_tags_v<QueryTag, tags_seq>, "requested Tag absent from DiscreteElement");
@@ -211,7 +218,7 @@ class DiscreteElement
     }
 
     template <class QueryTag>
-    inline constexpr value_type const& uid() const noexcept
+    KOKKOS_FUNCTION constexpr value_type const& uid() const noexcept
     {
         using namespace detail;
         static_assert(in_tags_v<QueryTag, tags_seq>, "requested Tag absent from DiscreteElement");
@@ -219,26 +226,26 @@ class DiscreteElement
     }
 
     template <std::size_t N = sizeof...(Tags)>
-    constexpr inline std::enable_if_t<N == 1, value_type&> uid() noexcept
+    KOKKOS_FUNCTION constexpr std::enable_if_t<N == 1, value_type&> uid() noexcept
     {
         return m_values[0];
     }
 
     template <std::size_t N = sizeof...(Tags)>
-    constexpr inline std::enable_if_t<N == 1, value_type const&> uid() const noexcept
+    KOKKOS_FUNCTION constexpr std::enable_if_t<N == 1, value_type const&> uid() const noexcept
     {
         return m_values[0];
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteElement& operator++()
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator++()
     {
         ++m_values[0];
         return *this;
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteElement operator++(int)
+    KOKKOS_FUNCTION constexpr DiscreteElement operator++(int)
     {
         DiscreteElement const tmp = *this;
         ++m_values[0];
@@ -246,14 +253,14 @@ class DiscreteElement
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteElement& operator--()
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator--()
     {
         ++m_values[0];
         return *this;
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteElement operator--(int)
+    KOKKOS_FUNCTION constexpr DiscreteElement operator--(int)
     {
         DiscreteElement const tmp = *this;
         ++m_values[0];
@@ -261,7 +268,7 @@ class DiscreteElement
     }
 
     template <class... OTags>
-    constexpr inline DiscreteElement& operator+=(DiscreteVector<OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator+=(DiscreteVector<OTags...> const& rhs)
     {
         static_assert(((type_seq_contains_v<detail::TypeSeq<OTags>, tags_seq>)&&...));
         ((m_values[type_seq_rank_v<OTags, tags_seq>] += rhs.template get<OTags>()), ...);
@@ -273,14 +280,14 @@ class DiscreteElement
             std::size_t N = sizeof...(Tags),
             class = std::enable_if_t<N == 1>,
             class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-    constexpr inline DiscreteElement& operator+=(IntegralType const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator+=(IntegralType const& rhs)
     {
         m_values[0] += rhs;
         return *this;
     }
 
     template <class... OTags>
-    constexpr inline DiscreteElement& operator-=(DiscreteVector<OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator-=(DiscreteVector<OTags...> const& rhs)
     {
         static_assert(((type_seq_contains_v<detail::TypeSeq<OTags>, tags_seq>)&&...));
         ((m_values[type_seq_rank_v<OTags, tags_seq>] -= rhs.template get<OTags>()), ...);
@@ -292,7 +299,7 @@ class DiscreteElement
             std::size_t N = sizeof...(Tags),
             class = std::enable_if_t<N == 1>,
             class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-    constexpr inline DiscreteElement& operator-=(IntegralType const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteElement& operator-=(IntegralType const& rhs)
     {
         m_values[0] -= rhs;
         return *this;
@@ -317,7 +324,7 @@ std::ostream& operator<<(std::ostream& out, DiscreteElement<Head, Tags...> const
 
 
 template <class... Tags, class... OTags>
-constexpr inline bool operator==(
+KOKKOS_FUNCTION constexpr bool operator==(
         DiscreteElement<Tags...> const& lhs,
         DiscreteElement<OTags...> const& rhs) noexcept
 {
@@ -325,7 +332,7 @@ constexpr inline bool operator==(
 }
 
 template <class... Tags, class... OTags>
-constexpr inline bool operator!=(
+KOKKOS_FUNCTION constexpr bool operator!=(
         DiscreteElement<Tags...> const& lhs,
         DiscreteElement<OTags...> const& rhs) noexcept
 {
@@ -333,25 +340,33 @@ constexpr inline bool operator!=(
 }
 
 template <class Tag>
-constexpr inline bool operator<(DiscreteElement<Tag> const& lhs, DiscreteElement<Tag> const& rhs)
+KOKKOS_FUNCTION constexpr bool operator<(
+        DiscreteElement<Tag> const& lhs,
+        DiscreteElement<Tag> const& rhs)
 {
     return lhs.uid() < rhs.uid();
 }
 
 template <class Tag>
-constexpr inline bool operator<=(DiscreteElement<Tag> const& lhs, DiscreteElement<Tag> const& rhs)
+KOKKOS_FUNCTION constexpr bool operator<=(
+        DiscreteElement<Tag> const& lhs,
+        DiscreteElement<Tag> const& rhs)
 {
     return lhs.uid() <= rhs.uid();
 }
 
 template <class Tag>
-constexpr inline bool operator>(DiscreteElement<Tag> const& lhs, DiscreteElement<Tag> const& rhs)
+KOKKOS_FUNCTION constexpr bool operator>(
+        DiscreteElement<Tag> const& lhs,
+        DiscreteElement<Tag> const& rhs)
 {
     return lhs.uid() > rhs.uid();
 }
 
 template <class Tag>
-constexpr inline bool operator>=(DiscreteElement<Tag> const& lhs, DiscreteElement<Tag> const& rhs)
+KOKKOS_FUNCTION constexpr bool operator>=(
+        DiscreteElement<Tag> const& lhs,
+        DiscreteElement<Tag> const& rhs)
 {
     return lhs.uid() >= rhs.uid();
 }
@@ -359,7 +374,7 @@ constexpr inline bool operator>=(DiscreteElement<Tag> const& lhs, DiscreteElemen
 /// right external binary operators: +, -
 
 template <class... Tags, class... OTags>
-constexpr inline DiscreteElement<Tags...> operator+(
+KOKKOS_FUNCTION constexpr DiscreteElement<Tags...> operator+(
         DiscreteElement<Tags...> const& lhs,
         DiscreteVector<OTags...> const& rhs)
 {
@@ -375,7 +390,7 @@ template <
         class IntegralType,
         class = std::enable_if_t<std::is_integral_v<IntegralType>>,
         class = std::enable_if_t<!is_discrete_vector_v<IntegralType>>>
-constexpr inline DiscreteElement<Tag> operator+(
+KOKKOS_FUNCTION constexpr DiscreteElement<Tag> operator+(
         DiscreteElement<Tag> const& lhs,
         IntegralType const& rhs)
 {
@@ -383,7 +398,7 @@ constexpr inline DiscreteElement<Tag> operator+(
 }
 
 template <class... Tags, class... OTags>
-constexpr inline DiscreteElement<Tags...> operator-(
+KOKKOS_FUNCTION constexpr DiscreteElement<Tags...> operator-(
         DiscreteElement<Tags...> const& lhs,
         DiscreteVector<OTags...> const& rhs)
 {
@@ -399,7 +414,7 @@ template <
         class IntegralType,
         class = std::enable_if_t<std::is_integral_v<IntegralType>>,
         class = std::enable_if_t<!is_discrete_vector_v<IntegralType>>>
-constexpr inline DiscreteElement<Tag> operator-(
+KOKKOS_FUNCTION constexpr DiscreteElement<Tag> operator-(
         DiscreteElement<Tag> const& lhs,
         IntegralType const& rhs)
 {
@@ -409,7 +424,7 @@ constexpr inline DiscreteElement<Tag> operator-(
 /// binary operator: -
 
 template <class... Tags, class... OTags>
-constexpr inline DiscreteVector<Tags...> operator-(
+KOKKOS_FUNCTION constexpr DiscreteVector<Tags...> operator-(
         DiscreteElement<Tags...> const& lhs,
         DiscreteElement<OTags...> const& rhs)
 {
diff --git a/include/ddc/discrete_space.hpp b/include/ddc/discrete_space.hpp
index d3276bbdf..54e0df3b4 100644
--- a/include/ddc/discrete_space.hpp
+++ b/include/ddc/discrete_space.hpp
@@ -75,19 +75,19 @@ class gpu_proxy
     alignas(T) Kokkos::Array<std::byte, sizeof(T)> m_data;
 
 public:
-    DDC_INLINE_FUNCTION
+    KOKKOS_FORCEINLINE_FUNCTION
     T* operator->()
     {
         return reinterpret_cast<T*>(m_data.data());
     }
 
-    DDC_INLINE_FUNCTION
+    KOKKOS_FORCEINLINE_FUNCTION
     T& operator*()
     {
         return *reinterpret_cast<T*>(m_data.data());
     }
 
-    DDC_INLINE_FUNCTION
+    KOKKOS_FORCEINLINE_FUNCTION
     T* data()
     {
         return reinterpret_cast<T*>(m_data.data());
@@ -187,7 +187,7 @@ std::enable_if_t<2 <= sizeof...(Args), std::tuple<Args...>> init_discrete_space(
 }
 
 template <class DDim, class MemorySpace = DDC_CURRENT_KOKKOS_SPACE>
-DDC_INLINE_FUNCTION detail::ddim_impl_t<DDim, MemorySpace> const& discrete_space()
+KOKKOS_FORCEINLINE_FUNCTION detail::ddim_impl_t<DDim, MemorySpace> const& discrete_space()
 {
     if constexpr (std::is_same_v<MemorySpace, Kokkos::HostSpace>) {
         return detail::g_discrete_space_dual<DDim>->get_host();
diff --git a/include/ddc/discrete_vector.hpp b/include/ddc/discrete_vector.hpp
index acea572ec..09ee7d898 100644
--- a/include/ddc/discrete_vector.hpp
+++ b/include/ddc/discrete_vector.hpp
@@ -7,6 +7,8 @@
 #include <ostream>
 #include <utility>
 
+#include <Kokkos_Core.hpp>
+
 #include "ddc/detail/type_seq.hpp"
 
 namespace ddc {
@@ -33,19 +35,20 @@ inline constexpr bool is_discrete_vector_v = IsDiscreteVector<T>::value;
 using DiscreteVectorElement = std::ptrdiff_t;
 
 template <class QueryTag, class... Tags>
-inline constexpr DiscreteVectorElement const& get(DiscreteVector<Tags...> const& tuple) noexcept
+KOKKOS_FUNCTION constexpr DiscreteVectorElement const& get(
+        DiscreteVector<Tags...> const& tuple) noexcept
 {
     return tuple.template get<QueryTag>();
 }
 
 template <class QueryTag, class... Tags>
-inline constexpr DiscreteVectorElement& get(DiscreteVector<Tags...>& tuple) noexcept
+KOKKOS_FUNCTION constexpr DiscreteVectorElement& get(DiscreteVector<Tags...>& tuple) noexcept
 {
     return tuple.template get<QueryTag>();
 }
 
 template <class QueryTag, class... Tags>
-inline constexpr DiscreteVectorElement const& get_or(
+KOKKOS_FUNCTION constexpr DiscreteVectorElement const& get_or(
         DiscreteVector<Tags...> const& tuple,
         DiscreteVectorElement const& default_value) noexcept
 {
@@ -55,13 +58,13 @@ inline constexpr DiscreteVectorElement const& get_or(
 /// Unary operators: +, -
 
 template <class... Tags>
-constexpr inline DiscreteVector<Tags...> operator+(DiscreteVector<Tags...> const& x)
+KOKKOS_FUNCTION constexpr DiscreteVector<Tags...> operator+(DiscreteVector<Tags...> const& x)
 {
     return x;
 }
 
 template <class... Tags>
-constexpr inline DiscreteVector<Tags...> operator-(DiscreteVector<Tags...> const& x)
+KOKKOS_FUNCTION constexpr DiscreteVector<Tags...> operator-(DiscreteVector<Tags...> const& x)
 {
     return DiscreteVector<Tags...>((-get<Tags>(x))...);
 }
@@ -69,7 +72,7 @@ constexpr inline DiscreteVector<Tags...> operator-(DiscreteVector<Tags...> const
 /// Internal binary operators: +, -
 
 template <class... Tags, class... OTags>
-constexpr inline auto operator+(
+KOKKOS_FUNCTION constexpr auto operator+(
         DiscreteVector<Tags...> const& lhs,
         DiscreteVector<OTags...> const& rhs)
 {
@@ -88,7 +91,7 @@ constexpr inline auto operator+(
 }
 
 template <class... Tags, class... OTags>
-constexpr inline auto operator-(
+KOKKOS_FUNCTION constexpr auto operator-(
         DiscreteVector<Tags...> const& lhs,
         DiscreteVector<OTags...> const& rhs)
 {
@@ -107,7 +110,7 @@ constexpr inline auto operator-(
 }
 
 template <class Tag, class IntegralType, class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-constexpr inline DiscreteVector<Tag> operator+(
+KOKKOS_FUNCTION constexpr DiscreteVector<Tag> operator+(
         DiscreteVector<Tag> const& lhs,
         IntegralType const& rhs)
 {
@@ -115,7 +118,7 @@ constexpr inline DiscreteVector<Tag> operator+(
 }
 
 template <class IntegralType, class Tag, class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-constexpr inline DiscreteVector<Tag> operator+(
+KOKKOS_FUNCTION constexpr DiscreteVector<Tag> operator+(
         IntegralType const& lhs,
         DiscreteVector<Tag> const& rhs)
 {
@@ -123,7 +126,7 @@ constexpr inline DiscreteVector<Tag> operator+(
 }
 
 template <class Tag, class IntegralType, class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-constexpr inline DiscreteVector<Tag> operator-(
+KOKKOS_FUNCTION constexpr DiscreteVector<Tag> operator-(
         DiscreteVector<Tag> const& lhs,
         IntegralType const& rhs)
 {
@@ -131,7 +134,7 @@ constexpr inline DiscreteVector<Tag> operator-(
 }
 
 template <class IntegralType, class Tag, class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-constexpr inline DiscreteVector<Tag> operator-(
+KOKKOS_FUNCTION constexpr DiscreteVector<Tag> operator-(
         IntegralType const& lhs,
         DiscreteVector<Tag> const& rhs)
 {
@@ -144,25 +147,29 @@ template <
         class IntegralType,
         class... Tags,
         class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-constexpr inline auto operator*(IntegralType const& lhs, DiscreteVector<Tags...> const& rhs)
+KOKKOS_FUNCTION constexpr auto operator*(
+        IntegralType const& lhs,
+        DiscreteVector<Tags...> const& rhs)
 {
     return DiscreteVector<Tags...>((lhs * get<Tags>(rhs))...);
 }
 
 template <class... QueryTags, class... Tags>
-inline constexpr DiscreteVector<QueryTags...> select(DiscreteVector<Tags...> const& arr) noexcept
+KOKKOS_FUNCTION constexpr DiscreteVector<QueryTags...> select(
+        DiscreteVector<Tags...> const& arr) noexcept
 {
     return DiscreteVector<QueryTags...>(arr);
 }
 
 template <class... QueryTags, class... Tags>
-inline constexpr DiscreteVector<QueryTags...> select(DiscreteVector<Tags...>&& arr) noexcept
+KOKKOS_FUNCTION constexpr DiscreteVector<QueryTags...> select(
+        DiscreteVector<Tags...>&& arr) noexcept
 {
     return DiscreteVector<QueryTags...>(std::move(arr));
 }
 
 template <class QueryTag, class HeadTag, class... TailTags>
-constexpr DiscreteVector<QueryTag> const& take(
+KOKKOS_FUNCTION constexpr DiscreteVector<QueryTag> const& take(
         DiscreteVector<HeadTag> const& head,
         DiscreteVector<TailTags> const&... tags)
 {
@@ -188,12 +195,12 @@ template <class Tag>
 class ConversionOperators<DiscreteVector<Tag>>
 {
 public:
-    constexpr inline operator DiscreteVectorElement const &() const noexcept
+    KOKKOS_FUNCTION constexpr operator DiscreteVectorElement const &() const noexcept
     {
         return static_cast<DiscreteVector<Tag> const*>(this)->m_values[0];
     }
 
-    constexpr inline operator DiscreteVectorElement&() noexcept
+    KOKKOS_FUNCTION constexpr operator DiscreteVectorElement&() noexcept
     {
         return static_cast<DiscreteVector<Tag>*>(this)->m_values[0];
     }
@@ -203,7 +210,7 @@ namespace detail {
 
 /// Returns a reference to the underlying `std::array`
 template <class... Tags>
-constexpr inline std::array<DiscreteVectorElement, sizeof...(Tags)>& array(
+KOKKOS_FUNCTION constexpr std::array<DiscreteVectorElement, sizeof...(Tags)>& array(
         DiscreteVector<Tags...>& v) noexcept
 {
     return v.m_values;
@@ -211,7 +218,7 @@ constexpr inline std::array<DiscreteVectorElement, sizeof...(Tags)>& array(
 
 /// Returns a reference to the underlying `std::array`
 template <class... Tags>
-constexpr inline std::array<DiscreteVectorElement, sizeof...(Tags)> const& array(
+KOKKOS_FUNCTION constexpr std::array<DiscreteVectorElement, sizeof...(Tags)> const& array(
         DiscreteVector<Tags...> const& v) noexcept
 {
     return v.m_values;
@@ -228,10 +235,10 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
 {
     friend class ConversionOperators<DiscreteVector<Tags...>>;
 
-    friend constexpr std::array<DiscreteVectorElement, sizeof...(Tags)>& detail::array<Tags...>(
-            DiscreteVector<Tags...>& v) noexcept;
-    friend constexpr std::array<DiscreteVectorElement, sizeof...(Tags)> const& detail::array<
-            Tags...>(DiscreteVector<Tags...> const& v) noexcept;
+    friend KOKKOS_FUNCTION constexpr std::array<DiscreteVectorElement, sizeof...(Tags)>& detail::
+            array<Tags...>(DiscreteVector<Tags...>& v) noexcept;
+    friend KOKKOS_FUNCTION constexpr std::array<DiscreteVectorElement, sizeof...(Tags)> const&
+    detail::array<Tags...>(DiscreteVector<Tags...> const& v) noexcept;
 
     using tags_seq = detail::TypeSeq<Tags...>;
 
@@ -239,26 +246,28 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
     std::array<DiscreteVectorElement, sizeof...(Tags)> m_values;
 
 public:
-    static constexpr std::size_t size() noexcept
+    static KOKKOS_FUNCTION constexpr std::size_t size() noexcept
     {
         return sizeof...(Tags);
     }
 
 public:
-    inline constexpr DiscreteVector() = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteVector() = default;
 
-    inline constexpr DiscreteVector(DiscreteVector const&) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteVector(DiscreteVector const&) = default;
 
-    inline constexpr DiscreteVector(DiscreteVector&&) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteVector(DiscreteVector&&) = default;
 
     template <class... OTags>
-    explicit inline constexpr DiscreteVector(DiscreteVector<OTags> const&... other) noexcept
+    explicit KOKKOS_FUNCTION constexpr DiscreteVector(
+            DiscreteVector<OTags> const&... other) noexcept
         : m_values {take<Tags>(other...).value()...}
     {
     }
 
     template <class... OTags>
-    explicit inline constexpr DiscreteVector(DiscreteVector<OTags...> const& other) noexcept
+    explicit KOKKOS_FUNCTION constexpr DiscreteVector(
+            DiscreteVector<OTags...> const& other) noexcept
         : m_values {other.template get<Tags>()...}
     {
     }
@@ -268,43 +277,47 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
             class = std::enable_if_t<(std::is_convertible_v<Params, DiscreteVectorElement> && ...)>,
             class = std::enable_if_t<(!is_discrete_vector_v<Params> && ...)>,
             class = std::enable_if_t<sizeof...(Params) == sizeof...(Tags)>>
-    explicit inline constexpr DiscreteVector(Params const&... params) noexcept
+    explicit KOKKOS_FUNCTION constexpr DiscreteVector(Params const&... params) noexcept
         : m_values {static_cast<DiscreteVectorElement>(params)...}
     {
     }
 
-    constexpr inline DiscreteVector& operator=(DiscreteVector const& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION ~DiscreteVector() = default;
+
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteVector& operator=(DiscreteVector const& other)
+            = default;
 
-    constexpr inline DiscreteVector& operator=(DiscreteVector&& other) = default;
+    KOKKOS_DEFAULTED_FUNCTION constexpr DiscreteVector& operator=(DiscreteVector&& other) = default;
 
     template <class... OTags>
-    constexpr inline DiscreteVector& operator=(DiscreteVector<OTags...> const& other) noexcept
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator=(
+            DiscreteVector<OTags...> const& other) noexcept
     {
         m_values = other.m_values;
         return *this;
     }
 
     template <class... OTags>
-    constexpr inline DiscreteVector& operator=(DiscreteVector<OTags...>&& other) noexcept
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator=(DiscreteVector<OTags...>&& other) noexcept
     {
         m_values = std::move(other.m_values);
         return *this;
     }
 
     template <class... OTags>
-    constexpr inline bool operator==(DiscreteVector<OTags...> const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr bool operator==(DiscreteVector<OTags...> const& rhs) const noexcept
     {
         return ((m_values[type_seq_rank_v<Tags, tags_seq>] == rhs.template get<Tags>()) && ...);
     }
 
     template <class... OTags>
-    constexpr inline bool operator!=(DiscreteVector<OTags...> const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr bool operator!=(DiscreteVector<OTags...> const& rhs) const noexcept
     {
         return !(*this == rhs);
     }
 
     template <class QueryTag>
-    inline constexpr DiscreteVectorElement& get() noexcept
+    KOKKOS_FUNCTION constexpr DiscreteVectorElement& get() noexcept
     {
         using namespace detail;
         static_assert(in_tags_v<QueryTag, tags_seq>, "requested Tag absent from DiscreteVector");
@@ -312,7 +325,7 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
     }
 
     template <class QueryTag>
-    inline constexpr DiscreteVectorElement const& get() const noexcept
+    KOKKOS_FUNCTION constexpr DiscreteVectorElement const& get() const noexcept
     {
         using namespace detail;
         static_assert(in_tags_v<QueryTag, tags_seq>, "requested Tag absent from DiscreteVector");
@@ -320,7 +333,8 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
     }
 
     template <class QueryTag>
-    DiscreteVectorElement const& get_or(DiscreteVectorElement const& default_value) const&
+    KOKKOS_FUNCTION constexpr DiscreteVectorElement const& get_or(
+            DiscreteVectorElement const& default_value) const&
     {
         DDC_IF_NVCC_THEN_PUSH_AND_SUPPRESS(implicit_return_from_non_void_function)
         if constexpr (in_tags_v<QueryTag, tags_seq>) {
@@ -332,20 +346,21 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
     }
 
     template <std::size_t N = sizeof...(Tags)>
-    constexpr inline std::enable_if_t<N == 1, DiscreteVectorElement const&> value() const noexcept
+    KOKKOS_FUNCTION constexpr std::enable_if_t<N == 1, DiscreteVectorElement const&> value()
+            const noexcept
     {
         return m_values[0];
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteVector& operator++()
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator++()
     {
         ++m_values[0];
         return *this;
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteVector operator++(int)
+    KOKKOS_FUNCTION constexpr DiscreteVector operator++(int)
     {
         DiscreteVector const tmp = *this;
         ++m_values[0];
@@ -353,14 +368,14 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteVector& operator--()
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator--()
     {
         ++m_values[0];
         return *this;
     }
 
     template <std::size_t N = sizeof...(Tags), class = std::enable_if_t<N == 1>>
-    constexpr inline DiscreteVector operator--(int)
+    KOKKOS_FUNCTION constexpr DiscreteVector operator--(int)
     {
         DiscreteVector const tmp = *this;
         ++m_values[0];
@@ -368,7 +383,7 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
     }
 
     template <class... OTags>
-    constexpr inline DiscreteVector& operator+=(DiscreteVector<OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator+=(DiscreteVector<OTags...> const& rhs)
     {
         static_assert(((type_seq_contains_v<detail::TypeSeq<OTags>, tags_seq>)&&...));
         ((m_values[type_seq_rank_v<OTags, tags_seq>] += rhs.template get<OTags>()), ...);
@@ -380,14 +395,14 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
             std::size_t N = sizeof...(Tags),
             class = std::enable_if_t<N == 1>,
             class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-    constexpr inline DiscreteVector& operator+=(IntegralType const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator+=(IntegralType const& rhs)
     {
         m_values[0] += rhs;
         return *this;
     }
 
     template <class... OTags>
-    constexpr inline DiscreteVector& operator-=(DiscreteVector<OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator-=(DiscreteVector<OTags...> const& rhs)
     {
         static_assert(((type_seq_contains_v<detail::TypeSeq<OTags>, tags_seq>)&&...));
         ((m_values[type_seq_rank_v<OTags, tags_seq>] -= rhs.template get<OTags>()), ...);
@@ -399,14 +414,14 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
             std::size_t N = sizeof...(Tags),
             class = std::enable_if_t<N == 1>,
             class = std::enable_if_t<std::is_integral_v<IntegralType>>>
-    constexpr inline DiscreteVector& operator-=(IntegralType const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator-=(IntegralType const& rhs)
     {
         m_values[0] -= rhs;
         return *this;
     }
 
     template <class... OTags>
-    constexpr inline DiscreteVector& operator*=(DiscreteVector<OTags...> const& rhs)
+    KOKKOS_FUNCTION constexpr DiscreteVector& operator*=(DiscreteVector<OTags...> const& rhs)
     {
         static_assert(((type_seq_contains_v<detail::TypeSeq<OTags>, tags_seq>)&&...));
         ((m_values[type_seq_rank_v<OTags, tags_seq>] *= rhs.template get<OTags>()), ...);
@@ -415,13 +430,15 @@ class DiscreteVector : public ConversionOperators<DiscreteVector<Tags...>>
 };
 
 template <class Tag>
-constexpr inline bool operator<(DiscreteVector<Tag> const& lhs, DiscreteVector<Tag> const& rhs)
+KOKKOS_FUNCTION constexpr bool operator<(
+        DiscreteVector<Tag> const& lhs,
+        DiscreteVector<Tag> const& rhs)
 {
     return lhs.value() < rhs.value();
 }
 
 template <class Tag, class IntegralType>
-constexpr inline bool operator<(DiscreteVector<Tag> const& lhs, IntegralType const& rhs)
+KOKKOS_FUNCTION constexpr bool operator<(DiscreteVector<Tag> const& lhs, IntegralType const& rhs)
 {
     return lhs.value() < rhs;
 }
@@ -441,4 +458,5 @@ std::ostream& operator<<(std::ostream& out, DiscreteVector<Head, Tags...> const&
     out << ")";
     return out;
 }
+
 } // namespace ddc
diff --git a/include/ddc/dual_discretization.hpp b/include/ddc/dual_discretization.hpp
index 0ca16a1be..265d9afeb 100644
--- a/include/ddc/dual_discretization.hpp
+++ b/include/ddc/dual_discretization.hpp
@@ -48,7 +48,7 @@ class DualDiscretization
     }
 
     template <class MemorySpace>
-    DDC_INLINE_FUNCTION typename DDim::template Impl<MemorySpace> const& get()
+    KOKKOS_FORCEINLINE_FUNCTION typename DDim::template Impl<MemorySpace> const& get()
     {
         if constexpr (std::is_same_v<MemorySpace, Kokkos::HostSpace>) {
             return m_host;
@@ -67,12 +67,12 @@ class DualDiscretization
         }
     }
 
-    DDC_INLINE_FUNCTION DDimImplHost const& get_host()
+    KOKKOS_FORCEINLINE_FUNCTION DDimImplHost const& get_host()
     {
         return m_host;
     }
 
-    DDC_INLINE_FUNCTION DDimImplDevice const& get_device()
+    KOKKOS_FORCEINLINE_FUNCTION DDimImplDevice const& get_device()
     {
 #if defined(__CUDACC__) || defined(__HIPCC__)
         return m_device_on_host;
diff --git a/include/ddc/kernels/fft.hpp b/include/ddc/kernels/fft.hpp
index 9807bd550..5db1e2664 100644
--- a/include/ddc/kernels/fft.hpp
+++ b/include/ddc/kernels/fft.hpp
@@ -91,13 +91,13 @@ constexpr bool is_complex_v = is_complex<T>::value;
 
 // LastSelector: returns a if Dim==Last, else b
 template <typename T, typename Dim, typename Last>
-constexpr T LastSelector(const T a, const T b)
+KOKKOS_FUNCTION constexpr T LastSelector(const T a, const T b)
 {
     return std::is_same_v<Dim, Last> ? a : b;
 }
 
 template <typename T, typename Dim, typename First, typename Second, typename... Tail>
-constexpr T LastSelector(const T a, const T b)
+KOKKOS_FUNCTION constexpr T LastSelector(const T a, const T b)
 {
     return LastSelector<T, Dim, Second, Tail...>(a, b);
 }
@@ -401,10 +401,10 @@ void core(
     else if constexpr (std::is_same_v<ExecSpace, Kokkos::OpenMP>) {
         if constexpr (std::is_same_v<real_type_t<Tin>, float>) {
             fftwf_init_threads();
-            fftwf_plan_with_nthreads(ExecSpace::concurrency());
+            fftwf_plan_with_nthreads(execSpace.concurrency());
         } else {
             fftw_init_threads();
-            fftw_plan_with_nthreads(ExecSpace::concurrency());
+            fftw_plan_with_nthreads(execSpace.concurrency());
         }
         _fftw_plan<Tin> plan = _fftw_plan_many_dft<Tin, Tout>(
                 kwargs.direction == ddc::FFT_Direction::FORWARD ? FFTW_FORWARD : FFTW_BACKWARD,
diff --git a/include/ddc/non_uniform_point_sampling.hpp b/include/ddc/non_uniform_point_sampling.hpp
index 2f1a849ca..648bd2102 100644
--- a/include/ddc/non_uniform_point_sampling.hpp
+++ b/include/ddc/non_uniform_point_sampling.hpp
@@ -60,7 +60,7 @@ class NonUniformPointSampling
 
         /// @brief Construct a `NonUniformPointSampling` using a C++20 "common range".
         template <class InputRange>
-        explicit inline constexpr Impl(InputRange const& points)
+        explicit Impl(InputRange const& points)
         {
             if constexpr (Kokkos::is_view<InputRange>::value) {
                 Kokkos::deep_copy(m_points, points);
@@ -75,7 +75,7 @@ class NonUniformPointSampling
 
         /// @brief Construct a `NonUniformPointSampling` using a pair of iterators.
         template <class InputIt>
-        inline constexpr Impl(InputIt points_begin, InputIt points_end)
+        Impl(InputIt points_begin, InputIt points_end)
         {
             std::vector<continuous_element_type> host_points(points_begin, points_end);
             Kokkos::View<continuous_element_type*, Kokkos::HostSpace>
@@ -96,14 +96,14 @@ class NonUniformPointSampling
 
         ~Impl() = default;
 
-        constexpr std::size_t size() const
+        KOKKOS_FUNCTION std::size_t size() const
         {
             return m_points.size();
         }
 
         /// @brief Convert a mesh index into a position in `CDim`
-        constexpr continuous_element_type coordinate(
-                discrete_element_type const& icoord) const noexcept
+        KOKKOS_FUNCTION continuous_element_type
+        coordinate(discrete_element_type const& icoord) const noexcept
         {
             return m_points(icoord.uid());
         }
@@ -134,40 +134,37 @@ std::ostream& operator<<(std::ostream& out, DDimImpl const& mesh)
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> coordinate(
-        DiscreteElement<NonUniformPointSampling<CDim>> const& c)
+KOKKOS_FUNCTION Coordinate<CDim> coordinate(DiscreteElement<NonUniformPointSampling<CDim>> const& c)
 {
     return discrete_space<NonUniformPointSampling<CDim>>().coordinate(c);
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> distance_at_left(
-        DiscreteElement<NonUniformPointSampling<CDim>> i)
+KOKKOS_FUNCTION Coordinate<CDim> distance_at_left(DiscreteElement<NonUniformPointSampling<CDim>> i)
 {
     return coordinate(i) - coordinate(i - 1);
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> distance_at_right(
-        DiscreteElement<NonUniformPointSampling<CDim>> i)
+KOKKOS_FUNCTION Coordinate<CDim> distance_at_right(DiscreteElement<NonUniformPointSampling<CDim>> i)
 {
     return coordinate(i + 1) - coordinate(i);
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rmin(DiscreteDomain<NonUniformPointSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rmin(DiscreteDomain<NonUniformPointSampling<CDim>> const& d)
 {
     return coordinate(d.front());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rmax(DiscreteDomain<NonUniformPointSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rmax(DiscreteDomain<NonUniformPointSampling<CDim>> const& d)
 {
     return coordinate(d.back());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rlength(DiscreteDomain<NonUniformPointSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rlength(DiscreteDomain<NonUniformPointSampling<CDim>> const& d)
 {
     return rmax(d) - rmin(d);
 }
diff --git a/include/ddc/periodic_sampling.hpp b/include/ddc/periodic_sampling.hpp
index 0335159e7..42dc02ad7 100644
--- a/include/ddc/periodic_sampling.hpp
+++ b/include/ddc/periodic_sampling.hpp
@@ -73,7 +73,7 @@ class PeriodicSampling
          * @param step   the real distance between two points of mesh distance 1
          * @param n_period   the number of steps in a period
          */
-        constexpr Impl(continuous_element_type origin, Real step, std::size_t n_period)
+        Impl(continuous_element_type origin, Real step, std::size_t n_period)
             : m_origin(origin)
             , m_step(step)
             , m_n_period(n_period)
@@ -85,32 +85,32 @@ class PeriodicSampling
         ~Impl() = default;
 
         /// @brief Lower bound index of the mesh
-        constexpr continuous_element_type origin() const noexcept
+        KOKKOS_FUNCTION continuous_element_type origin() const noexcept
         {
             return m_origin;
         }
 
         /// @brief Lower bound index of the mesh
-        constexpr discrete_element_type front() const noexcept
+        KOKKOS_FUNCTION discrete_element_type front() const noexcept
         {
             return discrete_element_type {0};
         }
 
         /// @brief Spacing step of the mesh
-        constexpr Real step() const
+        KOKKOS_FUNCTION Real step() const
         {
             return m_step;
         }
 
         /// @brief Number of steps in a period
-        constexpr std::size_t n_period() const
+        KOKKOS_FUNCTION std::size_t n_period() const
         {
             return m_n_period;
         }
 
         /// @brief Convert a mesh index into a position in `CDim`
-        constexpr continuous_element_type coordinate(
-                discrete_element_type const& icoord) const noexcept
+        KOKKOS_FUNCTION continuous_element_type
+        coordinate(discrete_element_type const& icoord) const noexcept
         {
             return m_origin
                    + continuous_element_type(
@@ -243,7 +243,7 @@ std::ostream& operator<<(std::ostream& out, DDimImpl const& mesh)
 
 /// @brief Lower bound index of the mesh
 template <class DDim>
-DDC_INLINE_FUNCTION std::
+KOKKOS_FUNCTION std::
         enable_if_t<is_periodic_sampling_v<DDim>, typename DDim::continuous_element_type>
         origin() noexcept
 {
@@ -252,53 +252,51 @@ DDC_INLINE_FUNCTION std::
 
 /// @brief Lower bound index of the mesh
 template <class DDim>
-DDC_INLINE_FUNCTION std::
-        enable_if_t<is_periodic_sampling_v<DDim>, typename DDim::discrete_element_type>
-        front() noexcept
+KOKKOS_FUNCTION std::enable_if_t<is_periodic_sampling_v<DDim>, typename DDim::discrete_element_type>
+front() noexcept
 {
     return discrete_space<DDim>().front();
 }
 
 /// @brief Spacing step of the mesh
 template <class DDim>
-DDC_INLINE_FUNCTION std::enable_if_t<is_periodic_sampling_v<DDim>, Real> step() noexcept
+KOKKOS_FUNCTION std::enable_if_t<is_periodic_sampling_v<DDim>, Real> step() noexcept
 {
     return discrete_space<DDim>().step();
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION constexpr Coordinate<CDim> coordinate(
-        DiscreteElement<PeriodicSampling<CDim>> const& c)
+KOKKOS_FUNCTION Coordinate<CDim> coordinate(DiscreteElement<PeriodicSampling<CDim>> const& c)
 {
     return discrete_space<PeriodicSampling<CDim>>().coordinate(c);
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> distance_at_left(DiscreteElement<PeriodicSampling<CDim>>)
+KOKKOS_FUNCTION Coordinate<CDim> distance_at_left(DiscreteElement<PeriodicSampling<CDim>>)
 {
     return Coordinate<CDim>(step<PeriodicSampling<CDim>>());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> distance_at_right(DiscreteElement<PeriodicSampling<CDim>>)
+KOKKOS_FUNCTION Coordinate<CDim> distance_at_right(DiscreteElement<PeriodicSampling<CDim>>)
 {
     return Coordinate<CDim>(step<PeriodicSampling<CDim>>());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rmin(DiscreteDomain<PeriodicSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rmin(DiscreteDomain<PeriodicSampling<CDim>> const& d)
 {
     return coordinate(d.front());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rmax(DiscreteDomain<PeriodicSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rmax(DiscreteDomain<PeriodicSampling<CDim>> const& d)
 {
     return coordinate(d.back());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rlength(DiscreteDomain<PeriodicSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rlength(DiscreteDomain<PeriodicSampling<CDim>> const& d)
 {
     return rmax(d) - rmin(d);
 }
diff --git a/include/ddc/reducer.hpp b/include/ddc/reducer.hpp
index b70e9d998..820459e5f 100644
--- a/include/ddc/reducer.hpp
+++ b/include/ddc/reducer.hpp
@@ -3,6 +3,8 @@
 #include <limits>
 #include <utility>
 
+#include <Kokkos_Core.hpp>
+
 namespace ddc::reducer {
 
 template <class T>
@@ -10,7 +12,8 @@ struct sum
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs + rhs;
     }
@@ -21,7 +24,8 @@ struct prod
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs * rhs;
     }
@@ -31,7 +35,8 @@ struct land
 {
     using value_type = bool;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs && rhs;
     }
@@ -41,7 +46,8 @@ struct lor
 {
     using value_type = bool;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs || rhs;
     }
@@ -52,7 +58,8 @@ struct band
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs & rhs;
     }
@@ -63,7 +70,8 @@ struct bor
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs | rhs;
     }
@@ -74,7 +82,8 @@ struct bxor
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs ^ rhs;
     }
@@ -85,7 +94,8 @@ struct min
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs < rhs ? lhs : rhs;
     }
@@ -96,7 +106,8 @@ struct max
 {
     using value_type = T;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return lhs > rhs ? lhs : rhs;
     }
@@ -107,7 +118,8 @@ struct minmax
 {
     using value_type = std::pair<T, T>;
 
-    constexpr value_type operator()(value_type const& lhs, value_type const& rhs) const noexcept
+    KOKKOS_FUNCTION constexpr value_type operator()(value_type const& lhs, value_type const& rhs)
+            const noexcept
     {
         return value_type(
                 lhs.first < rhs.first ? lhs.first : rhs.first,
diff --git a/include/ddc/transform_reduce.hpp b/include/ddc/transform_reduce.hpp
index aa9708126..3b0d358cd 100644
--- a/include/ddc/transform_reduce.hpp
+++ b/include/ddc/transform_reduce.hpp
@@ -133,7 +133,6 @@ class TransformReducerKokkosLambdaAdapter
     Functor functor;
 
 public:
-    DDC_FORCEINLINE_FUNCTION
     TransformReducerKokkosLambdaAdapter(Reducer const& r, Functor const& f) : reducer(r), functor(f)
     {
     }
diff --git a/include/ddc/uniform_point_sampling.hpp b/include/ddc/uniform_point_sampling.hpp
index 514ae1691..b5f21012f 100644
--- a/include/ddc/uniform_point_sampling.hpp
+++ b/include/ddc/uniform_point_sampling.hpp
@@ -69,7 +69,7 @@ class UniformPointSampling
          * @param origin the real coordinate of mesh coordinate 0
          * @param step   the real distance between two points of mesh distance 1
          */
-        constexpr Impl(continuous_element_type origin, Real step) : m_origin(origin), m_step(step)
+        Impl(continuous_element_type origin, Real step) : m_origin(origin), m_step(step)
         {
             assert(step > 0);
         }
@@ -77,26 +77,26 @@ class UniformPointSampling
         ~Impl() = default;
 
         /// @brief Lower bound index of the mesh
-        constexpr continuous_element_type origin() const noexcept
+        KOKKOS_FUNCTION continuous_element_type origin() const noexcept
         {
             return m_origin;
         }
 
         /// @brief Lower bound index of the mesh
-        constexpr discrete_element_type front() const noexcept
+        KOKKOS_FUNCTION discrete_element_type front() const noexcept
         {
             return discrete_element_type {0};
         }
 
         /// @brief Spacing step of the mesh
-        constexpr Real step() const
+        KOKKOS_FUNCTION Real step() const
         {
             return m_step;
         }
 
         /// @brief Convert a mesh index into a position in `CDim`
-        constexpr continuous_element_type coordinate(
-                discrete_element_type const& icoord) const noexcept
+        KOKKOS_FUNCTION continuous_element_type
+        coordinate(discrete_element_type const& icoord) const noexcept
         {
             return m_origin + continuous_element_type(icoord.uid()) * m_step;
         }
@@ -215,15 +215,14 @@ std::ostream& operator<<(std::ostream& out, DDimImpl const& mesh)
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> coordinate(
-        DiscreteElement<UniformPointSampling<CDim>> const& c)
+KOKKOS_FUNCTION Coordinate<CDim> coordinate(DiscreteElement<UniformPointSampling<CDim>> const& c)
 {
     return discrete_space<UniformPointSampling<CDim>>().coordinate(c);
 }
 
 /// @brief Lower bound index of the mesh
 template <class DDim>
-DDC_INLINE_FUNCTION std::
+KOKKOS_FUNCTION std::
         enable_if_t<is_uniform_sampling_v<DDim>, typename DDim::continuous_element_type>
         origin() noexcept
 {
@@ -232,46 +231,45 @@ DDC_INLINE_FUNCTION std::
 
 /// @brief Lower bound index of the mesh
 template <class DDim>
-DDC_INLINE_FUNCTION std::
-        enable_if_t<is_uniform_sampling_v<DDim>, typename DDim::discrete_element_type>
-        front() noexcept
+KOKKOS_FUNCTION std::enable_if_t<is_uniform_sampling_v<DDim>, typename DDim::discrete_element_type>
+front() noexcept
 {
     return discrete_space<DDim>().front();
 }
 
 /// @brief Spacing step of the mesh
 template <class DDim>
-DDC_INLINE_FUNCTION std::enable_if_t<is_uniform_sampling_v<DDim>, Real> step() noexcept
+KOKKOS_FUNCTION std::enable_if_t<is_uniform_sampling_v<DDim>, Real> step() noexcept
 {
     return discrete_space<DDim>().step();
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> distance_at_left(DiscreteElement<UniformPointSampling<CDim>>)
+KOKKOS_FUNCTION Coordinate<CDim> distance_at_left(DiscreteElement<UniformPointSampling<CDim>>)
 {
     return Coordinate<CDim>(step<UniformPointSampling<CDim>>());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> distance_at_right(DiscreteElement<UniformPointSampling<CDim>>)
+KOKKOS_FUNCTION Coordinate<CDim> distance_at_right(DiscreteElement<UniformPointSampling<CDim>>)
 {
     return Coordinate<CDim>(step<UniformPointSampling<CDim>>());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rmin(DiscreteDomain<UniformPointSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rmin(DiscreteDomain<UniformPointSampling<CDim>> const& d)
 {
     return coordinate(d.front());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rmax(DiscreteDomain<UniformPointSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rmax(DiscreteDomain<UniformPointSampling<CDim>> const& d)
 {
     return coordinate(d.back());
 }
 
 template <class CDim>
-DDC_INLINE_FUNCTION Coordinate<CDim> rlength(DiscreteDomain<UniformPointSampling<CDim>> const& d)
+KOKKOS_FUNCTION Coordinate<CDim> rlength(DiscreteDomain<UniformPointSampling<CDim>> const& d)
 {
     return rmax(d) - rmin(d);
 }
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 004468464..fb09849f9 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -1,5 +1,8 @@
 # SPDX-License-Identifier: MIT
 cmake_minimum_required(VERSION 3.15)
+
+add_subdirectory(discrete_space)
+
 add_executable(ddc_tests
     main.cpp
     aligned_allocator.cpp
diff --git a/tests/discrete_space/CMakeLists.txt b/tests/discrete_space/CMakeLists.txt
new file mode 100644
index 000000000..c9b0a6049
--- /dev/null
+++ b/tests/discrete_space/CMakeLists.txt
@@ -0,0 +1,16 @@
+add_library(discrete_space_tests_lib STATIC discrete_space.cpp)
+target_link_libraries(discrete_space_tests_lib
+    PUBLIC
+    GTest::gtest
+    DDC::DDC
+)
+
+add_executable(discrete_space_tests main.cpp)
+target_link_libraries(discrete_space_tests
+    PUBLIC
+        discrete_space_tests_lib
+        GTest::gtest
+        DDC::DDC
+)
+
+gtest_discover_tests(discrete_space_tests)
diff --git a/tests/discrete_space/discrete_space.cpp b/tests/discrete_space/discrete_space.cpp
new file mode 100644
index 000000000..a30e4c5b9
--- /dev/null
+++ b/tests/discrete_space/discrete_space.cpp
@@ -0,0 +1,21 @@
+// SPDX-License-Identifier: MIT
+
+#include <ddc/ddc.hpp>
+
+#include <gtest/gtest.h>
+
+#include "discrete_space.hpp"
+
+struct DimX;
+using DDimX = ddc::UniformPointSampling<DimX>;
+
+void do_not_optimize_away_discrete_space_tests() {}
+
+TEST(DiscreteSpace, Initialization)
+{
+    // ddc::init_discrete_space<QMeshX>(std::vector<ddc::Coordinate<QDimX>>());
+    ddc::init_discrete_space(DDimX::
+                                     init(ddc::Coordinate<DimX>(0),
+                                          ddc::Coordinate<DimX>(1),
+                                          ddc::DiscreteVector<DDimX>(2)));
+}
diff --git a/tests/discrete_space/discrete_space.hpp b/tests/discrete_space/discrete_space.hpp
new file mode 100644
index 000000000..26f61aece
--- /dev/null
+++ b/tests/discrete_space/discrete_space.hpp
@@ -0,0 +1,5 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+void do_not_optimize_away_discrete_space_tests();
diff --git a/tests/discrete_space/main.cpp b/tests/discrete_space/main.cpp
new file mode 100644
index 000000000..379b55161
--- /dev/null
+++ b/tests/discrete_space/main.cpp
@@ -0,0 +1,26 @@
+// SPDX-License-Identifier: MIT
+
+#include <ddc/ddc.hpp>
+
+#include <gtest/gtest.h>
+
+#include "discrete_space.hpp"
+
+namespace testing::internal {
+// accessing gtest internals is not very clean, but gtest provides no public access...
+extern bool g_help_flag;
+} // namespace testing::internal
+
+int main(int argc, char** argv)
+{
+    ::testing::InitGoogleTest(&argc, argv);
+    do_not_optimize_away_discrete_space_tests();
+    if (::testing::GTEST_FLAG(list_tests) || ::testing::internal::g_help_flag) {
+        // do not initialize DDC just to list tests or show help so as to be able to run
+        // a Cuda/Hip/etc. enabled DDC with no device available
+        return RUN_ALL_TESTS();
+    } else {
+        ddc::ScopeGuard scope(argc, argv);
+        return RUN_ALL_TESTS();
+    }
+}