From b1b35e9d1e9e7f329356ad61be39d13c3609e41a Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Tue, 2 Jul 2019 16:26:06 +0000 Subject: [PATCH] Fix a bug in CUDA matrix-free with adapted mesh when using Volta On Volta the code would deadlock. Refactor the code so that __syncthreads is not called inside a code that be executed differently for different threads. --- .../deal.II/matrix_free/cuda_fe_evaluation.h | 10 +- .../matrix_free/cuda_hanging_nodes_internal.h | 277 +++++++++--------- 2 files changed, 143 insertions(+), 144 deletions(-) diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index 8bb8b6c4ef..9635ba219f 100644 --- a/include/deal.II/matrix_free/cuda_fe_evaluation.h +++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h @@ -238,9 +238,8 @@ namespace CUDAWrappers // Use the read-only data cache. values[idx] = __ldg(&src[src_idx]); - if (constraint_mask) - internal::resolve_hanging_nodes(constraint_mask, - values); + internal::resolve_hanging_nodes(constraint_mask, + values); __syncthreads(); } @@ -258,9 +257,8 @@ namespace CUDAWrappers { static_assert(n_components_ == 1, "This function only supports FE with one \ components"); - if (constraint_mask) - internal::resolve_hanging_nodes(constraint_mask, - values); + internal::resolve_hanging_nodes(constraint_mask, + values); const unsigned int idx = (threadIdx.x % n_q_points_1d) + diff --git a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h index a07e1c2db1..9a9b828749 100644 --- a/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h +++ b/include/deal.II/matrix_free/cuda_hanging_nodes_internal.h @@ -758,7 +758,8 @@ namespace CUDAWrappers bool transpose, typename Number> __device__ inline void - interpolate_boundary_2d(const unsigned int constr, Number *values) + 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; @@ -766,72 +767,73 @@ namespace CUDAWrappers const unsigned int this_type = (direction == 0) ? internal::constr_type_x : internal::constr_type_y; - if (constr & (((direction == 0) ? internal::constr_face_y : 0) | - ((direction == 1) ? internal::constr_face_x : 0))) + 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) ? internal::constr_face_y : 0) | + ((direction == 1) ? internal::constr_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))); + if (constrained_face && constrained_dof) { - const unsigned int interp_idx = (direction == 0) ? x_idx : y_idx; + const bool type = constraint_mask & this_type; - __syncthreads(); - - Number t = 0; - const bool flag = - ((direction == 0) && - ((constr & internal::constr_type_y) ? (y_idx == 0) : - (y_idx == fe_degree))) || - ((direction == 1) && - ((constr & internal::constr_type_x) ? (x_idx == 0) : - (x_idx == fe_degree))); - - if (flag) + if (type) { - const bool type = constr & this_type; - - if (type) + for (unsigned int i = 0; i <= fe_degree; ++i) { - 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 ? - internal::constraint_weights[i * (fe_degree + 1) + - interp_idx] : - internal::constraint_weights[interp_idx * - (fe_degree + 1) + - i]; - t += w * values[real_idx]; - } + const unsigned int real_idx = + (direction == 0) ? index2(i, y_idx) : + index2(x_idx, i); + + const Number w = + transpose ? + internal::constraint_weights[i * (fe_degree + 1) + + interp_idx] : + internal::constraint_weights[interp_idx * + (fe_degree + 1) + + i]; + t += w * values[real_idx]; } - else + } + else + { + for (unsigned int i = 0; i <= fe_degree; ++i) { - 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 ? - 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]; - t += w * values[real_idx]; - } + const unsigned int real_idx = + (direction == 0) ? index2(i, y_idx) : + index2(x_idx, i); + + 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]; + t += w * values[real_idx]; } } + } - __syncthreads(); + __syncthreads(); - if (flag) - values[index2(x_idx, y_idx)] = t; - } + if (constrained_face && constrained_dof) + values[index2(x_idx, y_idx)] = t; } @@ -841,7 +843,8 @@ namespace CUDAWrappers bool transpose, typename Number> __device__ inline void - interpolate_boundary_3d(const unsigned int constr, Number *values) + 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; @@ -874,84 +877,81 @@ namespace CUDAWrappers (direction == 1) ? internal::constr_edge_zx : internal::constr_edge_xy; - - if (constr & (face1 | face2 | edge)) + 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)); + (direction == 0) ? z_idx : (direction == 1) ? x_idx : y_idx; + + if (constrained_face && constrained_dof) { - 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; - - __syncthreads(); - - Number t = 0; - const bool on_face1 = (constr & face1_type) ? - (face1_idx == 0) : - (face1_idx == fe_degree); - const bool on_face2 = (constr & face2_type) ? - (face2_idx == 0) : - (face2_idx == fe_degree); - const bool flag = (((constr & face1) && on_face1) || - ((constr & face2) && on_face2) || - ((constr & edge) && on_face1 && on_face2)); - - if (flag) + const bool type = constraint_mask & this_type; + if (type) { - const bool type = constr & this_type; - if (type) + for (unsigned int i = 0; i <= fe_degree; ++i) { - 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 ? - internal::constraint_weights[i * (fe_degree + 1) + - interp_idx] : - internal::constraint_weights[interp_idx * - (fe_degree + 1) + - i]; - t += w * values[real_idx]; - } + 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 ? + internal::constraint_weights[i * (fe_degree + 1) + + interp_idx] : + internal::constraint_weights[interp_idx * + (fe_degree + 1) + + i]; + t += w * values[real_idx]; } - else + } + else + { + for (unsigned int i = 0; i <= fe_degree; ++i) { - 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 ? - 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]; - t += w * values[real_idx]; - } + 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 ? + 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]; + t += w * values[real_idx]; } } + } - __syncthreads(); + __syncthreads(); - if (flag) - values[index3(x_idx, y_idx, z_idx)] = t; - } + if (constrained_face && constrained_dof) + values[index3(x_idx, y_idx, z_idx)] = t; } } // namespace internal @@ -967,26 +967,27 @@ namespace CUDAWrappers */ template __device__ void - resolve_hanging_nodes(const unsigned int constr, Number *values) + resolve_hanging_nodes(const unsigned int constraint_mask, Number *values) { if (dim == 2) { - internal::interpolate_boundary_2d(constr, - values); - internal::interpolate_boundary_2d(constr, - values); + internal::interpolate_boundary_2d( + constraint_mask, values); + + internal::interpolate_boundary_2d( + constraint_mask, values); } else if (dim == 3) { // Interpolate y and z faces (x-direction) - internal::interpolate_boundary_3d(constr, - values); + internal::interpolate_boundary_3d( + constraint_mask, values); // Interpolate x and z faces (y-direction) - internal::interpolate_boundary_3d(constr, - values); + internal::interpolate_boundary_3d( + constraint_mask, values); // Interpolate x and y faces (z-direction) - internal::interpolate_boundary_3d(constr, - values); + internal::interpolate_boundary_3d( + constraint_mask, values); } } } // namespace internal -- 2.39.5