From 1a98526a81ec0d3f3ab3837108ba5f1a29f89e48 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 10 Dec 2021 20:10:59 +0100 Subject: [PATCH] Add tests for MatrixFree + SHMEM --- .../deal.II/matrix_free/evaluation_kernels.h | 2 + tests/matrix_free/ecl.h | 320 ++++++++++++++++++ tests/matrix_free/ecl_01.cc | 288 +--------------- ...=> ecl_01.with_p4est=true.mpirun=1.output} | 0 tests/matrix_free/ecl_03.cc | 65 ++++ .../ecl_03.with_p4est=true.mpirun=5.output | 5 + tests/matrix_free/ecl_04.cc | 205 +++++++++++ .../ecl_04.with_p4est=true.mpirun=5.output | 7 + 8 files changed, 609 insertions(+), 283 deletions(-) create mode 100644 tests/matrix_free/ecl.h rename tests/matrix_free/{ecl_01.output => ecl_01.with_p4est=true.mpirun=1.output} (100%) create mode 100644 tests/matrix_free/ecl_03.cc create mode 100644 tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output create mode 100644 tests/matrix_free/ecl_04.cc create mode 100644 tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h index dc84931d63..81ff490f76 100644 --- a/include/deal.II/matrix_free/evaluation_kernels.h +++ b/include/deal.II/matrix_free/evaluation_kernels.h @@ -3996,6 +3996,8 @@ namespace internal proc.value( temp1[reorientate(v, i)][v], vector_ptrs[v][index_array_nodal[v][i]]); + else if (integrate == false) + temp1[i][v] = 0.0; } } } diff --git a/tests/matrix_free/ecl.h b/tests/matrix_free/ecl.h new file mode 100644 index 0000000000..01efcb6598 --- /dev/null +++ b/tests/matrix_free/ecl.h @@ -0,0 +1,320 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +template +class CosineFunction : public Function +{ +public: + Number + value(const Point &p, const unsigned int component = 0) const + { + (void)component; + return std::cos(2 * numbers::PI * p[0]); + } +}; + +template > +void +test(const unsigned int n_refinements = 1, + const bool print_vector = true, + const MPI_Comm comm = MPI_COMM_SELF) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + GridGenerator::hyper_cube(tria, 0.0, 1.0, true); + + if (false) + { + std::vector::cell_iterator>> + periodic_faces; + + if (dim >= 1) + dealii::GridTools::collect_periodic_faces( + tria, 0, 1, 0, periodic_faces); + + if (dim >= 2) + dealii::GridTools::collect_periodic_faces( + tria, 2, 3, 1, periodic_faces); + + if (dim >= 3) + dealii::GridTools::collect_periodic_faces( + tria, 4, 5, 2, periodic_faces); + + tria.add_periodicity(periodic_faces); + } + + tria.refine_global(n_refinements); + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + QGauss<1> quad(n_points); + + AffineConstraints constraint; + + using MF = MatrixFree; + + typename MF::AdditionalData additional_data; + additional_data.mapping_update_flags = update_values | update_gradients; + additional_data.mapping_update_flags_inner_faces = + update_values | update_gradients; + additional_data.mapping_update_flags_boundary_faces = + update_values | update_gradients; + additional_data.tasks_parallel_scheme = + MF::AdditionalData::TasksParallelScheme::none; + // additional_data.mapping_update_flags_faces_by_cells = update_values; + // additional_data.hold_all_faces_to_owned_cells = true; + additional_data.communicator_sm = comm; + + MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data); + + MF matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + FEEvaluation phi( + matrix_free); + FEFaceEvaluation + phi_m(matrix_free, true); + FEFaceEvaluation + phi_p(matrix_free, false); + + + CosineFunction function; + VectorTools::interpolate(dof_handler, function, src); + + dst = 0.0; + + /** + * Face-centric loop + */ + matrix_free.template loop( + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(false, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(false, true); + phi.set_dof_values(dst); + } + }, + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int face = range.first; face < range.second; ++face) + { + phi_m.reinit(face); + phi_p.reinit(face); + + phi_m.read_dof_values(src); + phi_m.evaluate(true, true); + phi_p.read_dof_values(src); + phi_p.evaluate(true, true); + VectorizedArrayType sigmaF = + (std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_m.get_normal_vector(0) * + phi_p.inverse_jacobian(0))[dim - 1])) * + (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = + (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; + VectorizedArrayType average_valgrad = + phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q); + average_valgrad = + average_value * 2. * sigmaF - average_valgrad * 0.5; + phi_m.submit_normal_derivative(-average_value, q); + phi_p.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + phi_p.submit_value(-average_valgrad, q); + } + phi_m.integrate(true, true); + phi_m.distribute_local_to_global(dst); + phi_p.integrate(true, true); + phi_p.distribute_local_to_global(dst); + } + }, + [&](const auto &, auto &dst, const auto &src, const auto face_range) { + for (unsigned int face = face_range.first; face < face_range.second; + face++) + { + phi_m.reinit(face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) * + Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = phi_m.get_value(q); + VectorizedArrayType average_valgrad = + -phi_m.get_normal_derivative(q); + average_valgrad += average_value * sigmaF * 2.0; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + + phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients); + phi_m.distribute_local_to_global(dst); + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::gradients, + MF::DataAccessOnFaces::gradients); + + deallog << dst.l2_norm() << std::endl; + + if (print_vector) + dst.print(deallog.get_file_stream()); + + dst = 0.0; + + /** + * Element-centric loop + */ + matrix_free.template loop_cell_centric( + [&](const auto &, auto &dst, const auto &src, const auto range) { + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.read_dof_values(src); + phi.evaluate(false, true, false); + for (unsigned int q = 0; q < phi.n_q_points; ++q) + phi.submit_gradient(phi.get_gradient(q), q); + phi.integrate(false, true); + + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + { + auto bids = + matrix_free.get_faces_by_cells_boundary_id(cell, face); + + if (bids[0] != numbers::internal_face_boundary_id) + { + phi_m.reinit(cell, face); + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + VectorizedArrayType sigmaF = + std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) * + Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = phi_m.get_value(q); + VectorizedArrayType average_valgrad = + -phi_m.get_normal_derivative(q); + average_valgrad += average_value * sigmaF * 2.0; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + + phi_m.integrate(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) + phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q]; + } + else + { + phi_m.reinit(cell, face); + phi_p.reinit(cell, face); + + phi_m.read_dof_values(src); + phi_m.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + phi_p.read_dof_values(src); + phi_p.evaluate(EvaluationFlags::values | + EvaluationFlags::gradients); + + VectorizedArrayType sigmaF = + (std::abs((phi_m.get_normal_vector(0) * + phi_m.inverse_jacobian(0))[dim - 1]) + + std::abs((phi_m.get_normal_vector(0) * + phi_p.inverse_jacobian(0))[dim - 1])) * + (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + VectorizedArrayType average_value = + (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; + VectorizedArrayType average_valgrad = + phi_m.get_normal_derivative(q) + + phi_p.get_normal_derivative(q); + average_valgrad = + average_value * 2. * sigmaF - average_valgrad * 0.5; + phi_m.submit_normal_derivative(-average_value, q); + phi_m.submit_value(average_valgrad, q); + } + phi_m.integrate(EvaluationFlags::values | + EvaluationFlags::gradients); + + for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) + phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q]; + } + } + + phi.set_dof_values(dst); + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::gradients); + + deallog << dst.l2_norm() << std::endl; + + if (print_vector) + dst.print(deallog.get_file_stream()); +} diff --git a/tests/matrix_free/ecl_01.cc b/tests/matrix_free/ecl_01.cc index ad85f1d5a2..a335c637d5 100644 --- a/tests/matrix_free/ecl_01.cc +++ b/tests/matrix_free/ecl_01.cc @@ -30,295 +30,17 @@ #include -#include "../tests.h" +#include "ecl.h" // Compare the evaluation of a SIP Laplace operator with face- and element- // centric loops on a structured grid. -template -class CosineFunction : public Function -{ -public: - Number - value(const Point &p, const unsigned int component = 0) const - { - (void)component; - return std::cos(2 * numbers::PI * p[0]); - } -}; - -template > -void -test(const unsigned int n_refinements = 1) -{ - using VectorType = LinearAlgebra::distributed::Vector; - - Triangulation tria; - GridGenerator::hyper_cube(tria, 0.0, 1.0, true); - - if (false) - { - std::vector::cell_iterator>> - periodic_faces; - - if (dim >= 1) - dealii::GridTools::collect_periodic_faces( - tria, 0, 1, 0, periodic_faces); - - if (dim >= 2) - dealii::GridTools::collect_periodic_faces( - tria, 2, 3, 1, periodic_faces); - - if (dim >= 3) - dealii::GridTools::collect_periodic_faces( - tria, 4, 5, 2, periodic_faces); - - tria.add_periodicity(periodic_faces); - } - - tria.refine_global(n_refinements); - - FE_DGQ fe(fe_degree); - DoFHandler dof_handler(tria); - dof_handler.distribute_dofs(fe); - - MappingQ mapping(1); - QGauss<1> quad(n_points); - - AffineConstraints constraint; - - using MF = MatrixFree; - - typename MF::AdditionalData additional_data; - additional_data.mapping_update_flags = update_values | update_gradients; - additional_data.mapping_update_flags_inner_faces = - update_values | update_gradients; - additional_data.mapping_update_flags_boundary_faces = - update_values | update_gradients; - additional_data.tasks_parallel_scheme = - MF::AdditionalData::TasksParallelScheme::none; - // additional_data.mapping_update_flags_faces_by_cells = update_values; - // additional_data.hold_all_faces_to_owned_cells = true; - - MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data); - - MF matrix_free; - matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); - - VectorType src, dst; - - matrix_free.initialize_dof_vector(src); - matrix_free.initialize_dof_vector(dst); - - FEEvaluation phi( - matrix_free); - FEFaceEvaluation - phi_m(matrix_free, true); - FEFaceEvaluation - phi_p(matrix_free, false); - - - CosineFunction function; - VectorTools::interpolate(dof_handler, function, src); - - dst = 0.0; - - /** - * Face-centric loop - */ - matrix_free.template loop( - [&](const auto &, auto &dst, const auto &src, const auto range) { - for (unsigned int cell = range.first; cell < range.second; ++cell) - { - phi.reinit(cell); - phi.read_dof_values(src); - phi.evaluate(false, true, false); - for (unsigned int q = 0; q < phi.n_q_points; ++q) - phi.submit_gradient(phi.get_gradient(q), q); - phi.integrate(false, true); - phi.set_dof_values(dst); - } - }, - [&](const auto &, auto &dst, const auto &src, const auto range) { - for (unsigned int face = range.first; face < range.second; ++face) - { - phi_m.reinit(face); - phi_p.reinit(face); - - phi_m.read_dof_values(src); - phi_m.evaluate(true, true); - phi_p.read_dof_values(src); - phi_p.evaluate(true, true); - VectorizedArrayType sigmaF = - (std::abs((phi_m.get_normal_vector(0) * - phi_m.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_m.get_normal_vector(0) * - phi_p.inverse_jacobian(0))[dim - 1])) * - (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); - - for (unsigned int q = 0; q < phi_m.n_q_points; ++q) - { - VectorizedArrayType average_value = - (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; - VectorizedArrayType average_valgrad = - phi_m.get_normal_derivative(q) + phi_p.get_normal_derivative(q); - average_valgrad = - average_value * 2. * sigmaF - average_valgrad * 0.5; - phi_m.submit_normal_derivative(-average_value, q); - phi_p.submit_normal_derivative(-average_value, q); - phi_m.submit_value(average_valgrad, q); - phi_p.submit_value(-average_valgrad, q); - } - phi_m.integrate(true, true); - phi_m.distribute_local_to_global(dst); - phi_p.integrate(true, true); - phi_p.distribute_local_to_global(dst); - } - }, - [&](const auto &, auto &dst, const auto &src, const auto face_range) { - for (unsigned int face = face_range.first; face < face_range.second; - face++) - { - phi_m.reinit(face); - phi_m.read_dof_values(src); - phi_m.evaluate(EvaluationFlags::values | EvaluationFlags::gradients); - VectorizedArrayType sigmaF = - std::abs((phi_m.get_normal_vector(0) * - phi_m.inverse_jacobian(0))[dim - 1]) * - Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; - - for (unsigned int q = 0; q < phi_m.n_q_points; ++q) - { - VectorizedArrayType average_value = phi_m.get_value(q); - VectorizedArrayType average_valgrad = - -phi_m.get_normal_derivative(q); - average_valgrad += average_value * sigmaF * 2.0; - phi_m.submit_normal_derivative(-average_value, q); - phi_m.submit_value(average_valgrad, q); - } - - phi_m.integrate(EvaluationFlags::values | EvaluationFlags::gradients); - phi_m.distribute_local_to_global(dst); - } - }, - dst, - src, - false, - MF::DataAccessOnFaces::gradients, - MF::DataAccessOnFaces::gradients); - - deallog << dst.l2_norm() << std::endl; - dst.print(deallog.get_file_stream()); - - dst = 0.0; - - /** - * Element-centric loop - */ - matrix_free.template loop_cell_centric( - [&](const auto &, auto &dst, const auto &src, const auto range) { - for (unsigned int cell = range.first; cell < range.second; ++cell) - { - phi.reinit(cell); - phi.read_dof_values(src); - phi.evaluate(false, true, false); - for (unsigned int q = 0; q < phi.n_q_points; ++q) - phi.submit_gradient(phi.get_gradient(q), q); - phi.integrate(false, true); - - for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; - ++face) - { - auto bids = - matrix_free.get_faces_by_cells_boundary_id(cell, face); - - if (bids[0] != numbers::internal_face_boundary_id) - { - phi_m.reinit(cell, face); - phi_m.read_dof_values(src); - phi_m.evaluate(EvaluationFlags::values | - EvaluationFlags::gradients); - VectorizedArrayType sigmaF = - std::abs((phi_m.get_normal_vector(0) * - phi_m.inverse_jacobian(0))[dim - 1]) * - Number(std::max(fe_degree, 1) * (fe_degree + 1.0)) * 2.0; - - for (unsigned int q = 0; q < phi_m.n_q_points; ++q) - { - VectorizedArrayType average_value = phi_m.get_value(q); - VectorizedArrayType average_valgrad = - -phi_m.get_normal_derivative(q); - average_valgrad += average_value * sigmaF * 2.0; - phi_m.submit_normal_derivative(-average_value, q); - phi_m.submit_value(average_valgrad, q); - } - - phi_m.integrate(EvaluationFlags::values | - EvaluationFlags::gradients); - - for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) - phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q]; - } - else - { - phi_m.reinit(cell, face); - phi_p.reinit(cell, face); - - phi_m.read_dof_values(src); - phi_m.evaluate(EvaluationFlags::values | - EvaluationFlags::gradients); - phi_p.read_dof_values(src); - phi_p.evaluate(EvaluationFlags::values | - EvaluationFlags::gradients); - - VectorizedArrayType sigmaF = - (std::abs((phi_m.get_normal_vector(0) * - phi_m.inverse_jacobian(0))[dim - 1]) + - std::abs((phi_m.get_normal_vector(0) * - phi_p.inverse_jacobian(0))[dim - 1])) * - (Number)(std::max(fe_degree, 1) * (fe_degree + 1.0)); - - for (unsigned int q = 0; q < phi_m.n_q_points; ++q) - { - VectorizedArrayType average_value = - (phi_m.get_value(q) - phi_p.get_value(q)) * 0.5; - VectorizedArrayType average_valgrad = - phi_m.get_normal_derivative(q) + - phi_p.get_normal_derivative(q); - average_valgrad = - average_value * 2. * sigmaF - average_valgrad * 0.5; - phi_m.submit_normal_derivative(-average_value, q); - phi_m.submit_value(average_valgrad, q); - } - phi_m.integrate(EvaluationFlags::values | - EvaluationFlags::gradients); - - for (unsigned int q = 0; q < phi.dofs_per_cell; ++q) - phi.begin_dof_values()[q] += phi_m.begin_dof_values()[q]; - } - } - - phi.set_dof_values(dst); - } - }, - dst, - src, - false, - MF::DataAccessOnFaces::gradients); - - deallog << dst.l2_norm() << std::endl; - dst.print(deallog.get_file_stream()); -} - int -main() +main(int argc, char **argv) { - initlog(); + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); test<2, 2, 3, double, VectorizedArray>(); } diff --git a/tests/matrix_free/ecl_01.output b/tests/matrix_free/ecl_01.with_p4est=true.mpirun=1.output similarity index 100% rename from tests/matrix_free/ecl_01.output rename to tests/matrix_free/ecl_01.with_p4est=true.mpirun=1.output diff --git a/tests/matrix_free/ecl_03.cc b/tests/matrix_free/ecl_03.cc new file mode 100644 index 0000000000..62ce5fbb8d --- /dev/null +++ b/tests/matrix_free/ecl_03.cc @@ -0,0 +1,65 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include "ecl.h" + + +// Like ecl_01.cc but running with shared-memory MPI enabled. + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + + // test without shared memory + { + test<2, 2, 3, double, VectorizedArray>(6, false, MPI_COMM_SELF); + } + + // test with shared memory + { + const auto rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + MPI_Comm subcommunicator; + MPI_Comm_split_type(MPI_COMM_WORLD, + MPI_COMM_TYPE_SHARED, + rank, + MPI_INFO_NULL, + &subcommunicator); + + test<2, 2, 3, double, VectorizedArray>(6, false, subcommunicator); + + MPI_Comm_free(&subcommunicator); + } +} diff --git a/tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output b/tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output new file mode 100644 index 0000000000..62b728caba --- /dev/null +++ b/tests/matrix_free/ecl_03.with_p4est=true.mpirun=5.output @@ -0,0 +1,5 @@ + +DEAL::210.386 +DEAL::210.386 +DEAL::210.386 +DEAL::210.386 diff --git a/tests/matrix_free/ecl_04.cc b/tests/matrix_free/ecl_04.cc new file mode 100644 index 0000000000..e746eca834 --- /dev/null +++ b/tests/matrix_free/ecl_04.cc @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2020 - 2021 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. +// +// --------------------------------------------------------------------- + +#include + +#include + +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + + +// Test FEFaceEvaluation::gather_evaluate() for values, ECL, shared-memory MPI, +// and different geometries (i.e., different orientations). + +template +class ExactSolution : public Function +{ +public: + ExactSolution() + {} + + virtual double + value(const Point &p, const unsigned int /*component*/ = 0) const + { + double result = p[0]; + + for (unsigned i = 1; i < dim; ++i) + result *= p[i]; + + return result; + } +}; + + +template > +void +test(const unsigned int geometry, const MPI_Comm comm = MPI_COMM_SELF) +{ + using VectorType = LinearAlgebra::distributed::Vector; + + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + + if (geometry == 0) + { + GridGenerator::hyper_cube(tria); + + if (dim == 2) + tria.refine_global(4); + else + tria.refine_global(3); + } + else if (geometry == 1) + { + GridGenerator::hyper_shell(tria, Point(), 0.5, 1.0); + + if (dim == 2) + tria.refine_global(4); + else + tria.refine_global(3); + } + else if (geometry == 2) + { + GridGenerator::hyper_ball(tria, Point(), 1.0); + + if (dim == 2) + tria.refine_global(4); + else + tria.refine_global(3); + } + + FE_DGQ fe(fe_degree); + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MappingQ mapping(1); + QGauss<1> quad(n_points); + + AffineConstraints constraint; + + using MF = MatrixFree; + + typename MF::AdditionalData additional_data; + additional_data.mapping_update_flags = update_values; + additional_data.mapping_update_flags_inner_faces = update_values; + additional_data.mapping_update_flags_boundary_faces = update_values; + additional_data.mapping_update_flags_faces_by_cells = update_values; + additional_data.hold_all_faces_to_owned_cells = true; + additional_data.communicator_sm = comm; + + MatrixFreeTools::categorize_by_boundary_ids(tria, additional_data); + + MF matrix_free; + matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data); + + VectorType src, dst; + + matrix_free.initialize_dof_vector(src); + matrix_free.initialize_dof_vector(dst); + + VectorTools::interpolate(dof_handler, ExactSolution(), src); + + dst = 0.0; + + /** + * Element-centric loop + */ + matrix_free.template loop_cell_centric( + [&](const auto &, auto &dst, const auto &src, const auto range) { + FEFaceEvaluation + phi_m(matrix_free, true); + FEFaceEvaluation + phi_p(matrix_free, false); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + face++) + { + auto bids = + matrix_free.get_faces_by_cells_boundary_id(cell, face); + + if (bids[0] != numbers::internal_face_boundary_id) + continue; + + phi_m.reinit(cell, face); + phi_p.reinit(cell, face); + + phi_m.gather_evaluate(src, EvaluationFlags::values); + phi_p.gather_evaluate(src, EvaluationFlags::values); + + for (unsigned int q = 0; q < phi_m.n_q_points; ++q) + { + const auto u_minus = phi_m.get_value(q); + const auto u_plus = phi_p.get_value(q); + + for (unsigned int v = 0; v < VectorizedArray::size(); + ++v) + { + Assert(std::abs(u_minus[v] - u_plus[v]) < 1e-10, + ExcMessage("Entries do not match!")); + } + } + } + } + }, + dst, + src, + false, + MF::DataAccessOnFaces::values); + + deallog << "OK" << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + mpi_initlog(); + + const auto rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + + MPI_Comm subcommunicator; + MPI_Comm_split_type(MPI_COMM_WORLD, + MPI_COMM_TYPE_SHARED, + rank, + MPI_INFO_NULL, + &subcommunicator); + + for (unsigned int i = 0; i < 3; ++i) + test<2, 3, 4, double>(i, subcommunicator); + + for (unsigned int i = 0; i < 3; ++i) + test<3, 3, 4, double>(i, subcommunicator); + + MPI_Comm_free(&subcommunicator); +} diff --git a/tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output b/tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output new file mode 100644 index 0000000000..0aa61ff573 --- /dev/null +++ b/tests/matrix_free/ecl_04.with_p4est=true.mpirun=5.output @@ -0,0 +1,7 @@ + +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK +DEAL::OK -- 2.39.5