From b14f35ebc2b8f477289a061dec5cb5312af92acf Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 11 Aug 2019 21:08:43 +0200 Subject: [PATCH] Add test --- tests/matrix_free/loop_boundary.cc | 235 ++++++++++++++++++ ...h_p4est=true.with_mpi=true.mpirun=1.output | 71 ++++++ ...h_p4est=true.with_mpi=true.mpirun=2.output | 71 ++++++ 3 files changed, 377 insertions(+) create mode 100644 tests/matrix_free/loop_boundary.cc create mode 100644 tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output create mode 100644 tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output diff --git a/tests/matrix_free/loop_boundary.cc b/tests/matrix_free/loop_boundary.cc new file mode 100644 index 0000000000..a06429f196 --- /dev/null +++ b/tests/matrix_free/loop_boundary.cc @@ -0,0 +1,235 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 matrix-free loop with continuous +// elements when both a cell and a boundary function are given. In an initial +// implementation, we used to miss to exchange some ghost entries upon the +// compress() functionality within the matrix-free loop + +#include + +#include + +#include + +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +do_test(const unsigned int n_refine, const bool overlap_communication) +{ + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria, 0.1, 0.96); + tria.refine_global(n_refine); + + FE_Q fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + AffineConstraints constraints; + constraints.close(); + + MatrixFree matrix_free; + typename MatrixFree::AdditionalData add_data; + add_data.mapping_update_flags = update_values | update_gradients | + update_JxW_values | update_quadrature_points; + add_data.mapping_update_flags_boundary_faces = + update_values | update_JxW_values | update_quadrature_points; + add_data.overlap_communication_computation = overlap_communication; + matrix_free.reinit(dof_handler, + constraints, + QGauss<1>(fe.degree + 1), + add_data); + + LinearAlgebra::distributed::Vector in, ref, test; + matrix_free.initialize_dof_vector(in); + for (unsigned int i = 0; i < in.get_partitioner()->local_size(); ++i) + in.local_element(i) = in.get_partitioner()->local_to_global(i); + matrix_free.initialize_dof_vector(ref); + matrix_free.initialize_dof_vector(test); + + std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + cell_func = [](const MatrixFree & data, + LinearAlgebra::distributed::Vector & out, + const LinearAlgebra::distributed::Vector &in, + const std::pair &range) -> void { + FEEvaluation eval(data); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + eval.reinit(cell); + eval.gather_evaluate(in, false, true); + for (unsigned int q = 0; q < eval.n_q_points; ++q) + { + eval.submit_gradient(eval.get_gradient(q), q); + eval.submit_value(eval.quadrature_point(q).square(), q); + } + eval.integrate_scatter(true, true, out); + } + }; + + std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + boundary_func = + [](const MatrixFree & data, + LinearAlgebra::distributed::Vector & out, + const LinearAlgebra::distributed::Vector &in, + const std::pair & range) -> void { + FEFaceEvaluation eval(data, true); + + for (unsigned int face = range.first; face < range.second; ++face) + { + eval.reinit(face); + eval.gather_evaluate(in, true, false); + for (unsigned int q = 0; q < eval.n_q_points; ++q) + { + eval.submit_value(eval.quadrature_point(q).square() - + 6. * eval.get_value(q), + q); + } + eval.integrate_scatter(true, false, out); + } + }; + + // compute reference result + in.update_ghost_values(); + cell_func(matrix_free, + ref, + in, + std::make_pair(0U, matrix_free.n_cell_batches())); + boundary_func(matrix_free, + ref, + in, + std::make_pair(matrix_free.n_inner_face_batches(), + matrix_free.n_inner_face_batches() + + matrix_free.n_boundary_face_batches())); + ref.compress(VectorOperation::add); + in.zero_out_ghosts(); + + std::function &, + LinearAlgebra::distributed::Vector &, + const LinearAlgebra::distributed::Vector &, + const std::pair &)> + inner_face_func; + + + // compute result through loop with various update settings. Note that we do + // no compute on inner faces, so MatrixFree::DataAccessOnFaces::none + // should be enough. + matrix_free.loop(cell_func, + inner_face_func, + boundary_func, + test, + in, + true, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + + deallog << "Number of dofs: " << dof_handler.n_dofs() + << (overlap_communication ? " with overlap" : " without overlap") + << std::endl; + + test -= ref; + deallog << "Error loop 1: " << test.linfty_norm() << std::endl; + matrix_free.loop(cell_func, + inner_face_func, + boundary_func, + test, + in, + true, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::values); + test -= ref; + deallog << "Error loop 2: " << test.linfty_norm() << std::endl; + matrix_free.loop(cell_func, + inner_face_func, + boundary_func, + test, + in, + true, + MatrixFree::DataAccessOnFaces::none, + MatrixFree::DataAccessOnFaces::none); + test -= ref; + deallog << "Error loop 3: " << test.linfty_norm() << std::endl; + matrix_free.loop(cell_func, + inner_face_func, + boundary_func, + test, + in, + true, + MatrixFree::DataAccessOnFaces::values, + MatrixFree::DataAccessOnFaces::none); + test -= ref; + deallog << "Error loop 4: " << test.linfty_norm() << std::endl; + matrix_free.loop(cell_func, + inner_face_func, + boundary_func, + test, + in, + true, + MatrixFree::DataAccessOnFaces::unspecified, + MatrixFree::DataAccessOnFaces::unspecified); + test -= ref; + deallog << "Error loop 5: " << test.linfty_norm() << std::endl; + + // compute again, now only cell loop + ref = 0; + in.update_ghost_values(); + cell_func(matrix_free, + ref, + in, + std::make_pair(0U, matrix_free.n_cell_batches())); + ref.compress(VectorOperation::add); + in.zero_out_ghosts(); + + matrix_free.cell_loop(cell_func, test, in, true); + test -= ref; + deallog << "Error cell loop: " << test.linfty_norm() << std::endl; +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + mpi_initlog(); + + do_test<2, 1>(2, false); + do_test<2, 1>(2, true); + do_test<2, 1>(3, false); + do_test<2, 1>(3, true); + do_test<2, 2>(2, false); + do_test<2, 2>(2, true); + do_test<2, 2>(3, false); + do_test<2, 2>(3, true); + do_test<3, 2>(2, false); + do_test<3, 2>(2, true); +} diff --git a/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output new file mode 100644 index 0000000000..ea287438cb --- /dev/null +++ b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=1.output @@ -0,0 +1,71 @@ + +DEAL::Number of dofs: 25 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 25 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 289 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 289 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 729 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 729 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 diff --git a/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output new file mode 100644 index 0000000000..ea287438cb --- /dev/null +++ b/tests/matrix_free/loop_boundary.with_p4est=true.with_mpi=true.mpirun=2.output @@ -0,0 +1,71 @@ + +DEAL::Number of dofs: 25 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 25 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 81 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 289 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 289 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 729 without overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 +DEAL::Number of dofs: 729 with overlap +DEAL::Error loop 1: 0.00000 +DEAL::Error loop 2: 0.00000 +DEAL::Error loop 3: 0.00000 +DEAL::Error loop 4: 0.00000 +DEAL::Error loop 5: 0.00000 +DEAL::Error cell loop: 0.00000 -- 2.39.5