Add helper functions to use varying coefficient and matrix-free

diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h
index 5c57033113bfdc706ab009e129560cb83fedd179..7bdec6767650e80cd4363bce1b6555b274892b64 100644 (file)
--- a/include/deal.II/matrix_free/cuda_matrix_free.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.h
@@ -210,7 +210,7 @@ namespace CUDAWrappers
 
     /**
      * This method runs the loop over all cells and apply the local operation on
-     * each element in parallel. @p func is a functor which is appplied on each color.
+     * each element in parallel. @p func is a functor which is applied on each color.
      */
     template <typename functor, typename VectorType>
     void
@@ -218,6 +218,25 @@ namespace CUDAWrappers
               const VectorType &src,
               VectorType &      dst) const;
 
+    /**
+     * This method runs the loop over all cells and apply the local operation on
+     * each element in parallel. This function is very similar to cell_loop()
+     * but it uses a simpler functor.
+     *
+     * @p func needs to define
+     * \code
+     *  __device__ void operator()(
+     *    const unsigned int                                          cell,
+     *    const typename CUDAWrappers::MatrixFree<dim, Number>::Data *gpu_data);
+     * static const unsigned int n_dofs_1d;
+     * static const unsigned int n_local_dofs;
+     * static const unsigned int n_q_points;
+     * \endcode
+     */
+    template <typename Functor>
+    void
+    evaluate_coefficients(Functor func) const;
+
     /**
      * Copy the values of the constrained entries from @p src to @p dst. This is
      * used to impose zero Dirichlet boundary condition.
@@ -482,6 +501,63 @@ namespace CUDAWrappers
                      1) : 1;
     /* clang-format on */
   }
+
+
+
+  /**
+   * Compute the quadrature point index in the local cell of a given thread.
+   *
+   * @relates MatrixFree
+   */
+  template <int dim>
+  __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));
+  }
+
+
+
+  /**
+   * Return the quadrature point index local of a given thread. The index is
+   * only unique for a given MPI process.
+   *
+   * @relates MatrixFree
+   */
+  template <int dim, typename Number>
+  __device__ inline unsigned int
+  local_q_point_id(
+    const unsigned int                                          cell,
+    const typename CUDAWrappers::MatrixFree<dim, Number>::Data *data,
+    const unsigned int                                          n_q_points_1d,
+    const unsigned int                                          n_q_points)
+  {
+    return (data->row_start / data->padding_length + cell) * n_q_points +
+           q_point_id_in_cell<dim>(n_q_points_1d);
+  }
+
+
+
+  /**
+   * Return the quadrature point associated with a given thread.
+   *
+   * @relates 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)
+  {
+    return *(data->q_points + data->padding_length * cell +
+             q_point_id_in_cell<dim>(n_q_points_1d));
+  }
 } // namespace CUDAWrappers
 
 DEAL_II_NAMESPACE_CLOSE

diff --git a/include/deal.II/matrix_free/cuda_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
index cb6e4081527f2bac1a32f4e7090e3cfab86c150f..ecfd74a09a9266b57a313c4a232564637fc74c89 100644 (file)
--- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
@@ -474,7 +474,7 @@ namespace CUDAWrappers
 
     template <int dim, typename Number, typename functor>
     __global__ void
-    apply_kernel_shmem(const functor &                              func,
+    apply_kernel_shmem(functor                                      func,
                        const typename MatrixFree<dim, Number>::Data gpu_data,
                        const Number *                               src,
                        Number *                                     dst)
@@ -504,6 +504,29 @@ namespace CUDAWrappers
       if (cell < gpu_data.n_cells)
         func(cell, &gpu_data, &shared_data, src, dst);
     }
+
+
+
+    template <int dim, typename Number, typename Functor>
+    __global__ void
+    evaluate_coeff(Functor                                      func,
+                   const typename MatrixFree<dim, Number>::Data gpu_data)
+    {
+      constexpr unsigned int cells_per_block =
+        cells_per_block_shmem(dim, Functor::n_dofs_1d - 1);
+
+      constexpr unsigned int n_dofs_per_block =
+        cells_per_block * Functor::n_local_dofs;
+      constexpr unsigned int n_q_points_per_block =
+        cells_per_block * Functor::n_q_points;
+
+      const unsigned int local_cell = threadIdx.x / Functor::n_dofs_1d;
+      const unsigned int cell =
+        local_cell + cells_per_block * (blockIdx.x + gridDim.x * blockIdx.y);
+
+      if (cell < gpu_data.n_cells)
+        func(cell, &gpu_data);
+    }
   } // namespace internal
 
 
@@ -732,6 +755,18 @@ namespace CUDAWrappers
 
 
 
+  template <int dim, typename Number>
+  template <typename Functor>
+  void
+  MatrixFree<dim, Number>::evaluate_coefficients(Functor func) const
+  {
+    for (unsigned int i = 0; i < n_colors; ++i)
+      internal::evaluate_coeff<dim, Number, Functor>
+        <<<grid_dim[i], block_dim[i]>>>(func, get_data(i));
+  }
+
+
+
   template <int dim, typename Number>
   std::size_t
   MatrixFree<dim, Number>::memory_consumption() const