From 5a0a7adf67e6505f841bda374a0a438f42a671e2 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Mon, 10 Dec 2018 20:29:47 +0000 Subject: [PATCH] Adapt current tests to be able to switch between constant and varying coefficients --- tests/cuda/matrix_free_matrix_vector_01.cu | 2 +- tests/cuda/matrix_free_matrix_vector_02.cu | 2 +- tests/cuda/matrix_free_matrix_vector_03.cu | 2 +- tests/cuda/matrix_free_matrix_vector_04.cu | 2 +- tests/cuda/matrix_free_matrix_vector_05.cu | 2 +- tests/cuda/matrix_free_matrix_vector_06.cu | 2 +- tests/cuda/matrix_free_matrix_vector_06a.cu | 2 +- tests/cuda/matrix_free_matrix_vector_09.cu | 2 +- tests/cuda/matrix_free_matrix_vector_10.cu | 4 +- tests/cuda/matrix_free_matrix_vector_19.cu | 4 +- tests/cuda/matrix_vector_common.h | 53 ++++++--- tests/cuda/matrix_vector_mf.h | 116 +++++++++++++++++++- 12 files changed, 165 insertions(+), 28 deletions(-) diff --git a/tests/cuda/matrix_free_matrix_vector_01.cu b/tests/cuda/matrix_free_matrix_vector_01.cu index efd928d119..ce08bb12f1 100644 --- a/tests/cuda/matrix_free_matrix_vector_01.cu +++ b/tests/cuda/matrix_free_matrix_vector_01.cu @@ -45,5 +45,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_02.cu b/tests/cuda/matrix_free_matrix_vector_02.cu index c6388d8024..166b73cf87 100644 --- a/tests/cuda/matrix_free_matrix_vector_02.cu +++ b/tests/cuda/matrix_free_matrix_vector_02.cu @@ -49,5 +49,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_03.cu b/tests/cuda/matrix_free_matrix_vector_03.cu index cd556e16ee..147febd3f5 100644 --- a/tests/cuda/matrix_free_matrix_vector_03.cu +++ b/tests/cuda/matrix_free_matrix_vector_03.cu @@ -72,5 +72,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_04.cu b/tests/cuda/matrix_free_matrix_vector_04.cu index 61dd8b491d..c4c1085670 100644 --- a/tests/cuda/matrix_free_matrix_vector_04.cu +++ b/tests/cuda/matrix_free_matrix_vector_04.cu @@ -48,5 +48,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_05.cu b/tests/cuda/matrix_free_matrix_vector_05.cu index 54ce17e032..48efc829d1 100644 --- a/tests/cuda/matrix_free_matrix_vector_05.cu +++ b/tests/cuda/matrix_free_matrix_vector_05.cu @@ -82,5 +82,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_06.cu b/tests/cuda/matrix_free_matrix_vector_06.cu index 55a4540e3d..704255a86e 100644 --- a/tests/cuda/matrix_free_matrix_vector_06.cu +++ b/tests/cuda/matrix_free_matrix_vector_06.cu @@ -78,5 +78,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_06a.cu b/tests/cuda/matrix_free_matrix_vector_06a.cu index c9c5d1cb42..7aa0b084fd 100644 --- a/tests/cuda/matrix_free_matrix_vector_06a.cu +++ b/tests/cuda/matrix_free_matrix_vector_06a.cu @@ -74,5 +74,5 @@ test() fe_degree, double, LinearAlgebra::distributed::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_09.cu b/tests/cuda/matrix_free_matrix_vector_09.cu index 1224ed4b68..24f36ed4f3 100644 --- a/tests/cuda/matrix_free_matrix_vector_09.cu +++ b/tests/cuda/matrix_free_matrix_vector_09.cu @@ -76,5 +76,5 @@ test() fe_degree, double, LinearAlgebra::CUDAWrappers::Vector, - fe_degree + 1>(dof, constraints); + fe_degree + 1>(dof, constraints, tria.n_active_cells()); } diff --git a/tests/cuda/matrix_free_matrix_vector_10.cu b/tests/cuda/matrix_free_matrix_vector_10.cu index fd26148e4e..63a887d5c0 100644 --- a/tests/cuda/matrix_free_matrix_vector_10.cu +++ b/tests/cuda/matrix_free_matrix_vector_10.cu @@ -115,11 +115,13 @@ test() mf_data.reinit( mapping, dof, constraints, quad, MPI_COMM_WORLD, additional_data); + const unsigned int coef_size = + tria.n_locally_owned_active_cells() * std::pow(fe_degree + 1, dim); MatrixFreeTest> - mf(mf_data); + mf(mf_data, coef_size); LinearAlgebra::distributed::Vector in_dev( owned_set, MPI_COMM_WORLD); LinearAlgebra::distributed::Vector out_dev( diff --git a/tests/cuda/matrix_free_matrix_vector_19.cu b/tests/cuda/matrix_free_matrix_vector_19.cu index 18c2f1b86b..c60eb40588 100644 --- a/tests/cuda/matrix_free_matrix_vector_19.cu +++ b/tests/cuda/matrix_free_matrix_vector_19.cu @@ -111,11 +111,13 @@ test() mf_data.reinit( mapping, dof, constraints, quad, MPI_COMM_WORLD, additional_data); + const unsigned int coef_size = + tria.n_locally_owned_active_cells() * std::pow(fe_degree + 1, dim); MatrixFreeTest> - mf(mf_data); + mf(mf_data, coef_size); LinearAlgebra::distributed::Vector in_dev( owned_set, MPI_COMM_WORLD); LinearAlgebra::distributed::Vector out_dev( diff --git a/tests/cuda/matrix_vector_common.h b/tests/cuda/matrix_vector_common.h index 662086f595..80dbb8aec5 100644 --- a/tests/cuda/matrix_vector_common.h +++ b/tests/cuda/matrix_vector_common.h @@ -19,6 +19,7 @@ // general, with and without hanging nodes). It also tests the multithreading #include +#include #include #include @@ -52,6 +53,20 @@ template void test(); +template +double +ConstantCoefficient(const Point & /*p*/) +{ + return 10.; +} + + +template +double +VaryingCoefficient(const Point &p) +{ + return 10. / (0.05 + 2. * p.square()); +} template void do_test(const DoFHandler & dof, - const AffineConstraints &constraints) + const AffineConstraints &constraints, + const unsigned int n_locally_owned_cells, + const bool constant_coefficient = true) { deallog << "Testing " << dof.get_fe().get_name() << std::endl; @@ -76,7 +93,10 @@ do_test(const DoFHandler & dof, mf_data.reinit(mapping, dof, constraints, quad, additional_data); const unsigned int n_dofs = dof.n_dofs(); - MatrixFreeTest mf(mf_data); + MatrixFreeTest mf( + mf_data, + n_locally_owned_cells * std::pow(n_q_points_1d, dim), + constant_coefficient); Vector in_host(n_dofs), out_host(n_dofs); LinearAlgebra::ReadWriteVector in(n_dofs), out(n_dofs); VectorType in_device(n_dofs); @@ -110,7 +130,7 @@ do_test(const DoFHandler & dof, dof.get_fe(), quadrature_formula, update_values | update_gradients | - update_JxW_values); + update_quadrature_points | update_JxW_values); const unsigned int dofs_per_cell = dof.get_fe().dofs_per_cell; const unsigned int n_q_points = quadrature_formula.size(); @@ -126,15 +146,22 @@ do_test(const DoFHandler & dof, fe_values.reinit(cell); for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) - for (unsigned int i = 0; i < dofs_per_cell; ++i) - { - for (unsigned int j = 0; j < dofs_per_cell; ++j) - cell_matrix(i, j) += ((fe_values.shape_grad(i, q_point) * - fe_values.shape_grad(j, q_point) + - 10. * fe_values.shape_value(i, q_point) * - fe_values.shape_value(j, q_point)) * - fe_values.JxW(q_point)); - } + { + const auto coef = + constant_coefficient ? + ConstantCoefficient(fe_values.quadrature_point(q_point)) : + VaryingCoefficient(fe_values.quadrature_point(q_point)); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + ((fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) + + coef * fe_values.shape_value(i, q_point) * + fe_values.shape_value(j, q_point)) * + fe_values.JxW(q_point)); + } + } cell->get_dof_indices(local_dof_indices); constraints.distribute_local_to_global(cell_matrix, @@ -169,7 +196,7 @@ main() { deallog.push("2d"); test<2, 1>(); - // test<2, 2>(); + // test<2, 2>(); test<2, 3>(); deallog.pop(); deallog.push("3d"); diff --git a/tests/cuda/matrix_vector_mf.h b/tests/cuda/matrix_vector_mf.h index 4deba79eae..4d62d78cec 100644 --- a/tests/cuda/matrix_vector_mf.h +++ b/tests/cuda/matrix_vector_mf.h @@ -26,14 +26,37 @@ #include "../tests.h" +__device__ unsigned int +compute_thread_id() +{ + const unsigned int thread_num_in_block = + threadIdx.x + blockDim.x * threadIdx.y + + blockDim.x * blockDim.y * threadIdx.z; + const unsigned int n_threads_per_block = blockDim.x * blockDim.y * blockDim.z; + const unsigned int blocks_num_in_grid = + blockIdx.x + gridDim.x * blockIdx.y + gridDim.x * gridDim.y * blockIdx.z; + + return thread_num_in_block + blocks_num_in_grid * n_threads_per_block; +} + + + template class HelmholtzOperatorQuad { public: + __device__ + HelmholtzOperatorQuad(Number coef) + : coef(coef) + {} + __device__ void operator()( CUDAWrappers::FEEvaluation * fe_eval, const unsigned int q_point) const; + +private: + Number coef; }; @@ -44,7 +67,7 @@ __device__ void HelmholtzOperatorQuad:: CUDAWrappers::FEEvaluation *fe_eval, const unsigned int q) const { - fe_eval->submit_value(Number(10) * fe_eval->get_value(q), q); + fe_eval->submit_value(coef * fe_eval->get_value(q), q); fe_eval->submit_gradient(fe_eval->get_gradient(q), q); } @@ -54,6 +77,10 @@ template class HelmholtzOperator { public: + HelmholtzOperator(Number *coefficient) + : coef(coefficient) + {} + __device__ void operator()( const unsigned int cell, @@ -67,6 +94,8 @@ public: dealii::Utilities::pow(fe_degree + 1, dim); static const unsigned int n_q_points = dealii::Utilities::pow(n_q_points_1d, dim); + + Number *coef; }; @@ -80,18 +109,68 @@ operator()(const unsigned int cell, const Number * src, Number * dst) const { + const unsigned int pos = CUDAWrappers::local_q_point_id( + cell, gpu_data, n_dofs_1d, n_q_points); + CUDAWrappers::FEEvaluation fe_eval( cell, gpu_data, shared_data); fe_eval.read_dof_values(src); fe_eval.evaluate(true, true); fe_eval.apply_quad_point_operations( - HelmholtzOperatorQuad()); + HelmholtzOperatorQuad(coef[pos])); fe_eval.integrate(true, true); fe_eval.distribute_local_to_global(dst); } +template +class VaryingCoefficientFunctor +{ +public: + VaryingCoefficientFunctor(Number *coefficient) + : coef(coefficient) + {} + + __device__ void + operator()( + const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data); + + static const unsigned int n_dofs_1d = fe_degree + 1; + static const unsigned int n_local_dofs = + dealii::Utilities::pow(fe_degree + 1, dim); + static const unsigned int n_q_points = + dealii::Utilities::pow(n_q_points_1d, dim); + +private: + Number *coef; +}; + + + +template +__device__ void +VaryingCoefficientFunctor:: +operator()(const unsigned int cell, + const typename CUDAWrappers::MatrixFree::Data *gpu_data) +{ + const unsigned int pos = CUDAWrappers::local_q_point_id( + cell, gpu_data, n_dofs_1d, n_q_points); + const auto q_point = + CUDAWrappers::get_quadrature_point(cell, gpu_data, n_dofs_1d); + + Number p_square = 0.; + for (unsigned int i = 0; i < dim; ++i) + { + Number coord = q_point[i]; + p_square += coord * coord; + } + coef[pos] = 10. / (0.05 + 2. * p_square); +} + + + template &data_in) - : data(data_in){}; + MatrixFreeTest(const CUDAWrappers::MatrixFree &data_in, + const unsigned int size, + const bool constant_coeff = true); void vmult(VectorType &dst, const VectorType &src); private: const CUDAWrappers::MatrixFree &data; + LinearAlgebra::CUDAWrappers::Vector coef; }; +template +MatrixFreeTest:: + MatrixFreeTest(const CUDAWrappers::MatrixFree &data_in, + const unsigned int size, + const bool constant_coeff) + : data(data_in) +{ + coef.reinit(size); + if (constant_coeff) + { + coef.add(10.); + } + else + { + VaryingCoefficientFunctor functor( + coef.get_values()); + data.evaluate_coefficients(functor); + } +} + template ::vmult( const VectorType &src) { dst = static_cast(0.); - HelmholtzOperator helmholtz_operator; + HelmholtzOperator helmholtz_operator( + coef.get_values()); data.cell_loop(helmholtz_operator, src, dst); data.copy_constrained_values(src, dst); } -- 2.39.5