From 3bf9bc3d839db2098ab79d278607348525c8478f Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Thu, 15 Jul 2021 09:28:48 +0200 Subject: [PATCH] Rename cuda_hanging_nodes_internal.h to hanging_nodes_internal.h --- .../deal.II/matrix_free/cuda_fe_evaluation.h | 2 +- .../matrix_free/cuda_matrix_free.templates.h | 4 +- ...es_internal.h => hanging_nodes_internal.h} | 431 ++++++++++-------- 3 files changed, 237 insertions(+), 200 deletions(-) rename include/deal.II/matrix_free/{cuda_hanging_nodes_internal.h => hanging_nodes_internal.h} (73%) diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index db88558e34..cf4195a9dc 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_matrix_free.templates.h b/include/deal.II/matrix_free/cuda_matrix_free.templates.h index dd34ac8838..ccfc52a9a1 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -31,7 +31,7 @@ # include # include -# include +# include # include # include @@ -242,7 +242,7 @@ namespace CUDAWrappers const unsigned int q_points_per_cell; const UpdateFlags & update_flags; const unsigned int padding_length; - HangingNodes hanging_nodes; + dealii::internal::MatrixFreeFunctions::HangingNodes hanging_nodes; }; diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/hanging_nodes_internal.h similarity index 73% rename from include/deal.II/matrix_free/cuda_hanging_nodes_internal.h rename to include/deal.II/matrix_free/hanging_nodes_internal.h index ab14ffdbb2..23a60d61c3 100644 --- a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/hanging_nodes_internal.h @@ -13,27 +13,72 @@ // // --------------------------------------------------------------------- -#ifndef dealii_cuda_hanging_nodes_internal_h -#define dealii_cuda_hanging_nodes_internal_h +#ifndef dealii_hanging_nodes_internal_h +#define dealii_hanging_nodes_internal_h #include -#ifdef DEAL_II_COMPILER_CUDA_AWARE - -# include -# include +#include +#include -# include -# include +#include +#include -# include -# include +#include +#include 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 { /** * This class creates the mask used in the treatment of hanging nodes in @@ -102,73 +147,38 @@ namespace CUDAWrappers const DoFHandler & dof_handler; }; - namespace internal + // Here is the system for how we store constraint types in a binary mask. + // This is not a complete contradiction-free system, i.e., there are + // invalid states that we just assume that we never get. + + // If the mask is zero, there are no constraints. Then, there are three + // different fields with one bit per dimension. The first field determines + // the type, or the position of an element along each direction. The + // second field determines if there is a constrained face with that + // direction as normal. The last field determines if there is a + // constrained edge of a given pair of coordinate planes, but where + // neither of the corresponding faces are constrained (only valid in 3D). + + // The element is placed in the 'first position' along *-axis. These also + // determine which face is constrained. For example, in 2D, if + // face_x and type are set, then x = 0 is constrained. + + namespace ConstraintTypes { - __constant__ double - constraint_weights[(mf_max_elem_degree + 1) * (mf_max_elem_degree + 1)]; - - // Here is the system for how we store constraint types in a binary mask. - // This is not a complete contradiction-free system, i.e., there are - // invalid states that we just assume that we never get. - - // If the mask is zero, there are no constraints. Then, there are three - // different fields with one bit per dimension. The first field determines - // the type, or the position of an element along each direction. The - // second field determines if there is a constrained face with that - // direction as normal. The last field determines if there is a - // constrained edge of a given pair of coordinate planes, but where - // neither of the corresponding faces are constrained (only valid in 3D). - - // The element is placed in the 'first position' along *-axis. These also - // determine which face is constrained. For example, in 2D, if - // constr_face_x and constr_type are set, then x = 0 is constrained. - constexpr unsigned int constr_type_x = 1 << 0; - constexpr unsigned int constr_type_y = 1 << 1; - constexpr unsigned int constr_type_z = 1 << 2; + constexpr unsigned int type_x = 1 << 0; + constexpr unsigned int type_y = 1 << 1; + constexpr unsigned int type_z = 1 << 2; // Element has as a constraint at * = 0 or * = fe_degree face - constexpr unsigned int constr_face_x = 1 << 3; - constexpr unsigned int constr_face_y = 1 << 4; - constexpr unsigned int constr_face_z = 1 << 5; + constexpr unsigned int face_x = 1 << 3; + constexpr unsigned int face_y = 1 << 4; + constexpr unsigned int face_z = 1 << 5; // Element has as a constraint at * = 0 or * = fe_degree edge - constexpr unsigned int constr_edge_xy = 1 << 6; - constexpr unsigned int constr_edge_yz = 1 << 7; - constexpr unsigned int constr_edge_zx = 1 << 8; - - 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(internal::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 - + constexpr unsigned int edge_xy = 1 << 6; + constexpr unsigned int edge_yz = 1 << 7; + constexpr unsigned int edge_zx = 1 << 8; + } // namespace ConstraintTypes template @@ -185,7 +195,14 @@ namespace CUDAWrappers // Set up line-to-cell mapping for edge constraints (only if dim = 3) setup_line_to_cell(); - internal::setup_constraint_weigths(fe_degree); +#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 } @@ -261,7 +278,7 @@ namespace CUDAWrappers child->line(neighbor_line)->index(); // Now add all active cells - for (const auto cl : line_to_cells[line_idx]) + for (const auto &cl : line_to_cells[line_idx]) line_to_cells[child_line_idx].push_back(cl); } } @@ -322,19 +339,19 @@ namespace CUDAWrappers { if (face < 2) { - mask |= internal::constr_face_x; + mask |= ConstraintTypes::face_x; if (face == 0) - mask |= internal::constr_type_x; + mask |= ConstraintTypes::type_x; if (subface == 0) - mask |= internal::constr_type_y; + mask |= ConstraintTypes::type_y; } else { - mask |= internal::constr_face_y; + mask |= ConstraintTypes::face_y; if (face == 2) - mask |= internal::constr_type_y; + mask |= ConstraintTypes::type_y; if (subface == 0) - mask |= internal::constr_type_x; + mask |= ConstraintTypes::type_x; } // Reorder neighbor_dofs and copy into faceth face of @@ -379,35 +396,35 @@ namespace CUDAWrappers // YZ-plane if (face < 2) { - mask |= internal::constr_face_x; + mask |= ConstraintTypes::face_x; if (face == 0) - mask |= internal::constr_type_x; + mask |= ConstraintTypes::type_x; if (subface % 2 == 0) - mask |= internal::constr_type_y; + mask |= ConstraintTypes::type_y; if (subface / 2 == 0) - mask |= internal::constr_type_z; + mask |= ConstraintTypes::type_z; } // XZ-plane else if (face < 4) { - mask |= internal::constr_face_y; + mask |= ConstraintTypes::face_y; if (face == 2) - mask |= internal::constr_type_y; + mask |= ConstraintTypes::type_y; if (subface % 2 == 0) - mask |= internal::constr_type_z; + mask |= ConstraintTypes::type_z; if (subface / 2 == 0) - mask |= internal::constr_type_x; + mask |= ConstraintTypes::type_x; } // XY-plane else { - mask |= internal::constr_face_z; + mask |= ConstraintTypes::face_z; if (face == 4) - mask |= internal::constr_type_z; + mask |= ConstraintTypes::type_z; if (subface % 2 == 0) - mask |= internal::constr_type_x; + mask |= ConstraintTypes::type_x; if (subface / 2 == 0) - mask |= internal::constr_type_y; + mask |= ConstraintTypes::type_y; } // Offset if upper/right/back face @@ -454,54 +471,54 @@ namespace CUDAWrappers // edge mask, what is the types of the faces it belong to, and what is // the type along the edge. const unsigned int line_to_edge[12][4] = { - {internal::constr_face_x | internal::constr_face_z, - internal::constr_edge_zx, - internal::constr_type_x | internal::constr_type_z, - internal::constr_type_y}, - {internal::constr_face_x | internal::constr_face_z, - internal::constr_edge_zx, - internal::constr_type_z, - internal::constr_type_y}, - {internal::constr_face_y | internal::constr_face_z, - internal::constr_edge_yz, - internal::constr_type_y | internal::constr_type_z, - internal::constr_type_x}, - {internal::constr_face_y | internal::constr_face_z, - internal::constr_edge_yz, - internal::constr_type_z, - internal::constr_type_x}, - {internal::constr_face_x | internal::constr_face_z, - internal::constr_edge_zx, - internal::constr_type_x, - internal::constr_type_y}, - {internal::constr_face_x | internal::constr_face_z, - internal::constr_edge_zx, + {ConstraintTypes::face_x | ConstraintTypes::face_z, + ConstraintTypes::edge_zx, + ConstraintTypes::type_x | ConstraintTypes::type_z, + ConstraintTypes::type_y}, + {ConstraintTypes::face_x | ConstraintTypes::face_z, + ConstraintTypes::edge_zx, + ConstraintTypes::type_z, + ConstraintTypes::type_y}, + {ConstraintTypes::face_y | ConstraintTypes::face_z, + ConstraintTypes::edge_yz, + ConstraintTypes::type_y | ConstraintTypes::type_z, + ConstraintTypes::type_x}, + {ConstraintTypes::face_y | ConstraintTypes::face_z, + ConstraintTypes::edge_yz, + ConstraintTypes::type_z, + ConstraintTypes::type_x}, + {ConstraintTypes::face_x | ConstraintTypes::face_z, + ConstraintTypes::edge_zx, + ConstraintTypes::type_x, + ConstraintTypes::type_y}, + {ConstraintTypes::face_x | ConstraintTypes::face_z, + ConstraintTypes::edge_zx, 0, - internal::constr_type_y}, - {internal::constr_face_y | internal::constr_face_z, - internal::constr_edge_yz, - internal::constr_type_y, - internal::constr_type_x}, - {internal::constr_face_y | internal::constr_face_z, - internal::constr_edge_yz, + ConstraintTypes::type_y}, + {ConstraintTypes::face_y | ConstraintTypes::face_z, + ConstraintTypes::edge_yz, + ConstraintTypes::type_y, + ConstraintTypes::type_x}, + {ConstraintTypes::face_y | ConstraintTypes::face_z, + ConstraintTypes::edge_yz, 0, - internal::constr_type_x}, - {internal::constr_face_x | internal::constr_face_y, - internal::constr_edge_xy, - internal::constr_type_x | internal::constr_type_y, - internal::constr_type_z}, - {internal::constr_face_x | internal::constr_face_y, - internal::constr_edge_xy, - internal::constr_type_y, - internal::constr_type_z}, - {internal::constr_face_x | internal::constr_face_y, - internal::constr_edge_xy, - internal::constr_type_x, - internal::constr_type_z}, - {internal::constr_face_x | internal::constr_face_y, - internal::constr_edge_xy, + ConstraintTypes::type_x}, + {ConstraintTypes::face_x | ConstraintTypes::face_y, + ConstraintTypes::edge_xy, + ConstraintTypes::type_x | ConstraintTypes::type_y, + ConstraintTypes::type_z}, + {ConstraintTypes::face_x | ConstraintTypes::face_y, + ConstraintTypes::edge_xy, + ConstraintTypes::type_y, + ConstraintTypes::type_z}, + {ConstraintTypes::face_x | ConstraintTypes::face_y, + ConstraintTypes::edge_xy, + ConstraintTypes::type_x, + ConstraintTypes::type_z}, + {ConstraintTypes::face_x | ConstraintTypes::face_y, + ConstraintTypes::edge_xy, 0, - internal::constr_type_z}}; + ConstraintTypes::type_z}}; for (unsigned int local_line = 0; local_line < GeometryInfo::lines_per_cell; @@ -725,12 +742,19 @@ namespace CUDAWrappers else if (subface == 2) subface = 1; } + } // namespace MatrixFreeFunctions +} // namespace internal - //-------------------------------------------------------------------------// - // Functions for resolving the hanging node constraints // - //-------------------------------------------------------------------------// +#ifdef DEAL_II_COMPILER_CUDA_AWARE +namespace CUDAWrappers +{ + namespace internal + { + //--------------------------------------------------------------------------// + // Functions for resolving the hanging node constraints on the GPU // + //--------------------------------------------------------------------------// namespace internal { template @@ -763,7 +787,9 @@ namespace CUDAWrappers const unsigned int y_idx = threadIdx.y; const unsigned int this_type = - (direction == 0) ? internal::constr_type_x : internal::constr_type_y; + (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; @@ -772,19 +798,27 @@ namespace CUDAWrappers // given face. const bool constrained_face = (constraint_mask & - (((direction == 0) ? internal::constr_face_y : 0) | - ((direction == 1) ? internal::constr_face_x : 0))); + (((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 & internal::constr_type_y) ? - (y_idx == 0) : - (y_idx == fe_degree))) || - ((direction == 1) && ((constraint_mask & internal::constr_type_x) ? - (x_idx == 0) : - (x_idx == fe_degree))); + ((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) { @@ -800,11 +834,8 @@ namespace CUDAWrappers const Number w = transpose ? - internal::constraint_weights[i * (fe_degree + 1) + - interp_idx] : - internal::constraint_weights[interp_idx * - (fe_degree + 1) + - i]; + constraint_weights[i * (fe_degree + 1) + interp_idx] : + constraint_weights[interp_idx * (fe_degree + 1) + i]; t += w * values[real_idx]; } } @@ -818,12 +849,11 @@ namespace CUDAWrappers const Number w = transpose ? - internal::constraint_weights[(fe_degree - i) * - (fe_degree + 1) + - fe_degree - interp_idx] : - internal::constraint_weights[(fe_degree - interp_idx) * - (fe_degree + 1) + - fe_degree - i]; + 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]; } } @@ -853,32 +883,44 @@ namespace CUDAWrappers const unsigned int z_idx = threadIdx.z; const unsigned int this_type = - (direction == 0) ? internal::constr_type_x : - (direction == 1) ? internal::constr_type_y : - internal::constr_type_z; + (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) ? internal::constr_type_y : - (direction == 1) ? internal::constr_type_z : - internal::constr_type_x; + (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) ? internal::constr_type_z : - (direction == 1) ? internal::constr_type_x : - internal::constr_type_y; - - // If computing in x-direction, need to match against constr_face_y or - // constr_face_z - const unsigned int face1 = (direction == 0) ? internal::constr_face_y : - (direction == 1) ? - internal::constr_face_z : - internal::constr_face_x; - const unsigned int face2 = (direction == 0) ? internal::constr_face_z : - (direction == 1) ? - internal::constr_face_x : - internal::constr_face_y; - const unsigned int edge = (direction == 0) ? internal::constr_edge_yz : - (direction == 1) ? - internal::constr_edge_zx : - internal::constr_edge_xy; + (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); @@ -917,11 +959,8 @@ namespace CUDAWrappers const Number w = transpose ? - internal::constraint_weights[i * (fe_degree + 1) + - interp_idx] : - internal::constraint_weights[interp_idx * - (fe_degree + 1) + - i]; + constraint_weights[i * (fe_degree + 1) + interp_idx] : + constraint_weights[interp_idx * (fe_degree + 1) + i]; t += w * values[real_idx]; } } @@ -938,12 +977,11 @@ namespace CUDAWrappers const Number w = transpose ? - internal::constraint_weights[(fe_degree - i) * - (fe_degree + 1) + - fe_degree - interp_idx] : - internal::constraint_weights[(fe_degree - interp_idx) * - (fe_degree + 1) + - fe_degree - i]; + 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]; } } @@ -997,9 +1035,8 @@ namespace CUDAWrappers } } // namespace internal } // namespace CUDAWrappers +#endif DEAL_II_NAMESPACE_CLOSE #endif - -#endif -- 2.39.5