From: Bruno Turcksin Date: Thu, 16 Sep 2021 20:53:31 +0000 (+0000) Subject: Refactor hanging_nodes_internals.h X-Git-Tag: v9.4.0-rc1~983^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=19105e31ecb29751496ec7c59da0580fdc1ba5d0;p=dealii.git Refactor hanging_nodes_internals.h --- diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index 965ff1c27e..d913082f36 100644 --- a/include/deal.II/matrix_free/cuda_fe_evaluation.h +++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h @@ -26,10 +26,10 @@ # include # include +# include # include # include # include -# include DEAL_II_NAMESPACE_OPEN diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h new file mode 100644 index 0000000000..d76f56fa25 --- /dev/null +++ b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h @@ -0,0 +1,349 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii_cuda_hanging_nodes_internal_h +#define dealii_cuda_hanging_nodes_internal_h + +#include + +#ifdef DEAL_II_COMPILER_CUDA_AWARE + +# include + +# include + +DEAL_II_NAMESPACE_OPEN +namespace CUDAWrappers +{ + 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 dealii::internal::MatrixFreeFunctions::ConstraintKinds + constraint_mask, + Number *values) + { + const unsigned int x_idx = threadIdx.x % (fe_degree + 1); + const unsigned int y_idx = threadIdx.y; + + const auto this_type = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::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::ConstraintKinds::face_y : + dealii::internal::MatrixFreeFunctions::ConstraintKinds:: + unconstrained) | + ((direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x : + dealii::internal::MatrixFreeFunctions::ConstraintKinds:: + unconstrained))) != + dealii::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained; + + // 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::ConstraintKinds::type_y) != + dealii::internal::MatrixFreeFunctions::ConstraintKinds:: + unconstrained) ? + (y_idx == 0) : + (y_idx == fe_degree))) || + ((direction == 1) && + (((constraint_mask & + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x) != + dealii::internal::MatrixFreeFunctions::ConstraintKinds:: + unconstrained) ? + (x_idx == 0) : + (x_idx == fe_degree))); + + if (constrained_face && constrained_dof) + { + const bool type = (constraint_mask & this_type) != + dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained; + + 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 dealii::internal::MatrixFreeFunctions::ConstraintKinds + 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 auto this_type = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x : + (direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z; + const auto face1_type = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y : + (direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x; + const auto face2_type = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z : + (direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y; + + // If computing in x-direction, need to match against face_y or + // face_z + const auto face1 = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y : + (direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x; + const auto face2 = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z : + (direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y; + const auto edge = + (direction == 0) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_yz : + (direction == 1) ? + dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_zx : + dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_xy; + const auto 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) != + dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) ? + (face1_idx == 0) : + (face1_idx == fe_degree); + const bool on_face2 = ((constraint_mask & face2_type) != + dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) ? + (face2_idx == 0) : + (face2_idx == fe_degree); + const bool constrained_dof = + ((((constraint_mask & face1) != dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) && + on_face1) || + (((constraint_mask & face2) != dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) && + on_face2) || + (((constraint_mask & edge) != dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) && + on_face1 && on_face2)); + + if ((constrained_face != dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) && + constrained_dof) + { + const bool type = (constraint_mask & this_type) != + dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained; + 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 != dealii::internal::MatrixFreeFunctions:: + ConstraintKinds::unconstrained) && + constrained_dof) + values[index3(x_idx, y_idx, z_idx)] = t; + + __syncthreads(); + } + + + + /** + * This function resolves the hanging nodes using tensor product. + * + * The implementation of this class is explained in Section 3 of + * @cite ljungkvist2017matrix and in Section 3.4 of + * @cite kronbichler2019multigrid. + */ + template + __device__ void + resolve_hanging_nodes( + const dealii::internal::MatrixFreeFunctions::ConstraintKinds + constraint_mask, + Number *values) + { + if (dim == 2) + { + interpolate_boundary_2d(constraint_mask, + values); + + interpolate_boundary_2d(constraint_mask, + values); + } + else if (dim == 3) + { + // Interpolate y and z faces (x-direction) + interpolate_boundary_3d(constraint_mask, + values); + // Interpolate x and z faces (y-direction) + interpolate_boundary_3d(constraint_mask, + values); + // Interpolate x and y faces (z-direction) + interpolate_boundary_3d(constraint_mask, + values); + } + } + } // namespace internal +} // namespace CUDAWrappers + +DEAL_II_NAMESPACE_CLOSE +#endif + +#endif diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h index 2fa34f40c1..31f322fa18 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.h @@ -38,11 +38,19 @@ # include # include -# include DEAL_II_NAMESPACE_OPEN +// Forward declaration +namespace internal +{ + namespace MatrixFreeFunctions + { + enum class ConstraintKinds : std::uint16_t; + } +} // namespace internal + namespace CUDAWrappers { // forward declaration 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 1bd8d75285..5039981471 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -23,6 +23,7 @@ #ifdef DEAL_II_COMPILER_CUDA_AWARE +# include # include # include @@ -31,6 +32,7 @@ # include # include +# include # include # include diff --git a/include/deal.II/matrix_free/hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h index f13aee5beb..b92fcf2316 100644 --- a/include/deal.II/matrix_free/hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/hanging_nodes_internal.h @@ -18,7 +18,6 @@ #include -#include #include #include @@ -855,327 +854,6 @@ namespace internal } // namespace internal - -#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)]; - - //------------------------------------------------------------------------// - // 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 dealii::internal::MatrixFreeFunctions::ConstraintKinds - constraint_mask, - Number *values) - { - const unsigned int x_idx = threadIdx.x % (fe_degree + 1); - const unsigned int y_idx = threadIdx.y; - - const auto this_type = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::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::ConstraintKinds::face_y : - dealii::internal::MatrixFreeFunctions::ConstraintKinds:: - unconstrained) | - ((direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x : - dealii::internal::MatrixFreeFunctions::ConstraintKinds:: - unconstrained))) != - dealii::internal::MatrixFreeFunctions::ConstraintKinds::unconstrained; - - // 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::ConstraintKinds::type_y) != - dealii::internal::MatrixFreeFunctions::ConstraintKinds:: - unconstrained) ? - (y_idx == 0) : - (y_idx == fe_degree))) || - ((direction == 1) && - (((constraint_mask & - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x) != - dealii::internal::MatrixFreeFunctions::ConstraintKinds:: - unconstrained) ? - (x_idx == 0) : - (x_idx == fe_degree))); - - if (constrained_face && constrained_dof) - { - const bool type = (constraint_mask & this_type) != - dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained; - - 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 dealii::internal::MatrixFreeFunctions::ConstraintKinds - 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 auto this_type = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x : - (direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z; - const auto face1_type = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y : - (direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x; - const auto face2_type = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_z : - (direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_x : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::type_y; - - // If computing in x-direction, need to match against face_y or - // face_z - const auto face1 = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y : - (direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x; - const auto face2 = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_z : - (direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_x : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::face_y; - const auto edge = - (direction == 0) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_yz : - (direction == 1) ? - dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_zx : - dealii::internal::MatrixFreeFunctions::ConstraintKinds::edge_xy; - const auto 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) != - dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) ? - (face1_idx == 0) : - (face1_idx == fe_degree); - const bool on_face2 = ((constraint_mask & face2_type) != - dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) ? - (face2_idx == 0) : - (face2_idx == fe_degree); - const bool constrained_dof = - ((((constraint_mask & face1) != dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) && - on_face1) || - (((constraint_mask & face2) != dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) && - on_face2) || - (((constraint_mask & edge) != dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) && - on_face1 && on_face2)); - - if ((constrained_face != dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) && - constrained_dof) - { - const bool type = (constraint_mask & this_type) != - dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained; - 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 != dealii::internal::MatrixFreeFunctions:: - ConstraintKinds::unconstrained) && - constrained_dof) - values[index3(x_idx, y_idx, z_idx)] = t; - - __syncthreads(); - } - - - - /** - * This function resolves the hanging nodes using tensor product. - * - * The implementation of this class is explained in Section 3 of - * @cite ljungkvist2017matrix and in Section 3.4 of - * @cite kronbichler2019multigrid. - */ - template - __device__ void - resolve_hanging_nodes( - const dealii::internal::MatrixFreeFunctions::ConstraintKinds - constraint_mask, - Number *values) - { - if (dim == 2) - { - interpolate_boundary_2d(constraint_mask, - values); - - interpolate_boundary_2d(constraint_mask, - values); - } - else if (dim == 3) - { - // Interpolate y and z faces (x-direction) - interpolate_boundary_3d(constraint_mask, - values); - // Interpolate x and z faces (y-direction) - interpolate_boundary_3d(constraint_mask, - values); - // Interpolate x and y faces (z-direction) - interpolate_boundary_3d(constraint_mask, - values); - } - } - } // namespace internal -} // namespace CUDAWrappers -#endif - DEAL_II_NAMESPACE_CLOSE #endif