Topology build speed-ups and (including threading)

cpp/dolfinx/io/VTKHDF.h

@@ -468,6 +468,6 @@ mesh::Mesh<U> read_mesh(MPI_Comm comm, const std::filesystem::path& filename,
                                                         cells_local.end());
   return mesh::create_mesh(comm, comm, cells_span, coordinate_elements, comm,
                            points_pruned, {(std::size_t)x_shape[0], gdim}, part,
-                           max_facet_to_cell_links);
+                           max_facet_to_cell_links, 1);
 }
 } // namespace dolfinx::io::VTKHDF

cpp/dolfinx/mesh/Topology.cpp

@@ -10,7 +10,10 @@
 #include "topologycomputation.h"
 #include "utils.h"
 #include <algorithm>
+#include <array>
+#include <boost/sort/sort.hpp>
 #include <dolfinx/common/IndexMap.h>
+#include <dolfinx/common/local_range.h>
 #include <dolfinx/common/log.h>
 #include <dolfinx/common/sort.h>
 #include <dolfinx/graph/AdjacencyList.h>
@@ -19,6 +22,7 @@
 #include <numeric>
 #include <random>
 #include <set>
+#include <thread>
 
 using namespace dolfinx;
 using namespace dolfinx::mesh;
@@ -40,7 +44,8 @@
 /// @return Map from global index to sharing ranks for each index in
 /// indices. The owner rank is the first as the first in the of ranks.
 graph::AdjacencyList<int>
-determine_sharing_ranks(MPI_Comm comm, std::span<const std::int64_t> indices)
+determine_sharing_ranks(MPI_Comm comm, std::span<const std::int64_t> indices,
+                        int num_threads)
 {
   common::Timer timer("Topology: determine shared index ownership");
 
@@ -63,7 +68,13 @@
       int dest = dolfinx::MPI::index_owner(size, idx, global_range);
       dest_to_index.push_back({dest, static_cast<int>(dest_to_index.size())});
     }
-    std::ranges::sort(dest_to_index);
+    if (num_threads > 1)
+    {
+      boost::sort::block_indirect_sort(dest_to_index.begin(),
+                                       dest_to_index.end(), num_threads);
+    }
+    else
+      std::ranges::sort(dest_to_index);
   }
 
   // Build list of neighbour dest ranks and count number of indices to
@@ -135,10 +146,13 @@
 
   // Build {global index, pos, src} list
   std::vector<std::array<std::int64_t, 3>> indices_list;
-  for (std::size_t p = 0; p < recv_disp0.size() - 1; ++p)
-    for (std::int32_t i = recv_disp0[p]; i < recv_disp0[p + 1]; ++i)
-      indices_list.push_back({recv_buffer0[i], i, int(p)});
-  std::ranges::sort(indices_list);
+  {
+    common::Timer timer("Topology: build and sort transposed index list");
+    for (std::size_t p = 0; p < recv_disp0.size() - 1; ++p)
+      for (std::int32_t i = recv_disp0[p]; i < recv_disp0[p + 1]; ++i)
+        indices_list.push_back({recv_buffer0[i], i, int(p)});
+    std::ranges::sort(indices_list);
+  }
 
   // Find which ranks have each index
   std::vector<std::int32_t> num_items_per_dest1(recv_disp0.size() - 1, 0);
@@ -286,7 +300,7 @@
 std::array<std::vector<std::int64_t>, 2> vertex_ownership_groups(
     const std::vector<std::span<const std::int64_t>>& cells_owned,
     const std::vector<std::span<const std::int64_t>>& cells_ghost,
-    std::span<const std::int64_t> boundary_vertices)
+    std::span<const std::int64_t> boundary_vertices, int num_threads)
 {
   common::Timer timer("Topology: determine vertex ownership groups (owned, "
                       "undetermined, unowned)");
@@ -300,10 +314,18 @@
     local_vertex_set.insert(local_vertex_set.end(), c.begin(), c.end());
 
   {
-    dolfinx::radix_sort(local_vertex_set);
+    if (num_threads > 1)
+    {
+      boost::sort::block_indirect_sort(local_vertex_set.begin(),
+                                       local_vertex_set.end(), num_threads);
+    }
+    else
+      dolfinx::radix_sort(local_vertex_set);
+
     auto [unique_end, range_end] = std::ranges::unique(local_vertex_set);
     local_vertex_set.erase(unique_end, range_end);
   }
+
   // Build set of ghost cell vertices (attached to a ghost cell)
   std::vector<std::int64_t> ghost_vertex_set;
   ghost_vertex_set.reserve(
@@ -313,10 +335,18 @@
     ghost_vertex_set.insert(ghost_vertex_set.end(), c.begin(), c.end());
 
   {
-    dolfinx::radix_sort(ghost_vertex_set);
+    if (num_threads > 1)
+    {
+      boost::sort::block_indirect_sort(ghost_vertex_set.begin(),
+                                       ghost_vertex_set.end(), num_threads);
+    }
+    else
+      dolfinx::radix_sort(ghost_vertex_set);
+
     auto [unique_end, range_end] = std::ranges::unique(ghost_vertex_set);
     ghost_vertex_set.erase(unique_end, range_end);
   }
+
   // Build difference 1: Vertices attached only to owned cells, and
   // therefore owned by this rank
   std::vector<std::int64_t> owned_vertices;
@@ -329,9 +359,9 @@
   std::ranges::set_difference(ghost_vertex_set, local_vertex_set,
                               std::back_inserter(unowned_vertices));
 
-  // TODO Check this in debug mode only?
-  // Sanity check
-  // No vertices in unowned should also be in boundary...
+#ifndef NDEBUG
+  // Sanity check: no vertices in unowned should also be in boundary.
+  // Test in DEBUG mode only because of cost.
   std::vector<std::int64_t> unowned_vertices_in_error;
   std::ranges::set_intersection(unowned_vertices, boundary_vertices,
                                 std::back_inserter(unowned_vertices_in_error));
@@ -341,6 +371,7 @@
     throw std::runtime_error(
         "Adding boundary vertices in ghost cells not allowed.");
   }
+#endif
 
   return {std::move(owned_vertices), std::move(unowned_vertices)};
 }
@@ -1093,17 +1124,17 @@
 
   // Create sets of owned and unowned vertices from the cell ownership
   // and the list of boundary vertices
-  auto [owned_vertices, unowned_vertices]
-      = vertex_ownership_groups(owned_cells, ghost_cells, boundary_vertices);
+  auto [owned_vertices, unowned_vertices] = vertex_ownership_groups(
+      owned_cells, ghost_cells, boundary_vertices, num_threads);
+
+  timer1.stop();
+  timer1.flush();
 
   // For each vertex whose ownership needs determining, find the sharing
   // ranks. The first index in the list of ranks for a vertex is the
   // owner (as determined by determine_sharing_ranks).
   const graph::AdjacencyList<int> global_vertex_to_ranks
-      = determine_sharing_ranks(comm, boundary_vertices);
-
-  timer1.stop();
-  timer1.flush();
+      = determine_sharing_ranks(comm, boundary_vertices, num_threads);
 
   // Iterate over vertices that have 'unknown' ownership, and if flagged
   // as owned by determine_sharing_ranks update ownership status
@@ -1138,7 +1169,7 @@
   // Number all owned vertices, iterating over vertices cell-wise
   std::vector<std::int32_t> local_vertex_indices(owned_vertices.size(), -1);
   {
-    common::Timer timer3("Topology: 3 ");
+    common::Timer timer3("Topology: number owned vertices");
     std::int32_t v = 0;
     for (std::span<const std::int64_t> cells_t : cells)
     {

cpp/dolfinx/mesh/generation.h

@@ -285,20 +285,20 @@ create_interval(MPI_Comm comm, std::int64_t n, std::array<T, 2> p,
     if (gdim == 1)
     {
       return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF,
-                         x1d, {npts, 1}, partitioner, 2, reorder_fn);
+                         x1d, {npts, 1}, partitioner, 2, 1, reorder_fn);
     }
     std::vector<T> x(npts * gdim, T(0));
     for (std::size_t i = 0; i < npts; i++)
       x[i * gdim] = x1d[i];
     return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x,
                        {npts, static_cast<std::size_t>(gdim)}, partitioner, 2,
-                       reorder_fn);
+                       1, reorder_fn);
   }
   else
   {
     return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL,
                        std::vector<T>{}, {0, static_cast<std::size_t>(gdim)},
-                       partitioner, 2, reorder_fn);
+                       partitioner, 2, 1, reorder_fn);
   }
 }
 
@@ -422,7 +422,7 @@ Mesh<T> build_tet(MPI_Comm comm, MPI_Comm subcomm,
   }
 
   return create_mesh(comm, subcomm, cells, element, subcomm, x,
-                     {x.size() / 3, 3}, partitioner, 2, reorder_fn);
+                     {x.size() / 3, 3}, partitioner, 2, 1, reorder_fn);
 }
 
 template <std::floating_point T>
@@ -467,7 +467,7 @@ build_hex(MPI_Comm comm, MPI_Comm subcomm, std::array<std::array<T, 3>, 2> p,
   }
 
   return create_mesh(comm, subcomm, cells, element, subcomm, x,
-                     {x.size() / 3, 3}, partitioner, 2, reorder_fn);
+                     {x.size() / 3, 3}, partitioner, 2, 1, reorder_fn);
 }
 
 template <std::floating_point T>
@@ -515,7 +515,7 @@ Mesh<T> build_prism(MPI_Comm comm, MPI_Comm subcomm,
   }
 
   return create_mesh(comm, subcomm, cells, element, subcomm, x,
-                     {x.size() / 3, 3}, partitioner, 2, reorder_fn);
+                     {x.size() / 3, 3}, partitioner, 2, 1, reorder_fn);
 }
 
 template <std::floating_point T>
@@ -672,7 +672,7 @@ Mesh<T> build_tri(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
     if (gdim == 2)
     {
       return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x,
-                         {npts, 2}, partitioner, 2, reorder_fn);
+                         {npts, 2}, partitioner, 2, 1, reorder_fn);
     }
     std::vector<T> xg(npts * gdim, T(0));
     for (std::size_t i = 0; i < npts; i++)
@@ -682,13 +682,13 @@ Mesh<T> build_tri(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
     }
     return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, xg,
                        {npts, static_cast<std::size_t>(gdim)}, partitioner, 2,
-                       reorder_fn);
+                       1, reorder_fn);
   }
   else
   {
     return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL,
                        std::vector<T>{}, {0, static_cast<std::size_t>(gdim)},
-                       partitioner, 2, reorder_fn);
+                       partitioner, 2, 1, reorder_fn);
   }
 }
 
@@ -738,7 +738,7 @@ Mesh<T> build_quad(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
     if (gdim == 2)
     {
       return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, x,
-                         {npts, 2}, partitioner, 2, reorder_fn);
+                         {npts, 2}, partitioner, 2, 1, reorder_fn);
     }
     std::vector<T> xg(npts * gdim, T(0));
     for (std::size_t i = 0; i < npts; i++)
@@ -748,13 +748,13 @@ Mesh<T> build_quad(MPI_Comm comm, std::array<std::array<T, 2>, 2> p,
     }
     return create_mesh(comm, MPI_COMM_SELF, cells, element, MPI_COMM_SELF, xg,
                        {npts, static_cast<std::size_t>(gdim)}, partitioner, 2,
-                       reorder_fn);
+                       1, reorder_fn);
   }
   else
   {
     return create_mesh(comm, MPI_COMM_NULL, {}, element, MPI_COMM_NULL,
                        std::vector<T>{}, {0, static_cast<std::size_t>(gdim)},
-                       partitioner, 2, reorder_fn);
+                       partitioner, 2, 1, reorder_fn);
   }
 }
 } // namespace impl