From 8934fff64c55ee40ed5cbb7c50f36d5cd342b068 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sat, 18 Jun 2022 20:59:12 +0200 Subject: [PATCH] MatrixFree: allow categorization of ghost cells --- source/matrix_free/task_info.cc | 103 +++++++-- tests/matrix_free/categorize_ghost_01.cc | 207 ++++++++++++++++++ ...gorize_ghost_01.with_p4est.mpirun=2.output | 9 + 3 files changed, 299 insertions(+), 20 deletions(-) create mode 100644 tests/matrix_free/categorize_ghost_01.cc create mode 100644 tests/matrix_free/categorize_ghost_01.with_p4est.mpirun=2.output diff --git a/source/matrix_free/task_info.cc b/source/matrix_free/task_info.cc index 30fb250604..04b46ea1fc 100644 --- a/source/matrix_free/task_info.cc +++ b/source/matrix_free/task_info.cc @@ -1014,10 +1014,62 @@ namespace internal blocks = {0, comm_begin, comm_end, end}; } - // Step 7: Fill in the data by batches for the locally owned cells. + // Step 7: sort ghost cells according to the category + std::vector> tight_category_map_ghost; + + if (cell_vectorization_categories.empty() == false) + { + tight_category_map_ghost.reserve(n_ghost_cells); + + std::set used_categories; + for (unsigned int i = 0; i < n_ghost_cells; ++i) + used_categories.insert( + cell_vectorization_categories[i + n_active_cells]); + + std::vector used_categories_vector( + used_categories.size()); + n_categories = 0; + for (const auto &it : used_categories) + used_categories_vector[n_categories++] = it; + + std::vector counters(n_categories, 0); + + for (unsigned int i = 0; i < n_ghost_cells; ++i) + { + const unsigned int index = + std::lower_bound( + used_categories_vector.begin(), + used_categories_vector.end(), + cell_vectorization_categories[i + n_active_cells]) - + used_categories_vector.begin(); + AssertIndexRange(index, used_categories_vector.size()); + tight_category_map_ghost.emplace_back( + std::array{{index, i}}); + + // account for padding in the hp and strict case + if (categories_are_hp || cell_vectorization_categories_strict) + counters[index]++; + } + + // insert padding + for (unsigned int i = 0; i < counters.size(); ++i) + if (counters[i] % n_lanes != 0) + for (unsigned int j = counters[i] % n_lanes; j < n_lanes; ++j) + tight_category_map_ghost.emplace_back( + std::array{ + {i, numbers::invalid_unsigned_int}}); + + std::sort(tight_category_map_ghost.begin(), + tight_category_map_ghost.end()); + } + + // Step 8: Fill in the data by batches for the locally owned cells. const unsigned int n_cell_batches = batch_order.size(); const unsigned int n_ghost_batches = - (n_ghost_cells + n_lanes - 1) / n_lanes; + ((tight_category_map_ghost.empty() ? n_ghost_cells : + tight_category_map_ghost.size()) + + n_lanes - 1) / + n_lanes; incompletely_filled_vectorization.resize(n_cell_batches + n_ghost_batches); @@ -1054,26 +1106,37 @@ namespace internal } AssertDimension(counter, n_active_cells); - // Step 8: Treat the ghost cells - for (unsigned int cell = n_active_cells; - cell < n_active_cells + n_ghost_cells; - ++cell) + // Step 9: Treat the ghost cells + if (tight_category_map_ghost.empty()) + { + for (unsigned int cell = 0; cell < n_ghost_cells; ++cell) + renumbering[n_active_cells + cell] = n_active_cells + cell; + + if ((n_ghost_cells % n_lanes) != 0u) + incompletely_filled_vectorization.back() = n_ghost_cells % n_lanes; + } + else { - if (!cell_vectorization_categories.empty()) - Assert( - cell_vectorization_categories[cell] == - cell_vectorization_categories[n_active_cells], - ExcMessage( - "Currently, all ghost cells need to have the same category, " - "but got " + - std::to_string(cell_vectorization_categories[cell]) + " and " + - std::to_string(cell_vectorization_categories[n_active_cells]) + - ". Please check MatrixFree::AdditionalData::cell_vectorization_category " - "or the active FE index of the ghost cells in the hp case!")); - renumbering[cell] = cell; + for (unsigned int k = 0, ptr = 0; k < n_ghost_batches; + ++k, ptr += n_lanes) + { + unsigned int j = 0; + + for (; + j < n_lanes && (ptr + j < tight_category_map_ghost.size()) && + (tight_category_map_ghost[ptr + j][1] != + numbers::invalid_unsigned_int); + ++j) + renumbering[counter++] = + n_active_cells + tight_category_map_ghost[ptr + j][1]; + + if (j < n_lanes) + incompletely_filled_vectorization[n_cell_batches + k] = j; + } + + AssertDimension(counter, n_active_cells + n_ghost_cells); } - if ((n_ghost_cells % n_lanes) != 0u) - incompletely_filled_vectorization.back() = n_ghost_cells % n_lanes; + cell_partition_data.push_back(n_cell_batches + n_ghost_batches); partition_row_index.back() = cell_partition_data.size() - 1; diff --git a/tests/matrix_free/categorize_ghost_01.cc b/tests/matrix_free/categorize_ghost_01.cc new file mode 100644 index 0000000000..2f41019ad7 --- /dev/null +++ b/tests/matrix_free/categorize_ghost_01.cc @@ -0,0 +1,207 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + + + +// tests categorization of ghost cells + +#include + +#include + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +test() +{ + using Number = double; + using VectorizedArrayType = VectorizedArray; + + parallel::shared::Triangulation tria( + MPI_COMM_WORLD, + ::Triangulation::none, + true, + parallel::shared::Triangulation::partition_custom_signal); + + tria.signals.create.connect([&]() { + for (const auto &cell : tria.active_cell_iterators()) + if (cell->center()[1] < 0.5) + cell->set_subdomain_id(0); + else + cell->set_subdomain_id(1); + }); + + GridGenerator::subdivided_hyper_rectangle(tria, + {10, 10}, + {0.0, 0.0}, + {1.0, 1.0}); + + // caterorization - not strict + { + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(FE_Q(1)); + + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_quadrature_points; + data.mapping_update_flags_boundary_faces = update_quadrature_points; + data.mapping_update_flags_inner_faces = update_quadrature_points; + data.hold_all_faces_to_owned_cells = true; + + std::vector cell_vectorization_category( + tria.n_active_cells()); + for (unsigned int i = 0; i < cell_vectorization_category.size(); ++i) + cell_vectorization_category[i] = i % 10; + + data.cell_vectorization_category = cell_vectorization_category; + data.cell_vectorization_categories_strict = false; + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handler, + AffineConstraints(), + QGauss(2), + data); + + for (unsigned int cell = matrix_free.n_cell_batches(); + cell < + matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches(); + ++cell) + { + unsigned int category = 0; + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + category = std::max( + category, + cell_vectorization_category[matrix_free.get_cell_iterator(cell, v) + ->active_cell_index()]); + + AssertDimension(category, matrix_free.get_cell_category(cell)); + } + } + deallog << "OK" << std::endl; + + // caterorization - strict + { + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(FE_Q(1)); + + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_quadrature_points; + data.mapping_update_flags_boundary_faces = update_quadrature_points; + data.mapping_update_flags_inner_faces = update_quadrature_points; + data.hold_all_faces_to_owned_cells = true; + + std::vector cell_vectorization_category( + tria.n_active_cells()); + for (unsigned int i = 0; i < cell_vectorization_category.size(); ++i) + cell_vectorization_category[i] = i % 10; + + data.cell_vectorization_category = cell_vectorization_category; + data.cell_vectorization_categories_strict = true; + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handler, + AffineConstraints(), + QGauss(2), + data); + + for (unsigned int cell = matrix_free.n_cell_batches(); + cell < + matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches(); + ++cell) + { + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + AssertDimension( + cell_vectorization_category[matrix_free.get_cell_iterator(cell, v) + ->active_cell_index()], + matrix_free.get_cell_category(cell)); + } + } + deallog << "OK" << std::endl; + + + // hp + { + DoFHandler dof_handler(tria); + + unsigned int counter = 0; + for (const auto &cell : dof_handler.active_cell_iterators()) + if (cell->is_locally_owned()) + cell->set_active_fe_index((counter++) % 10); + + hp::FECollection fe_collection; + + for (unsigned int i = 0; i < 10; ++i) + fe_collection.push_back(FE_Q(1)); + + dof_handler.distribute_dofs(fe_collection); + + typename MatrixFree::AdditionalData data; + data.tasks_parallel_scheme = MatrixFree::AdditionalData::none; + data.mapping_update_flags = update_quadrature_points; + data.mapping_update_flags_boundary_faces = update_quadrature_points; + data.mapping_update_flags_inner_faces = update_quadrature_points; + data.hold_all_faces_to_owned_cells = true; + + MatrixFree matrix_free; + matrix_free.reinit(MappingQ1(), + dof_handler, + AffineConstraints(), + QGauss(2), + data); + + for (unsigned int cell = matrix_free.n_cell_batches(); + cell < + matrix_free.n_cell_batches() + matrix_free.n_ghost_cell_batches(); + ++cell) + { + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + AssertDimension( + matrix_free.get_cell_iterator(cell, v)->active_fe_index(), + matrix_free.get_cell_category(cell)); + } + } + deallog << "OK" << std::endl; +} + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + + MPILogInitAll log; + + AssertDimension(Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD), 2); + + test<2>(); +} diff --git a/tests/matrix_free/categorize_ghost_01.with_p4est.mpirun=2.output b/tests/matrix_free/categorize_ghost_01.with_p4est.mpirun=2.output new file mode 100644 index 0000000000..5693ef953a --- /dev/null +++ b/tests/matrix_free/categorize_ghost_01.with_p4est.mpirun=2.output @@ -0,0 +1,9 @@ + +DEAL:0::OK +DEAL:0::OK +DEAL:0::OK + +DEAL:1::OK +DEAL:1::OK +DEAL:1::OK + -- 2.39.5