From: Bruno Turcksin Date: Thu, 16 Mar 2023 18:12:15 +0000 (-0400) Subject: Move JxW to Kokkos::View X-Git-Tag: v9.5.0-rc1~349^2~12 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=77d7cb35cc325834dc47a3c2698c2248f2f9c0b0;p=dealii.git Move JxW to Kokkos::View --- diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index fd603d9ad9..541c0ab852 100644 --- a/include/deal.II/matrix_free/cuda_fe_evaluation.h +++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h @@ -18,6 +18,8 @@ #include +#include "deal.II/base/memory_space.h" + #ifdef DEAL_II_WITH_CUDA # include @@ -31,6 +33,8 @@ # include # include +# include + DEAL_II_NAMESPACE_OPEN /** @@ -247,7 +251,15 @@ namespace CUDAWrappers const bool use_coloring; Number *inv_jac; - Number *JxW; + // We would like to use + // Kokkos::Subview, int, decltype(Kokkos::ALL)> + // but we get error: incomplete type is not allowed. I cannot reproduce + // outside of deal.II. Need to investigate more. + Kokkos::Subview, + int, + Kokkos::pair> + JxW; // Internal buffer Number *values; @@ -271,11 +283,14 @@ namespace CUDAWrappers , mf_object_id(data->id) , constraint_mask(data->constraint_mask[cell_id]) , use_coloring(data->use_coloring) + , JxW(Kokkos::subview( + data->JxW, + cell_id, + Kokkos::pair(0, Utilities::pow(n_q_points_1d, dim)))) , values(shdata->values) { local_to_global = data->local_to_global + padding_length * cell_id; inv_jac = data->inv_jacobian + padding_length * cell_id; - JxW = data->JxW + padding_length * cell_id; for (unsigned int i = 0; i < dim; ++i) gradients[i] = shdata->gradients[i]; diff --git a/include/deal.II/matrix_free/cuda_matrix_free.h b/include/deal.II/matrix_free/cuda_matrix_free.h index 078798384c..74d6955cd1 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.h @@ -19,9 +19,12 @@ #include +#include "deal.II/base/memory_space.h" + #ifdef DEAL_II_WITH_CUDA # include +# include # include # include # include @@ -118,15 +121,12 @@ namespace CUDAWrappers /** * Constructor. */ - AdditionalData( - const ParallelizationScheme parallelization_scheme = parallel_in_elem, - const UpdateFlags mapping_update_flags = update_gradients | - update_JxW_values | - update_quadrature_points, - const bool use_coloring = false, - const bool overlap_communication_computation = false) - : parallelization_scheme(parallelization_scheme) - , mapping_update_flags(mapping_update_flags) + AdditionalData(const UpdateFlags mapping_update_flags = + update_gradients | update_JxW_values | + update_quadrature_points, + const bool use_coloring = false, + const bool overlap_communication_computation = false) + : mapping_update_flags(mapping_update_flags) , use_coloring(use_coloring) , overlap_communication_computation(overlap_communication_computation) { @@ -142,11 +142,6 @@ namespace CUDAWrappers ExcMessage( "Overlapping communication and coloring are incompatible options. Only one of them can be enabled.")); } - /** - * Parallelization scheme used, parallelization over degrees of freedom or - * over cells. - */ - ParallelizationScheme parallelization_scheme; /** * This flag is used to determine which quantities should be cached. This * class can cache data needed for gradient computations (inverse @@ -195,9 +190,9 @@ namespace CUDAWrappers Number *inv_jacobian; /** - * Pointer to the Jacobian times the weights. + * Kokkos::View of the Jacobian times the weights. */ - Number *JxW; + Kokkos::View JxW; /** * ID of the associated MatrixFree object. @@ -458,12 +453,6 @@ namespace CUDAWrappers */ int my_id; - /** - * Parallelization scheme used, parallelization over degrees of freedom or - * over cells. - */ - ParallelizationScheme parallelization_scheme; - /** * If true, use graph coloring. Otherwise, use atomic operations. Graph * coloring ensures bitwise reproducibility but is slower on Pascal and @@ -531,10 +520,11 @@ namespace CUDAWrappers std::vector inv_jacobian; /** - * Vector of pointer to the Jacobian times the weights associated to the - * cells of each color. + * Vector of Kokkos::View to the Jacobian times the weights associated to + * the cells of each color. */ - std::vector JxW; + std::vector> + JxW; /** * Pointer to the constrained degrees of freedom. @@ -830,7 +820,7 @@ namespace CUDAWrappers if (update_flags & update_JxW_values) { data_host.JxW.resize(n_elements); - Utilities::CUDA::copy_to_host(data.JxW, data_host.JxW); + Utilities::CUDA::copy_to_host(data.JxW.data(), data_host.JxW); } data_host.constraint_mask.resize(data_host.n_cells); 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 9df3cfee63..69c9f885b6 100644 --- a/include/deal.II/matrix_free/cuda_matrix_free.templates.h +++ b/include/deal.II/matrix_free/cuda_matrix_free.templates.h @@ -19,8 +19,12 @@ #include +#include "deal.II/base/memory_space.h" + #include +#include + #ifdef DEAL_II_WITH_CUDA # include @@ -231,8 +235,8 @@ namespace CUDAWrappers // Host data std::vector local_to_global_host; std::vector> q_points_host; - std::vector JxW_host; - std::vector inv_jacobian_host; + Kokkos::View JxW; + std::vector inv_jacobian_host; std::vector constraint_mask_host; // Local buffer @@ -334,20 +338,13 @@ namespace CUDAWrappers // TODO this should be a templated parameter. const unsigned int n_dofs_1d = fe_degree + 1; - if (data->parallelization_scheme == - MatrixFree::parallel_in_elem) - { - if (dim == 1) - data->block_dim[color] = dim3(n_dofs_1d * cells_per_block); - else if (dim == 2) - data->block_dim[color] = - dim3(n_dofs_1d * cells_per_block, n_dofs_1d); - else - data->block_dim[color] = - dim3(n_dofs_1d * cells_per_block, n_dofs_1d, n_dofs_1d); - } + if (dim == 1) + data->block_dim[color] = dim3(n_dofs_1d * cells_per_block); + else if (dim == 2) + data->block_dim[color] = dim3(n_dofs_1d * cells_per_block, n_dofs_1d); else - data->block_dim[color] = dim3(cells_per_block); + data->block_dim[color] = + dim3(n_dofs_1d * cells_per_block, n_dofs_1d, n_dofs_1d); local_to_global_host.resize(n_cells * padding_length); @@ -355,7 +352,11 @@ namespace CUDAWrappers q_points_host.resize(n_cells * padding_length); if (update_flags & update_JxW_values) - JxW_host.resize(n_cells * padding_length); + JxW = Kokkos::View( + Kokkos::view_alloc("JxW_" + std::to_string(color), + Kokkos::WithoutInitializing), + n_cells, + dofs_per_cell); if (update_flags & update_gradients) inv_jacobian_host.resize(n_cells * padding_length * dim * dim); @@ -411,10 +412,13 @@ namespace CUDAWrappers if (update_flags & update_JxW_values) { + // FIXME too many deep_copy + auto JxW_host = Kokkos::create_mirror_view_and_copy( + MemorySpace::Host::kokkos_space{}, JxW); std::vector JxW_values_double = fe_values.get_JxW_values(); - const unsigned int offset = cell_id * padding_length; for (unsigned int i = 0; i < q_points_per_cell; ++i) - JxW_host[i + offset] = static_cast(JxW_values_double[i]); + JxW_host(cell_id, i) = static_cast(JxW_values_double[i]); + Kokkos::deep_copy(JxW, JxW_host); } if (update_flags & update_gradients) @@ -438,10 +442,6 @@ namespace CUDAWrappers const unsigned int n_cells = data->n_cells[color]; // Local-to-global mapping - if (data->parallelization_scheme == - MatrixFree::parallel_over_elem) - transpose_in_place(local_to_global_host, n_cells, padding_length); - alloc_and_copy( &data->local_to_global[color], ArrayView(local_to_global_host.data(), @@ -451,10 +451,6 @@ namespace CUDAWrappers // Quadrature points if (update_flags & update_quadrature_points) { - if (data->parallelization_scheme == - MatrixFree::parallel_over_elem) - transpose_in_place(q_points_host, n_cells, padding_length); - alloc_and_copy(&data->q_points[color], ArrayView>( q_points_host.data(), q_points_host.size()), @@ -464,14 +460,7 @@ namespace CUDAWrappers // Jacobian determinants/quadrature weights if (update_flags & update_JxW_values) { - if (data->parallelization_scheme == - MatrixFree::parallel_over_elem) - transpose_in_place(JxW_host, n_cells, padding_length); - - alloc_and_copy(&data->JxW[color], - ArrayView(JxW_host.data(), - JxW_host.size()), - n_cells * padding_length); + data->JxW[color] = JxW; } // Inverse jacobians @@ -485,15 +474,6 @@ namespace CUDAWrappers padding_length * n_cells, dim * dim); - // Transpose second time means we get the following index order: - // q*n_cells*dim*dim + i*n_cells + cell_id which is good for an - // element-level parallelization - if (data->parallelization_scheme == - MatrixFree::parallel_over_elem) - transpose_in_place(inv_jacobian_host, - n_cells * dim * dim, - padding_length); - alloc_and_copy(&data->inv_jacobian[color], ArrayView(inv_jacobian_host.data(), inv_jacobian_host.size()), @@ -740,10 +720,6 @@ namespace CUDAWrappers Utilities::CUDA::free(inv_jacobian_color_ptr); inv_jacobian.clear(); - for (auto &JxW_color_ptr : JxW) - Utilities::CUDA::free(JxW_color_ptr); - JxW.clear(); - for (auto &constraint_mask_color_ptr : constraint_mask) Utilities::CUDA::free(constraint_mask_color_ptr); constraint_mask.clear(); @@ -892,6 +868,7 @@ namespace CUDAWrappers n_constrained_dofs * sizeof(unsigned int); // For each color, add local_to_global, inv_jacobian, JxW, and q_points. + // FIXME for (unsigned int i = 0; i < n_colors; ++i) { bytes += n_cells[i] * padding_length * sizeof(unsigned int) + @@ -927,12 +904,7 @@ namespace CUDAWrappers if (update_flags & update_gradients) update_flags |= update_JxW_values; - if (additional_data.parallelization_scheme != parallel_over_elem && - additional_data.parallelization_scheme != parallel_in_elem) - AssertThrow(false, ExcMessage("Invalid parallelization scheme.")); - - this->parallelization_scheme = additional_data.parallelization_scheme; - this->use_coloring = additional_data.use_coloring; + this->use_coloring = additional_data.use_coloring; this->overlap_communication_computation = additional_data.overlap_communication_computation;