]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Let AffineConstraints::make_consistent_in_parallel() iterate until converged.
authorMarc Fehling <mafehling.git@gmail.com>
Wed, 26 Jun 2024 11:09:43 +0000 (13:09 +0200)
committerMarc Fehling <mafehling.git@gmail.com>
Mon, 1 Jul 2024 16:02:33 +0000 (18:02 +0200)
doc/news/changes/minor/20240620BangerthFehlingKronbichlerMunch [new file with mode: 0644]
include/deal.II/lac/affine_constraints.h
include/deal.II/lac/affine_constraints.templates.h

diff --git a/doc/news/changes/minor/20240620BangerthFehlingKronbichlerMunch b/doc/news/changes/minor/20240620BangerthFehlingKronbichlerMunch
new file mode 100644 (file)
index 0000000..88f8cf1
--- /dev/null
@@ -0,0 +1,19 @@
+Fixed: AffineConstraints::make_consistent_in_parallel() now iteratively
+updates the constraints on user-specified DoFs until the set of constraints
+globally converged on all subdomains.
+<br>
+The function now also updates the locally stored constraints of the
+underlying AffineConstraints object. After using this function, it might be
+necessary to use the IndexSet retrieved by
+AffineConstraints::get_local_lines() when initializing data structures,
+for example:
+<code>
+DynamicSparsityPattern dsp(constraints.get_local_lines());
+LinearAlgebra::distributed::Vector<double> solution(
+  locally_owned_dofs, constraints.get_local_lines(), mpi_communicator);
+</code>
+This used to be an issue in parallel hp-adaptive applications, when
+finite elements of different types have constraints on faces between
+locally relevant and artificial cells.
+<br>
+(Wolfgang Bangerth, Marc Fehling, Martin Kronbichler, Peter Munch, 2024/06/20)
index 30a058e4de859e46268c6e73eacdc1fcaeb6d32d..f5cbecf8096f8a492644058ba966819ff20a2fad 100644 (file)
@@ -1990,13 +1990,38 @@ public:
                             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().
+   * Make the constraints that are locally stored on the current object
+   * consistent on all MPI processes in a distributed computation.
+   *
+   * The IndexSet @p locally_owned_dofs conforms to all DoFs that are locally
+   * owned on the current process, i.e., DoFHandler::locally_owned_dofs().
+   *
+   * The IndexSet @p constraints_to_make_consistent shall contain all DoF
+   * indices for which we want to store the complete constraints for. Which DoFs
+   * these might be depends on your actual use case: For most applications in
+   * which you are using continuous Galerkin methods, you would want to use
+   * locally active DoFs here, obtained via
+   * DoFTools::extract_locally_active_dofs(). However in discontinuous Galerkin
+   * methods, you might need consistent information on all locally relevant
+   * DoFs due to the need to compute face integrals against ghosted cells,
+   * for which you need DoFTools::extract_locally_relevant_dofs().
+   *
+   * This function updates the locally stored constraints (or local_lines) of
+   * this object, whenever a DoF is constrained against another that we have not
+   * known of before. After calling this function, it might be necessary to use
+   * the IndexSet retrieved by get_local_lines() when initializing data
+   * structures built on top of the AffineConstraints object, for example:
+   * @code
+   * DynamicSparsityPattern dsp(constraints.get_local_lines());
+   * LinearAlgebra::distributed::Vector<double> solution(
+   *   locally_owned_dofs, constraints.get_local_lines(), mpi_communicator);
+   * @endcode
+   *
+   * @note One should call this function before calling close().
    */
   void
   make_consistent_in_parallel(const IndexSet &locally_owned_dofs,
-                              const IndexSet &locally_relevant_dofs,
+                              const IndexSet &constraints_to_make_consistent,
                               const MPI_Comm  mpi_communicator);
 
   /**
index 902ea9ad83c6ac73700a7698b6145bc96b6aba20..98c1d654413b3a62f4b8ae5e6b393c241266a84a 100644 (file)
@@ -597,36 +597,97 @@ template <typename number>
 void
 AffineConstraints<number>::make_consistent_in_parallel(
   const IndexSet &locally_owned_dofs,
-  const IndexSet &locally_relevant_dofs,
+  const IndexSet &constraints_to_make_consistent_,
   const MPI_Comm  mpi_communicator)
 {
   if (Utilities::MPI::n_mpi_processes(mpi_communicator) == 1)
-    return; // nothing to do, since serial
+    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();
+  Assert(this->local_lines.is_empty() == false,
+         ExcMessage(
+           "This functionality requires that the AffineConstraints object "
+           "knows for which degrees of freedom it can store constraints. "
+           "Please initialize this object with the corresponding index sets."));
+
+  // Container for indices that are constrained or that other indices are
+  // constrained against. Generously reserve memory for it, assuming the worst
+  // case that each constraint line involves two distinct DoF indices.
+  std::vector<types::global_dof_index> constrained_indices;
+  constrained_indices.reserve(2 * n_constraints());
+
+  // This IndexSet keeps track of the locally stored constraints on this
+  // AffineConstraints object: If we have received constrained indices that we
+  // don't know of, we need to expand our playing field.
+  IndexSet locally_stored_constraints;
+
+  // This IndexSet contains those DoFs about which we want to know all
+  // constraints. We will receive constraints for these DoFs against other DoFs,
+  // which might be constrained themselves. To successfully resolve all chains
+  // of constraints, we need to know the constraints of all these DoFs, and we
+  // keep track of the relevant DoFs here.
+  IndexSet constraints_to_make_consistent = constraints_to_make_consistent_;
+
+  // This IndexSet stores DoFs from the previous iteration. If the old and new
+  // index sets match, we have converged.
+  IndexSet constraints_made_consistent;
+
+  const unsigned int max_iterations  = 10;
+  unsigned int       iteration_count = 0;
+  for (; iteration_count < max_iterations; ++iteration_count)
+    {
+      // 1) Get all locally relevant constraints ("temporal constraint matrix").
+      const auto temporal_constraint_matrix =
+        internal::compute_locally_relevant_constraints(
+          *this,
+          locally_owned_dofs,
+          constraints_to_make_consistent,
+          mpi_communicator);
+
+      // 2) Add untracked DoFs to the index sets.
+      constrained_indices.clear();
+      for (const auto &line : temporal_constraint_matrix)
+        {
+          constrained_indices.push_back(line.index);
+          for (const auto &entry : line.entries)
+            constrained_indices.push_back(entry.first);
+        }
+      std::sort(constrained_indices.begin(), constrained_indices.end());
+
+      locally_stored_constraints = this->local_lines;
+      locally_stored_constraints.add_indices(constrained_indices.begin(),
+                                             constrained_indices.end());
+
+      constraints_made_consistent = constraints_to_make_consistent;
+      constraints_to_make_consistent.add_indices(constrained_indices.begin(),
+                                                 constrained_indices.end());
+
+      // 3) Clear and refill this constraint matrix.
+      this->reinit(locally_owned_dofs, locally_stored_constraints);
+      for (const auto &line : temporal_constraint_matrix)
+        this->add_constraint(line.index, line.entries, line.inhomogeneity);
+
+      // 4) Stop loop if converged.
+      const auto constraints_converged =
+        Utilities::MPI::min(static_cast<unsigned int>(
+                              constraints_to_make_consistent ==
+                              constraints_made_consistent),
+                            mpi_communicator);
+      if (constraints_converged)
+        break;
+    }
 
-  // 3) refill this constraint matrix
-  for (const auto &line : temporal_constraint_matrix)
-    this->add_constraint(line.index, line.entries, line.inhomogeneity);
+  AssertThrow(iteration_count < max_iterations,
+              ExcMessage(
+                "make_consistent_in_parallel() did not converge after " +
+                Utilities::to_string(max_iterations) + " iterations."));
 
-#ifdef DEBUG
   Assert(this->is_consistent_in_parallel(
            Utilities::MPI::all_gather(mpi_communicator, locally_owned_dofs),
-           locally_relevant_dofs,
+           constraints_to_make_consistent,
            mpi_communicator),
          ExcInternalError());
-#endif
 }
 
 

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.