]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fill CUDAWrappers::internal::constraint_weights inside of ReinitHelper 12571/head
authorPeter Munch <peterrmuench@gmail.com>
Fri, 16 Jul 2021 16:21:01 +0000 (18:21 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 19 Jul 2021 16:56:41 +0000 (18:56 +0200)
include/deal.II/matrix_free/cuda_matrix_free.templates.h
include/deal.II/matrix_free/hanging_nodes_internal.h

index ccfc52a9a1dccfe965c7bf529182973d4b9c4da0..c1db385d10a67682a4dcd6d63790e156d50f0323 100644 (file)
@@ -271,6 +271,11 @@ namespace CUDAWrappers
       , padding_length(data->get_padding_length())
       , hanging_nodes(fe_degree, dof_handler, lexicographic_inv)
     {
+      cudaError_t error_code = cudaMemcpyToSymbol(
+        constraint_weights,
+        shape_info.data.front().subface_interpolation_matrix.data(),
+        sizeof(double) * fe.n_dofs_per_face(0) * fe.n_dofs_per_face(0));
+
       local_dof_indices.resize(data->dofs_per_cell);
       lexicographic_dof_indices.resize(dofs_per_cell);
     }
index 23a60d61c347ffab7da7eeb726b2ddc96e8896b1..0dc758cee96691e5724efd1df2a22db8ce631579 100644 (file)
 
 DEAL_II_NAMESPACE_OPEN
 
-#ifdef DEAL_II_COMPILER_CUDA_AWARE
-namespace CUDAWrappers
-{
-  namespace internal
-  {
-    __constant__ double
-      constraint_weights[(CUDAWrappers::mf_max_elem_degree + 1) *
-                         (CUDAWrappers::mf_max_elem_degree + 1)];
-
-    template <int dim>
-    void
-    setup_constraint_weigths(unsigned int fe_degree)
-    {
-      const unsigned int face_no =
-        0; // we assume that all faces have the same number of dofs
-
-      FE_Q<2>            fe_q(fe_degree);
-      FullMatrix<double> interpolation_matrix(fe_q.n_dofs_per_face(face_no),
-                                              fe_q.n_dofs_per_face(face_no));
-      fe_q.get_subface_interpolation_matrix(fe_q,
-                                            0,
-                                            interpolation_matrix,
-                                            face_no);
-
-      std::vector<unsigned int> mapping =
-        FETools::lexicographic_to_hierarchic_numbering<1>(fe_degree);
-
-      FullMatrix<double> mapped_matrix(fe_q.n_dofs_per_face(face_no),
-                                       fe_q.n_dofs_per_face(face_no));
-      for (unsigned int i = 0; i < fe_q.n_dofs_per_face(face_no); ++i)
-        for (unsigned int j = 0; j < fe_q.n_dofs_per_face(face_no); ++j)
-          mapped_matrix(i, j) = interpolation_matrix(mapping[i], mapping[j]);
-
-      cudaError_t error_code =
-        cudaMemcpyToSymbol(constraint_weights,
-                           &mapped_matrix[0][0],
-                           sizeof(double) * fe_q.n_dofs_per_face(face_no) *
-                             fe_q.n_dofs_per_face(face_no));
-      AssertCuda(error_code);
-    }
-  } // namespace internal
-} // namespace CUDAWrappers
-
-#endif
-
-
-
 namespace internal
 {
   namespace MatrixFreeFunctions
@@ -194,15 +147,6 @@ namespace internal
     {
       // Set up line-to-cell mapping for edge constraints (only if dim = 3)
       setup_line_to_cell();
-
-#ifndef DEAL_II_COMPILER_CUDA_AWARE
-      Assert(false,
-             ExcMessage(
-               "This function needs DEAL_II_COMPILER_CUDA_AWARE to enabled!"));
-      (void)fe_degree;
-#else
-      CUDAWrappers::internal::setup_constraint_weigths<dim>(fe_degree);
-#endif
     }
 
 
@@ -752,251 +696,248 @@ namespace CUDAWrappers
 {
   namespace internal
   {
-    //--------------------------------------------------------------------------//
-    // Functions for resolving the hanging node constraints on the GPU //
-    //--------------------------------------------------------------------------//
-    namespace internal
+    __constant__ double
+      constraint_weights[(CUDAWrappers::mf_max_elem_degree + 1) *
+                         (CUDAWrappers::mf_max_elem_degree + 1)];
+
+    //------------------------------------------------------------------------//
+    // Functions for resolving the hanging node constraints on the GPU        //
+    //------------------------------------------------------------------------//
+    template <unsigned int size>
+    __device__ inline unsigned int
+    index2(unsigned int i, unsigned int j)
+    {
+      return i + size * j;
+    }
+
+
+
+    template <unsigned int size>
+    __device__ inline unsigned int
+    index3(unsigned int i, unsigned int j, unsigned int k)
+    {
+      return i + size * j + size * size * k;
+    }
+
+
+
+    template <unsigned int fe_degree,
+              unsigned int direction,
+              bool         transpose,
+              typename Number>
+    __device__ inline void
+    interpolate_boundary_2d(const unsigned int constraint_mask, Number *values)
+    {
+      const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
+      const unsigned int y_idx = threadIdx.y;
+
+      const unsigned int this_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y;
+
+      const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
+
+      Number t = 0;
+      // Flag is true if dof is constrained for the given direction and the
+      // given face.
+      const bool constrained_face =
+        (constraint_mask &
+         (((direction == 0) ?
+             dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y :
+             0) |
+          ((direction == 1) ?
+             dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x :
+             0)));
+
+      // Flag is true if for the given direction, the dof is constrained with
+      // the right type and is on the correct side (left (= 0) or right (=
+      // fe_degree))
+      const bool constrained_dof =
+        ((direction == 0) &&
+         ((constraint_mask &
+           dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y) ?
+            (y_idx == 0) :
+            (y_idx == fe_degree))) ||
+        ((direction == 1) &&
+         ((constraint_mask &
+           dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x) ?
+            (x_idx == 0) :
+            (x_idx == fe_degree)));
+
+      if (constrained_face && constrained_dof)
+        {
+          const bool type = constraint_mask & this_type;
+
+          if (type)
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
+                                       index2<fe_degree + 1>(x_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[i * (fe_degree + 1) + interp_idx] :
+                      constraint_weights[interp_idx * (fe_degree + 1) + i];
+                  t += w * values[real_idx];
+                }
+            }
+          else
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
+                                       index2<fe_degree + 1>(x_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[(fe_degree - i) * (fe_degree + 1) +
+                                         fe_degree - interp_idx] :
+                      constraint_weights[(fe_degree - interp_idx) *
+                                           (fe_degree + 1) +
+                                         fe_degree - i];
+                  t += w * values[real_idx];
+                }
+            }
+        }
+
+      // The synchronization is done for all the threads in one block with
+      // each block being assigned to one element.
+      __syncthreads();
+      if (constrained_face && constrained_dof)
+        values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
+
+      __syncthreads();
+    }
+
+
+
+    template <unsigned int fe_degree,
+              unsigned int direction,
+              bool         transpose,
+              typename Number>
+    __device__ inline void
+    interpolate_boundary_3d(const unsigned int constraint_mask, Number *values)
     {
-      template <unsigned int size>
-      __device__ inline unsigned int
-      index2(unsigned int i, unsigned int j)
-      {
-        return i + size * j;
-      }
-
-
-
-      template <unsigned int size>
-      __device__ inline unsigned int
-      index3(unsigned int i, unsigned int j, unsigned int k)
-      {
-        return i + size * j + size * size * k;
-      }
-
-
-
-      template <unsigned int fe_degree,
-                unsigned int direction,
-                bool         transpose,
-                typename Number>
-      __device__ inline void
-      interpolate_boundary_2d(const unsigned int constraint_mask,
-                              Number *           values)
-      {
-        const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
-        const unsigned int y_idx = threadIdx.y;
-
-        const unsigned int this_type =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y;
-
-        const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx;
-
-        Number t = 0;
-        // Flag is true if dof is constrained for the given direction and the
-        // given face.
-        const bool constrained_face =
-          (constraint_mask &
-           (((direction == 0) ?
-               dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y :
-               0) |
-            ((direction == 1) ?
-               dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x :
-               0)));
-
-        // Flag is true if for the given direction, the dof is constrained with
-        // the right type and is on the correct side (left (= 0) or right (=
-        // fe_degree))
-        const bool constrained_dof =
-          ((direction == 0) &&
-           ((constraint_mask &
-             dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y) ?
-              (y_idx == 0) :
-              (y_idx == fe_degree))) ||
-          ((direction == 1) &&
-           ((constraint_mask &
-             dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x) ?
-              (x_idx == 0) :
-              (x_idx == fe_degree)));
-
-        if (constrained_face && constrained_dof)
-          {
-            const bool type = constraint_mask & this_type;
-
-            if (type)
-              {
-                for (unsigned int i = 0; i <= fe_degree; ++i)
-                  {
-                    const unsigned int real_idx =
-                      (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
-                                         index2<fe_degree + 1>(x_idx, i);
-
-                    const Number w =
-                      transpose ?
-                        constraint_weights[i * (fe_degree + 1) + interp_idx] :
-                        constraint_weights[interp_idx * (fe_degree + 1) + i];
-                    t += w * values[real_idx];
-                  }
-              }
-            else
-              {
-                for (unsigned int i = 0; i <= fe_degree; ++i)
-                  {
-                    const unsigned int real_idx =
-                      (direction == 0) ? index2<fe_degree + 1>(i, y_idx) :
-                                         index2<fe_degree + 1>(x_idx, i);
-
-                    const Number w =
-                      transpose ?
-                        constraint_weights[(fe_degree - i) * (fe_degree + 1) +
-                                           fe_degree - interp_idx] :
-                        constraint_weights[(fe_degree - interp_idx) *
-                                             (fe_degree + 1) +
-                                           fe_degree - i];
-                    t += w * values[real_idx];
-                  }
-              }
-          }
-
-        // The synchronization is done for all the threads in one block with
-        // each block being assigned to one element.
-        __syncthreads();
-        if (constrained_face && constrained_dof)
-          values[index2<fe_degree + 1>(x_idx, y_idx)] = t;
-
-        __syncthreads();
-      }
-
-
-
-      template <unsigned int fe_degree,
-                unsigned int direction,
-                bool         transpose,
-                typename Number>
-      __device__ inline void
-      interpolate_boundary_3d(const unsigned int constraint_mask,
-                              Number *           values)
-      {
-        const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
-        const unsigned int y_idx = threadIdx.y;
-        const unsigned int z_idx = threadIdx.z;
-
-        const unsigned int this_type =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x :
-            (direction == 1) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_z;
-        const unsigned int face1_type =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y :
-            (direction == 1) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_z :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x;
-        const unsigned int face2_type =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_z :
-            (direction == 1) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y;
-
-        // If computing in x-direction, need to match against face_y or
-        // face_z
-        const unsigned int face1 =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y :
-            (direction == 1) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_z :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x;
-        const unsigned int face2 =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_z :
-            (direction == 1) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y;
-        const unsigned int edge =
-          (direction == 0) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::edge_yz :
-            (direction == 1) ?
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::edge_zx :
-            dealii::internal::MatrixFreeFunctions::ConstraintTypes::edge_xy;
-        const unsigned int constrained_face =
-          constraint_mask & (face1 | face2 | edge);
-
-        const unsigned int interp_idx =
-          (direction == 0) ? x_idx : (direction == 1) ? y_idx : z_idx;
-        const unsigned int face1_idx =
-          (direction == 0) ? y_idx : (direction == 1) ? z_idx : x_idx;
-        const unsigned int face2_idx =
-          (direction == 0) ? z_idx : (direction == 1) ? x_idx : y_idx;
-
-        Number     t        = 0;
-        const bool on_face1 = (constraint_mask & face1_type) ?
-                                (face1_idx == 0) :
-                                (face1_idx == fe_degree);
-        const bool on_face2 = (constraint_mask & face2_type) ?
-                                (face2_idx == 0) :
-                                (face2_idx == fe_degree);
-        const bool constrained_dof =
-          (((constraint_mask & face1) && on_face1) ||
-           ((constraint_mask & face2) && on_face2) ||
-           ((constraint_mask & edge) && on_face1 && on_face2));
-
-        if (constrained_face && constrained_dof)
-          {
-            const bool type = constraint_mask & this_type;
-            if (type)
-              {
-                for (unsigned int i = 0; i <= fe_degree; ++i)
-                  {
-                    const unsigned int real_idx =
-                      (direction == 0) ?
-                        index3<fe_degree + 1>(i, y_idx, z_idx) :
-                        (direction == 1) ?
-                        index3<fe_degree + 1>(x_idx, i, z_idx) :
-                        index3<fe_degree + 1>(x_idx, y_idx, i);
-
-                    const Number w =
-                      transpose ?
-                        constraint_weights[i * (fe_degree + 1) + interp_idx] :
-                        constraint_weights[interp_idx * (fe_degree + 1) + i];
-                    t += w * values[real_idx];
-                  }
-              }
-            else
-              {
-                for (unsigned int i = 0; i <= fe_degree; ++i)
-                  {
-                    const unsigned int real_idx =
-                      (direction == 0) ?
-                        index3<fe_degree + 1>(i, y_idx, z_idx) :
-                        (direction == 1) ?
-                        index3<fe_degree + 1>(x_idx, i, z_idx) :
-                        index3<fe_degree + 1>(x_idx, y_idx, i);
-
-                    const Number w =
-                      transpose ?
-                        constraint_weights[(fe_degree - i) * (fe_degree + 1) +
-                                           fe_degree - interp_idx] :
-                        constraint_weights[(fe_degree - interp_idx) *
-                                             (fe_degree + 1) +
-                                           fe_degree - i];
-                    t += w * values[real_idx];
-                  }
-              }
-          }
-
-        // The synchronization is done for all the threads in one block with
-        // each block being assigned to one element.
-        __syncthreads();
-
-        if (constrained_face && constrained_dof)
-          values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
-
-        __syncthreads();
-      }
-    } // namespace internal
+      const unsigned int x_idx = threadIdx.x % (fe_degree + 1);
+      const unsigned int y_idx = threadIdx.y;
+      const unsigned int z_idx = threadIdx.z;
+
+      const unsigned int this_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x :
+          (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_z;
+      const unsigned int face1_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y :
+          (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_z :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x;
+      const unsigned int face2_type =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_z :
+          (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_x :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::type_y;
+
+      // If computing in x-direction, need to match against face_y or
+      // face_z
+      const unsigned int face1 =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y :
+          (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_z :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x;
+      const unsigned int face2 =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_z :
+          (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_x :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::face_y;
+      const unsigned int edge =
+        (direction == 0) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::edge_yz :
+          (direction == 1) ?
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::edge_zx :
+          dealii::internal::MatrixFreeFunctions::ConstraintTypes::edge_xy;
+      const unsigned int constrained_face =
+        constraint_mask & (face1 | face2 | edge);
+
+      const unsigned int interp_idx =
+        (direction == 0) ? x_idx : (direction == 1) ? y_idx : z_idx;
+      const unsigned int face1_idx =
+        (direction == 0) ? y_idx : (direction == 1) ? z_idx : x_idx;
+      const unsigned int face2_idx =
+        (direction == 0) ? z_idx : (direction == 1) ? x_idx : y_idx;
+
+      Number     t        = 0;
+      const bool on_face1 = (constraint_mask & face1_type) ?
+                              (face1_idx == 0) :
+                              (face1_idx == fe_degree);
+      const bool on_face2 = (constraint_mask & face2_type) ?
+                              (face2_idx == 0) :
+                              (face2_idx == fe_degree);
+      const bool constrained_dof =
+        (((constraint_mask & face1) && on_face1) ||
+         ((constraint_mask & face2) && on_face2) ||
+         ((constraint_mask & edge) && on_face1 && on_face2));
+
+      if (constrained_face && constrained_dof)
+        {
+          const bool type = constraint_mask & this_type;
+          if (type)
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
+                                       (direction == 1) ?
+                                       index3<fe_degree + 1>(x_idx, i, z_idx) :
+                                       index3<fe_degree + 1>(x_idx, y_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[i * (fe_degree + 1) + interp_idx] :
+                      constraint_weights[interp_idx * (fe_degree + 1) + i];
+                  t += w * values[real_idx];
+                }
+            }
+          else
+            {
+              for (unsigned int i = 0; i <= fe_degree; ++i)
+                {
+                  const unsigned int real_idx =
+                    (direction == 0) ? index3<fe_degree + 1>(i, y_idx, z_idx) :
+                                       (direction == 1) ?
+                                       index3<fe_degree + 1>(x_idx, i, z_idx) :
+                                       index3<fe_degree + 1>(x_idx, y_idx, i);
+
+                  const Number w =
+                    transpose ?
+                      constraint_weights[(fe_degree - i) * (fe_degree + 1) +
+                                         fe_degree - interp_idx] :
+                      constraint_weights[(fe_degree - interp_idx) *
+                                           (fe_degree + 1) +
+                                         fe_degree - i];
+                  t += w * values[real_idx];
+                }
+            }
+        }
+
+      // The synchronization is done for all the threads in one block with
+      // each block being assigned to one element.
+      __syncthreads();
+
+      if (constrained_face && constrained_dof)
+        values[index3<fe_degree + 1>(x_idx, y_idx, z_idx)] = t;
+
+      __syncthreads();
+    }
 
 
 
@@ -1014,23 +955,23 @@ namespace CUDAWrappers
     {
       if (dim == 2)
         {
-          internal::interpolate_boundary_2d<fe_degree, 0, transpose>(
-            constraint_mask, values);
+          interpolate_boundary_2d<fe_degree, 0, transpose>(constraint_mask,
+                                                           values);
 
-          internal::interpolate_boundary_2d<fe_degree, 1, transpose>(
-            constraint_mask, values);
+          interpolate_boundary_2d<fe_degree, 1, transpose>(constraint_mask,
+                                                           values);
         }
       else if (dim == 3)
         {
           // Interpolate y and z faces (x-direction)
-          internal::interpolate_boundary_3d<fe_degree, 0, transpose>(
-            constraint_mask, values);
+          interpolate_boundary_3d<fe_degree, 0, transpose>(constraint_mask,
+                                                           values);
           // Interpolate x and z faces (y-direction)
-          internal::interpolate_boundary_3d<fe_degree, 1, transpose>(
-            constraint_mask, values);
+          interpolate_boundary_3d<fe_degree, 1, transpose>(constraint_mask,
+                                                           values);
           // Interpolate x and y faces (z-direction)
-          internal::interpolate_boundary_3d<fe_degree, 2, transpose>(
-            constraint_mask, values);
+          interpolate_boundary_3d<fe_degree, 2, transpose>(constraint_mask,
+                                                           values);
         }
     }
   } // namespace internal

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.