From: Peter Munch Date: Fri, 9 Dec 2022 08:21:45 +0000 (+0100) Subject: MatrixFree: disable overlap of communication/computation for FE_Nothing X-Git-Tag: v9.6.0-rc1~146^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14553%2Fhead;p=dealii.git MatrixFree: disable overlap of communication/computation for FE_Nothing --- diff --git a/include/deal.II/matrix_free/matrix_free.templates.h b/include/deal.II/matrix_free/matrix_free.templates.h index 4cca89c2fe..2f0611a0d0 100644 --- a/include/deal.II/matrix_free/matrix_free.templates.h +++ b/include/deal.II/matrix_free/matrix_free.templates.h @@ -1882,6 +1882,24 @@ MatrixFree::initialize_indices( b); } + bool overlap_communication_computation = + additional_data.overlap_communication_computation; + + if (overlap_communication_computation) + { + for (unsigned int no = 0; no < dof_handlers.size(); ++no) + for (unsigned int fe_no = 0; + fe_no < dof_handlers[no]->get_fe_collection().size(); + ++fe_no) + if ((additional_data.mapping_update_flags_inner_faces != + update_default) && + (dof_handlers[no]->get_fe(fe_no).n_dofs_per_cell() == 0)) + // disable overlapping of communication and computation if + // FE_Nothing is used and face integrals are performed + // see #14342 and #14553. + overlap_communication_computation = false; + } + const unsigned int n_lanes = VectorizedArrayType::size(); task_info.vectorization_length = n_lanes; internal::MatrixFreeFunctions::ConstraintValues constraint_values; @@ -1897,7 +1915,7 @@ MatrixFree::initialize_indices( additional_data.cell_vectorization_categories_strict, do_face_integrals, additional_data.mapping_update_flags_inner_faces != update_default, - additional_data.overlap_communication_computation, + overlap_communication_computation, task_info, cell_level_index, dof_info, diff --git a/tests/matrix_free/fe_nothing_02.cc b/tests/matrix_free/fe_nothing_02.cc new file mode 100644 index 0000000000..1d2b8e2739 --- /dev/null +++ b/tests/matrix_free/fe_nothing_02.cc @@ -0,0 +1,268 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 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 + +#include "../tests.h" + +// In https://github.com/dealii/dealii/pull/14342, we have observed that +// overlapping comunication/computation does not work properly since +// ghost faces might be assigned to the wrong partition if +// FE_Nothing is used. This is a test which shows this problem by +// running the same problem on a large mesh where 50% of the +// cells are deactivated and on a smaller mesh. + +using namespace dealii; + +template +class Fu : public Function +{ +public: + double + value(const Point &point, const unsigned int = 0) const + { + return std::sin(point[0] * 2 * numbers::PI) * + std::sin(point[1] * 2 * numbers::PI); + } +}; + +template +void +test(const unsigned int n_refinements) +{ + using Number = double; + using VectorizedArrayType = VectorizedArray; + using VectorType = LinearAlgebra::distributed::Vector; + using FECellIntegrator = + FEEvaluation; + using FEFaceIntegrator = + FEFaceEvaluation; + + for (unsigned int i = 0; i < 2; ++i) + { + parallel::distributed::Triangulation tria(MPI_COMM_WORLD); + if (i == 0) + GridGenerator::subdivided_hyper_rectangle( + tria, {2, 2}, {0.0, 0.0}, {1.0, 1.0}, true); + else + GridGenerator::subdivided_hyper_rectangle( + tria, {2, 1}, {0.0, 0.0}, {1.0, 0.5}, true); + + tria.refine_global(n_refinements); + + // repartition such that interface between ranks intersects the + // FE_Q-FE_Nothing interface + tria.signals.weight.connect([](const auto &cell, const auto /*status*/) { + if (cell->center()[1] < 0.5) + return 10u; + else + return 8u; + }); + + tria.repartition(); + + DoFHandler dof_handler(tria); + + for (const auto &cell : + filter_iterators(dof_handler.active_cell_iterators(), + IteratorFilters::LocallyOwnedCell())) + { + if (cell->center()[1] < 0.5) + cell->set_active_fe_index(0); + else + cell->set_active_fe_index(1); + } + + dof_handler.distribute_dofs( + hp::FECollection(FE_Q(1), FE_Nothing(1))); + + typename MatrixFree::AdditionalData + data; + data.tasks_parallel_scheme = + MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_gradients; + data.mapping_update_flags_boundary_faces = update_gradients; + data.mapping_update_flags_inner_faces = update_gradients; + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handler, + AffineConstraints(), + QGauss(2), + data); + + VectorType dst, src; + matrix_free.initialize_dof_vector(dst); + matrix_free.initialize_dof_vector(src); + + VectorTools::interpolate(dof_handler, Fu(), src); + + const auto cell_operation = [&](const auto &matrix_free, + auto &dst, + const auto &src, + const auto range) { + if (matrix_free.get_cell_range_category(range) != 0) + return; + FECellIntegrator phi(matrix_free); + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + phi.reinit(cell); + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (const auto q : phi.quadrature_point_indices()) + phi.submit_gradient(phi.get_gradient(q), q); + + phi.integrate_scatter(EvaluationFlags::gradients, dst); + } + }; + + const auto face_operation = [&](const auto &matrix_free, + auto &dst, + const auto &src, + const auto range) { + const auto category = matrix_free.get_face_range_category(range); + + const unsigned int type = + static_cast(category.first == 0) + + static_cast(category.second == 0); + + if (type != 1) + return; + + const bool is_interior_face = category.first == 0; + FEFaceIntegrator phi(matrix_free, is_interior_face); + + for (unsigned int face = range.first; face < range.second; ++face) + { + phi.reinit(face); + + const double scaling = + phi.at_boundary() ? (phi.boundary_id() + 1.0) : (2.0 * dim); + + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (const auto q : phi.quadrature_point_indices()) + { + auto gradient = phi.get_gradient(q); + auto normal = phi.get_normal_vector(q); + + if (is_interior_face == false) // fix sign! + normal *= -1.0; + + phi.submit_value(gradient * normal * scaling, q); + } + + phi.integrate_scatter(EvaluationFlags::values, dst); + } + }; + + const auto boundary_operation = + [&](const auto &matrix_free, auto &, const auto &, const auto range) { + const auto category = matrix_free.get_face_range_category(range); + const bool is_interior_face = category.first == 0; + + if (!is_interior_face) + return; + FEFaceIntegrator phi(matrix_free, is_interior_face); + + for (unsigned int face = range.first; face < range.second; ++face) + { + phi.reinit(face); + + const double scaling = + phi.at_boundary() ? (phi.boundary_id() + 1.0) : (2.0 * dim); + + phi.gather_evaluate(src, EvaluationFlags::gradients); + + for (const auto q : phi.quadrature_point_indices()) + { + auto gradient = phi.get_gradient(q); + auto normal = phi.get_normal_vector(q); + + if (is_interior_face == false) // fix sign! + normal *= -1.0; + + phi.submit_value(gradient * normal * scaling, q); + } + + phi.integrate_scatter(EvaluationFlags::values, dst); + } + }; + + matrix_free.template cell_loop(cell_operation, + dst, + src, + true); + + deallog << dst.l2_norm() << std::endl; + + matrix_free.template loop( + cell_operation, face_operation, boundary_operation, dst, src, true); + + deallog << dst.l2_norm() << std::endl; + + + // optional output visualization + if (true) + { + DataOut data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(src, "src"); + data_out.add_data_vector(dst, "dst"); + + dealii::Vector partitioning(tria.n_active_cells()); + for (unsigned int j = 0; j < partitioning.size(); ++j) + { + partitioning(j) = tria.locally_owned_subdomain(); + } + data_out.add_data_vector(partitioning, "partitioning"); + + data_out.build_patches(); + + DataOutBase::VtkFlags output_flags; + output_flags.cycle = i; + data_out.set_flags(output_flags); + + data_out.write_vtu_in_parallel("solution-" + Utilities::to_string(i) + + ".vtu", + MPI_COMM_WORLD); + } + } +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1); + mpi_initlog(); + + test<2>(3); +} diff --git a/tests/matrix_free/fe_nothing_02.with_p4est=true.mpirun=3.output b/tests/matrix_free/fe_nothing_02.with_p4est=true.mpirun=3.output new file mode 100644 index 0000000000..9b6a91ce94 --- /dev/null +++ b/tests/matrix_free/fe_nothing_02.with_p4est=true.mpirun=3.output @@ -0,0 +1,5 @@ + +DEAL::2.44839 +DEAL::7.42519 +DEAL::2.44839 +DEAL::7.42519