From bd54f0870a072a6b38c6212337897fa79f5c3d2b Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Tue, 27 Aug 2024 20:49:23 -0600 Subject: [PATCH] Don't merge, just drop constraints. --- include/deal.II/lac/affine_constraints.h | 43 +++- .../lac/affine_constraints.templates.h | 208 ++++++------------ 2 files changed, 103 insertions(+), 148 deletions(-) diff --git a/include/deal.II/lac/affine_constraints.h b/include/deal.II/lac/affine_constraints.h index 76438d9b5d..b2332a8db4 100644 --- a/include/deal.II/lac/affine_constraints.h +++ b/include/deal.II/lac/affine_constraints.h @@ -1994,17 +1994,44 @@ public: * 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(). + * owned on the current process, that is typically what you get from + * 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 + * these might be depends on your actual use case; in general, these are + * the degrees of freedoms for which you compute contributions to the matrix + * and right hand side of a linear system. 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(). + * DoFTools::extract_locally_active_dofs(). This is because you compute + * contributions for the linear system only on the locally owned cells, + * and the degrees of freedom located there -- i.e., exactly the locally + * active degrees of freedom (see + * @ref GlossLocallyActiveDof "this glossary entry"). On the other hand, + * for discontinuous Galerkin methods, you might need consistent information + * on all locally relevant DoFs due to the need to compute face integrals + * between locally owned and ghost cells, with contributions computed also + * for matrix and right hand side entries that correspond to degrees of + * freedom on ghost cells, i.e., DoFs that are owned by other processes. + * In that case, you need to know about constraints also for degrees of + * freedom on locally relevant DoFs, and in that case you want to pass + * as second argument to this function what you get from + * DoFTools::extract_locally_relevant_dofs(). It is worth noting that the + * same happens if you use a continuous Galerkin method but if your + * bilinear form contains terms computed on faces -- for example to + * penalize certain terms when solving the biharmonic equation (see, + * for example, step-47). + * + * At the end of this function, the current object + * only stores constraints for degrees of freedom whose indices are + * listed in the second argument. In other words, while + * DoFTools::make_hanging_node_constraints() may compute constraints + * also for degrees of freedom at the interface between two ghost cells + * (these DoFs are locally *relevant*), if you pass as second argument only + * the set of locally *active* DoFs, this function will at the end have + * forgotten constraints about all DoFs that are locally relevant but + * not locally active. * * This function updates the locally stored constraints (or local_lines) of * this object, whenever a DoF is constrained against another that we have not @@ -2017,7 +2044,9 @@ public: * locally_owned_dofs, constraints.get_local_lines(), mpi_communicator); * @endcode * - * @note One should call this function before calling close(). + * @note This function internally calls close(). As a consequence, whether + * or not you have called close() before, the current object will be + * closed after this function. */ void make_consistent_in_parallel(const IndexSet &locally_owned_dofs, diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index 16a9f1b689..2a7131a21c 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -320,149 +320,68 @@ namespace internal const unsigned int my_rank = Utilities::MPI::this_mpi_process(mpi_communicator); - // helper function - const auto sort_and_make_unique = [&constraints_in, &locally_owned_dofs]( - std::vector - &locally_relevant_constraints) { - (void)constraints_in; - if (locally_relevant_constraints.empty()) - return; - - for (ConstraintLine &entry : locally_relevant_constraints) - std::sort(entry.entries.begin(), - entry.entries.end(), - [](const auto &l1, const auto &l2) { - return l1.first < l2.first; + // First define a helper function that sorts and normalizes the constraints + // provided through the function's argument: + const auto sort_and_make_unique = + [](std::vector &locally_relevant_constraints) { + if (locally_relevant_constraints.empty()) + return; + + // First sort the array of constraints by their index: + std::sort(locally_relevant_constraints.begin(), + locally_relevant_constraints.end(), + [](const ConstraintLine &l1, const ConstraintLine &l2) { + return l1.index < l2.index; }); - std::sort(locally_relevant_constraints.begin(), - locally_relevant_constraints.end(), - [](const ConstraintLine &l1, const ConstraintLine &l2) { - return l1.index < l2.index; - }); - - const auto equal_with_tol = [](const number d, const number e) { - if (std::abs(std::real(d - e)) < - std::sqrt( - std::numeric_limits< - typename numbers::NumberTraits::real_type>::epsilon())) - return true; - else - return false; + // It is possible that two processes have computed constraints + // for the same DoF differently (for example, in parallel because + // they have visited different faces of the same cell because + // on different processes, different neighbor cells are artificial. + // + // This is ok, as long as after resolution of all chains of + // constraints we end up with the same constraint. But at this + // point, we just don't know yet. We deal with this by simply + // dropping constraints for DoFs for which there is a previous + // constraint in the list: + locally_relevant_constraints.erase( + std::unique(locally_relevant_constraints.begin(), + locally_relevant_constraints.end(), + [](const ConstraintLine &a, const ConstraintLine &b) { + return (a.index == b.index); + }), + locally_relevant_constraints.end()); + + // For those constraints that survive, sort the right hand side arrays: + for (ConstraintLine &entry : locally_relevant_constraints) + std::sort(entry.entries.begin(), + entry.entries.end(), + [](const auto &l1, const auto &l2) { + return l1.first < l2.first; + }); }; - typename std::vector::iterator read_ptr = - locally_relevant_constraints.begin(); - typename std::vector::iterator write_ptr = - locally_relevant_constraints.begin(); - // go through sorted locally relevant constraints and look out for - // duplicates (same constrained index, same entries) or cases that need - // to be augmented (same index, different entries) - while (++read_ptr != locally_relevant_constraints.end()) - { - if (read_ptr->index != write_ptr->index) - { - // if we have a different index, use it here - if (++write_ptr != read_ptr) - *write_ptr = std::move(*read_ptr); - } - else // equal global dof index - { - ConstraintLine &a = *write_ptr; - const ConstraintLine &b = *read_ptr; - Assert(a.index == b.index, ExcInternalError()); - if (!equal_with_tol(a.inhomogeneity, b.inhomogeneity) && - locally_owned_dofs.is_element(b.index)) - { - Assert(equal_with_tol(b.inhomogeneity, - constraints_in.get_inhomogeneity( - b.index)), - ExcInternalError()); - a.inhomogeneity = b.inhomogeneity; - } - - typename ConstraintLine::Entries &av = a.entries; - const typename ConstraintLine::Entries &bv = b.entries; - // check if entries vectors are equal - bool vectors_are_equal = (av.size() == bv.size()); - for (unsigned int i = 0; vectors_are_equal && i < av.size(); ++i) - { - if (av[i].first != bv[i].first) - vectors_are_equal = false; - else - Assert(equal_with_tol(av[i].second, bv[i].second), - ExcInternalError()); - } - - // merge entries vectors if different, otherwise ignore the - // second entry - if (!vectors_are_equal) - { - typename ConstraintLine::Entries cv; - cv.reserve(av.size() + bv.size()); - - unsigned int i = 0; - unsigned int j = 0; - - for (; (i < av.size()) && (j < bv.size());) - { - if (av[i].first == bv[j].first) - { - Assert(equal_with_tol(av[i].second, bv[j].second), - ExcInternalError()); - - cv.push_back(av[i]); - ++i; - ++j; - } - else if (av[i].first < bv[j].first) - { - cv.push_back(av[i]); - ++i; - } - else - { - cv.push_back(bv[j]); - ++j; - } - } - - for (; i < av.size(); ++i) - cv.push_back(av[i]); - - for (; j < bv.size(); ++j) - cv.push_back(bv[j]); - - std::swap(av, cv); - } - } - } - ++write_ptr; - locally_relevant_constraints.erase(write_ptr, - locally_relevant_constraints.end()); - }; - // step 0: Collect the indices of constrained DoFs we know of: - IndexSet constrained_indices(locally_owned_dofs.size()); - - std::vector constrained_indices_temp; - constrained_indices_temp.reserve(constraints_in.n_constraints()); - for (const auto &line : constraints_in.get_lines()) - constrained_indices_temp.push_back(line.index); - std::sort(constrained_indices_temp.begin(), constrained_indices_temp.end()); - constrained_indices.add_indices(constrained_indices_temp.begin(), - constrained_indices_temp.end()); + IndexSet my_constraint_indices(locally_owned_dofs.size()); + { + std::vector indices; + indices.reserve(constraints_in.n_constraints()); + for (const auto &line : constraints_in.get_lines()) + indices.push_back(line.index); + std::sort(indices.begin(), indices.end()); + my_constraint_indices.add_indices(indices.begin(), indices.end()); + } // step 1: Identify the owners of DoFs we know to be constrained. // Of course, we are likely to own a good share of the DoFs // we know to be constrained -- and that is ok. - std::vector constrained_indices_owners( - constrained_indices.n_elements()); + std::vector owners_of_my_constraints( + my_constraint_indices.n_elements()); Utilities::MPI::internal::ComputeIndexOwner::ConsensusAlgorithmsPayload constrained_indices_process(locally_owned_dofs, - constrained_indices, + my_constraint_indices, mpi_communicator, - constrained_indices_owners, + owners_of_my_constraints, true); Utilities::MPI::ConsensusAlgorithms::Selector< @@ -471,20 +390,20 @@ namespace internal consensus_algorithm; consensus_algorithm.run(constrained_indices_process, mpi_communicator); - // step 2: collect all locally owned constraints + // step 2: Collect all locally owned constraints into a data structure + // that we can send to other processes that want to know + // about them. const std::map constrained_indices_by_ranks = constrained_indices_process.get_requesters(); - { - // ... collect data and sort according to owner std::map> send_data; - for (unsigned int i = 0; i < constrained_indices_owners.size(); ++i) + for (unsigned int i = 0; i < owners_of_my_constraints.size(); ++i) { ConstraintLine entry; const types::global_dof_index index = - constrained_indices.nth_index_in_set(i); + my_constraint_indices.nth_index_in_set(i); entry.index = index; @@ -495,10 +414,10 @@ namespace internal *constraints = constraints_in.get_constraint_entries(index)) entry.entries = *constraints; - if (constrained_indices_owners[i] == my_rank) + if (owners_of_my_constraints[i] == my_rank) locally_relevant_constraints.push_back(entry); else - send_data[constrained_indices_owners[i]].push_back(entry); + send_data[owners_of_my_constraints[i]].push_back(entry); } // Now exchange this data between processes: @@ -515,7 +434,7 @@ namespace internal } // step 3: communicate constraints so that each process knows how the - // locally locally relevant dofs are constrained + // locally relevant dofs are constrained { // ... determine owners of locally relevant dofs IndexSet locally_relevant_dofs_non_local = locally_relevant_dofs; @@ -593,11 +512,18 @@ AffineConstraints::make_consistent_in_parallel( const IndexSet &constraints_to_make_consistent_, const MPI_Comm mpi_communicator) { + // We need to resolve chains of constraints in this function, and that + // is exactly what close() does. So call this function at the beginning + // to make sure we start working with constraints that are already + // resolved to the degree possible before exchanging information. + // We do this here so that it is also done if we choose the early + // exit in the following statement. + if (sorted == false) + close(); + if (Utilities::MPI::n_mpi_processes(mpi_communicator) == 1) return; // Nothing to do, since serial. - Assert(sorted == false, ExcMatrixIsClosed()); - Assert(this->local_lines.size() > 0, ExcMessage( "This functionality requires that the AffineConstraints object " -- 2.39.5