From: Marc Fehling Date: Wed, 26 Jun 2024 11:09:43 +0000 (+0200) Subject: Let AffineConstraints::make_consistent_in_parallel() iterate until converged. X-Git-Tag: v9.6.0-rc1~148^2~1 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3bae6142c1af5ba4d6137c00aef64eff7b878339;p=dealii.git Let AffineConstraints::make_consistent_in_parallel() iterate until converged. --- diff --git a/doc/news/changes/minor/20240620BangerthFehlingKronbichlerMunch b/doc/news/changes/minor/20240620BangerthFehlingKronbichlerMunch new file mode 100644 index 0000000000..88f8cf192c --- /dev/null +++ b/doc/news/changes/minor/20240620BangerthFehlingKronbichlerMunch @@ -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. +
+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: + +DynamicSparsityPattern dsp(constraints.get_local_lines()); +LinearAlgebra::distributed::Vector solution( + locally_owned_dofs, constraints.get_local_lines(), mpi_communicator); + +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. +
+(Wolfgang Bangerth, Marc Fehling, Martin Kronbichler, Peter Munch, 2024/06/20) diff --git a/include/deal.II/lac/affine_constraints.h b/include/deal.II/lac/affine_constraints.h index 30a058e4de..f5cbecf809 100644 --- a/include/deal.II/lac/affine_constraints.h +++ b/include/deal.II/lac/affine_constraints.h @@ -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 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); /** diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index 902ea9ad83..98c1d65441 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -597,36 +597,97 @@ template void AffineConstraints::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 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( + 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 }