]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement AffineConstraints::make_consistent_in_parallel() 12430/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 9 Jun 2021 23:53:53 +0000 (01:53 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 17 Jun 2021 13:21:04 +0000 (15:21 +0200)
include/deal.II/base/mpi_tags.h
include/deal.II/lac/affine_constraints.h
include/deal.II/lac/affine_constraints.templates.h
tests/lac/constraints_make_consistent_in_parallel_01.cc [new file with mode: 0644]
tests/lac/constraints_make_consistent_in_parallel_01.with_p4est=true.mpirun=4.output [new file with mode: 0644]

index a59594d59c2e0639df944b78a1ea7452f3d9e82f..2aca35e75becd93be30895086d60e122721f9b59 100644 (file)
@@ -149,6 +149,10 @@ namespace Utilities
           // GridTools::internal::distributed_compute_point_locations
           distributed_compute_point_locations,
 
+          // AffineConstraints::make_consistent_in_parallel()
+          affine_constraints_make_consistent_in_parallel_0,
+          affine_constraints_make_consistent_in_parallel_1,
+
         };
       } // namespace Tags
     }   // namespace internal
index d8025e47d014733f3d95da366c068ee07c5caaa7..f324dd3930378866ae05dd34e70ab3ed412b983a 100644 (file)
@@ -1711,6 +1711,16 @@ public:
                             const MPI_Comm &             mpi_communicator,
                             const bool                   verbose = false) const;
 
+  /**
+   * Make the current object consistent on all processors
+   * in a distributed computation. One should call this function before
+   * calling close().
+   */
+  void
+  make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
+                              const IndexSet &locally_relevant_dofs,
+                              const MPI_Comm  mpi_communicator);
+
   /**
    * Exception
    *
index 5d30fd2e71ff672ac6165206000cc531965b755a..4b113e64376c14a0da9fb01313b42e0f27a3ecc5 100644 (file)
@@ -20,6 +20,7 @@
 
 #include <deal.II/base/cuda_size.h>
 #include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/mpi_compute_index_owner_internal.h>
 #include <deal.II/base/table.h>
 #include <deal.II/base/thread_local_storage.h>
 
@@ -190,6 +191,361 @@ AffineConstraints<number>::is_consistent_in_parallel(
 }
 
 
+namespace internal
+{
+  template <typename number>
+  std::vector<
+    std::tuple<types::global_dof_index,
+               number,
+               std::vector<std::pair<types::global_dof_index, number>>>>
+  compute_locally_relevant_constraints(
+    const dealii::AffineConstraints<number> &constraints_in,
+    const IndexSet &                         locally_owned_dofs,
+    const IndexSet &                         locally_relevant_dofs,
+    const MPI_Comm                           mpi_communicator)
+  {
+    using ConstraintType =
+      std::tuple<types::global_dof_index,
+                 number,
+                 std::vector<std::pair<types::global_dof_index, number>>>;
+
+    // The result vector filled step by step.
+    std::vector<ConstraintType> locally_relevant_constraints;
+
+#ifndef DEAL_II_WITH_MPI
+    AssertThrow(false, ExcNotImplemented()); // one should not come here
+    (void)constraints_in;
+    (void)locally_owned_dofs;
+    (void)locally_relevant_dofs;
+    (void)mpi_communicator;
+#else
+
+    const unsigned int my_rank =
+      Utilities::MPI::this_mpi_process(mpi_communicator);
+
+    // helper function
+    const auto sort_constraints = [&]() {
+      std::sort(locally_relevant_constraints.begin(),
+                locally_relevant_constraints.end(),
+                [](const auto &a, const auto &b) {
+                  return std::get<0>(a) < std::get<0>(b);
+                });
+
+      locally_relevant_constraints.erase(
+        std::unique(locally_relevant_constraints.begin(),
+                    locally_relevant_constraints.end(),
+                    [](const auto &a, const auto &b) {
+                      return std::get<0>(a) == std::get<0>(b);
+                    }),
+        locally_relevant_constraints.end());
+    };
+
+    // 0) collect constrained indices of the current object
+    IndexSet constrained_indices(locally_owned_dofs.size());
+
+    std::vector<types::global_dof_index> constrained_indices_temp;
+    for (const auto &line : constraints_in.get_lines())
+      constrained_indices_temp.push_back(line.index);
+
+    constrained_indices.add_indices(constrained_indices_temp.begin(),
+                                    constrained_indices_temp.end());
+
+    // step 1: identify owners of constraints
+    std::vector<unsigned int> constrained_indices_owners(
+      constrained_indices.n_elements());
+    Utilities::MPI::internal::ComputeIndexOwner::ConsensusAlgorithmsPayload
+      constrained_indices_process(locally_owned_dofs,
+                                  constrained_indices,
+                                  mpi_communicator,
+                                  constrained_indices_owners,
+                                  true);
+
+    Utilities::MPI::ConsensusAlgorithms::Selector<
+      std::pair<types::global_dof_index, types::global_dof_index>,
+      unsigned int>(constrained_indices_process, mpi_communicator)
+      .run();
+
+    // step 2: collect all locally owned constraints
+    const auto constrained_indices_by_ranks =
+      constrained_indices_process.get_requesters();
+
+    {
+      const unsigned int tag = Utilities::MPI::internal::Tags::
+        affine_constraints_make_consistent_in_parallel_0;
+
+      // ... collect data and sort according to owner
+      std::map<unsigned int, std::vector<ConstraintType>> send_data_temp;
+
+      for (unsigned int i = 0; i < constrained_indices_owners.size(); ++i)
+        {
+          ConstraintType entry;
+
+          const types::global_dof_index index =
+            constrained_indices.nth_index_in_set(i);
+
+          std::get<0>(entry) = index;
+
+          if (constraints_in.is_inhomogeneously_constrained(index))
+            std::get<1>(entry) = constraints_in.get_inhomogeneity(index);
+
+          const auto constraints = constraints_in.get_constraint_entries(index);
+
+          if (constraints)
+            for (const auto &i : *constraints)
+              std::get<2>(entry).push_back(i);
+
+          if (constrained_indices_owners[i] == my_rank)
+            locally_relevant_constraints.push_back(entry);
+          else
+            send_data_temp[constrained_indices_owners[i]].push_back(entry);
+        }
+
+      std::map<unsigned int, std::vector<char>> send_data;
+
+      for (const auto &i : send_data_temp)
+        send_data[i.first] = Utilities::pack(i.second, false);
+
+      std::vector<MPI_Request> requests;
+      requests.reserve(send_data.size());
+
+      // ... send data
+      for (const auto &i : send_data)
+        {
+          if (i.first == my_rank)
+            continue;
+
+          requests.push_back({});
+
+          const int ierr = MPI_Isend(i.second.data(),
+                                     i.second.size(),
+                                     MPI_CHAR,
+                                     i.first,
+                                     tag,
+                                     mpi_communicator,
+                                     &requests.back());
+          AssertThrowMPI(ierr);
+        }
+
+      // ... receive data
+      unsigned int n_rec_ranks = 0;
+
+      for (const auto &i : constrained_indices_by_ranks)
+        if (i.first != my_rank)
+          n_rec_ranks++;
+
+      for (unsigned int i = 0; i < n_rec_ranks; ++i)
+        {
+          MPI_Status status;
+          int ierr = MPI_Probe(MPI_ANY_SOURCE, tag, mpi_communicator, &status);
+          AssertThrowMPI(ierr);
+
+          int message_length;
+          ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
+          AssertThrowMPI(ierr);
+
+          std::vector<char> buffer(message_length);
+
+          ierr = MPI_Recv(buffer.data(),
+                          buffer.size(),
+                          MPI_CHAR,
+                          status.MPI_SOURCE,
+                          tag,
+                          mpi_communicator,
+                          MPI_STATUS_IGNORE);
+          AssertThrowMPI(ierr);
+
+          const auto data =
+            Utilities::unpack<std::vector<ConstraintType>>(buffer, false);
+
+          for (const auto &i : data)
+            locally_relevant_constraints.push_back(i);
+        }
+
+      const int ierr =
+        MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+      AssertThrowMPI(ierr);
+
+      sort_constraints();
+    }
+
+    // step 3: communicate constraints so that each process know how the
+    // locally locally relevant dofs are constrained
+    {
+      const unsigned int tag = Utilities::MPI::internal::Tags::
+        affine_constraints_make_consistent_in_parallel_1;
+
+      // ... determine owners of locally relevant dofs
+      IndexSet locally_relevant_dofs_non_local = locally_relevant_dofs;
+      locally_relevant_dofs_non_local.subtract_set(locally_owned_dofs);
+
+      std::vector<unsigned int> locally_relevant_dofs_owners(
+        locally_relevant_dofs_non_local.n_elements());
+      Utilities::MPI::internal::ComputeIndexOwner::ConsensusAlgorithmsPayload
+        locally_relevant_dofs_process(locally_owned_dofs,
+                                      locally_relevant_dofs_non_local,
+                                      mpi_communicator,
+                                      locally_relevant_dofs_owners,
+                                      true);
+
+      Utilities::MPI::ConsensusAlgorithms::Selector<
+        std::pair<types::global_dof_index, types::global_dof_index>,
+        unsigned int>(locally_relevant_dofs_process, mpi_communicator)
+        .run();
+
+      const auto locally_relevant_dofs_by_ranks =
+        locally_relevant_dofs_process.get_requesters();
+
+      std::map<unsigned int, std::vector<char>> send_data;
+
+      std::vector<MPI_Request> requests;
+      requests.reserve(send_data.size());
+
+      // ... send data
+      for (const auto &rank_and_indices : locally_relevant_dofs_by_ranks)
+        {
+          Assert(rank_and_indices.first != my_rank, ExcInternalError());
+
+          std::vector<ConstraintType> data;
+
+          for (const auto index : rank_and_indices.second)
+            {
+              // note: at this stage locally_relevant_constraints still
+              // contains only locally owned constraints
+              const auto prt =
+                std::find_if(locally_relevant_constraints.begin(),
+                             locally_relevant_constraints.end(),
+                             [index](const auto &a) {
+                               return std::get<0>(a) == index;
+                             });
+              if (prt != locally_relevant_constraints.end())
+                data.push_back(*prt);
+            }
+
+          send_data[rank_and_indices.first] = Utilities::pack(data, false);
+
+          requests.push_back({});
+
+          const int ierr = MPI_Isend(send_data[rank_and_indices.first].data(),
+                                     send_data[rank_and_indices.first].size(),
+                                     MPI_CHAR,
+                                     rank_and_indices.first,
+                                     tag,
+                                     mpi_communicator,
+                                     &requests.back());
+          AssertThrowMPI(ierr);
+        }
+
+      // ... receive data
+      const unsigned int n_rec_ranks = [&]() {
+        // count number of ranks from where data will be received from
+        // by looping locally_relevant_dofs_owners and identifying unique
+        // rank (ignoring the current rank)
+
+        std::set<unsigned int> rec_ranks;
+
+        for (const unsigned int rank : locally_relevant_dofs_owners)
+          if (rank != my_rank)
+            rec_ranks.insert(rank);
+
+        return rec_ranks.size();
+      }();
+
+      for (unsigned int counter = 0; counter < n_rec_ranks; ++counter)
+        {
+          MPI_Status status;
+          int ierr = MPI_Probe(MPI_ANY_SOURCE, tag, mpi_communicator, &status);
+          AssertThrowMPI(ierr);
+
+          int message_length;
+          ierr = MPI_Get_count(&status, MPI_CHAR, &message_length);
+          AssertThrowMPI(ierr);
+
+          std::vector<char> buffer(message_length);
+
+          ierr = MPI_Recv(buffer.data(),
+                          buffer.size(),
+                          MPI_CHAR,
+                          status.MPI_SOURCE,
+                          tag,
+                          mpi_communicator,
+                          MPI_STATUS_IGNORE);
+          AssertThrowMPI(ierr);
+
+          const auto received_locally_relevant_constrain =
+            Utilities::unpack<std::vector<ConstraintType>>(buffer, false);
+
+          for (const auto &locally_relevant_constrain :
+               received_locally_relevant_constrain)
+            locally_relevant_constraints.push_back(locally_relevant_constrain);
+        }
+
+      const int ierr =
+        MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
+      AssertThrowMPI(ierr);
+
+      sort_constraints();
+    }
+
+#endif
+
+    return locally_relevant_constraints;
+  }
+} // namespace internal
+
+
+
+template <typename number>
+void
+AffineConstraints<number>::make_consistent_in_parallel(
+  const IndexSet &locally_owned_dofs,
+  const IndexSet &locally_relevant_dofs,
+  const MPI_Comm  mpi_communicator)
+{
+  if (Utilities::MPI::job_supports_mpi() == false ||
+      Utilities::MPI::n_mpi_processes(mpi_communicator) == 1)
+    return; // nothing to do, since serial
+
+  Assert(sorted == false, ExcMatrixIsClosed());
+
+  // 1) get all locally relevant constraints ("temporal constraint matrix")
+  const auto temporal_constraint_matrix =
+    internal::compute_locally_relevant_constraints(*this,
+                                                   locally_owned_dofs,
+                                                   locally_relevant_dofs,
+                                                   mpi_communicator);
+
+  // 2) clear the content of this constraint matrix
+  lines.clear();
+  lines_cache.clear();
+
+  // 3) refill this constraint matrix
+  for (const auto &line : temporal_constraint_matrix)
+    {
+      const types::global_dof_index index = std::get<0>(line);
+
+      // ... line
+      this->add_line(index);
+
+      // ... inhomogeneity
+      if (std::get<1>(line) != number())
+        this->set_inhomogeneity(index, std::get<1>(line));
+
+      // ... entries
+      if (std::get<2>(line).size() > 0)
+        for (const auto &j : std::get<2>(line))
+          this->add_entry(index, j.first, j.second);
+    }
+
+#ifdef DEBUG
+  Assert(this->is_consistent_in_parallel(
+           Utilities::MPI::all_gather(mpi_communicator, locally_owned_dofs),
+           locally_relevant_dofs,
+           mpi_communicator),
+         ExcInternalError());
+#endif
+}
+
+
 
 template <typename number>
 void
diff --git a/tests/lac/constraints_make_consistent_in_parallel_01.cc b/tests/lac/constraints_make_consistent_in_parallel_01.cc
new file mode 100644 (file)
index 0000000..d92dd22
--- /dev/null
@@ -0,0 +1,121 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test AffineConstraints::make_consistent_in_parallel().
+
+#include <deal.II/base/conditional_ostream.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/affine_constraints.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+/**
+ * mpirun -np 4 ./constraints_01
+ *
+ * 6-----7-----8
+ * | 0/2 | 1/3 |
+ * 2-----3-----5
+ * | 0/0 | 0/1 |  ... with dof-index, material-id/rank
+ * 0-----1-----4
+ *
+ * Expected output: {3}, {3, 5}, {3, 7}, {3, 5, 7}
+ * Actual output:   {},  {3, 5}, {3, 7}, {3, 5, 7}
+ *
+ * ... code to reproduce https://github.com/dealii/dealii/issues/11725
+ */
+template <typename Number>
+IndexSet
+collect_lines(const AffineConstraints<Number> &constraints,
+              const unsigned int               size)
+{
+  IndexSet lines_local(size);
+  for (const auto &line : constraints.get_lines())
+    lines_local.add_index(line.index);
+  return lines_local;
+}
+
+template <int dim, int spacedim>
+void
+test(const DoFHandler<dim, spacedim> &dof_handler,
+     const IndexSet &                 locally_relevant_dofs)
+{
+  AffineConstraints<double> constraints;
+
+  std::vector<types::global_dof_index> dof_indices(
+    dof_handler.get_fe().n_dofs_per_face());
+
+  for (const auto &cell : dof_handler.active_cell_iterators())
+    for (const auto face : cell->face_indices())
+      if (cell->is_locally_owned() && !cell->at_boundary(face) &&
+          cell->material_id() != cell->neighbor(face)->material_id())
+        {
+          cell->face(face)->get_dof_indices(dof_indices);
+          for (const auto i : dof_indices)
+            constraints.add_line(i);
+        }
+
+  const auto a = collect_lines(constraints, dof_handler.n_dofs());
+  a.print(deallog.get_file_stream());
+
+  constraints.make_consistent_in_parallel(dof_handler.locally_owned_dofs(),
+                                          locally_relevant_dofs,
+                                          dof_handler.get_communicator());
+
+  const auto b = collect_lines(constraints, dof_handler.n_dofs());
+  b.print(deallog.get_file_stream());
+}
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi(argc, argv, 1);
+  MPILogInitAll                    all;
+
+  const int dim = 2;
+
+  ConditionalOStream pcout(std::cout,
+                           Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) ==
+                             0);
+
+  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  GridGenerator::subdivided_hyper_cube(tria, 2);
+
+  DoFHandler dof_handler(tria);
+  dof_handler.distribute_dofs(FE_Q<dim>(1));
+
+  for (const auto &cell : tria.active_cell_iterators())
+    if (cell->center()[0] > 0.5 && cell->center()[1] > 0.5)
+      cell->set_material_id(1);
+
+  IndexSet locally_active_dofs;
+  DoFTools::extract_locally_active_dofs(dof_handler, locally_active_dofs);
+  test(dof_handler, locally_active_dofs);
+
+  IndexSet locally_relevant_dofs;
+  DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
+  test(dof_handler, locally_relevant_dofs);
+}
diff --git a/tests/lac/constraints_make_consistent_in_parallel_01.with_p4est=true.mpirun=4.output b/tests/lac/constraints_make_consistent_in_parallel_01.with_p4est=true.mpirun=4.output
new file mode 100644 (file)
index 0000000..e1392e6
--- /dev/null
@@ -0,0 +1,23 @@
+
+{}
+{3}
+{}
+{3, 5, 7}
+
+{3, 5}
+{3, 5}
+{3, 5}
+{3, 5, 7}
+
+
+{3, 7}
+{3, 7}
+{3, 7}
+{3, 5, 7}
+
+
+{3, 5, 7}
+{3, 5, 7}
+{3, 5, 7}
+{3, 5, 7}
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.