Move q_points to Kokkos::View

diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h
index 069e98e74a52f62c6e75f044d978f39a2910293e..d6588813f6966176ca43e674c970b04c17cef16d 100644 (file)
--- a/include/deal.II/matrix_free/cuda_matrix_free.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.h
@@ -20,6 +20,7 @@
 #include <deal.II/base/config.h>
 
 #include "deal.II/base/memory_space.h"
+#include "deal.II/base/utilities.h"
 
 #ifdef DEAL_II_WITH_CUDA
 
@@ -174,9 +175,9 @@ namespace CUDAWrappers
     struct Data
     {
       /**
-       * Pointer to the quadrature points.
+       * Kokkos::View of the quadrature points.
        */
-      point_type *q_points;
+      Kokkos::View<point_type **, MemorySpace::Default::kokkos_space> q_points;
 
       /**
        * Map the position in the local vector to the position in the global
@@ -503,10 +504,11 @@ namespace CUDAWrappers
     std::vector<unsigned int> n_cells;
 
     /**
-     * Vector of pointers to the quadrature points associated to the cells of
-     * each color.
+     * Vector of Kokkos::View to the quadrature points associated to the cells
+     * of each color.
      */
-    std::vector<point_type *> q_points;
+    std::vector<Kokkos::View<point_type **, MemorySpace::Default::kokkos_space>>
+      q_points;
 
     /**
      * Map the position in the local vector to the position in the global
@@ -714,8 +716,7 @@ namespace CUDAWrappers
       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));
+    return data->q_points(cell, q_point_id_in_cell<dim>(n_q_points_1d));
   }
 
   /**
@@ -726,9 +727,11 @@ namespace CUDAWrappers
   struct DataHost
   {
     /**
-     * Vector of quadrature points.
+     * Kokkos::View of quadrature points on the host.
      */
-    std::vector<Point<dim, Number>> q_points;
+    typename Kokkos::View<Point<dim, Number> **,
+                          MemorySpace::Default::kokkos_space>::HostMirror
+      q_points;
 
     /**
      * Map the position in the local vector to the position in the global
@@ -737,14 +740,17 @@ namespace CUDAWrappers
     std::vector<types::global_dof_index> local_to_global;
 
     /**
-     * Vector of inverse Jacobians.
+     * Kokkos::View of inverse Jacobians on the host.
      */
-    std::vector<Number> inv_jacobian;
+    typename Kokkos::View<Number **[dim][dim],
+                          MemorySpace::Default::kokkos_space>::HostMirror
+      inv_jacobian;
 
     /**
-     * Vector of Jacobian times the weights.
+     * Kokkos::View of Jacobian times the weights on the host.
      */
-    std::vector<Number> JxW;
+    typename Kokkos::View<Number **,
+                          MemorySpace::Default::kokkos_space>::HostMirror JxW;
 
     /**
      * ID of the associated MatrixFree object.
@@ -805,8 +811,8 @@ namespace CUDAWrappers
       data_host.n_cells * data_host.padding_length;
     if (update_flags & update_quadrature_points)
       {
-        data_host.q_points.resize(n_elements);
-        Utilities::CUDA::copy_to_host(data.q_points, data_host.q_points);
+        data_host.q_points = Kokkos::create_mirror(data.q_points);
+        Kokkos::deep_copy(data_host.q_points, data.q_points);
       }
 
     data_host.local_to_global.resize(n_elements);
@@ -815,15 +821,14 @@ namespace CUDAWrappers
 
     if (update_flags & update_gradients)
       {
-        data_host.inv_jacobian.resize(n_elements * dim * dim);
-        Utilities::CUDA::copy_to_host(data.inv_jacobian.data(),
-                                      data_host.inv_jacobian);
+        data_host.inv_jacobian = Kokkos::create_mirror(data.inv_jacobian);
+        Kokkos::deep_copy(data_host.inv_jacobian, data.inv_jacobian);
       }
 
     if (update_flags & update_JxW_values)
       {
-        data_host.JxW.resize(n_elements);
-        Utilities::CUDA::copy_to_host(data.JxW.data(), data_host.JxW);
+        data_host.JxW = Kokkos::create_mirror(data.JxW);
+        Kokkos::deep_copy(data_host.JxW, data.JxW);
       }
 
     data_host.constraint_mask.resize(data_host.n_cells);
@@ -865,7 +870,7 @@ namespace CUDAWrappers
                             const DataHost<dim, Number> &data,
                             const unsigned int           i)
   {
-    return data.q_points[data.padding_length * cell + i];
+    return data.q_points(cell, i);
   }
 
 

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 c50828c1808813251cb6f49183900b61910dcec5..d23c1b82bc10e42da251a9eb2283d9a51512e438 100644 (file)
--- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h
+++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h
@@ -234,7 +234,8 @@ namespace CUDAWrappers
       MatrixFree<dim, Number> *data;
       // Host data
       std::vector<types::global_dof_index> local_to_global_host;
-      std::vector<Point<dim, Number>>      q_points_host;
+      Kokkos::View<Point<dim, Number> **, MemorySpace::Default::kokkos_space>
+                                                                  q_points;
       Kokkos::View<Number **, MemorySpace::Default::kokkos_space> JxW;
       Kokkos::View<Number **[dim][dim], MemorySpace::Default::kokkos_space>
         inv_jacobian;
@@ -350,7 +351,12 @@ namespace CUDAWrappers
       local_to_global_host.resize(n_cells * padding_length);
 
       if (update_flags & update_quadrature_points)
-        q_points_host.resize(n_cells * padding_length);
+        q_points = Kokkos::View<Point<dim, Number> **,
+                                MemorySpace::Default::kokkos_space>(
+          Kokkos::view_alloc("q_points_" + std::to_string(color),
+                             Kokkos::WithoutInitializing),
+          n_cells,
+          q_points_per_cell);
 
       if (update_flags & update_JxW_values)
         JxW = Kokkos::View<Number **, MemorySpace::Default::kokkos_space>(
@@ -409,11 +415,14 @@ namespace CUDAWrappers
       // Quadrature points
       if (update_flags & update_quadrature_points)
         {
-          const std::vector<Point<dim>> &q_points =
+          // FIXME too many deep_copy
+          auto q_points_host = Kokkos::create_mirror_view_and_copy(
+            MemorySpace::Host::kokkos_space{}, q_points);
+          const std::vector<Point<dim>> &q_points_vec =
             fe_values.get_quadrature_points();
-          std::copy(q_points.begin(),
-                    q_points.end(),
-                    &q_points_host[cell_id * padding_length]);
+          for (unsigned int i = 0; i < q_points_per_cell; ++i)
+            q_points_host(cell_id, i) = q_points_vec[i];
+          Kokkos::deep_copy(q_points, q_points_host);
         }
 
       if (update_flags & update_JxW_values)
@@ -460,10 +469,7 @@ namespace CUDAWrappers
       // Quadrature points
       if (update_flags & update_quadrature_points)
         {
-          alloc_and_copy(&data->q_points[color],
-                         ArrayView<const Point<dim, Number>>(
-                           q_points_host.data(), q_points_host.size()),
-                         n_cells * padding_length);
+          data->q_points[color] = q_points;
         }
 
       // Jacobian determinants/quadrature weights
@@ -706,10 +712,6 @@ namespace CUDAWrappers
   void
   MatrixFree<dim, Number>::free()
   {
-    for (auto &q_points_color_ptr : q_points)
-      Utilities::CUDA::free(q_points_color_ptr);
-    q_points.clear();
-
     for (auto &local_to_global_color_ptr : local_to_global)
       Utilities::CUDA::free(local_to_global_color_ptr);
     local_to_global.clear();