From: Maximilian Bergbauer Date: Thu, 16 Mar 2023 16:56:39 +0000 (+0100) Subject: Add consistent strategy to merge constraints in parallel X-Git-Tag: v9.6.0-rc1~115^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c64cd1ad9e0df23cbc043344350f2ef00aa4e925;p=dealii.git Add consistent strategy to merge constraints in parallel --- diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index 98c1d65441..5bc85a5ab7 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -321,23 +321,125 @@ namespace internal Utilities::MPI::this_mpi_process(mpi_communicator); // helper function - const auto sort_and_make_unique = - [](std::vector &constraints) { - std::sort(constraints.begin(), - constraints.end(), - [](const ConstraintType &l1, const ConstraintType &l2) { - return l1.index < l2.index; + const auto sort_and_make_unique = [&constraints_in, &locally_owned_dofs]( + std::vector + &locally_relevant_constraints) { + if (locally_relevant_constraints.empty()) + return; + + for (auto &entry : locally_relevant_constraints) + std::sort(entry.entries.begin(), + entry.entries.end(), + [](const auto &l1, const auto &l2) { + return l1.first < l2.first; }); - constraints.erase(std::unique(constraints.begin(), - constraints.end(), - [](const ConstraintType &l1, - const ConstraintType &l2) { - return l1.index == l2.index; - }), - constraints.end()); + std::sort(locally_relevant_constraints.begin(), + locally_relevant_constraints.end(), + [](const ConstraintType &l1, const ConstraintType &l2) { + return l1.index < l2.index; + }); + + 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; }; + auto read_ptr = locally_relevant_constraints.begin(); + auto 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 + { + auto &a = *write_ptr; + const auto &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; + } + + auto &av = a.entries; + const auto &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 dealii::AffineConstraints< + number>::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()); + }; + // 0) collect constrained indices of the current object IndexSet constrained_indices(locally_owned_dofs.size()); diff --git a/tests/lac/constraints_make_consistent_in_parallel_03.cc b/tests/lac/constraints_make_consistent_in_parallel_03.cc new file mode 100644 index 0000000000..37fa665648 --- /dev/null +++ b/tests/lac/constraints_make_consistent_in_parallel_03.cc @@ -0,0 +1,75 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2018 - 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + + +// test AffineConstraints::make_consistent_in_parallel for case +// where constraints need to be combined between different ranks + +#include + +#include + +#include "../tests.h" + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv); + MPILogInitAll all; + + const unsigned int my_proc = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + const unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + + AssertThrow(n_procs == 3, ExcMessage("Only working on 3 ranks")); + + const unsigned int n_indices = 5; + IndexSet owned_elements(n_indices); + IndexSet relevant_elements(n_indices); + AffineConstraints constraints; + if (my_proc == 0) + { + owned_elements.add_index(0); + relevant_elements.add_range(0, 3); + constraints.reinit(owned_elements, relevant_elements); + } + else if (my_proc == 1) + { + owned_elements.add_range(1, 3); + relevant_elements.add_range(0, n_indices); + constraints.reinit(owned_elements, relevant_elements); + constraints.add_line(0); + constraints.add_entry(0, 1, 0.5); + constraints.add_entry(0, 2, 0.5); + constraints.add_line(2); + constraints.add_entry(2, 3, 1.5); + } + else if (my_proc == 2) + { + owned_elements.add_range(3, n_indices); + relevant_elements.add_range(0, n_indices); + constraints.reinit(owned_elements, relevant_elements); + constraints.add_line(0); + constraints.add_entry(0, 1, 0.5); + constraints.add_entry(0, 3, 0.5); + } + + constraints.make_consistent_in_parallel(owned_elements, + relevant_elements, + MPI_COMM_WORLD); + constraints.close(); + constraints.print(deallog.get_file_stream()); + + return 0; +} diff --git a/tests/lac/constraints_make_consistent_in_parallel_03.mpirun=3.output b/tests/lac/constraints_make_consistent_in_parallel_03.mpirun=3.output new file mode 100644 index 0000000000..452ceebf8d --- /dev/null +++ b/tests/lac/constraints_make_consistent_in_parallel_03.mpirun=3.output @@ -0,0 +1,14 @@ + + 0 1: 0.500000 + 0 3: 1.25000 + 2 3: 1.50000 + + 0 1: 0.500000 + 0 3: 1.25000 + 2 3: 1.50000 + + + 0 1: 0.500000 + 0 3: 1.25000 + 2 3: 1.50000 +