MatrixFree: allow categorization of ghost cells

diff --git a/source/matrix_free/task_info.cc b/source/matrix_free/task_info.cc
index 30fb2506049321761be2234a0d271ad1fe053fba..04b46ea1fc0bdf8f71a72f2fbfa6ba8c57dbd634 100644 (file)
--- 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<std::array<unsigned int, 2>> tight_category_map_ghost;
+
+      if (cell_vectorization_categories.empty() == false)
+        {
+          tight_category_map_ghost.reserve(n_ghost_cells);
+
+          std::set<unsigned int> used_categories;
+          for (unsigned int i = 0; i < n_ghost_cells; ++i)
+            used_categories.insert(
+              cell_vectorization_categories[i + n_active_cells]);
+
+          std::vector<unsigned int> used_categories_vector(
+            used_categories.size());
+          n_categories = 0;
+          for (const auto &it : used_categories)
+            used_categories_vector[n_categories++] = it;
+
+          std::vector<unsigned int> 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<unsigned int, 2>{{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<unsigned int, 2>{
+                    {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 (file)
index 0000000..2f41019
--- /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 <deal.II/distributed/shared_tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+test()
+{
+  using Number              = double;
+  using VectorizedArrayType = VectorizedArray<Number>;
+
+  parallel::shared::Triangulation<dim> tria(
+    MPI_COMM_WORLD,
+    ::Triangulation<dim>::none,
+    true,
+    parallel::shared::Triangulation<dim>::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<dim> dof_handler(tria);
+    dof_handler.distribute_dofs(FE_Q<dim>(1));
+
+    typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, double>::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<unsigned int> 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<dim, Number, VectorizedArrayType> matrix_free;
+    matrix_free.reinit(MappingQ1<dim>(),
+                       dof_handler,
+                       AffineConstraints<Number>(),
+                       QGauss<dim>(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<dim> dof_handler(tria);
+    dof_handler.distribute_dofs(FE_Q<dim>(1));
+
+    typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, double>::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<unsigned int> 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<dim, Number, VectorizedArrayType> matrix_free;
+    matrix_free.reinit(MappingQ1<dim>(),
+                       dof_handler,
+                       AffineConstraints<Number>(),
+                       QGauss<dim>(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<dim> 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<dim> fe_collection;
+
+    for (unsigned int i = 0; i < 10; ++i)
+      fe_collection.push_back(FE_Q<dim>(1));
+
+    dof_handler.distribute_dofs(fe_collection);
+
+    typename MatrixFree<dim, Number, VectorizedArrayType>::AdditionalData data;
+    data.tasks_parallel_scheme = MatrixFree<dim, double>::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<dim, Number, VectorizedArrayType> matrix_free;
+    matrix_free.reinit(MappingQ1<dim>(),
+                       dof_handler,
+                       AffineConstraints<Number>(),
+                       QGauss<dim>(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 (file)
index 0000000..5693ef9
--- /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
+