From a97d948ed502c2486497cae745bd0b87488e5fb6 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 16 Jul 2021 18:21:01 +0200 Subject: [PATCH] Fill CUDAWrappers::internal::constraint_weights inside of ReinitHelper --- .../matrix_free/cuda_matrix_free.templates.h | 5 + .../matrix_free/hanging_nodes_internal.h | 561 ++++++++---------- 2 files changed, 256 insertions(+), 310 deletions(-) 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 ccfc52a9a1..c1db385d10 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -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); } diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h index 23a60d61c3..0dc758cee9 100644 --- a/include/deal.II/matrix_free/hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/hanging_nodes_internal.h @@ -29,53 +29,6 @@ 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 - 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 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 mapping = - FETools::lexicographic_to_hierarchic_numbering<1>(fe_degree); - - FullMatrix 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(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 + __device__ inline unsigned int + index2(unsigned int i, unsigned int j) + { + return i + size * j; + } + + + + template + __device__ inline unsigned int + index3(unsigned int i, unsigned int j, unsigned int k) + { + return i + size * j + size * size * k; + } + + + + template + __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(i, y_idx) : + index2(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(i, y_idx) : + index2(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(x_idx, y_idx)] = t; + + __syncthreads(); + } + + + + template + __device__ inline void + interpolate_boundary_3d(const unsigned int constraint_mask, Number *values) { - template - __device__ inline unsigned int - index2(unsigned int i, unsigned int j) - { - return i + size * j; - } - - - - template - __device__ inline unsigned int - index3(unsigned int i, unsigned int j, unsigned int k) - { - return i + size * j + size * size * k; - } - - - - template - __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(i, y_idx) : - index2(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(i, y_idx) : - index2(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(x_idx, y_idx)] = t; - - __syncthreads(); - } - - - - template - __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(i, y_idx, z_idx) : - (direction == 1) ? - index3(x_idx, i, z_idx) : - index3(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(i, y_idx, z_idx) : - (direction == 1) ? - index3(x_idx, i, z_idx) : - index3(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(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(i, y_idx, z_idx) : + (direction == 1) ? + index3(x_idx, i, z_idx) : + index3(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(i, y_idx, z_idx) : + (direction == 1) ? + index3(x_idx, i, z_idx) : + index3(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(x_idx, y_idx, z_idx)] = t; + + __syncthreads(); + } @@ -1014,23 +955,23 @@ namespace CUDAWrappers { if (dim == 2) { - internal::interpolate_boundary_2d( - constraint_mask, values); + interpolate_boundary_2d(constraint_mask, + values); - internal::interpolate_boundary_2d( - constraint_mask, values); + interpolate_boundary_2d(constraint_mask, + values); } else if (dim == 3) { // Interpolate y and z faces (x-direction) - internal::interpolate_boundary_3d( - constraint_mask, values); + interpolate_boundary_3d(constraint_mask, + values); // Interpolate x and z faces (y-direction) - internal::interpolate_boundary_3d( - constraint_mask, values); + interpolate_boundary_3d(constraint_mask, + values); // Interpolate x and y faces (z-direction) - internal::interpolate_boundary_3d( - constraint_mask, values); + interpolate_boundary_3d(constraint_mask, + values); } } } // namespace internal -- 2.39.5