From f5890276c3012ecd55f4e148517561c5504029d3 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 23 Sep 2019 14:34:41 +0200 Subject: [PATCH] Eliminate q_point from functor and other GPU functions --- doc/news/changes/minor/20190924PeterMunch | 5 + .../deal.II/matrix_free/cuda_fe_evaluation.h | 244 ++++++++++++++++-- tests/cuda/matrix_free_no_index_initialize.cu | 13 +- tests/cuda/matrix_vector_mf.h | 14 +- 4 files changed, 235 insertions(+), 41 deletions(-) create mode 100644 doc/news/changes/minor/20190924PeterMunch diff --git a/doc/news/changes/minor/20190924PeterMunch b/doc/news/changes/minor/20190924PeterMunch new file mode 100644 index 0000000000..9232cc6c4d --- /dev/null +++ b/doc/news/changes/minor/20190924PeterMunch @@ -0,0 +1,5 @@ +New: Provide new methods in CUDAWrappers::FEEvaluation that do neither take the +local dof index nor the quadrature point but recompute them based on the +thread id. +
+(Peter Munch, 2019/09/24) diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index b2bb935ccf..4c37b94981 100644 --- a/include/deal.II/matrix_free/cuda_fe_evaluation.h +++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h @@ -38,6 +38,25 @@ DEAL_II_NAMESPACE_OPEN */ namespace CUDAWrappers { + namespace internal + { + /** + * Compute the dof/quad index for a given thread id, dimension, and + * number of points in each space dimensions. + */ + template + __device__ inline unsigned int + compute_index() + { + return (dim == 1 ? + threadIdx.x % n_points_1d : + dim == 2 ? + threadIdx.x % n_points_1d + n_points_1d * threadIdx.y : + threadIdx.x % n_points_1d + + n_points_1d * (threadIdx.y + n_points_1d * threadIdx.z)); + } + } // namespace internal + /** * This class provides all the functions necessary to evaluate functions at * quadrature points and cell integrations. In functionality, this class is @@ -134,48 +153,103 @@ namespace CUDAWrappers /** * Return the value of a finite element function at quadrature point * number @p q_point after a call to @p evaluate(true,...). + * + * @deprecated Use the version without parameters instead. + */ + DEAL_II_DEPRECATED __device__ value_type + get_value(const unsigned int q_point) const; + + /** + * Same as above, except that the quadrature point is computed from thread + * id. */ __device__ value_type - get_value(const unsigned int q_point) const; + get_value() const; /** * Return the value of a finite element function at degree of freedom * @p dof after a call to integrate() or before a call to evaluate(). + * + * @deprecated Use the version without parameters instead. + */ + DEAL_II_DEPRECATED __device__ value_type + get_dof_value(const unsigned int dof) const; + + /** + * Same as above, except that the local dof index is computed from the + * thread id. */ __device__ value_type - get_dof_value(const unsigned int dof) const; + get_dof_value() const; /** * Write a value to the field containing the values on quadrature points * with component @p q_point. Access to the same fields as through @p * get_value(). This specifies the value which is tested by all basis * function on the current cell and integrated over. + * + * @deprecated Use the version without parameters instead. + */ + DEAL_II_DEPRECATED __device__ void + submit_value(const value_type &val_in, const unsigned int q_point); + + /** + * Same as above, except that the quadrature point is computed from the + * thread id. */ __device__ void - submit_value(const value_type &val_in, const unsigned int q_point); + submit_value(const value_type &val_in); /** * Write a value to the field containing the values for the degree of * freedom with index @p dof after a call to integrate() or before * calling evaluate(). Access through the same fields as through * get_dof_value(). + * + * @deprecated Use the version without parameters instead. + */ + DEAL_II_DEPRECATED __device__ void + submit_dof_value(const value_type &val_in, const unsigned int dof); + + /** + * Same as above, except that the local dof index is computed from the + * thread id. */ __device__ void - submit_dof_value(const value_type &val_in, const unsigned int dof); + submit_dof_value(const value_type &val_in); /** * Return the gradient of a finite element function at quadrature point * number @p q_point after a call to @p evaluate(...,true). + * + * @deprecated Use the version without parameters instead. + */ + DEAL_II_DEPRECATED __device__ gradient_type + get_gradient(const unsigned int q_point) const; + + /** + * Same as above, except that the quadrature point is computed from the + * thread id. */ __device__ gradient_type - get_gradient(const unsigned int q_point) const; + get_gradient() const; /** * Write a contribution that is tested by the gradient to the field - * containing the values on quadrature points with component @p q_point + * containing the values on quadrature points with component @p q_point. + * + * @deprecated Use the version without parameters instead. + */ + DEAL_II_DEPRECATED __device__ void + submit_gradient(const gradient_type &grad_in, const unsigned int q_point); + + + /** + * Same as above, except that the quadrature point is computed from the + * thread id. */ __device__ void - submit_gradient(const gradient_type &grad_in, const unsigned int q_point); + submit_gradient(const gradient_type &grad_in); // clang-format off /** @@ -187,11 +261,29 @@ namespace CUDAWrappers * CUDAWrappers::FEEvaluation *fe_eval, * const unsigned int q_point) const; * \endcode + * + * @deprecated Use apply_for_each_quad_point() instead. + */ + // clang-format on + template + DEAL_II_DEPRECATED __device__ void + apply_quad_point_operations(const Functor &func); + + // clang-format off + /** + * Same as above, except that the functor @func only takes a single input + * argument (fe_eval) and computes the quadrature point from the thread id. + * + * @p func needs to define + * \code + * __device__ void operator()( + * CUDAWrappers::FEEvaluation *fe_eval) const; + * \endcode */ // clang-format on template __device__ void - apply_quad_point_operations(const Functor &func); + apply_for_each_quad_point(const Functor &func); private: types::global_dof_index *local_to_global; @@ -249,10 +341,7 @@ namespace CUDAWrappers { static_assert(n_components_ == 1, "This function only supports FE with one \ components"); - const unsigned int idx = - (threadIdx.x % n_q_points_1d) + - (dim > 1 ? threadIdx.y : 0) * n_q_points_1d + - (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d; + const unsigned int idx = internal::compute_index(); const types::global_dof_index src_idx = local_to_global[idx]; // Use the read-only data cache. @@ -280,10 +369,8 @@ namespace CUDAWrappers internal::resolve_hanging_nodes(constraint_mask, values); - const unsigned int idx = - (threadIdx.x % n_q_points_1d) + - (dim > 1 ? threadIdx.y : 0) * n_q_points_1d + - (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d; + const unsigned int idx = internal::compute_index(); + const types::global_dof_index destination_idx = local_to_global[idx]; if (use_coloring) @@ -384,6 +471,24 @@ namespace CUDAWrappers + template + __device__ typename FEEvaluation::value_type + FEEvaluation:: + get_value() const + { + return get_value(internal::compute_index()); + } + + + template + __device__ typename FEEvaluation::value_type + FEEvaluation:: + get_dof_value() const + { + return get_dof_value(internal::compute_index()); + } + + + template + __device__ void + FEEvaluation:: + submit_value(const value_type &val_in) + { + submit_value(val_in, internal::compute_index()); + } + + + template + __device__ void + FEEvaluation:: + submit_dof_value(const value_type &val_in) + { + submit_dof_value(val_in, internal::compute_index()); + } + + + template + __device__ typename FEEvaluation::gradient_type + FEEvaluation:: + get_gradient() const + { + static_assert(n_components_ == 1, "This function only supports FE with one \ + components"); + + return get_gradient(internal::compute_index()); + } + + + template + __device__ void + FEEvaluation:: + submit_gradient(const gradient_type &grad_in) + { + submit_gradient(grad_in, internal::compute_index()); + } + + + template :: apply_quad_point_operations(const Functor &func) { - const unsigned int q_point = - (dim == 1 ? - threadIdx.x % n_q_points_1d : - dim == 2 ? - threadIdx.x % n_q_points_1d + n_q_points_1d * threadIdx.y : - threadIdx.x % n_q_points_1d + - n_q_points_1d * (threadIdx.y + n_q_points_1d * threadIdx.z)); - func(this, q_point); + func(this, internal::compute_index()); + + __syncthreads(); + } + + + + template + template + __device__ void + FEEvaluation:: + apply_for_each_quad_point(const Functor &func) + { + func(this); __syncthreads(); } diff --git a/tests/cuda/matrix_free_no_index_initialize.cu b/tests/cuda/matrix_free_no_index_initialize.cu index 85b2a780ac..a66a64fbd9 100644 --- a/tests/cuda/matrix_free_no_index_initialize.cu +++ b/tests/cuda/matrix_free_no_index_initialize.cu @@ -49,25 +49,20 @@ public: CUDAWrappers::FEEvaluation fe_eval(cell, gpu_data, shared_data); - const unsigned int tid = - (threadIdx.x % n_q_points_1d) + - (dim > 1 ? threadIdx.y : 0) * n_q_points_1d + - (dim > 2 ? threadIdx.z : 0) * n_q_points_1d * n_q_points_1d; - // set to unit vector - fe_eval.submit_dof_value(1., tid); + fe_eval.submit_dof_value(1.); __syncthreads(); fe_eval.evaluate(/*evaluate_values =*/true, /*evaluate_gradients=*/true); #ifndef __APPLE__ // values should evaluate to one, derivatives to zero - assert(fe_eval.get_value(tid) == 1.); + assert(fe_eval.get_value() == 1.); for (unsigned int e = 0; e < dim; ++e) - assert(fe_eval.get_gradient(tid)[e] == 0.); + assert(fe_eval.get_gradient()[e] == 0.); fe_eval.integrate(/*integrate_values = */ true, /*integrate_gradients=*/true); - assert(fe_eval.get_dof_value(tid) == 1.); + assert(fe_eval.get_dof_value() == 1.); #endif } diff --git a/tests/cuda/matrix_vector_mf.h b/tests/cuda/matrix_vector_mf.h index ef4a0c3d4f..e27fdc1d0e 100644 --- a/tests/cuda/matrix_vector_mf.h +++ b/tests/cuda/matrix_vector_mf.h @@ -37,8 +37,7 @@ public: __device__ void operator()( CUDAWrappers::FEEvaluation - * fe_eval, - const unsigned int q_point) const; + *fe_eval) const; private: Number coef; @@ -48,12 +47,11 @@ private: template __device__ void HelmholtzOperatorQuad:: - operator()( - CUDAWrappers::FEEvaluation *fe_eval, - const unsigned int q) const + operator()(CUDAWrappers::FEEvaluation + *fe_eval) const { - fe_eval->submit_value(coef * fe_eval->get_value(q), q); - fe_eval->submit_gradient(fe_eval->get_gradient(q), q); + fe_eval->submit_value(coef * fe_eval->get_value()); + fe_eval->submit_gradient(fe_eval->get_gradient()); } @@ -101,7 +99,7 @@ operator()(const unsigned int cell, cell, gpu_data, shared_data); fe_eval.read_dof_values(src); fe_eval.evaluate(true, true); - fe_eval.apply_quad_point_operations( + fe_eval.apply_for_each_quad_point( HelmholtzOperatorQuad(coef[pos])); fe_eval.integrate(true, true); fe_eval.distribute_local_to_global(dst); -- 2.39.5