From 0323bac571d2f195af22f106cdfffec25513cfac Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 20 Sep 2018 19:12:34 +0200 Subject: [PATCH] Fix zeroing vector without any non-constrained dofs. --- .../deal.II/matrix_free/dof_info.templates.h | 4 + tests/matrix_free/matrix_vector_24.cc | 194 ++++++++++++++++++ tests/matrix_free/matrix_vector_24.output | 7 + 3 files changed, 205 insertions(+) create mode 100644 tests/matrix_free/matrix_vector_24.cc create mode 100644 tests/matrix_free/matrix_vector_24.output diff --git a/include/deal.II/matrix_free/dof_info.templates.h b/include/deal.II/matrix_free/dof_info.templates.h index 58ceba9879..b46c0e91e4 100644 --- a/include/deal.II/matrix_free/dof_info.templates.h +++ b/include/deal.II/matrix_free/dof_info.templates.h @@ -1197,6 +1197,10 @@ namespace internal } } } + // ensure that all indices are touched at least during the last round + for (auto &i : touched_by) + if (i == numbers::invalid_unsigned_int) + i = task_info.cell_partition_data.back() - 1; vector_zero_range_list_index.resize( 1 + task_info diff --git a/tests/matrix_free/matrix_vector_24.cc b/tests/matrix_free/matrix_vector_24.cc new file mode 100644 index 0000000000..44275e5b20 --- /dev/null +++ b/tests/matrix_free/matrix_vector_24.cc @@ -0,0 +1,194 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2016 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. +// +// --------------------------------------------------------------------- + + + +// this function tests the correctness of the implementation of matrix free +// matrix-vector products for a special case of linear elements where all DoFs +// are subject to constraints + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" +#include "matrix_vector_mf.h" + + + +template , + int n_q_points_1d = fe_degree + 1> +class MatrixFreeVariant +{ +public: + typedef VectorizedArray vector_t; + + MatrixFreeVariant(const MatrixFree &data_in) + : data(data_in) + {} + + void + vmult(VectorType &dst, const VectorType &src) const + { + const std::function< + void(const MatrixFree &, + VectorType &, + const VectorType &, + const std::pair &)> + wrap = helmholtz_operator; + data.cell_loop(wrap, dst, src, true); + for (auto i : data.get_constrained_dofs()) + dst(i) += src(i); + } + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + + FE_Q fe(fe_degree); + DoFHandler dof(tria); + dof.distribute_dofs(fe); + AffineConstraints constraints; + VectorTools::interpolate_boundary_values(dof, + 0, + ZeroFunction(), + constraints); + constraints.close(); + + deallog << "Testing " << dof.get_fe().get_name() << std::endl; + + using number = double; + + MatrixFree mf_data; + { + const QGauss<1> quad(fe_degree + 1); + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + mf_data.reinit(dof, constraints, quad, data); + } + + MatrixFreeVariant, + fe_degree + 1> + mf(mf_data); + LinearAlgebra::distributed::Vector in(dof.n_dofs()), + out(dof.n_dofs()); + LinearAlgebra::distributed::Vector out_dist(in); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + const double entry = random_value(); + in(i) = entry; + } + + out = in; + mf.vmult(out, in); + + // assemble sparse matrix with (\nabla v, \nabla u) + (v, 10 * u) + SparsityPattern sparsity; + { + DynamicSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern(dof, csp, constraints, true); + sparsity.copy_from(csp); + } + SparseMatrix sparse_matrix(sparsity); + { + QGauss quadrature_formula(fe_degree + 1); + + FEValues fe_values(dof.get_fe(), + quadrature_formula, + update_values | update_gradients | + 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::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) + 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)); + } + + cell->get_dof_indices(local_dof_indices); + constraints.distribute_local_to_global(cell_matrix, + local_dof_indices, + sparse_matrix); + } + } + // set matrix entries to constrained rows to 1 for consistency with + // matrix-free + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + if (constraints.is_constrained(i)) + sparse_matrix.set(i, i, 1); + + sparse_matrix.vmult(out_dist, in); + out -= out_dist; + const double diff_norm = out.linfty_norm() / out_dist.linfty_norm(); + + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} + + +int +main() +{ + initlog(); + deallog << std::setprecision(3); + + test<2, 1>(); + test<3, 1>(); +} diff --git a/tests/matrix_free/matrix_vector_24.output b/tests/matrix_free/matrix_vector_24.output new file mode 100644 index 0000000000..92d8da9226 --- /dev/null +++ b/tests/matrix_free/matrix_vector_24.output @@ -0,0 +1,7 @@ + +DEAL::Testing FE_Q<2>(1) +DEAL::Norm of difference: 0.00 +DEAL:: +DEAL::Testing FE_Q<3>(1) +DEAL::Norm of difference: 0.00 +DEAL:: -- 2.39.5