Eliminate q_point from functor and other GPU functions

diff --git a/doc/news/changes/minor/20190924PeterMunch b/doc/news/changes/minor/20190924PeterMunch
new file mode 100644 (file)
index 0000000..9232cc6
--- /dev/null
+++ b/doc/news/changes/minor/20190924PeterMunch
@@ -0,0 +1,5 @@
+New: Provide new methods in CUDAWrappers::FEEvaluation that do neither take the 
+local dof index nor the quadrature point but recompute them based on the 
+thread id.
+<br>
+(Peter Munch, 2019/09/24)

diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h
index b2bb935ccfa457352bee0c4447edd6e00d16a1e7..4c37b94981c39c20ca9d5a4cf5d38085bfd3e046 100644 (file)
--- a/include/deal.II/matrix_free/cuda_fe_evaluation.h
+++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h
@@ -38,6 +38,25 @@ DEAL_II_NAMESPACE_OPEN
  */
 namespace CUDAWrappers
 {
+  namespace internal
+  {
+    /**
+     * Compute the dof/quad index for a given thread id, dimension, and
+     * number of points in each space dimensions.
+     */
+    template <int dim, int n_points_1d>
+    __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));
+    }
+  } // namespace internal
+
   /**
    * This class provides all the functions necessary to evaluate functions at
    * quadrature points and cell integrations. In functionality, this class is
@@ -134,48 +153,103 @@ namespace CUDAWrappers
     /**
      * Return the value of a finite element function at quadrature point
      * number @p q_point after a call to @p evaluate(true,...).
+     *
+     * @deprecated Use the version without parameters instead.
+     */
+    DEAL_II_DEPRECATED __device__ value_type
+                                  get_value(const unsigned int q_point) const;
+
+    /**
+     * Same as above, except that the quadrature point is computed from thread
+     * id.
      */
     __device__ value_type
-               get_value(const unsigned int q_point) const;
+               get_value() const;
 
     /**
      * Return the value of a finite element function at degree of freedom
      * @p dof after a call to integrate() or before a call to evaluate().
+     *
+     * @deprecated Use the version without parameters instead.
+     */
+    DEAL_II_DEPRECATED __device__ value_type
+                                  get_dof_value(const unsigned int dof) const;
+
+    /**
+     * Same as above, except that the local dof index is computed from the
+     * thread id.
      */
     __device__ value_type
-               get_dof_value(const unsigned int dof) const;
+               get_dof_value() const;
 
     /**
      * Write a value to the field containing the values on quadrature points
      * with component @p q_point. Access to the same fields as through @p
      * get_value(). This specifies the value which is tested by all basis
      * function on the current cell and integrated over.
+     *
+     * @deprecated Use the version without parameters instead.
+     */
+    DEAL_II_DEPRECATED __device__ void
+                       submit_value(const value_type &val_in, const unsigned int q_point);
+
+    /**
+     * Same as above, except that the quadrature point is computed from the
+     * thread id.
      */
     __device__ void
-    submit_value(const value_type &val_in, const unsigned int q_point);
+    submit_value(const value_type &val_in);
 
     /**
      * Write a value to the field containing the values for the degree of
      * freedom with index @p dof after a call to integrate() or before
      * calling evaluate(). Access through the same fields as through
      * get_dof_value().
+     *
+     * @deprecated Use the version without parameters instead.
+     */
+    DEAL_II_DEPRECATED __device__ void
+                       submit_dof_value(const value_type &val_in, const unsigned int dof);
+
+    /**
+     * Same as above, except that the local dof index is computed from the
+     * thread id.
      */
     __device__ void
-    submit_dof_value(const value_type &val_in, const unsigned int dof);
+    submit_dof_value(const value_type &val_in);
 
     /**
      * Return the gradient of a finite element function at quadrature point
      * number @p q_point after a call to @p evaluate(...,true).
+     *
+     * @deprecated Use the version without parameters instead.
+     */
+    DEAL_II_DEPRECATED __device__ gradient_type
+                                  get_gradient(const unsigned int q_point) const;
+
+    /**
+     * Same as above, except that the quadrature point is computed from the
+     * thread id.
      */
     __device__ gradient_type
-               get_gradient(const unsigned int q_point) const;
+               get_gradient() const;
 
     /**
      * Write a contribution that is tested by the gradient to the field
-     * containing the values on quadrature points with component @p q_point
+     * containing the values on quadrature points with component @p q_point.
+     *
+     * @deprecated Use the version without parameters instead.
+     */
+    DEAL_II_DEPRECATED __device__ void
+                       submit_gradient(const gradient_type &grad_in, const unsigned int q_point);
+
+
+    /**
+     * Same as above, except that the quadrature point is computed from the
+     * thread id.
      */
     __device__ void
-    submit_gradient(const gradient_type &grad_in, const unsigned int q_point);
+    submit_gradient(const gradient_type &grad_in);
 
     // clang-format off
     /**
@@ -187,11 +261,29 @@ namespace CUDAWrappers
      *   CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, n_components, Number> *fe_eval,
      *   const unsigned int                                                               q_point) const;
      * \endcode
+     * 
+     * @deprecated Use apply_for_each_quad_point() instead.
+     */
+    // clang-format on
+    template <typename Functor>
+    DEAL_II_DEPRECATED __device__ void
+                       apply_quad_point_operations(const Functor &func);
+
+    // clang-format off
+    /**
+     * Same as above, except that the functor @func only takes a single input 
+     * argument (fe_eval) and computes the quadrature point from the thread id.
+     *
+     * @p func needs to define
+     * \code
+     * __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
-    apply_quad_point_operations(const Functor &func);
+    apply_for_each_quad_point(const Functor &func);
 
   private:
     types::global_dof_index *local_to_global;
@@ -249,10 +341,7 @@ namespace CUDAWrappers
   {
     static_assert(n_components_ == 1, "This function only supports FE with one \
                   components");
-    const unsigned int idx =
-      (threadIdx.x % n_q_points_1d) +
-      (dim > 1 ? threadIdx.y : 0) * n_q_points_1d +
-      (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d;
+    const unsigned int idx = internal::compute_index<dim, n_q_points_1d>();
 
     const types::global_dof_index src_idx = local_to_global[idx];
     // Use the read-only data cache.
@@ -280,10 +369,8 @@ namespace CUDAWrappers
     internal::resolve_hanging_nodes<dim, fe_degree, true>(constraint_mask,
                                                           values);
 
-    const unsigned int idx =
-      (threadIdx.x % n_q_points_1d) +
-      (dim > 1 ? threadIdx.y : 0) * n_q_points_1d +
-      (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d;
+    const unsigned int idx = internal::compute_index<dim, n_q_points_1d>();
+
     const types::global_dof_index destination_idx = local_to_global[idx];
 
     if (use_coloring)
@@ -384,6 +471,24 @@ namespace CUDAWrappers
 
 
 
+  template <int dim,
+            int fe_degree,
+            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
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    get_value() const
+  {
+    return get_value(internal::compute_index<dim, n_q_points_1d>());
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -402,6 +507,24 @@ namespace CUDAWrappers
 
 
 
+  template <int dim,
+            int fe_degree,
+            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
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    get_dof_value() const
+  {
+    return get_dof_value(internal::compute_index<dim, fe_degree + 1>());
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -416,6 +539,20 @@ namespace CUDAWrappers
 
 
 
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components_,
+            typename Number>
+  __device__ void
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    submit_value(const value_type &val_in)
+  {
+    submit_value(val_in, internal::compute_index<dim, n_q_points_1d>());
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -430,6 +567,20 @@ namespace CUDAWrappers
 
 
 
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components_,
+            typename Number>
+  __device__ void
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    submit_dof_value(const value_type &val_in)
+  {
+    submit_dof_value(val_in, internal::compute_index<dim, fe_degree + 1>());
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -462,6 +613,27 @@ namespace CUDAWrappers
 
 
 
+  template <int dim,
+            int fe_degree,
+            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
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    get_gradient() const
+  {
+    static_assert(n_components_ == 1, "This function only supports FE with one \
+                  components");
+
+    return get_gradient(internal::compute_index<dim, n_q_points_1d>());
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -485,6 +657,20 @@ namespace CUDAWrappers
 
 
 
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components_,
+            typename Number>
+  __device__ void
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    submit_gradient(const gradient_type &grad_in)
+  {
+    submit_gradient(grad_in, internal::compute_index<dim, n_q_points_1d>());
+  }
+
+
+
   template <int dim,
             int fe_degree,
             int n_q_points_1d,
@@ -495,14 +681,24 @@ namespace CUDAWrappers
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     apply_quad_point_operations(const Functor &func)
   {
-    const unsigned int q_point =
-      (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));
-    func(this, q_point);
+    func(this, internal::compute_index<dim, n_q_points_1d>());
+
+    __syncthreads();
+  }
+
+
+
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components_,
+            typename Number>
+  template <typename Functor>
+  __device__ void
+  FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
+    apply_for_each_quad_point(const Functor &func)
+  {
+    func(this);
 
     __syncthreads();
   }

diff --git a/tests/cuda/matrix_free_no_index_initialize.cu b/tests/cuda/matrix_free_no_index_initialize.cu
index 85b2a780acffa2acefd44c916ee73a48b5ef538e..a66a64fbd9ccb87ba2ba47b60aaa6cd88c6f3fbd 100644 (file)
--- a/tests/cuda/matrix_free_no_index_initialize.cu
+++ b/tests/cuda/matrix_free_no_index_initialize.cu
@@ -49,25 +49,20 @@ public:
     CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number>
       fe_eval(cell, gpu_data, shared_data);
 
-    const unsigned int tid =
-      (threadIdx.x % n_q_points_1d) +
-      (dim > 1 ? threadIdx.y : 0) * n_q_points_1d +
-      (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d;
-
     // set to unit vector
-    fe_eval.submit_dof_value(1., tid);
+    fe_eval.submit_dof_value(1.);
     __syncthreads();
     fe_eval.evaluate(/*evaluate_values =*/true, /*evaluate_gradients=*/true);
 
 #ifndef __APPLE__
     // values should evaluate to one, derivatives to zero
-    assert(fe_eval.get_value(tid) == 1.);
+    assert(fe_eval.get_value() == 1.);
     for (unsigned int e = 0; e < dim; ++e)
-      assert(fe_eval.get_gradient(tid)[e] == 0.);
+      assert(fe_eval.get_gradient()[e] == 0.);
 
     fe_eval.integrate(/*integrate_values = */ true,
                       /*integrate_gradients=*/true);
-    assert(fe_eval.get_dof_value(tid) == 1.);
+    assert(fe_eval.get_dof_value() == 1.);
 #endif
   }
 

diff --git a/tests/cuda/matrix_vector_mf.h b/tests/cuda/matrix_vector_mf.h
index ef4a0c3d4f1ab53f6a08cb841bd738eb5774ca3b..e27fdc1d0ee7447f23d08d63c67d5655e5c3f0e1 100644 (file)
--- a/tests/cuda/matrix_vector_mf.h
+++ b/tests/cuda/matrix_vector_mf.h
@@ -37,8 +37,7 @@ public:
 
   __device__ void operator()(
     CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number>
-      *                fe_eval,
-    const unsigned int q_point) const;
+      *fe_eval) const;
 
 private:
   Number coef;
@@ -48,12 +47,11 @@ private:
 
 template <int dim, int fe_degree, typename Number, int n_q_points_1d>
 __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 unsigned int                                                    q) const
+                operator()(CUDAWrappers::FEEvaluation<dim, fe_degree, n_q_points_1d, 1, Number>
+             *fe_eval) const
 {
-  fe_eval->submit_value(coef * fe_eval->get_value(q), q);
-  fe_eval->submit_gradient(fe_eval->get_gradient(q), q);
+  fe_eval->submit_value(coef * fe_eval->get_value());
+  fe_eval->submit_gradient(fe_eval->get_gradient());
 }
 
 
@@ -101,7 +99,7 @@ operator()(const unsigned int                                          cell,
     cell, gpu_data, shared_data);
   fe_eval.read_dof_values(src);
   fe_eval.evaluate(true, true);
-  fe_eval.apply_quad_point_operations(
+  fe_eval.apply_for_each_quad_point(
     HelmholtzOperatorQuad<dim, fe_degree, Number, n_q_points_1d>(coef[pos]));
   fe_eval.integrate(true, true);
   fe_eval.distribute_local_to_global(dst);