Introduction of transfer operator assembly for geometric multigrid

cpp/doc/source/api.rst

@@ -14,6 +14,7 @@ generated documentation is `here <doxygen>`_.
    io
    la
    mesh
+   multigrid
    refinement
 
 * :ref:`genindex`

cpp/doc/source/multigrid.rst

@@ -0,0 +1,5 @@
+Multigrid (``dolfinx::multigrid``)
+====================================
+
+.. doxygennamespace:: dolfinx::multigrid
+   :project: DOLFINx

cpp/dolfinx/CMakeLists.txt

@@ -17,6 +17,7 @@ set(DOLFINX_DIRS
     io
     la
     mesh
+    multigrid
     nls
     refinement
 )

cpp/dolfinx/graph/partitioners.cpp

@@ -35,21 +35,8 @@ extern "C"
 
 using namespace dolfinx;
 
-namespace
-{
-/// @todo Is it un-documented that the owning rank must come first in
-/// reach list of edges?
-///
-/// @param[in] comm The communicator
-/// @param[in] graph Graph, using global indices for graph edges
-/// @param[in] node_disp The distribution of graph nodes across MPI
-/// ranks. The global index `gidx` of local index `lidx` is `lidx +
-/// node_disp[my_rank]`.
-/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]`
-/// is the destination of the node with local index `i`.
-/// @return Destination ranks for each local node.
 template <typename T>
-graph::AdjacencyList<int> compute_destination_ranks(
+graph::AdjacencyList<int> dolfinx::graph::compute_destination_ranks(
     MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
     const std::vector<T>& node_disp, const std::vector<T>& part)
 {
@@ -202,6 +189,18 @@ graph::AdjacencyList<int> compute_destination_ranks(
 
   return g;
 }
+
+/// @cond
+template graph::AdjacencyList<int> dolfinx::graph::compute_destination_ranks(
+    MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
+    const std::vector<int>& node_disp, const std::vector<int>& part);
+
+template graph::AdjacencyList<int> dolfinx::graph::compute_destination_ranks(
+    MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
+    const std::vector<unsigned long long>& node_disp,
+    const std::vector<unsigned long long>& part);
+/// @endcond
+
 //-----------------------------------------------------------------------------
 #ifdef HAS_PARMETIS
 template <typename T>
@@ -304,7 +303,6 @@ std::vector<int> refine(MPI_Comm comm, const graph::AdjacencyList<T>& adj_graph)
   //-----------------------------------------------------------------------------
 }
 #endif
-} // namespace
 
 //-----------------------------------------------------------------------------
 #ifdef HAS_PTSCOTCH
@@ -587,7 +585,7 @@ graph::partition_fn graph::parmetis::partitioner(double imbalance,
     {
       // FIXME: Is it implicit that the first entry is the owner?
       graph::AdjacencyList<int> dest
-          = compute_destination_ranks(pcomm, graph, node_disp, part);
+          = graph::compute_destination_ranks(pcomm, graph, node_disp, part);
       if (split_comm)
         MPI_Comm_free(&pcomm);
       return dest;
@@ -653,7 +651,7 @@ graph::partition_fn graph::kahip::partitioner(int mode, int seed,
     timer2.stop();
 
     if (ghosting)
-      return compute_destination_ranks(comm, graph, node_disp, part);
+      return graph::compute_destination_ranks(comm, graph, node_disp, part);
     else
     {
       return regular_adjacency_list(std::vector<int>(part.begin(), part.end()),

cpp/dolfinx/graph/partitioners.h

@@ -11,6 +11,23 @@
 
 namespace dolfinx::graph
 {
+
+/// @todo Is it un-documented that the owning rank must come first in
+/// reach list of edges?
+///
+/// @param[in] comm The communicator
+/// @param[in] graph Graph, using global indices for graph edges
+/// @param[in] node_disp The distribution of graph nodes across MPI
+/// ranks. The global index `gidx` of local index `lidx` is `lidx +
+/// node_disp[my_rank]`.
+/// @param[in] part The destination rank for owned nodes, i.e. `dest[i]`
+/// is the destination of the node with local index `i`.
+/// @return Destination ranks for each local node.
+template <typename T>
+graph::AdjacencyList<int> compute_destination_ranks(
+    MPI_Comm comm, const graph::AdjacencyList<std::int64_t>& graph,
+    const std::vector<T>& node_disp, const std::vector<T>& part);
+
 namespace scotch
 {
 #ifdef HAS_PTSCOTCH

cpp/dolfinx/multigrid/CMakeLists.txt

@@ -0,0 +1,5 @@
+set(HEADERS_multigrid
+    ${CMAKE_CURRENT_SOURCE_DIR}/inclusion.h
+    ${CMAKE_CURRENT_SOURCE_DIR}/transfer_matrix.h
+    PARENT_SCOPE
+)

cpp/dolfinx/multigrid/dolfinx_multigrid.h

@@ -0,0 +1,12 @@
+#pragma once
+
+/// @brief Multigrid algorithms.
+///
+namespace dolfinx::multigrid
+{
+}
+
+// DOLFINx multigrid interface
+
+#include <dolfinx/multigrid/inclusion.h>
+#include <dolfinx/multigrid/transfer_matrix.h>

cpp/dolfinx/multigrid/inclusion.h

@@ -0,0 +1,186 @@
+// Copyright (C) 2024 Paul T. Kühner
+//
+// This file is part of DOLFINX (https://www.fenicsproject.org)
+//
+// SPDX-License-Identifier:    LGPL-3.0-or-later
+
+#pragma once
+
+#include <algorithm>
+#include <concepts>
+#include <cstdint>
+#include <iterator>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+
+#include <mpi.h>
+
+#include "dolfinx/common/IndexMap.h"
+#include "dolfinx/la/SparsityPattern.h"
+#include "dolfinx/mesh/Mesh.h"
+
+namespace dolfinx::multigrid
+{
+
+/**
+ * @brief Gathers a global vector from combination of local data.
+ * @note Performs an all-to-all communication.
+ *
+ * @param local local data
+ * @param global_size number of global data entries
+ * @param comm MPI communicator
+ * @return std::vector<T> on communicator gathered global data.
+ */
+template <std::floating_point T>
+std::vector<T> gather_global(std::span<const T> local, std::int64_t global_size,
+                             MPI_Comm comm)
+{
+  // 1) exchange local sizes
+  std::vector<std::int32_t> local_sizes(dolfinx::MPI::size(comm));
+  {
+    std::array<std::int32_t, 1> tmp{local.size()};
+    MPI_Allgather(&tmp, 1, MPI_INT32_T, local_sizes.data(), 1, MPI_INT32_T,
+                  comm);
+  }
+
+  // 2) compute displacement vector
+  std::vector<std::int32_t> displacements(local_sizes.size() + 1, 0);
+  std::partial_sum(local_sizes.begin(), local_sizes.end(),
+                   displacements.begin() + 1);
+
+  // 3) Allgather global vector
+  std::vector<T> global(global_size);
+  MPI_Allgatherv(local.data(), local.size(), dolfinx::MPI::mpi_t<T>,
+                 global.data(), local_sizes.data(), displacements.data(),
+                 dolfinx::MPI::mpi_t<T>, comm);
+
+  return global;
+}
+
+/**
+ * @brief Computes an inclusion map, i.e. local list of global vertex indices of
+ * another mesh, between to meshes.
+ *
+ *
+ * @param mesh_from Coarser mesh (domain of the inclusion map)
+ * @param mesh_to Finer mesh (range of the inclusion map)
+ * @param allow_all_to_all If the vertices of `mesh_from` are not equally
+ * spatially parallelized as `mesh_to` an all-to-all gathering of all vertices
+ * in `mesh_to` is performed. If true, performs all-to-all gathering, otherwise
+ * throws an exception if this becomes necessary.
+ * @return std::vector<std::int64_t> Map from local vertex index in `mesh_from`
+ * to global vertex index in `mesh_to`, i.e. `mesh_from.geometry.x()[i:i+3] ==
+ * mesh_to.geometry.x()[map[i]:map[i]+3]`.
+ */
+template <std::floating_point T>
+std::vector<std::int64_t>
+inclusion_mapping(const dolfinx::mesh::Mesh<T>& mesh_from,
+                  const dolfinx::mesh::Mesh<T>& mesh_to,
+                  bool allow_all_to_all = false)
+{
+  {
+    // Check comms equal
+    int result;
+    MPI_Comm_compare(mesh_from.comm(), mesh_to.comm(), &result);
+    assert(result == MPI_CONGRUENT);
+  }
+
+  const common::IndexMap& im_from = *mesh_from.topology()->index_map(0);
+  const common::IndexMap& im_to = *mesh_to.topology()->index_map(0);
+
+  std::vector<std::int64_t> map(im_from.size_local() + im_from.num_ghosts(),
+                                -1);
+
+  std::span<const T> x_from = mesh_from.geometry().x();
+  std::span<const T> x_to = mesh_to.geometry().x();
+
+  auto build_global_to_local = [&](const auto& im)
+  {
+    return [&](std::int32_t idx)
+    {
+      std::array<std::int64_t, 1> tmp;
+      im.local_to_global(std::vector<std::int32_t>{idx}, tmp);
+      return tmp[0];
+    };
+  };
+
+  auto to_global_to = build_global_to_local(im_to);
+  auto to_global_from = build_global_to_local(im_from);
+
+  for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++)
+  {
+    std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i),
+                                      std::next(x_from.begin(), 3 * (i + 1)));
+    for (std::int64_t j = 0; j < im_to.size_local() + im_to.num_ghosts(); j++)
+    {
+      std::ranges::subrange vertex_to(std::next(x_to.begin(), 3 * j),
+                                      std::next(x_to.begin(), 3 * (j + 1)));
+
+      if (std::ranges::equal(
+              vertex_from, vertex_to, [](T a, T b)
+              { return std::abs(a - b) <= std::numeric_limits<T>::epsilon(); }))
+      {
+        assert(map[i] == -1);
+        map[i] = to_global_to(j);
+        break;
+      }
+    }
+  }
+
+  if (dolfinx::MPI::size(mesh_to.comm()) == 1)
+  {
+    // no communication required
+    assert(std::ranges::all_of(map, [](auto e) { return e >= 0; }));
+    return map;
+  }
+
+  bool all_found = std::ranges::all_of(map, [](auto e) { return e >= 0; });
+  MPI_Allreduce(MPI_IN_PLACE, &all_found, 1, MPI_CXX_BOOL, MPI_LAND,
+                mesh_to.comm());
+  if (all_found)
+    return map;
+
+  if (!allow_all_to_all)
+  {
+    throw std::runtime_error(
+        "Parallelization of mesh requires all to all communication to compute "
+        "inclusion map, but allow_all_to_all is not set.");
+  }
+
+  // Build global to vertex list
+  std::vector<T> global_x_to
+      = gather_global(mesh_to.geometry().x().subspan(0, im_to.size_local() * 3),
+                      im_to.size_global() * 3, mesh_to.comm());
+
+  // Recheck indices on global data structure
+  for (std::int32_t i = 0; i < im_from.size_local() + im_from.num_ghosts(); i++)
+  {
+    // TODO:
+    // if (map[i] >= 0)
+    //   continue;
+
+    std::ranges::subrange vertex_from(std::next(x_from.begin(), 3 * i),
+                                      std::next(x_from.begin(), 3 * (i + 1)));
+    for (std::int64_t j = 0; j < im_to.size_global(); j++)
+    {
+      std::ranges::subrange vertex_to(
+          std::next(global_x_to.begin(), 3 * j),
+          std::next(global_x_to.begin(), 3 * (j + 1)));
+
+      if (std::ranges::equal(
+              vertex_from, vertex_to, [](T a, T b)
+              { return std::abs(a - b) <= std::numeric_limits<T>::epsilon(); }))
+      {
+        map[i] = j;
+        break;
+      }
+    }
+  }
+
+  assert(std::ranges::all_of(map, [&](auto e)
+                             { return e >= 0 && e < im_to.size_global(); }));
+  return map;
+}
+
+} // namespace dolfinx::multigrid