]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Use DEAL_II_HOST_DEVICE in CUDA matrix-free framework
authorBruno Turcksin <bruno.turcksin@gmail.com>
Fri, 2 Dec 2022 20:35:32 +0000 (15:35 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Thu, 6 Apr 2023 13:10:39 +0000 (13:10 +0000)
include/deal.II/matrix_free/cuda_fe_evaluation.h
include/deal.II/matrix_free/cuda_hanging_nodes_internal.h
include/deal.II/matrix_free/cuda_matrix_free.h
include/deal.II/matrix_free/cuda_matrix_free.templates.h
include/deal.II/matrix_free/cuda_tensor_product_kernels.h
tests/cuda/coefficient_eval.cc
tests/cuda/matrix_free_no_index_initialize.cc
tests/cuda/matrix_vector_mf.h

index 39a6214f0f3cd5251d567336f9800946428fc58d..fd603d9ad9b24378ae9f792ab89712f8ec33f9a9 100644 (file)
@@ -45,14 +45,17 @@ namespace CUDAWrappers
      * number of points in each space dimensions.
      */
     template <int dim, int n_points_1d>
-    __device__ inline unsigned int
+    DEAL_II_HOST_DEVICE inline unsigned int
     compute_index()
     {
-      return (dim == 1 ? threadIdx.x % n_points_1d :
-              dim == 2 ? threadIdx.x % n_points_1d + n_points_1d * threadIdx.y :
-                         threadIdx.x % n_points_1d +
-                           n_points_1d *
-                             (threadIdx.y + n_points_1d * threadIdx.z));
+      KOKKOS_IF_ON_DEVICE(
+        return (dim == 1 ?
+                  threadIdx.x % n_points_1d :
+                dim == 2 ?
+                  threadIdx.x % n_points_1d + n_points_1d * threadIdx.y :
+                  threadIdx.x % n_points_1d +
+                    n_points_1d * (threadIdx.y + n_points_1d * threadIdx.z));)
+      KOKKOS_IF_ON_HOST(return 0;)
     }
   } // namespace internal
 
@@ -129,7 +132,7 @@ namespace CUDAWrappers
     /**
      * Constructor.
      */
-    __device__
+    DEAL_II_HOST_DEVICE
     FEEvaluation(const unsigned int       cell_id,
                  const data_type *        data,
                  SharedData<dim, Number> *shdata);
@@ -142,7 +145,7 @@ namespace CUDAWrappers
      * nodes, so once can see it as a similar function to
      * AffineConstraints::read_dof_valuess as well.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     read_dof_values(const Number *src);
 
     /**
@@ -151,7 +154,7 @@ namespace CUDAWrappers
      * during the write operation. The functionality is hence similar to the
      * function AffineConstraints::distribute_local_to_global.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     distribute_local_to_global(Number *dst) const;
 
     /**
@@ -161,7 +164,7 @@ namespace CUDAWrappers
      * computed. This function needs to be called before the functions
      * @p get_value() or @p get_gradient() give useful information.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     evaluate(const bool evaluate_val, const bool evaluate_grad);
 
     /**
@@ -171,49 +174,49 @@ namespace CUDAWrappers
      * @p integrate_val and @p integrate_grad are used to enable/disable some
      * of the values or the gradients.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     integrate(const bool integrate_val, const bool integrate_grad);
 
     /**
      * Same as above, except that the quadrature point is computed from thread
      * id.
      */
-    __device__ value_type
+    DEAL_II_HOST_DEVICE value_type
     get_value() const;
 
     /**
      * Same as above, except that the local dof index is computed from the
      * thread id.
      */
-    __device__ value_type
+    DEAL_II_HOST_DEVICE value_type
     get_dof_value() const;
 
     /**
      * Same as above, except that the quadrature point is computed from the
      * thread id.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     submit_value(const value_type &val_in);
 
     /**
      * Same as above, except that the local dof index is computed from the
      * thread id.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     submit_dof_value(const value_type &val_in);
 
     /**
      * Same as above, except that the quadrature point is computed from the
      * thread id.
      */
-    __device__ gradient_type
+    DEAL_II_HOST_DEVICE gradient_type
     get_gradient() const;
 
     /**
      * Same as above, except that the quadrature point is computed from the
      * thread id.
      */
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     submit_gradient(const gradient_type &grad_in);
 
     // clang-format off
@@ -223,13 +226,13 @@ namespace CUDAWrappers
      *
      * @p func needs to define
      * \code
-     * __device__ void operator()(
+     * DEAL_II_HOST_DEVICE void operator()(
      *   CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> *fe_eval) const;
      * \endcode
      */
     // clang-format on
     template <typename Functor>
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     apply_for_each_quad_point(const Functor &func);
 
   private:
@@ -258,7 +261,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__
+  DEAL_II_HOST_DEVICE
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     FEEvaluation(const unsigned int       cell_id,
                  const data_type *        data,
@@ -285,7 +288,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     read_dof_values(const Number *src)
   {
@@ -295,9 +298,8 @@ namespace CUDAWrappers
 
     const types::global_dof_index src_idx = local_to_global[idx];
     // Use the read-only data cache.
-    values[idx] = __ldg(&src[src_idx]);
-
-    __syncthreads();
+    KOKKOS_IF_ON_DEVICE(values[idx] = __ldg(&src[src_idx]); __syncthreads();)
+    KOKKOS_IF_ON_HOST(values[idx] = src[src_idx];)
 
     internal::resolve_hanging_nodes<dim, fe_degree, false>(constraint_mask,
                                                            values);
@@ -310,7 +312,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     distribute_local_to_global(Number *dst) const
   {
@@ -336,7 +338,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::evaluate(
     const bool evaluate_val,
     const bool evaluate_grad)
@@ -354,17 +356,17 @@ namespace CUDAWrappers
       {
         evaluator_tensor_product.value_and_gradient_at_quad_pts(values,
                                                                 gradients);
-        __syncthreads();
+        KOKKOS_IF_ON_DEVICE(__syncthreads();)
       }
     else if (evaluate_grad == true)
       {
         evaluator_tensor_product.gradient_at_quad_pts(values, gradients);
-        __syncthreads();
+        KOKKOS_IF_ON_DEVICE(__syncthreads();)
       }
     else if (evaluate_val == true)
       {
         evaluator_tensor_product.value_at_quad_pts(values);
-        __syncthreads();
+        KOKKOS_IF_ON_DEVICE(__syncthreads();)
       }
   }
 
@@ -375,7 +377,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::integrate(
     const bool integrate_val,
     const bool integrate_grad)
@@ -395,13 +397,13 @@ namespace CUDAWrappers
     else if (integrate_val == true)
       {
         evaluator_tensor_product.integrate_value(values);
-        __syncthreads();
+        KOKKOS_IF_ON_DEVICE(__syncthreads();)
       }
     else if (integrate_grad == true)
       {
         evaluator_tensor_product.template integrate_gradient<false>(values,
                                                                     gradients);
-        __syncthreads();
+        KOKKOS_IF_ON_DEVICE(__syncthreads();)
       }
   }
 
@@ -412,11 +414,11 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ typename FEEvaluation<dim,
-                                   fe_degree,
-                                   n_q_points_1d,
-                                   n_components_,
-                                   Number>::value_type
+  DEAL_II_HOST_DEVICE typename FEEvaluation<dim,
+                                            fe_degree,
+                                            n_q_points_1d,
+                                            n_components_,
+                                            Number>::value_type
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     get_value() const
   {
@@ -431,11 +433,11 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ typename FEEvaluation<dim,
-                                   fe_degree,
-                                   n_q_points_1d,
-                                   n_components_,
-                                   Number>::value_type
+  DEAL_II_HOST_DEVICE typename FEEvaluation<dim,
+                                            fe_degree,
+                                            n_q_points_1d,
+                                            n_components_,
+                                            Number>::value_type
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     get_dof_value() const
   {
@@ -450,7 +452,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_value(const value_type &val_in)
   {
@@ -465,7 +467,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_dof_value(const value_type &val_in)
   {
@@ -480,11 +482,11 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ typename FEEvaluation<dim,
-                                   fe_degree,
-                                   n_q_points_1d,
-                                   n_components_,
-                                   Number>::gradient_type
+  DEAL_II_HOST_DEVICE typename FEEvaluation<dim,
+                                            fe_degree,
+                                            n_q_points_1d,
+                                            n_components_,
+                                            Number>::gradient_type
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     get_gradient() const
   {
@@ -514,7 +516,7 @@ namespace CUDAWrappers
             int n_q_points_1d,
             int n_components_,
             typename Number>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_gradient(const gradient_type &grad_in)
   {
@@ -539,13 +541,13 @@ namespace CUDAWrappers
             int n_components_,
             typename Number>
   template <typename Functor>
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     apply_for_each_quad_point(const Functor &func)
   {
     func(this);
 
-    __syncthreads();
+    KOKKOS_IF_ON_DEVICE(__syncthreads();)
   }
 } // namespace CUDAWrappers
 
index e15070c511d9e68dd5e7e9baa10c4a89c2d8876a..c354e0b2b716f58b3806a119dc28e495b7101345 100644 (file)
@@ -24,6 +24,8 @@
 
 #  include <deal.II/matrix_free/hanging_nodes_internal.h>
 
+#  include <Kokkos_Macros.hpp>
+
 DEAL_II_NAMESPACE_OPEN
 namespace CUDAWrappers
 {
@@ -37,7 +39,7 @@ namespace CUDAWrappers
     // Functions for resolving the hanging node constraints on the GPU        //
     //------------------------------------------------------------------------//
     template <unsigned int size>
-    __device__ inline unsigned int
+    DEAL_II_HOST_DEVICE inline unsigned int
     index2(unsigned int i, unsigned int j)
     {
       return i + size * j;
@@ -46,7 +48,7 @@ namespace CUDAWrappers
 
 
     template <unsigned int size>
-    __device__ inline unsigned int
+    DEAL_II_HOST_DEVICE inline unsigned int
     index3(unsigned int i, unsigned int j, unsigned int k)
     {
       return i + size * j + size * size * k;
@@ -58,7 +60,7 @@ namespace CUDAWrappers
               unsigned int direction,
               bool         transpose,
               typename Number>
-    __device__ inline void
+    DEAL_II_HOST_DEVICE inline void
     interpolate_boundary_2d(
       const dealii::internal::MatrixFreeFunctions::ConstraintKinds
               constraint_mask,
@@ -151,11 +153,11 @@ namespace CUDAWrappers
 
       // The synchronization is done for all the threads in one block with
       // each block being assigned to one element.
-      __syncthreads();
+      KOKKOS_IF_ON_DEVICE(__syncthreads();)
       if (constrained_face && constrained_dof)
         values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
 
-      __syncthreads();
+      KOKKOS_IF_ON_DEVICE(__syncthreads();)
     }
 
 
@@ -164,7 +166,7 @@ namespace CUDAWrappers
               unsigned int direction,
               bool         transpose,
               typename Number>
-    __device__ inline void
+    DEAL_II_HOST_DEVICE inline void
     interpolate_boundary_3d(
       const dealii::internal::MatrixFreeFunctions::ConstraintKinds
               constraint_mask,
@@ -293,14 +295,14 @@ namespace CUDAWrappers
 
       // The synchronization is done for all the threads in one block with
       // each block being assigned to one element.
-      __syncthreads();
+      KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
       if ((constrained_face != dealii::internal::MatrixFreeFunctions::
                                  ConstraintKinds::unconstrained) &&
           constrained_dof)
         values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
 
-      __syncthreads();
+      KOKKOS_IF_ON_DEVICE(__syncthreads();)
     }
 
 
@@ -313,7 +315,7 @@ namespace CUDAWrappers
      * @cite kronbichler2019multigrid.
      */
     template <int dim, int fe_degree, bool transpose, typename Number>
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     resolve_hanging_nodes(
       const dealii::internal::MatrixFreeFunctions::ConstraintKinds
               constraint_mask,
index cf2c0a2c4c15d6fd21e01bf296c8ecad2539010e..fca3ad059bf92ce6c5875b387c3061179ecce976 100644 (file)
@@ -38,6 +38,8 @@
 #  include <deal.II/lac/cuda_vector.h>
 #  include <deal.II/lac/la_parallel_vector.h>
 
+#  include <Kokkos_Core.hpp>
+
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -296,7 +298,7 @@ namespace CUDAWrappers
      *
      * @p func needs to define
      * \code
-     * __device__ void operator()(
+     * DEAL_II_HOST_DEVICE void operator()(
      *   const unsigned int                                          cell,
      *   const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data,
      *   CUDAWrappers::SharedData<dim, Number> *                     shared_data,
@@ -321,7 +323,7 @@ namespace CUDAWrappers
      *
      * @p func needs to define
      * \code
-     *  __device__ void operator()(
+     *  DEAL_II_HOST_DEVICE void operator()(
      *    const unsigned int                                          cell,
      *    const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data);
      * static const unsigned int n_dofs_1d;
@@ -674,7 +676,7 @@ namespace CUDAWrappers
     /**
      * Constructor.
      */
-    __device__
+    DEAL_II_HOST_DEVICE
     SharedData(Number *vd, Number *gq[dim])
       : values(vd)
     {
@@ -700,7 +702,7 @@ namespace CUDAWrappers
   // This function determines the number of cells per block, possibly at compile
   // time (by virtue of being 'constexpr')
   // TODO this function should be rewritten using meta-programming
-  __host__ __device__ constexpr unsigned int
+  DEAL_II_HOST_DEVICE constexpr unsigned int
   cells_per_block_shmem(int dim, int fe_degree)
   {
     /* clang-format off */
@@ -727,14 +729,18 @@ namespace CUDAWrappers
    * @relates CUDAWrappers::MatrixFree
    */
   template <int dim>
-  __device__ inline unsigned int
+  DEAL_II_HOST_DEVICE inline unsigned int
   q_point_id_in_cell(const unsigned int n_q_points_1d)
   {
-    return (
-      dim == 1 ? threadIdx.x % n_q_points_1d :
-      dim == 2 ? threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
-                 threadIdx.x % n_q_points_1d +
-                   n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));
+    KOKKOS_IF_ON_DEVICE(
+      return (dim == 1 ?
+                threadIdx.x % n_q_points_1d :
+              dim == 2 ?
+                threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y :
+                threadIdx.x % n_q_points_1d +
+                  n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z));)
+
+    KOKKOS_IF_ON_HOST(AssertThrow(false, ExcInternalError()); return 0;)
   }
 
 
@@ -746,7 +752,7 @@ namespace CUDAWrappers
    * @relates CUDAWrappers::MatrixFree
    */
   template <int dim, typename Number>
-  __device__ inline unsigned int
+  DEAL_II_HOST_DEVICE inline unsigned int
   local_q_point_id(
     const unsigned int                                          cell,
     const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
@@ -765,11 +771,12 @@ namespace CUDAWrappers
    * @relates CUDAWrappers::MatrixFree
    */
   template <int dim, typename Number>
-  __device__ inline typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
-  get_quadrature_point(
-    const unsigned int                                          cell,
-    const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
-    const unsigned int                                          n_q_points_1d)
+  DEAL_II_HOST_DEVICE inline
+    typename CUDAWrappers::MatrixFree<dim, Number>::point_type &
+    get_quadrature_point(
+      const unsigned int                                          cell,
+      const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
+      const unsigned int                                          n_q_points_1d)
   {
     return *(data->q_points + data->padding_length * cell +
              q_point_id_in_cell<dim>(n_q_points_1d));
index 1fde44beb053bf5c32f87248f741617083d1e812..7ab5b105f89e502896a270318eaa802bd0b15d0a 100644 (file)
@@ -75,56 +75,56 @@ namespace CUDAWrappers
                                                   [data_array_size];
 
     template <typename Number>
-    __host__ __device__ inline DataArray<Number> &
-             get_global_shape_values(unsigned int i);
+    DEAL_II_HOST_DEVICE inline DataArray<Number> &
+    get_global_shape_values(unsigned int i);
 
     template <>
-    __host__ __device__ inline DataArray<double> &
-             get_global_shape_values<double>(unsigned int i)
+    DEAL_II_HOST_DEVICE inline DataArray<double> &
+    get_global_shape_values<double>(unsigned int i)
     {
       return global_shape_values_d[i];
     }
 
     template <>
-    __host__ __device__ inline DataArray<float> &
-             get_global_shape_values<float>(unsigned int i)
+    DEAL_II_HOST_DEVICE inline DataArray<float> &
+    get_global_shape_values<float>(unsigned int i)
     {
       return global_shape_values_f[i];
     }
 
     template <typename Number>
-    __host__ __device__ inline DataArray<Number> &
-             get_global_shape_gradients(unsigned int i);
+    DEAL_II_HOST_DEVICE inline DataArray<Number> &
+    get_global_shape_gradients(unsigned int i);
 
     template <>
-    __host__ __device__ inline DataArray<double> &
-             get_global_shape_gradients<double>(unsigned int i)
+    DEAL_II_HOST_DEVICE inline DataArray<double> &
+    get_global_shape_gradients<double>(unsigned int i)
     {
       return global_shape_gradients_d[i];
     }
 
     template <>
-    __host__ __device__ inline DataArray<float> &
-             get_global_shape_gradients<float>(unsigned int i)
+    DEAL_II_HOST_DEVICE inline DataArray<float> &
+    get_global_shape_gradients<float>(unsigned int i)
     {
       return global_shape_gradients_f[i];
     }
 
     // for collocation methods
     template <typename Number>
-    __host__ __device__ inline DataArray<Number> &
-             get_global_co_shape_gradients(unsigned int i);
+    DEAL_II_HOST_DEVICE inline DataArray<Number> &
+    get_global_co_shape_gradients(unsigned int i);
 
     template <>
-    __host__ __device__ inline DataArray<double> &
-             get_global_co_shape_gradients<double>(unsigned int i)
+    DEAL_II_HOST_DEVICE inline DataArray<double> &
+    get_global_co_shape_gradients<double>(unsigned int i)
     {
       return global_co_shape_gradients_d[i];
     }
 
     template <>
-    __host__ __device__ inline DataArray<float> &
-             get_global_co_shape_gradients<float>(unsigned int i)
+    DEAL_II_HOST_DEVICE inline DataArray<float> &
+    get_global_co_shape_gradients<float>(unsigned int i)
     {
       return global_co_shape_gradients_f[i];
     }
index c423b48dd608e7fa45ae69356eea2d27d14281e0..9a93c939d4c1517a7b6bcbd57b63b28c813a2d15 100644 (file)
@@ -83,7 +83,7 @@ namespace CUDAWrappers
       static constexpr unsigned int n_q_points =
         Utilities::pow(n_q_points_1d, dim);
 
-      __device__
+      DEAL_II_HOST_DEVICE
       EvaluatorTensorProduct(int mf_object_id);
 
       /**
@@ -91,7 +91,7 @@ namespace CUDAWrappers
        * points.
        */
       template <int direction, bool dof_to_quad, bool add, bool in_place>
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       values(Number shape_values[], const Number *in, Number *out) const;
 
       /**
@@ -99,40 +99,40 @@ namespace CUDAWrappers
        * points for a given @p direction.
        */
       template <int direction, bool dof_to_quad, bool add, bool in_place>
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       gradients(Number shape_gradients[], const Number *in, Number *out) const;
 
       /**
        * Helper function for values() and gradients().
        */
       template <int direction, bool dof_to_quad, bool add, bool in_place>
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       apply(Number shape_data[], const Number *in, Number *out) const;
 
       /**
        * Evaluate the finite element function at the quadrature points.
        */
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       value_at_quad_pts(Number *u);
 
       /**
        * Helper function for integrate(). Integrate the finite element function.
        */
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       integrate_value(Number *u);
 
       /**
        * Evaluate the gradients of the finite element function at the quadrature
        * points.
        */
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       gradient_at_quad_pts(const Number *const u, Number *grad_u[dim]);
 
       /**
        * Evaluate the values and the gradients of the finite element function at
        * the quadrature points.
        */
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       value_and_gradient_at_quad_pts(Number *const u, Number *grad_u[dim]);
 
       /**
@@ -140,14 +140,14 @@ namespace CUDAWrappers
        * element function.
        */
       template <bool add>
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       integrate_gradient(Number *u, Number *grad_u[dim]);
 
       /**
        * Helper function for integrate(). Integrate the values and the gradients
        * of the finite element function.
        */
-      __device__ void
+      DEAL_II_HOST_DEVICE void
       integrate_value_and_gradient(Number *u, Number *grad_u[dim]);
 
       const int mf_object_id;
@@ -156,7 +156,7 @@ namespace CUDAWrappers
 
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    __device__
+    DEAL_II_HOST_DEVICE
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -169,7 +169,7 @@ namespace CUDAWrappers
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
     template <int direction, bool dof_to_quad, bool add, bool in_place>
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -185,7 +185,7 @@ namespace CUDAWrappers
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
     template <int direction, bool dof_to_quad, bool add, bool in_place>
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -201,7 +201,7 @@ namespace CUDAWrappers
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
     template <int direction, bool dof_to_quad, bool add, bool in_place>
-    __device__ void
+    DEAL_II_HOST_DEVICE void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -210,17 +210,17 @@ namespace CUDAWrappers
                                           const Number *in,
                                           Number *      out) const
     {
-      const unsigned int i = (dim == 1) ? 0 : threadIdx.x % n_q_points_1d;
-      const unsigned int j = (dim == 3) ? threadIdx.y : 0;
-      const unsigned int q = (dim == 1) ? (threadIdx.x % n_q_points_1d) :
-                             (dim == 2) ? threadIdx.y :
-                                          threadIdx.z;
-
-      // This loop simply multiply the shape function at the quadrature point by
-      // the value finite element coefficient.
-      Number t = 0;
-      for (int k = 0; k < n_q_points_1d; ++k)
-        {
+      KOKKOS_IF_ON_DEVICE(
+        const unsigned int i = (dim == 1) ? 0 : threadIdx.x % n_q_points_1d;
+        const unsigned int j = (dim == 3) ? threadIdx.y : 0;
+        const unsigned int q = (dim == 1) ? (threadIdx.x % n_q_points_1d) :
+                               (dim == 2) ? threadIdx.y :
+                                            threadIdx.z;
+
+        // This loop simply multiply the shape function at the quadrature point
+        // by the value finite element coefficient.
+        Number t = 0;
+        for (int k = 0; k < n_q_points_1d; ++k) {
           const unsigned int shape_idx =
             dof_to_quad ? (q + k * n_q_points_1d) : (k + q * n_q_points_1d);
           const unsigned int source_idx =
@@ -231,24 +231,21 @@ namespace CUDAWrappers
                (in_place ? out[source_idx] : in[source_idx]);
         }
 
-      if (in_place)
-        __syncthreads();
+        if (in_place) __syncthreads();
 
-      const unsigned int destination_idx =
-        (direction == 0) ? (q + n_q_points_1d * (i + n_q_points_1d * j)) :
-        (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
-                           (i + n_q_points_1d * (j + n_q_points_1d * q));
+        const unsigned int destination_idx =
+          (direction == 0) ? (q + n_q_points_1d * (i + n_q_points_1d * j)) :
+          (direction == 1) ? (i + n_q_points_1d * (q + n_q_points_1d * j)) :
+                             (i + n_q_points_1d * (j + n_q_points_1d * q));
 
-      if (add)
-        out[destination_idx] += t;
-      else
-        out[destination_idx] = t;
+        if (add) out[destination_idx] += t;
+        else out[destination_idx] = t;)
     }
 
 
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    inline __device__ void
+    DEAL_II_HOST_DEVICE inline void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -268,7 +265,7 @@ namespace CUDAWrappers
             {
               values<0, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
 
@@ -278,10 +275,10 @@ namespace CUDAWrappers
             {
               values<0, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<2, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
 
@@ -298,7 +295,7 @@ namespace CUDAWrappers
 
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    inline __device__ void
+    DEAL_II_HOST_DEVICE inline void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -318,7 +315,7 @@ namespace CUDAWrappers
             {
               values<0, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
 
@@ -328,10 +325,10 @@ namespace CUDAWrappers
             {
               values<0, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<2, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
 
@@ -348,7 +345,7 @@ namespace CUDAWrappers
 
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    inline __device__ void
+    DEAL_II_HOST_DEVICE inline void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -372,7 +369,7 @@ namespace CUDAWrappers
               values<0, true, false, false>(
                 get_global_shape_values<Number>(mf_object_id), u, grad_u[1]);
 
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<1, true, false, true>(get_global_shape_values<Number>(
                                              mf_object_id),
@@ -394,7 +391,7 @@ namespace CUDAWrappers
               values<0, true, false, false>(
                 get_global_shape_values<Number>(mf_object_id), u, grad_u[2]);
 
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<1, true, false, true>(get_global_shape_values<Number>(
                                              mf_object_id),
@@ -409,7 +406,7 @@ namespace CUDAWrappers
                                            grad_u[2],
                                            grad_u[2]);
 
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<2, true, false, true>(get_global_shape_values<Number>(
                                              mf_object_id),
@@ -437,7 +434,7 @@ namespace CUDAWrappers
 
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    inline __device__ void
+    DEAL_II_HOST_DEVICE inline void
     EvaluatorTensorProduct<
       evaluate_general,
       dim,
@@ -452,7 +449,7 @@ namespace CUDAWrappers
             {
               values<0, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               gradients<0, true, false, false>(
                 get_global_co_shape_gradients<Number>(mf_object_id),
@@ -465,10 +462,10 @@ namespace CUDAWrappers
             {
               values<0, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               gradients<0, true, false, false>(
                 get_global_co_shape_gradients<Number>(mf_object_id),
@@ -485,13 +482,13 @@ namespace CUDAWrappers
             {
               values<0, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<2, true, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               gradients<0, true, false, false>(
                 get_global_co_shape_gradients<Number>(mf_object_id),
@@ -520,7 +517,7 @@ namespace CUDAWrappers
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
     template <bool add>
-    inline __device__ void
+    DEAL_II_HOST_DEVICE inline void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -550,11 +547,11 @@ namespace CUDAWrappers
                                             grad_u[1],
                                             grad_u[1]);
 
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<1, false, add, false>(
                 get_global_shape_values<Number>(mf_object_id), grad_u[0], u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               gradients<1, false, true, false>(
                 get_global_shape_gradients<Number>(mf_object_id), grad_u[1], u);
 
@@ -575,7 +572,7 @@ namespace CUDAWrappers
                                             grad_u[2],
                                             grad_u[2]);
 
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<1, false, false, true>(get_global_shape_values<Number>(
                                               mf_object_id),
@@ -590,14 +587,14 @@ namespace CUDAWrappers
                                             grad_u[2],
                                             grad_u[2]);
 
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<2, false, add, false>(
                 get_global_shape_values<Number>(mf_object_id), grad_u[0], u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<2, false, true, false>(
                 get_global_shape_values<Number>(mf_object_id), grad_u[1], u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               gradients<2, false, true, false>(
                 get_global_shape_gradients<Number>(mf_object_id), grad_u[2], u);
 
@@ -614,7 +611,7 @@ namespace CUDAWrappers
 
 
     template <int dim, int fe_degree, int n_q_points_1d, typename Number>
-    inline __device__ void
+    DEAL_II_HOST_DEVICE inline void
     EvaluatorTensorProduct<evaluate_general,
                            dim,
                            fe_degree,
@@ -631,7 +628,7 @@ namespace CUDAWrappers
                 get_global_co_shape_gradients<Number>(mf_object_id),
                 grad_u[0],
                 u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<0, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
@@ -644,19 +641,19 @@ namespace CUDAWrappers
                 get_global_co_shape_gradients<Number>(mf_object_id),
                 grad_u[1],
                 u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               gradients<0, false, true, false>(
                 get_global_co_shape_gradients<Number>(mf_object_id),
                 grad_u[0],
                 u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<1, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<0, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               break;
             }
@@ -666,27 +663,27 @@ namespace CUDAWrappers
                 get_global_co_shape_gradients<Number>(mf_object_id),
                 grad_u[2],
                 u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               gradients<1, false, true, false>(
                 get_global_co_shape_gradients<Number>(mf_object_id),
                 grad_u[1],
                 u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               gradients<0, false, true, false>(
                 get_global_co_shape_gradients<Number>(mf_object_id),
                 grad_u[0],
                 u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               values<2, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<1, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
               values<0, false, false, true>(
                 get_global_shape_values<Number>(mf_object_id), u, u);
-              __syncthreads();
+              KOKKOS_IF_ON_DEVICE(__syncthreads();)
 
               break;
             }
index 38da36089baec046dfabaa631efbdf856323e6fd..8db177d609abbf2ca4421c1ec2c200f50d92ae1d 100644 (file)
@@ -35,7 +35,7 @@ class DummyOperator
 public:
   DummyOperator() = default;
 
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   operator()(
     const unsigned int                                          cell,
     const typename CUDAWrappers::MatrixFree<dim, double>::Data *gpu_data,
@@ -53,7 +53,7 @@ public:
 
 
 template <int dim, int fe_degree>
-__device__ void
+DEAL_II_HOST_DEVICE void
 DummyOperator<dim, fe_degree>::operator()(
   const unsigned int                                          cell,
   const typename CUDAWrappers::MatrixFree<dim, double>::Data *gpu_data,
index 0adf975f14dee1171763bde062e7379e437a5cea..204658a41aac33c0796bd1d67bf9f24e27654bed 100644 (file)
@@ -38,7 +38,7 @@ public:
   MatrixFreeTest(const CUDAWrappers::MatrixFree<dim, Number> &data_in)
     : data(data_in){};
 
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   operator()(
     const unsigned int                                          cell,
     const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data,
@@ -51,7 +51,7 @@ public:
 
     // set to unit vector
     fe_eval.submit_dof_value(1.);
-    __syncthreads();
+    KOKKOS_IF_ON_DEVICE(__syncthreads();)
     fe_eval.evaluate(/*evaluate_values =*/true, /*evaluate_gradients=*/true);
 
 #ifndef __APPLE__
index baeb5d6e6f1fc0bf22431ed6d2d500145a52a740..20af64558eabcd6ef41d474d3dc39d6c503c2d8e 100644 (file)
@@ -30,12 +30,12 @@ template <int dim, int fe_degree, typename Number, int n_q_points_1d>
 class HelmholtzOperatorQuad
 {
 public:
-  __device__
+  DEAL_II_HOST_DEVICE
   HelmholtzOperatorQuad(Number coef)
     : coef(coef)
   {}
 
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   operator()(
     CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number>
       *fe_eval) const;
@@ -47,7 +47,7 @@ private:
 
 
 template <int dim, int fe_degree, typename Number, int n_q_points_1d>
-__device__ void
+DEAL_II_HOST_DEVICE void
 HelmholtzOperatorQuad<dim, fe_degree, Number, n_q_points_1d>::operator()(
   CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number> *fe_eval)
   const
@@ -66,7 +66,7 @@ public:
     : coef(coefficient)
   {}
 
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   operator()(
     const unsigned int                                          cell,
     const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data,
@@ -86,7 +86,7 @@ public:
 
 
 template <int dim, int fe_degree, typename Number, int n_q_points_1d>
-__device__ void
+DEAL_II_HOST_DEVICE void
 HelmholtzOperator<dim, fe_degree, Number, n_q_points_1d>::operator()(
   const unsigned int                                          cell,
   const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data,
@@ -117,7 +117,7 @@ public:
     : coef(coefficient)
   {}
 
-  __device__ void
+  DEAL_II_HOST_DEVICE void
   operator()(
     const unsigned int                                          cell,
     const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data);
@@ -135,7 +135,7 @@ private:
 
 
 template <int dim, int fe_degree, typename Number, int n_q_points_1d>
-__device__ void
+DEAL_II_HOST_DEVICE void
 VaryingCoefficientFunctor<dim, fe_degree, Number, n_q_points_1d>::operator()(
   const unsigned int                                          cell,
   const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data)

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.