From 4bd35833cc6d97e88f1327e6ba8e307b3c38769a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 20 Nov 2020 09:05:37 +0100 Subject: [PATCH] Add test --- .../matrix_vector_hp_no_template.cc | 227 ++++++++++++++++++ .../matrix_vector_hp_no_template.output | 5 + tests/matrix_free/matrix_vector_mf.h | 32 +++ 3 files changed, 264 insertions(+) create mode 100644 tests/matrix_free/matrix_vector_hp_no_template.cc create mode 100644 tests/matrix_free/matrix_vector_hp_no_template.output diff --git a/tests/matrix_free/matrix_vector_hp_no_template.cc b/tests/matrix_free/matrix_vector_hp_no_template.cc new file mode 100644 index 0000000000..4fad45acc6 --- /dev/null +++ b/tests/matrix_free/matrix_vector_hp_no_template.cc @@ -0,0 +1,227 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2013 - 2020 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 by comparing with the result of deal.II sparse +// matrix for hp DoFHandler on a hyperball mesh with hanging nodes and finite +// elements orders distributed randomly. + +#include + +#include +#include + +#include "../tests.h" + +#include "matrix_vector_common.h" + + + +template +class MatrixFreeTestHP +{ +public: + MatrixFreeTestHP(const MatrixFree &data_in) + : data(data_in){}; + + void + local_apply(const MatrixFree & data, + Vector & dst, + const Vector & src, + const std::pair &cell_range) const + { + for (unsigned int i = 1; i <= 9 - 2 * dim; ++i) // skip FE_Nothing + { + const auto cell_subrange = + data.create_cell_subrange_hp_by_index(cell_range, i); + + if (cell_subrange.second <= cell_subrange.first) + continue; + + helmholtz_operator_no_template(data, dst, src, cell_subrange, i, i); + } + } + + void + vmult(Vector &dst, const Vector &src) const + { + dst = 0; + data.cell_loop(&MatrixFreeTestHP::local_apply, this, dst, src); + }; + +private: + const MatrixFree &data; +}; + + + +template +void +test() +{ + if (fe_degree > 1) + return; + + typedef double number; + const SphericalManifold manifold; + Triangulation tria; + GridGenerator::hyper_ball(tria); + typename Triangulation::active_cell_iterator cell = tria.begin_active(), + endc = tria.end(); + for (; cell != endc; ++cell) + for (const unsigned int f : GeometryInfo::face_indices()) + if (cell->at_boundary(f)) + cell->face(f)->set_all_manifold_ids(0); + tria.set_manifold(0, manifold); + tria.refine_global(1); + + // refine a few cells + for (unsigned int i = 0; i < 11 - 3 * dim; ++i) + { + typename Triangulation::active_cell_iterator cell = + tria.begin_active(), + endc = tria.end(); + unsigned int counter = 0; + for (; cell != endc; ++cell, ++counter) + if (counter % (7 - i) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement(); + } + + const unsigned int max_degree = 9 - 2 * dim; + + hp::FECollection fe_collection; + hp::QCollection quadrature_collection; + hp::QCollection<1> quadrature_collection_mf; + + fe_collection.push_back(FE_Nothing()); + quadrature_collection.push_back(QGauss(1)); + quadrature_collection_mf.push_back(QGauss<1>(1)); + for (unsigned int deg = 1; deg <= max_degree; ++deg) + { + fe_collection.push_back(FE_Q(QGaussLobatto<1>(deg + 1))); + quadrature_collection.push_back(QGauss(deg + 1)); + quadrature_collection_mf.push_back(QGauss<1>(deg + 1)); + } + + hp::DoFHandler dof(tria); + // set the active FE index in a random order + { + typename hp::DoFHandler::active_cell_iterator cell = + dof.begin_active(), + endc = dof.end(); + for (; cell != endc; ++cell) + { + const unsigned int fe_index = Testing::rand() % max_degree; + cell->set_active_fe_index(fe_index + 1); + } + // We cannot set random cells to FE_Nothing. We get the following error + // The violated condition was: dominating_fe.n_dofs_per_face(face) <= + // subface_fe.n_dofs_per_face(face) + cell = dof.begin_active(); + cell->set_active_fe_index(0); + } + + // setup DoFs + dof.distribute_dofs(fe_collection); + AffineConstraints constraints; + DoFTools::make_hanging_node_constraints(dof, constraints); + VectorTools::interpolate_boundary_values(dof, + 0, + Functions::ZeroFunction(), + constraints); + constraints.close(); + DynamicSparsityPattern csp(dof.n_dofs(), dof.n_dofs()); + DoFTools::make_sparsity_pattern(dof, csp, constraints, false); + SparsityPattern sparsity; + sparsity.copy_from(csp); + SparseMatrix system_matrix(sparsity); + + // std::cout << "Number of cells: " << + // dof.get_triangulation().n_active_cells() << std::endl; std::cout << "Number + // of degrees of freedom: " << dof.n_dofs() << std::endl; std::cout << "Number + // of constraints: " << constraints.n_constraints() << std::endl; + + // set up MatrixFree + MatrixFree mf_data; + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + mf_data.reinit(dof, constraints, quadrature_collection_mf, data); + MatrixFreeTestHP mf(mf_data); + + // assemble sparse matrix with (\nabla v, + // \nabla u) + (v, 10 * u) + { + hp::FEValues hp_fe_values(fe_collection, + quadrature_collection, + update_values | update_gradients | + update_JxW_values); + FullMatrix cell_matrix; + std::vector local_dof_indices; + + typename hp::DoFHandler::active_cell_iterator cell = + dof.begin_active(), + endc = dof.end(); + for (; cell != endc; ++cell) + { + const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + + cell_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_matrix = 0; + hp_fe_values.reinit(cell); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + for (unsigned int q_point = 0; q_point < fe_values.n_quadrature_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)); + } + local_dof_indices.resize(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + + constraints.distribute_local_to_global(cell_matrix, + local_dof_indices, + system_matrix); + } + } + + // fill a right hand side vector with random + // numbers in unconstrained degrees of freedom + Vector src(dof.n_dofs()); + Vector result_spmv(src), result_mf(src); + + for (unsigned int i = 0; i < dof.n_dofs(); ++i) + { + if (constraints.is_constrained(i) == false) + src(i) = random_value(); + } + + // now perform matrix-vector product and check + // its correctness + system_matrix.vmult(result_spmv, src); + mf.vmult(result_mf, src); + + result_mf -= result_spmv; + const double diff_norm = result_mf.linfty_norm(); + deallog << "Norm of difference: " << diff_norm << std::endl << std::endl; +} diff --git a/tests/matrix_free/matrix_vector_hp_no_template.output b/tests/matrix_free/matrix_vector_hp_no_template.output new file mode 100644 index 0000000000..7252835322 --- /dev/null +++ b/tests/matrix_free/matrix_vector_hp_no_template.output @@ -0,0 +1,5 @@ + +DEAL:2d::Norm of difference: 0 +DEAL:2d:: +DEAL:3d::Norm of difference: 0 +DEAL:3d:: diff --git a/tests/matrix_free/matrix_vector_mf.h b/tests/matrix_free/matrix_vector_mf.h index afc1ebf6bf..8a77fccac9 100644 --- a/tests/matrix_free/matrix_vector_mf.h +++ b/tests/matrix_free/matrix_vector_mf.h @@ -56,6 +56,38 @@ helmholtz_operator(const MatrixFree &data, +template +void +helmholtz_operator_no_template( + const MatrixFree &data, + VectorType & dst, + const VectorType & src, + const std::pair & cell_range, + const unsigned int active_fe_index, + const unsigned int active_quad_index) +{ + typedef typename VectorType::value_type Number; + FEEvaluation fe_eval( + data, 0, 0, 0, active_fe_index, active_quad_index); + const unsigned int n_q_points = fe_eval.n_q_points; + + for (unsigned int cell = cell_range.first; cell < cell_range.second; ++cell) + { + fe_eval.reinit(cell); + fe_eval.read_dof_values(src); + fe_eval.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + for (unsigned int q = 0; q < n_q_points; ++q) + { + fe_eval.submit_value(Number(10) * fe_eval.get_value(q), q); + fe_eval.submit_gradient(fe_eval.get_gradient(q), q); + } + fe_eval.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + fe_eval.distribute_local_to_global(dst); + } +} + + + template