From: Peter Munch Date: Sun, 19 Jun 2022 08:23:31 +0000 (+0200) Subject: Make get_cell_range_category() more consistent with get_cell_category() X-Git-Tag: v9.4.0-rc2~1^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14028%2Fhead;p=dealii.git Make get_cell_range_category() more consistent with get_cell_category() --- diff --git a/include/deal.II/matrix_free/matrix_free.h b/include/deal.II/matrix_free/matrix_free.h index 68fb429be8..81cc9506a7 100644 --- a/include/deal.II/matrix_free/matrix_free.h +++ b/include/deal.II/matrix_free/matrix_free.h @@ -1887,8 +1887,14 @@ public: /** * Return the category the current batch range of cells was assigned to. * Categories run between the given values in the field - * AdditionalData::cell_vectorization_category for non-hp-DoFHandler types + * AdditionalData::cell_vectorization_category for the non-hp case * and return the active FE index in the hp-adaptive case. + * + * @note Following the behaviour of get_cell_category(), we return the + * maximum category of any cell batch. In the hp case, it is + * guaranteed that all cells and as a consequence all cell batches in a range + * have the same category. Otherwise, there may be different categories in + * different cell batches. */ unsigned int get_cell_range_category( @@ -1905,8 +1911,13 @@ public: /** * Return the category the current batch of cells was assigned to. Categories * run between the given values in the field - * AdditionalData::cell_vectorization_category for non-hp-DoFHandler types + * AdditionalData::cell_vectorization_category for the non-hp case * and return the active FE index in the hp-adaptive case. + * + * @note In the non-hp case, a category of a cell batch is given + * as the maximum category of any of its cell. In the hp case or the case that + * MatrixFree::AdditionalData::cell_vectorization_categories_strict was + * enabled, it is guaranteed that all cells have the same category. */ unsigned int get_cell_category(const unsigned int cell_batch_index) const; @@ -2800,14 +2811,10 @@ inline unsigned int MatrixFree::get_cell_range_category( const std::pair range) const { - const auto result = get_cell_category(range.first); + auto result = get_cell_category(range.first); -# ifdef DEBUG for (unsigned int i = range.first; i < range.second; ++i) - Assert(result == get_cell_category(i), - ExcMessage( - "The cell batches of the range have different categories!")); -# endif + result = std::max(result, get_cell_category(i)); return result; } @@ -2819,14 +2826,13 @@ inline std::pair MatrixFree::get_face_range_category( const std::pair range) const { - const auto result = get_face_category(range.first); + auto result = get_face_category(range.first); -# ifdef DEBUG for (unsigned int i = range.first; i < range.second; ++i) - Assert(result == get_face_category(i), - ExcMessage( - "The face batches of the range have different categories!")); -# endif + { + result.first = std::max(result.first, get_face_category(i).first); + result.second = std::max(result.second, get_face_category(i).second); + } return result; } diff --git a/tests/matrix_free/get_cell_catergory_01.cc b/tests/matrix_free/get_cell_catergory_01.cc new file mode 100644 index 0000000000..fa82e9a058 --- /dev/null +++ b/tests/matrix_free/get_cell_catergory_01.cc @@ -0,0 +1,246 @@ +// --------------------------------------------------------------------- +// +// 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 MatrixFree::get_cell_category() and +// MatrixFree::get_cell_range_category() + +#include + +#include + +#include +#include + +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +test() +{ + using Number = double; + using VectorizedArrayType = VectorizedArray; + + Triangulation tria; + + GridGenerator::hyper_cube(tria); + tria.refine_global(3); + + // 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; + + 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 / 7; + + 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); + + double dummy = 0; + + matrix_free.template cell_loop( + [&](const auto &matrix_free, auto &, const auto &, const auto range) { + unsigned int max_category_range = 0; + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + unsigned int max_category = 0; + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + { + const auto category = cell_vectorization_category + [matrix_free.get_cell_iterator(cell, v)->active_cell_index()]; + + // note: here and below we are using std::cout to output + // information that is useful for debugging; however, since the + // output depends on the hardware and might change if we + // internally modify the way we loop over cells, we don't write + // the data to the log but instead only check that the + // properties that are described in the documentation are + // fulfilled + std::cout << category << " "; + max_category = std::max(max_category, category); + } + + std::cout << std::endl; + + AssertDimension(max_category, + matrix_free.get_cell_category(cell)); // test + + max_category_range = std::max(max_category_range, max_category); + } + + AssertDimension(max_category_range, + matrix_free.get_cell_range_category(range)); // test + std::cout << std::endl; + }, + dummy, + dummy); + } + 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 / 7; + + 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); + + double dummy = 0; + + matrix_free.template cell_loop( + [&](const auto &matrix_free, auto &, const auto &, const auto range) { + unsigned int max_category_range = 0; + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + { + const auto category = cell_vectorization_category + [matrix_free.get_cell_iterator(cell, v)->active_cell_index()]; + std::cout << category << " "; + AssertDimension(category, + matrix_free.get_cell_category(cell)); // test + } + + std::cout << std::endl; + + max_category_range = + std::max(max_category_range, matrix_free.get_cell_category(cell)); + } + + AssertDimension(max_category_range, + matrix_free.get_cell_range_category(range)); // test + std::cout << std::endl; + }, + dummy, + dummy); + } + 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++) / 7); + + hp::FECollection fe_collection; + + for (unsigned int i = 0; i < (tria.n_cells() + 6) / 7; ++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); + + + double dummy = 0; + + matrix_free.template cell_loop( + [&](const auto &matrix_free, auto &, const auto &, const auto range) { + unsigned int max_category_range = 0; + + for (unsigned int cell = range.first; cell < range.second; ++cell) + { + for (unsigned int v = 0; + v < matrix_free.n_active_entries_per_cell_batch(cell); + ++v) + { + const auto category = + matrix_free.get_cell_iterator(cell, v)->active_fe_index(); + std::cout << category << " "; + AssertDimension(category, + matrix_free.get_cell_category(cell)); // test + AssertDimension( + category, matrix_free.get_cell_range_category(range)); // test + } + + std::cout << std::endl; + } + std::cout << std::endl; + }, + dummy, + dummy); + } + deallog << "OK" << std::endl; +} + +int +main(int argc, char **argv) +{ + initlog(); + + test<2>(); +} diff --git a/tests/matrix_free/get_cell_catergory_01.output b/tests/matrix_free/get_cell_catergory_01.output new file mode 100644 index 0000000000..fb71de2867 --- /dev/null +++ b/tests/matrix_free/get_cell_catergory_01.output @@ -0,0 +1,4 @@ + +DEAL::OK +DEAL::OK +DEAL::OK