From b600783b94ad2a53d85de15b0cf52391a6eca3f9 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Tue, 9 Jul 2019 18:32:52 -0400 Subject: [PATCH] Implement CUDAWrappers::FEEvaluation::get_dof_values/submit_dof_values --- doc/news/changes/minor/20190709DanielArndt | 5 + .../deal.II/matrix_free/cuda_fe_evaluation.h | 52 +++++- tests/cuda/matrix_free_no_index_initialize.cu | 158 ++++++++++++++++++ .../matrix_free_no_index_initialize.output | 3 + 4 files changed, 216 insertions(+), 2 deletions(-) create mode 100644 doc/news/changes/minor/20190709DanielArndt create mode 100644 tests/cuda/matrix_free_no_index_initialize.cu create mode 100644 tests/cuda/matrix_free_no_index_initialize.output diff --git a/doc/news/changes/minor/20190709DanielArndt b/doc/news/changes/minor/20190709DanielArndt new file mode 100644 index 0000000000..3b0a6588ae --- /dev/null +++ b/doc/news/changes/minor/20190709DanielArndt @@ -0,0 +1,5 @@ +New: CUDAWrappers::FEEvaluation::get_dof_values() and +CUDAWrappers::FEEvaluation::submit_dof_values() provide access to the values +stored for the degrees of freedom. +
+(Daniel Arndt, 2019/07/09) diff --git a/include/deal.II/matrix_free/cuda_fe_evaluation.h b/include/deal.II/matrix_free/cuda_fe_evaluation.h index 9635ba219f..755b6fdcda 100644 --- a/include/deal.II/matrix_free/cuda_fe_evaluation.h +++ b/include/deal.II/matrix_free/cuda_fe_evaluation.h @@ -132,20 +132,36 @@ namespace CUDAWrappers /** * Return the value of a finite element function at quadrature point - * number @p q_point after a call to @p evalue(true,...). + * number @p q_point after a call to @p evaluate(true,...). */ __device__ value_type get_value(const unsigned int q_point) 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(). + */ + __device__ value_type + get_dof_value(const unsigned int dof) 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 + * get_value(). This specifies the value which is tested by all basis * function on the current cell and integrated over. */ __device__ void submit_value(const value_type &val_in, const unsigned int q_point); + /** + * 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(). + */ + __device__ void + submit_dof_value(const value_type &val_in, const unsigned int dof); + /** * Return the gradient of a finite element function at quadrature point * number @p q_point after a call to @p evaluate(...,true). @@ -360,6 +376,24 @@ namespace CUDAWrappers + template + __device__ typename FEEvaluation::value_type + FEEvaluation:: + get_dof_value(const unsigned int dof) const + { + return values[dof]; + } + + + template + __device__ void + FEEvaluation:: + submit_dof_value(const value_type &val_in, const unsigned int dof) + { + values[dof] = val_in; + } + + + template +#include + +#include + +#include "../tests.h" + +// forward declare this function +template +void +test(); + +template +class MatrixFreeTest +{ +public: + static const unsigned int n_dofs_1d = fe_degree + 1; + static const unsigned int n_local_dofs = Utilities::pow(n_dofs_1d, dim); + static const unsigned int n_q_points = Utilities::pow(n_q_points_1d, dim); + + MatrixFreeTest(const CUDAWrappers::MatrixFree &data_in) + : data(data_in){}; + + __device__ void + operator()( + const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data, + CUDAWrappers::SharedData * shared_data, + const Number * src, + Number * dst) const + { + 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; + + bool const evaluate_values = true; + bool const evaluate_gradients = true; + + // set to unit vector + fe_eval.submit_dof_value(1., tid); + __syncthreads(); + fe_eval.evaluate(evaluate_values, evaluate_gradients); + + // values should evaluate to one, derivatives to zero + assert(fe_eval.get_dof_value(tid) == 1.); + assert(fe_eval.get_value(tid) == 1.); + for (unsigned int e = 0; e < dim; ++e) + assert(fe_eval.get_gradient(tid)[e] == 0.); + } + + + + void + test() const + { + LinearAlgebra::distributed::Vector dst_dummy; + LinearAlgebra::distributed::Vector src_dummy; + + data.cell_loop(*this, src_dummy, dst_dummy); + + // Check that the kernel was launched correctly + AssertCuda(cudaGetLastError()); + // Check that there was no problem during the execution of the kernel + AssertCuda(cudaDeviceSynchronize()); + deallog << "OK" << std::endl; + }; + +protected: + const CUDAWrappers::MatrixFree &data; +}; + + + +template +void +do_test(const DoFHandler & dof, + const AffineConstraints &constraints) +{ + CUDAWrappers::MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename CUDAWrappers::MatrixFree::AdditionalData data; + data.mapping_update_flags = update_values | update_gradients | + update_JxW_values | update_quadrature_points; + mf_data.reinit(dof, constraints, quad, data); + } + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + MatrixFreeTest mf(mf_data); + mf.test(); +} + + +int +main() +{ + initlog(); + deal_II_exceptions::disable_abort_on_exception(); + + test<2, 1>(); +} + + +template +void +test() +{ + const SphericalManifold manifold; + Triangulation tria; + GridGenerator::hyper_ball(tria); + for (const auto &cell : tria.active_cell_iterators()) + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + if (cell->at_boundary(f)) + cell->face(f)->set_all_manifold_ids(0); + tria.set_manifold(0, manifold); + + // refine first and last cell + tria.begin(tria.n_levels() - 1)->set_refine_flag(); + tria.last()->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + tria.refine_global(4 - dim); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + constraints.close(); + + do_test(dof, constraints); +} diff --git a/tests/cuda/matrix_free_no_index_initialize.output b/tests/cuda/matrix_free_no_index_initialize.output new file mode 100644 index 0000000000..09b5b737bc --- /dev/null +++ b/tests/cuda/matrix_free_no_index_initialize.output @@ -0,0 +1,3 @@ + +DEAL::Testing FE_Q<2>(1) +DEAL::OK -- 2.39.5