]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Support multiple components in Portable::MatrixFree. 17525/head
authorPeter Munch <peterrmuench@gmail.com>
Tue, 13 Aug 2024 14:46:05 +0000 (16:46 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 15 Aug 2024 19:29:33 +0000 (21:29 +0200)
doc/news/changes/minor/20240814Munch [new file with mode: 0644]
include/deal.II/matrix_free/portable_fe_evaluation.h
include/deal.II/matrix_free/portable_hanging_nodes_internal.h
include/deal.II/matrix_free/portable_matrix_free.h
include/deal.II/matrix_free/portable_matrix_free.templates.h
tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc [new file with mode: 0644]
tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20240814Munch b/doc/news/changes/minor/20240814Munch
new file mode 100644 (file)
index 0000000..59da529
--- /dev/null
@@ -0,0 +1,4 @@
+Improved: Portable::MatrixFree now supports FESystem(FE_Q, n_components)
+with arbitrary number of components.
+<br>
+(Peter Munch, Daniel Arndt, 2024/08/14)
index f94adf548d8a47c14db32177db4219fffb0a256b..d2db2377ca5b642684b1f594b2c563c0a7973e9c 100644 (file)
@@ -72,12 +72,19 @@ namespace Portable
     /**
      * An alias for scalar quantities.
      */
-    using value_type = Number;
+    using value_type = std::conditional_t<(n_components_ == 1),
+                                          Number,
+                                          Tensor<1, n_components_, Number>>;
 
     /**
      * An alias for vectorial quantities.
      */
-    using gradient_type = Tensor<1, dim, Number>;
+    using gradient_type = std::conditional_t<
+      n_components_ == 1,
+      Tensor<1, dim, Number>,
+      std::conditional_t<n_components_ == dim,
+                         Tensor<2, dim, Number>,
+                         Tensor<1, n_components_, Tensor<1, dim, Number>>>>;
 
     /**
      * An alias to kernel specific information.
@@ -265,22 +272,24 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     read_dof_values(const Number *src)
   {
-    static_assert(n_components_ == 1, "This function only supports FE with one \
-                  components");
     // Populate the scratch memory
-    Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
-                                                 n_q_points),
-                         [&](const int &i) {
-                           shared_data->values(i) =
-                             src[data->local_to_global(cell_id, i)];
-                         });
+    Kokkos::parallel_for(
+      Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
+      [&](const int &i) {
+        for (unsigned int c = 0; c < n_components_; ++c)
+          shared_data->values(i, c) =
+            src[data->local_to_global(cell_id, i + tensor_dofs_per_cell * c)];
+      });
     shared_data->team_member.team_barrier();
 
-    internal::resolve_hanging_nodes<dim, fe_degree, false>(
-      shared_data->team_member,
-      data->constraint_weights,
-      data->constraint_mask(cell_id),
-      shared_data->values);
+    for (unsigned int c = 0; c < n_components_; ++c)
+      {
+        internal::resolve_hanging_nodes<dim, fe_degree, false, Number>(
+          shared_data->team_member,
+          data->constraint_weights,
+          data->constraint_mask(cell_id),
+          Kokkos::subview(shared_data->values, Kokkos::ALL, c));
+      }
   }
 
 
@@ -294,31 +303,35 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     distribute_local_to_global(Number *dst) const
   {
-    static_assert(n_components_ == 1, "This function only supports FE with one \
-                  components");
-
-    internal::resolve_hanging_nodes<dim, fe_degree, true>(
-      shared_data->team_member,
-      data->constraint_weights,
-      data->constraint_mask(cell_id),
-      shared_data->values);
+    for (unsigned int c = 0; c < n_components_; ++c)
+      {
+        internal::resolve_hanging_nodes<dim, fe_degree, true, Number>(
+          shared_data->team_member,
+          data->constraint_weights,
+          data->constraint_mask(cell_id),
+          Kokkos::subview(shared_data->values, Kokkos::ALL, c));
+      }
 
     if (data->use_coloring)
       {
-        Kokkos::parallel_for(Kokkos::TeamThreadRange(shared_data->team_member,
-                                                     n_q_points),
-                             [&](const int &i) {
-                               dst[data->local_to_global(cell_id, i)] +=
-                                 shared_data->values(i);
-                             });
+        Kokkos::parallel_for(
+          Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
+          [&](const int &i) {
+            for (unsigned int c = 0; c < n_components_; ++c)
+              dst[data->local_to_global(cell_id,
+                                        i + tensor_dofs_per_cell * c)] +=
+                shared_data->values(i, c);
+          });
       }
     else
       {
         Kokkos::parallel_for(
           Kokkos::TeamThreadRange(shared_data->team_member, n_q_points),
           [&](const int &i) {
-            Kokkos::atomic_add(&dst[data->local_to_global(cell_id, i)],
-                               shared_data->values(i));
+            for (unsigned int c = 0; c < n_components_; ++c)
+              Kokkos::atomic_add(&dst[data->local_to_global(
+                                   cell_id, i + tensor_dofs_per_cell * c)],
+                                 shared_data->values(i, c));
           });
       }
   }
@@ -347,23 +360,31 @@ namespace Portable
                                data->shape_gradients,
                                data->co_shape_gradients);
 
-    if ((evaluate_flag & EvaluationFlags::values) &&
-        (evaluate_flag & EvaluationFlags::gradients))
-      {
-        evaluator_tensor_product.evaluate_values_and_gradients(
-          shared_data->values, shared_data->gradients);
-        shared_data->team_member.team_barrier();
-      }
-    else if (evaluate_flag & EvaluationFlags::gradients)
-      {
-        evaluator_tensor_product.evaluate_gradients(shared_data->values,
-                                                    shared_data->gradients);
-        shared_data->team_member.team_barrier();
-      }
-    else if (evaluate_flag & EvaluationFlags::values)
+    for (unsigned int c = 0; c < n_components_; ++c)
       {
-        evaluator_tensor_product.evaluate_values(shared_data->values);
-        shared_data->team_member.team_barrier();
+        if ((evaluate_flag & EvaluationFlags::values) &&
+            (evaluate_flag & EvaluationFlags::gradients))
+          {
+            evaluator_tensor_product.evaluate_values_and_gradients(
+              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
+              Kokkos::subview(
+                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
+            shared_data->team_member.team_barrier();
+          }
+        else if (evaluate_flag & EvaluationFlags::gradients)
+          {
+            evaluator_tensor_product.evaluate_gradients(
+              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
+              Kokkos::subview(
+                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
+            shared_data->team_member.team_barrier();
+          }
+        else if (evaluate_flag & EvaluationFlags::values)
+          {
+            evaluator_tensor_product.evaluate_values(
+              Kokkos::subview(shared_data->values, Kokkos::ALL, c));
+            shared_data->team_member.team_barrier();
+          }
       }
   }
 
@@ -406,22 +427,31 @@ namespace Portable
                                data->shape_gradients,
                                data->co_shape_gradients);
 
-    if ((integration_flag & EvaluationFlags::values) &&
-        (integration_flag & EvaluationFlags::gradients))
-      {
-        evaluator_tensor_product.integrate_values_and_gradients(
-          shared_data->values, shared_data->gradients);
-      }
-    else if (integration_flag & EvaluationFlags::values)
-      {
-        evaluator_tensor_product.integrate_values(shared_data->values);
-        shared_data->team_member.team_barrier();
-      }
-    else if (integration_flag & EvaluationFlags::gradients)
+
+    for (unsigned int c = 0; c < n_components_; ++c)
       {
-        evaluator_tensor_product.template integrate_gradients<false>(
-          shared_data->values, shared_data->gradients);
-        shared_data->team_member.team_barrier();
+        if ((integration_flag & EvaluationFlags::values) &&
+            (integration_flag & EvaluationFlags::gradients))
+          {
+            evaluator_tensor_product.integrate_values_and_gradients(
+              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
+              Kokkos::subview(
+                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
+          }
+        else if (integration_flag & EvaluationFlags::values)
+          {
+            evaluator_tensor_product.integrate_values(
+              Kokkos::subview(shared_data->values, Kokkos::ALL, c));
+            shared_data->team_member.team_barrier();
+          }
+        else if (integration_flag & EvaluationFlags::gradients)
+          {
+            evaluator_tensor_product.template integrate_gradients<false>(
+              Kokkos::subview(shared_data->values, Kokkos::ALL, c),
+              Kokkos::subview(
+                shared_data->gradients, Kokkos::ALL, Kokkos::ALL, c));
+            shared_data->team_member.team_barrier();
+          }
       }
   }
 
@@ -457,7 +487,17 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::get_value(
     int q_point) const
   {
-    return shared_data->values(q_point);
+    if constexpr (n_components_ == 1)
+      {
+        return shared_data->values(q_point, 0);
+      }
+    else
+      {
+        value_type result;
+        for (unsigned int c = 0; c < n_components; ++c)
+          result[c] = shared_data->values(q_point, c);
+        return result;
+      }
   }
 
 
@@ -475,7 +515,17 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     get_dof_value(int q_point) const
   {
-    return shared_data->values(q_point);
+    if constexpr (n_components_ == 1)
+      {
+        return shared_data->values(q_point, 0);
+      }
+    else
+      {
+        value_type result;
+        for (unsigned int c = 0; c < n_components; ++c)
+          result[c] = shared_data->values(q_point, c);
+        return result;
+      }
   }
 
 
@@ -489,7 +539,16 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_value(const value_type &val_in, int q_point)
   {
-    shared_data->values(q_point) = val_in * data->JxW(cell_id, q_point);
+    if constexpr (n_components_ == 1)
+      {
+        shared_data->values(q_point, 0) = val_in * data->JxW(cell_id, q_point);
+      }
+    else
+      {
+        for (unsigned int c = 0; c < n_components; ++c)
+          shared_data->values(q_point, c) =
+            val_in[c] * data->JxW(cell_id, q_point);
+      }
   }
 
 
@@ -503,7 +562,15 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_dof_value(const value_type &val_in, int q_point)
   {
-    shared_data->values(q_point) = val_in;
+    if constexpr (n_components_ == 1)
+      {
+        shared_data->values(q_point, 0) = val_in;
+      }
+    else
+      {
+        for (unsigned int c = 0; c < n_components; ++c)
+          shared_data->values(q_point, c) = val_in[c];
+      }
   }
 
 
@@ -521,17 +588,30 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     get_gradient(int q_point) const
   {
-    static_assert(n_components_ == 1, "This function only supports FE with one \
-                  components");
-
     gradient_type grad;
-    for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
+
+    if constexpr (n_components_ == 1)
       {
-        Number tmp = 0.;
-        for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
-          tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) *
-                 shared_data->gradients(q_point, d_2);
-        grad[d_1] = tmp;
+        for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
+          {
+            Number tmp = 0.;
+            for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
+              tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) *
+                     shared_data->gradients(q_point, d_2, 0);
+            grad[d_1] = tmp;
+          }
+      }
+    else
+      {
+        for (unsigned int c = 0; c < n_components; ++c)
+          for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
+            {
+              Number tmp = 0.;
+              for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
+                tmp += data->inv_jacobian(cell_id, q_point, d_2, d_1) *
+                       shared_data->gradients(q_point, d_2, c);
+              grad[c][d_1] = tmp;
+            }
       }
 
     return grad;
@@ -548,13 +628,30 @@ namespace Portable
   FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     submit_gradient(const gradient_type &grad_in, int q_point)
   {
-    for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
+    if constexpr (n_components_ == 1)
+      {
+        for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
+          {
+            Number tmp = 0.;
+            for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
+              tmp +=
+                data->inv_jacobian(cell_id, q_point, d_1, d_2) * grad_in[d_2];
+            shared_data->gradients(q_point, d_1, 0) =
+              tmp * data->JxW(cell_id, q_point);
+          }
+      }
+    else
       {
-        Number tmp = 0.;
-        for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
-          tmp += data->inv_jacobian(cell_id, q_point, d_1, d_2) * grad_in[d_2];
-        shared_data->gradients(q_point, d_1) =
-          tmp * data->JxW(cell_id, q_point);
+        for (unsigned int c = 0; c < n_components; ++c)
+          for (unsigned int d_1 = 0; d_1 < dim; ++d_1)
+            {
+              Number tmp = 0.;
+              for (unsigned int d_2 = 0; d_2 < dim; ++d_2)
+                tmp += data->inv_jacobian(cell_id, q_point, d_1, d_2) *
+                       grad_in[c][d_2];
+              shared_data->gradients(q_point, d_1, c) =
+                tmp * data->JxW(cell_id, q_point);
+            }
       }
   }
 
index 3139473b3ff856d9c3b4c698fcd8e40544ae2653..d54dc72862b08da5e6dc10e67add6867dac25e2b 100644 (file)
@@ -121,7 +121,8 @@ namespace Portable
     template <unsigned int fe_degree,
               unsigned int direction,
               bool         transpose,
-              typename Number>
+              typename Number,
+              typename ViewType>
     DEAL_II_HOST_DEVICE inline void
     interpolate_boundary_2d(
       const Kokkos::TeamPolicy<
@@ -130,11 +131,8 @@ namespace Portable
       Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
         constraint_weights,
       const dealii::internal::MatrixFreeFunctions::ConstraintKinds
-                                                           &constraint_mask,
-      Kokkos::View<Number *,
-                   MemorySpace::Default::kokkos_space::execution_space::
-                     scratch_memory_space,
-                   Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
+              &constraint_mask,
+      ViewType values)
     {
       constexpr unsigned int n_q_points_1d = fe_degree + 1;
       constexpr unsigned int n_q_points    = Utilities::pow(n_q_points_1d, 2);
@@ -242,7 +240,8 @@ namespace Portable
     template <unsigned int fe_degree,
               unsigned int direction,
               bool         transpose,
-              typename Number>
+              typename Number,
+              typename ViewType>
     DEAL_II_HOST_DEVICE inline void
     interpolate_boundary_3d(
       const Kokkos::TeamPolicy<
@@ -251,11 +250,8 @@ namespace Portable
       Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
         constraint_weights,
       const dealii::internal::MatrixFreeFunctions::ConstraintKinds
-                                                            constraint_mask,
-      Kokkos::View<Number *,
-                   MemorySpace::Default::kokkos_space::execution_space::
-                     scratch_memory_space,
-                   Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
+               constraint_mask,
+      ViewType values)
     {
       constexpr unsigned int n_q_points_1d = fe_degree + 1;
       constexpr unsigned int n_q_points    = Utilities::pow(n_q_points_1d, 3);
@@ -410,7 +406,11 @@ namespace Portable
      * @cite ljungkvist2017matrix and in Section 3.4 of
      * @cite kronbichler2019multigrid.
      */
-    template <int dim, int fe_degree, bool transpose, typename Number>
+    template <int  dim,
+              int  fe_degree,
+              bool transpose,
+              typename Number,
+              typename ViewType>
     DEAL_II_HOST_DEVICE void
     resolve_hanging_nodes(
       const Kokkos::TeamPolicy<
@@ -419,11 +419,8 @@ namespace Portable
       Kokkos::View<Number *, MemorySpace::Default::kokkos_space>
         constraint_weights,
       const dealii::internal::MatrixFreeFunctions::ConstraintKinds
-                                                            constraint_mask,
-      Kokkos::View<Number *,
-                   MemorySpace::Default::kokkos_space::execution_space::
-                     scratch_memory_space,
-                   Kokkos::MemoryTraits<Kokkos::Unmanaged>> values)
+               constraint_mask,
+      ViewType values)
     {
       if constexpr (dim == 2)
         {
index 2b25b1ef3ff551f8e32f706d0017d4cdf8dd57e6..73551b2d3180de5b209eeccca029282a9e4d7e1c 100644 (file)
@@ -217,6 +217,11 @@ namespace Portable
        */
       unsigned int n_cells;
 
+      /**
+       * Number of components.
+       */
+      unsigned int n_components;
+
       /**
        * Length of the padding.
        */
@@ -504,6 +509,16 @@ namespace Portable
      */
     unsigned int fe_degree;
 
+    /**
+     * Number of components.
+     */
+    unsigned int n_components;
+
+    /**
+     * Number of scalar degrees of freedom per cell.
+     */
+    unsigned int scalar_dofs_per_cell;
+
     /**
      * Number of degrees of freedom per cell.
      */
@@ -634,19 +649,19 @@ namespace Portable
     using TeamHandle = Kokkos::TeamPolicy<
       MemorySpace::Default::kokkos_space::execution_space>::member_type;
 
-    using SharedView1D = Kokkos::View<
-      Number *,
+    using SharedViewValues = Kokkos::View<
+      Number **,
       MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
       Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
-    using SharedView2D = Kokkos::View<
-      Number *[dim],
+    using SharedViewGradients = Kokkos::View<
+      Number ***,
       MemorySpace::Default::kokkos_space::execution_space::scratch_memory_space,
       Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
 
     DEAL_II_HOST_DEVICE
-    SharedData(const TeamHandle   &team_member,
-               const SharedView1D &values,
-               const SharedView2D &gradients)
+    SharedData(const TeamHandle          &team_member,
+               const SharedViewValues    &values,
+               const SharedViewGradients &gradients)
       : team_member(team_member)
       , values(values)
       , gradients(gradients)
@@ -660,12 +675,12 @@ namespace Portable
     /**
      * Memory for dof and quad values.
      */
-    SharedView1D values;
+    SharedViewValues values;
 
     /**
      * Memory for computed gradients in reference coordinate system.
      */
-    SharedView2D gradients;
+    SharedViewGradients gradients;
   };
 
 
index 4c8a3db6d11c3db020da95f1e1747617e4e36f66..ebbaeef14d337b6141d53c7c9bc86c7a5ce731f2 100644 (file)
@@ -82,6 +82,8 @@ namespace Portable
       const std::vector<unsigned int>     &lexicographic_inv;
       std::vector<types::global_dof_index> lexicographic_dof_indices;
       const unsigned int                   fe_degree;
+      const unsigned int                   n_components;
+      const unsigned int                   scalar_dofs_per_cell;
       const unsigned int                   dofs_per_cell;
       const unsigned int                   q_points_per_cell;
       const UpdateFlags                   &update_flags;
@@ -109,6 +111,8 @@ namespace Portable
                     update_values | update_gradients | update_JxW_values)
       , lexicographic_inv(shape_info.lexicographic_numbering)
       , fe_degree(data->fe_degree)
+      , n_components(data->n_components)
+      , scalar_dofs_per_cell(data->scalar_dofs_per_cell)
       , dofs_per_cell(data->dofs_per_cell)
       , q_points_per_cell(data->q_points_per_cell)
       , update_flags(update_flags)
@@ -180,7 +184,7 @@ namespace Portable
             Kokkos::view_alloc("JxW_" + std::to_string(color),
                                Kokkos::WithoutInitializing),
             n_cells,
-            dofs_per_cell);
+            scalar_dofs_per_cell);
 
       if (update_flags & update_gradients)
         data->inv_jacobian[color] =
@@ -188,7 +192,7 @@ namespace Portable
             Kokkos::view_alloc("inv_jacobian_" + std::to_string(color),
                                Kokkos::WithoutInitializing),
             n_cells,
-            dofs_per_cell);
+            scalar_dofs_per_cell);
 
       // Initialize to zero, i.e., unconstrained cell
       data->constraint_mask[color] =
@@ -354,13 +358,13 @@ namespace Portable
     {
       using TeamHandle = Kokkos::TeamPolicy<
         MemorySpace::Default::kokkos_space::execution_space>::member_type;
-      using SharedView1D =
-        Kokkos::View<Number *,
+      using SharedViewValues =
+        Kokkos::View<Number **,
                      MemorySpace::Default::kokkos_space::execution_space::
                        scratch_memory_space,
                      Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
-      using SharedView2D =
-        Kokkos::View<Number *[dim],
+      using SharedViewGradients =
+        Kokkos::View<Number ***,
                      MemorySpace::Default::kokkos_space::execution_space::
                        scratch_memory_space,
                      Kokkos::MemoryTraits<Kokkos::Unmanaged>>;
@@ -386,8 +390,11 @@ namespace Portable
       size_t
       team_shmem_size(int /*team_size*/) const
       {
-        return SharedView1D::shmem_size(Functor::n_local_dofs) +
-               SharedView2D::shmem_size(Functor::n_local_dofs);
+        return SharedViewValues::shmem_size(Functor::n_local_dofs,
+                                            gpu_data.n_components) +
+               SharedViewGradients::shmem_size(Functor::n_local_dofs,
+                                               dim,
+                                               gpu_data.n_components);
       }
 
 
@@ -396,8 +403,13 @@ namespace Portable
       operator()(const TeamHandle &team_member) const
       {
         // Get the scratch memory
-        SharedView1D values(team_member.team_shmem(), Functor::n_local_dofs);
-        SharedView2D gradients(team_member.team_shmem(), Functor::n_local_dofs);
+        SharedViewValues    values(team_member.team_shmem(),
+                                Functor::n_local_dofs,
+                                gpu_data.n_components);
+        SharedViewGradients gradients(team_member.team_shmem(),
+                                      Functor::n_local_dofs,
+                                      dim,
+                                      gpu_data.n_components);
 
         SharedData<dim, Number> shared_data(team_member, values, gradients);
         func(team_member.league_rank(), &gpu_data, &shared_data, src, dst);
@@ -504,6 +516,7 @@ namespace Portable
     data_copy.co_shape_gradients = co_shape_gradients;
     data_copy.constraint_weights = constraint_weights;
     data_copy.n_cells            = n_cells[color];
+    data_copy.n_components       = n_components;
     data_copy.padding_length     = padding_length;
     data_copy.row_start          = row_start[color];
     data_copy.use_coloring       = use_coloring;
@@ -716,8 +729,10 @@ namespace Portable
     padding_length = 1 << static_cast<unsigned int>(
                        std::ceil(dim * std::log2(fe_degree + 1.)));
 
-    dofs_per_cell     = fe.n_dofs_per_cell();
-    q_points_per_cell = Utilities::fixed_power<dim>(n_q_points_1d);
+    dofs_per_cell        = fe.n_dofs_per_cell();
+    n_components         = fe.n_components();
+    scalar_dofs_per_cell = dofs_per_cell / n_components;
+    q_points_per_cell    = Utilities::fixed_power<dim>(n_q_points_1d);
 
     ::dealii::internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info(quad,
                                                                           fe);
diff --git a/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc b/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.cc
new file mode 100644 (file)
index 0000000..10639ff
--- /dev/null
@@ -0,0 +1,383 @@
+// ------------------------------------------------------------------------
+//
+// SPDX-License-Identifier: LGPL-2.1-or-later
+// Copyright (C) 2019 - 2024 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// Part of the source code is dual licensed under Apache-2.0 WITH
+// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
+// governing the source code and code contributions can be found in
+// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
+//
+// ------------------------------------------------------------------------
+
+
+
+// Test vector valued FEM for Portable::MatrixFree and Porable::FEEvaluation.
+
+#include <deal.II/base/convergence_table.h>
+
+#include <deal.II/distributed/fully_distributed_tria.h>
+#include <deal.II/distributed/repartitioning_policy_tools.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+#include <deal.II/matrix_free/portable_fe_evaluation.h>
+#include <deal.II/matrix_free/portable_matrix_free.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+#include "Kokkos_Core.hpp"
+
+template <int dim,
+          int fe_degree,
+          int n_components,
+          typename Number,
+          typename MemorySpace>
+class LaplaceOperator;
+
+template <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperator<dim, fe_degree, n_components, Number, MemorySpace::Host>
+{
+public:
+  using VectorType =
+    LinearAlgebra::distributed::Vector<Number, MemorySpace::Host>;
+
+  LaplaceOperator(const unsigned int version)
+    : version(version)
+  {}
+
+  void
+  reinit(const Mapping<dim>              &mapping,
+         const DoFHandler<dim>           &dof_handler,
+         const AffineConstraints<Number> &constraints,
+         const Quadrature<1>             &quadrature)
+  {
+    typename MatrixFree<dim, Number>::AdditionalData additional_data;
+    additional_data.mapping_update_flags = update_gradients;
+
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, quadrature, additional_data);
+  }
+
+  void
+  initialize_dof_vector(VectorType &vec) const
+  {
+    matrix_free.initialize_dof_vector(vec);
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    matrix_free.cell_loop(&LaplaceOperator::local_apply, this, dst, src, true);
+  }
+
+private:
+  void
+  local_apply(const MatrixFree<dim, Number>               &data,
+              VectorType                                  &dst,
+              const VectorType                            &src,
+              const std::pair<unsigned int, unsigned int> &cell_range) const
+  {
+    FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, Number> phi(data);
+    for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell)
+      {
+        phi.reinit(cell);
+        phi.read_dof_values_plain(src);
+
+        if (version == 0)
+          phi.evaluate(EvaluationFlags::values);
+        else if (version == 1)
+          phi.evaluate(EvaluationFlags::gradients);
+        else if (version == 2)
+          phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+
+        for (unsigned int q = 0; q < phi.n_q_points; ++q)
+          {
+            if (version == 0 || version == 2)
+              phi.submit_value(phi.get_value(q), q);
+            if (version == 1 || version == 2)
+              phi.submit_gradient(phi.get_gradient(q), q);
+          }
+
+        if (version == 0)
+          phi.integrate(EvaluationFlags::values);
+        else if (version == 1)
+          phi.integrate(EvaluationFlags::gradients);
+        else if (version == 2)
+          phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+
+        phi.distribute_local_to_global(dst);
+      }
+  }
+
+  const unsigned int      version;
+  MatrixFree<dim, Number> matrix_free;
+};
+
+
+
+template <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperatorQuad
+{
+public:
+  DEAL_II_HOST_DEVICE
+  LaplaceOperatorQuad(const unsigned int version)
+    : version(version)
+  {}
+
+  DEAL_II_HOST_DEVICE void
+  operator()(
+    Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, Number>
+             *phi,
+    const int q) const
+  {
+    if (version == 0 || version == 2)
+      phi->submit_value(phi->get_value(q), q);
+    if (version == 1 || version == 2)
+      phi->submit_gradient(phi->get_gradient(q), q);
+  }
+
+private:
+  const unsigned int version;
+};
+
+template <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperatorLocal
+{
+public:
+  LaplaceOperatorLocal(const unsigned int version)
+    : version(version)
+  {}
+
+  DEAL_II_HOST_DEVICE void
+  operator()(const unsigned int                                      cell,
+             const typename Portable::MatrixFree<dim, Number>::Data *gpu_data,
+             Portable::SharedData<dim, Number> *shared_data,
+             const Number                      *src,
+             Number                            *dst) const
+  {
+    (void)cell; // TODO?
+
+    Portable::FEEvaluation<dim, fe_degree, fe_degree + 1, n_components, Number>
+      phi(
+        /*cell,*/ gpu_data, shared_data);
+    phi.read_dof_values(src);
+
+    if (version == 0)
+      phi.evaluate(EvaluationFlags::values);
+    else if (version == 1)
+      phi.evaluate(EvaluationFlags::gradients);
+    else if (version == 2)
+      phi.evaluate(EvaluationFlags::values | EvaluationFlags::gradients);
+
+    phi.apply_for_each_quad_point(
+      LaplaceOperatorQuad<dim, fe_degree, n_components, Number>(version));
+
+    if (version == 0)
+      phi.integrate(EvaluationFlags::values);
+    else if (version == 1)
+      phi.integrate(EvaluationFlags::gradients);
+    else if (version == 2)
+      phi.integrate(EvaluationFlags::values | EvaluationFlags::gradients);
+
+    phi.distribute_local_to_global(dst);
+  }
+  static const unsigned int n_dofs_1d    = fe_degree + 1;
+  static const unsigned int n_local_dofs = Utilities::pow(fe_degree + 1, dim);
+  static const unsigned int n_q_points   = Utilities::pow(fe_degree + 1, dim);
+
+private:
+  const unsigned int version;
+};
+
+template <int dim, int fe_degree, int n_components, typename Number>
+class LaplaceOperator<dim,
+                      fe_degree,
+                      n_components,
+                      Number,
+                      MemorySpace::Default>
+{
+public:
+  using VectorType =
+    LinearAlgebra::distributed::Vector<Number, MemorySpace::Default>;
+
+  LaplaceOperator(const unsigned int version)
+    : version(version)
+  {}
+
+  void
+  reinit(const Mapping<dim>              &mapping,
+         const DoFHandler<dim>           &dof_handler,
+         const AffineConstraints<Number> &constraints,
+         const Quadrature<1>             &quadrature)
+  {
+    typename Portable::MatrixFree<dim, Number>::AdditionalData additional_data;
+    additional_data.mapping_update_flags = update_JxW_values | update_gradients;
+
+    matrix_free.reinit(
+      mapping, dof_handler, constraints, quadrature, additional_data);
+  }
+
+  void
+  initialize_dof_vector(VectorType &vec) const
+  {
+    matrix_free.initialize_dof_vector(vec);
+  }
+
+  void
+  vmult(VectorType &dst, const VectorType &src) const
+  {
+    dst = 0.0; // TODO: annoying
+    LaplaceOperatorLocal<dim, fe_degree, n_components, Number> local_operator(
+      version);
+    matrix_free.cell_loop(local_operator, src, dst);
+    matrix_free.copy_constrained_values(src, dst); // TODO: annoying
+  }
+
+private:
+  const unsigned int                version;
+  Portable::MatrixFree<dim, Number> matrix_free;
+};
+
+
+
+template <int dim, typename T>
+class AnalyticalFunction : public Function<dim, T>
+{
+public:
+  AnalyticalFunction(const unsigned int n_components)
+    : Function<dim, T>(n_components)
+  {}
+
+  virtual T
+  value(const Point<dim, T> &p, const unsigned int component = 0) const override
+  {
+    double temp = 0.0;
+
+    for (unsigned int d = 0; d < dim; ++d)
+      temp += std::sin(p[d]);
+
+    return temp * (1.0 + component);
+  }
+};
+
+
+
+template <unsigned int dim,
+          const int    degree,
+          int          n_components,
+          typename MemorySpace>
+void
+run(const unsigned int n_refinements, ConvergenceTable &table)
+{
+  using Number     = double;
+  using VectorType = LinearAlgebra::distributed::Vector<Number, MemorySpace>;
+
+  Triangulation<dim> tria;
+
+  GridGenerator::hyper_cube(tria);
+  tria.refine_global(n_refinements);
+
+  const MappingQ1<dim> mapping;
+  const FE_Q<dim>      fe_q(degree);
+  const FESystem<dim>  fe(fe_q, n_components);
+  const QGauss<dim>    quadrature(degree + 1);
+
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  AffineConstraints<Number> constraints;
+  DoFTools::make_zero_boundary_constraints(dof_handler, constraints);
+  constraints.close();
+
+  for (unsigned int i = 0; i < 3; ++i)
+    {
+      LaplaceOperator<dim, degree, n_components, Number, MemorySpace>
+        laplace_operator(i);
+
+      laplace_operator.reinit(mapping,
+                              dof_handler,
+                              constraints,
+                              quadrature.get_tensor_basis()[0]);
+
+      VectorType src, dst;
+
+      laplace_operator.initialize_dof_vector(src);
+      laplace_operator.initialize_dof_vector(dst);
+
+      {
+        LinearAlgebra::distributed::Vector<Number> src_host(
+          src.get_partitioner());
+
+        VectorTools::create_right_hand_side<dim, dim>(
+          mapping,
+          dof_handler,
+          quadrature,
+          AnalyticalFunction<dim, Number>(n_components),
+          src_host,
+          constraints);
+
+        LinearAlgebra::ReadWriteVector<Number> rw_vector(
+          src.get_partitioner()->locally_owned_range());
+        rw_vector.import(src_host, VectorOperation::insert);
+        src.import(rw_vector, VectorOperation::insert);
+
+        dst = 0.0;
+      }
+
+      laplace_operator.vmult(dst, src);
+
+      table.add_value("fe_degree", degree);
+      table.add_value("n_refinements", n_refinements);
+      table.add_value("n_components", n_components);
+      table.add_value("n_dofs", dof_handler.n_dofs());
+
+      if (std::is_same_v<MemorySpace, dealii::MemorySpace::Host>)
+        table.add_value("version", "host");
+      else
+        table.add_value("version", "default");
+
+      table.add_value("norm", dst.l2_norm());
+      table.set_scientific("norm", true);
+    }
+}
+
+int
+main()
+{
+  initlog();
+  Kokkos::initialize();
+
+  const unsigned int dim           = 2;
+  const unsigned int fe_degree     = 3;
+  unsigned int       n_refinements = 3;
+
+  ConvergenceTable table;
+
+  run<dim, fe_degree, 1, MemorySpace::Host>(n_refinements, table);
+  run<dim, fe_degree, 1, MemorySpace::Default>(n_refinements, table);
+  run<dim, fe_degree, dim, MemorySpace::Host>(n_refinements, table);
+  run<dim, fe_degree, dim, MemorySpace::Default>(n_refinements, table);
+  run<dim, fe_degree, dim + 1, MemorySpace::Host>(n_refinements, table);
+  run<dim, fe_degree, dim + 1, MemorySpace::Default>(n_refinements, table);
+
+  table.write_text(deallog.get_file_stream());
+
+  Kokkos::finalize();
+}
diff --git a/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output b/tests/matrix_free_kokkos/matrix_vector_host_device_multi_component.output
new file mode 100644 (file)
index 0000000..d06c47d
--- /dev/null
@@ -0,0 +1,20 @@
+
+fe_degree n_refinements n_components n_dofs version    norm    
+3         3             1            625    host    1.1706e-04 
+3         3             1            625    host    1.2457e-01 
+3         3             1            625    host    1.2465e-01 
+3         3             1            625    default 1.1706e-04 
+3         3             1            625    default 1.2457e-01 
+3         3             1            625    default 1.2465e-01 
+3         3             2            1250   host    2.6174e-04 
+3         3             2            1250   host    2.7854e-01 
+3         3             2            1250   host    2.7872e-01 
+3         3             2            1250   default 2.6174e-04 
+3         3             2            1250   default 2.7854e-01 
+3         3             2            1250   default 2.7872e-01 
+3         3             3            1875   host    4.3798e-04 
+3         3             3            1875   host    4.6609e-01 
+3         3             3            1875   host    4.6639e-01 
+3         3             3            1875   default 4.3798e-04 
+3         3             3            1875   default 4.6609e-01 
+3         3             3            1875   default 4.6639e-01 

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.