From 09c91ba2886e2a37386931ff0a2e554107f6d392 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 29 May 2020 18:52:54 +0000 Subject: [PATCH] Add test reproducing a bug Bug: only one MatrixFree object is valid at a given time --- tests/cuda/matrix_free_multiple_objects.cu | 216 ++++++++++++++++++ .../cuda/matrix_free_multiple_objects.output | 5 + 2 files changed, 221 insertions(+) create mode 100644 tests/cuda/matrix_free_multiple_objects.cu create mode 100644 tests/cuda/matrix_free_multiple_objects.output diff --git a/tests/cuda/matrix_free_multiple_objects.cu b/tests/cuda/matrix_free_multiple_objects.cu new file mode 100644 index 0000000000..5ce71310c0 --- /dev/null +++ b/tests/cuda/matrix_free_multiple_objects.cu @@ -0,0 +1,216 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Reproduce a bug where only one matrix-free object is valid + +#include +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +#include "matrix_vector_mf.h" + +template +void +do_test(const DoFHandler<2> & dof, + MatrixFreeTest<2, + fe_degree, + double, + LinearAlgebra::CUDAWrappers::Vector, + n_q_points_1d> & mf, + unsigned int n_dofs, + MappingQGeneric<2> & mapping, + const AffineConstraints &constraints) +{ + Vector in_host(n_dofs), out_host(n_dofs); + LinearAlgebra::ReadWriteVector in(n_dofs), out(n_dofs); + LinearAlgebra::CUDAWrappers::Vector in_device(n_dofs); + LinearAlgebra::CUDAWrappers::Vector out_device(n_dofs); + + for (unsigned int i = 0; i < n_dofs; ++i) + { + if (constraints.is_constrained(i)) + continue; + const double entry = Testing::rand() / (double)RAND_MAX; + in(i) = entry; + } + + in_device.import(in, VectorOperation::insert); + mf.vmult(out_device, in_device); + cudaDeviceSynchronize(); + out.import(out_device, VectorOperation::insert); + + // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u) + SparsityPattern sparsity; + { + DynamicSparsityPattern csp(n_dofs, n_dofs); + DoFTools::make_sparsity_pattern(dof, csp, constraints, true); + sparsity.copy_from(csp); + } + SparseMatrix sparse_matrix(sparsity); + { + QGauss<2> quadrature_formula(n_q_points_1d); + + FEValues<2> fe_values(mapping, + dof.get_fe(), + quadrature_formula, + update_values | update_gradients | + 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(); + + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + + typename DoFHandler<2>::active_cell_iterator cell = dof.begin_active(), + endc = dof.end(); + for (; cell != endc; ++cell) + { + cell_matrix = 0; + fe_values.reinit(cell); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const auto coef = 10.; + 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, + local_dof_indices, + sparse_matrix); + } + } + for (unsigned i = 0; i < n_dofs; ++i) + in_host[i] = in[i]; + sparse_matrix.vmult(out_host, in_host); + + double out_norm = 0.; + for (unsigned i = 0; i < n_dofs; ++i) + out_norm += std::pow(out[i] - out_host[i], 2); + const double diff_norm = std::sqrt(out_norm) / out_host.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + +int +main() +{ + initlog(); + deallog.depth_console(0); + + deallog << std::setprecision(3); + + init_cuda(); + + Triangulation<2> tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(5 - 2); + AffineConstraints constraints; + constraints.close(); + bool constant_coefficient = true; + + // Create the first MatrixFree object + constexpr unsigned int fe_degree_1 = 1; + constexpr unsigned int n_q_points_1d_1 = fe_degree_1 + 1; + FE_Q<2> fe_1(fe_degree_1); + DoFHandler<2> dof_1(tria); + dof_1.distribute_dofs(fe_1); + MappingQGeneric<2> mapping_1(fe_degree_1); + CUDAWrappers::MatrixFree<2, double> mf_data_1; + CUDAWrappers::MatrixFree<2, double>::AdditionalData additional_data_1; + additional_data_1.mapping_update_flags = update_values | update_gradients | + update_JxW_values | + update_quadrature_points; + const QGauss<1> quad_1(n_q_points_1d_1); + mf_data_1.reinit(mapping_1, dof_1, constraints, quad_1, additional_data_1); + const unsigned int n_dofs_1 = dof_1.n_dofs(); + MatrixFreeTest<2, + fe_degree_1, + double, + LinearAlgebra::CUDAWrappers::Vector, + n_q_points_1d_1> + mf_1(mf_data_1, + n_dofs_1 * std::pow(n_q_points_1d_1, 2), + constant_coefficient); + + // Create the second MatrixFree object + constexpr unsigned int fe_degree_2 = 2; + constexpr unsigned int n_q_points_1d_2 = fe_degree_2 + 1; + FE_Q<2> fe_2(fe_degree_2); + DoFHandler<2> dof_2(tria); + dof_2.distribute_dofs(fe_2); + MappingQGeneric<2> mapping_2(fe_degree_2); + CUDAWrappers::MatrixFree<2, double> mf_data_2; + CUDAWrappers::MatrixFree<2, double>::AdditionalData additional_data_2; + additional_data_2.mapping_update_flags = update_values | update_gradients | + update_JxW_values | + update_quadrature_points; + const QGauss<1> quad_2(n_q_points_1d_2); + mf_data_2.reinit(mapping_2, dof_2, constraints, quad_2, additional_data_2); + const unsigned int n_dofs_2 = dof_2.n_dofs(); + MatrixFreeTest<2, + fe_degree_2, + double, + LinearAlgebra::CUDAWrappers::Vector, + n_q_points_1d_2> + mf_2(mf_data_2, + n_dofs_2 * std::pow(n_q_points_1d_2, 2), + constant_coefficient); + + // Perform MV with the first object + do_test( + dof_1, mf_1, n_dofs_1, mapping_1, constraints); + + // Perform MV with the second object + do_test( + dof_2, mf_2, n_dofs_2, mapping_2, constraints); + + return 0; +} diff --git a/tests/cuda/matrix_free_multiple_objects.output b/tests/cuda/matrix_free_multiple_objects.output new file mode 100644 index 0000000000..0f8bb9f28a --- /dev/null +++ b/tests/cuda/matrix_free_multiple_objects.output @@ -0,0 +1,5 @@ + +DEAL::Norm of difference: 1.57e-15 +DEAL:: +DEAL::Norm of difference: 2.37e-15 +DEAL:: -- 2.39.5