From: Sebastian Kinnewig Date: Fri, 12 Jan 2024 08:38:54 +0000 (+0100) Subject: add nonlocal_vector X-Git-Tag: relicensing~148^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=729542dc8fdfaad52a5f6cfae5ee3416ac7b058f;p=dealii.git add nonlocal_vector --- diff --git a/include/deal.II/lac/affine_constraints.templates.h b/include/deal.II/lac/affine_constraints.templates.h index b99edfb2fe..a5048ba3a8 100644 --- a/include/deal.II/lac/affine_constraints.templates.h +++ b/include/deal.II/lac/affine_constraints.templates.h @@ -24,6 +24,7 @@ #include #include #include +#include #include #include @@ -2666,12 +2667,10 @@ namespace internal IndexSet parallel_partitioner = locally_owned_elements; parallel_partitioner.add_indices(needed_elements); - const MPI_Comm mpi_comm = - Trilinos::teuchos_comm_to_mpi_comm(vec.trilinos_rcp()->getMap()->getComm()); + const MPI_Comm mpi_comm = Utilities::Trilinos::teuchos_comm_to_mpi_comm( + vec.trilinos_rcp()->getMap()->getComm()); - output.reinit(locally_owned_elements, - needed_elements, - mpi_comm); + output.reinit(locally_owned_elements, needed_elements, mpi_comm); output = vec; } diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h index 856bb7456e..a34e7d87ba 100644 --- a/include/deal.II/lac/trilinos_tpetra_vector.h +++ b/include/deal.II/lac/trilinos_tpetra_vector.h @@ -288,6 +288,7 @@ namespace LinearAlgebra Vector(const Teuchos::RCP V); /** + * TODO: This is not used * This constructor takes an IndexSet that defines how to distribute the * individual components among the MPI processors. Since it also * includes information about the size of the vector, this is all we @@ -296,6 +297,25 @@ namespace LinearAlgebra explicit Vector(const IndexSet ¶llel_partitioner, const MPI_Comm communicator); + /** + * In addition to just specifying one index set as in all the other + * methods above, this method allows to supply an additional set of + * ghost entries. + * + * Depending on whether the @p locally_relevant_or_ghost_entries argument uniquely + * subdivides elements among processors or not, the resulting vector may + * or may not have ghost elements. See the general documentation of this + * class for more information. + * + * @see + * @ref GlossGhostedVector "vectors with ghost elements" + + */ + explicit Vector(const IndexSet &locally_owned_entries, + const IndexSet &ghost_entries, + const MPI_Comm communicator, + const bool vector_writable = false); + /** * Reinit functionality. This function destroys the old vector content * and generates a new one based on the input partitioning. The flag @@ -307,7 +327,6 @@ namespace LinearAlgebra const MPI_Comm communicator, const bool omit_zeroing_entries = false); - /** * Reinit functionality. This function destroys the old vector content * and generates a new one based on the input partitioning. In addition @@ -325,7 +344,8 @@ namespace LinearAlgebra void reinit(const IndexSet &locally_owned_entries, const IndexSet &locally_relevant_or_ghost_entries, - const MPI_Comm communicator = MPI_COMM_WORLD); + const MPI_Comm communicator = MPI_COMM_WORLD, + const bool vector_writable = false); /** * Change the dimension to that of the vector @p V. The elements of @p V are not @@ -842,6 +862,12 @@ namespace LinearAlgebra create_tpetra_comm_pattern(const IndexSet &source_index_set, const MPI_Comm mpi_comm); + /** + * A boolean variable to hold information on whether the vector is + * compressed or not. + */ + bool compressed; + /** * Store whether the vector has ghost elements or not. * @@ -858,17 +884,11 @@ namespace LinearAlgebra Teuchos::RCP vector; /** - * Pointer to the Tpetra::Map storing the owned elemenets - * per rank. This map has a one-to-one mapping. - */ - Teuchos::RCP locally_owned_map; - - /** - * Pointer to the Tpetra::Map storing the relevant elements - * per rank. There are entries in this map, that belong to multiple - * ranks. + * A vector object in Trilinos to be used for collecting the non-local + * elements if the vector was constructed with an additional IndexSet + * describing ghost elements. */ - Teuchos::RCP locally_relevant_map; + Teuchos::RCP nonlocal_vector; /** * IndexSet of the elements of the last imported vector. @@ -902,7 +922,7 @@ namespace LinearAlgebra inline bool Vector::is_compressed() const { - return !has_ghost; + return compressed; } @@ -928,15 +948,24 @@ namespace LinearAlgebra const Number *values) { # if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) - auto vector_2d = vector->template getLocalView( + auto vector_2d_local = vector->template getLocalView( Tpetra::Access::ReadWrite); + auto vector_2d_nonlocal = + nonlocal_vector->template getLocalView( + Tpetra::Access::ReadWrite); # else vector->template sync(); - auto vector_2d = vector->template getLocalView(); + auto vector_2d_local = vector->template getLocalView(); + auto vector_2d_nonlocal = + nonlocal_vector->template getLocalView(); # endif - auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0); + + auto vector_1d_local = Kokkos::subview(vector_2d_local, Kokkos::ALL(), 0); + auto vector_1d_nonlocal = + Kokkos::subview(vector_2d_nonlocal, Kokkos::ALL(), 0); # if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) vector->template modify(); + nonlocal_vector->template modify(); # endif for (size_type i = 0; i < n_elements; ++i) @@ -945,6 +974,15 @@ namespace LinearAlgebra TrilinosWrappers::types::int_type local_row = vector->getMap()->getLocalElement(row); + // check if the index is in the non local set + bool nonlocal = false; + if (local_row == Teuchos::OrdinalTraits::invalid()) + { + local_row = nonlocal_vector->getMap()->getLocalElement(row); + nonlocal = true; + compressed = false; + } + # if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0) Assert( local_row != Teuchos::OrdinalTraits::invalid(), @@ -963,13 +1001,19 @@ namespace LinearAlgebra # endif if (local_row != Teuchos::OrdinalTraits::invalid()) - vector_1d(local_row) += values[i]; + if (nonlocal) + vector_1d_nonlocal(local_row) += values[i]; + else + vector_1d_local(local_row) += values[i]; } # if !DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) vector->template sync< typename Tpetra::Vector:: device_type::memory_space>(); + nonlocal_vector->template sync< + typename Tpetra::Vector:: + device_type::memory_space>(); # endif } @@ -981,6 +1025,10 @@ namespace LinearAlgebra const size_type *indices, const Number *values) { + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); + # if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) auto vector_2d = vector->template getLocalView( Tpetra::Access::ReadWrite); diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h index 2a7f11edb6..e119b4f5ad 100644 --- a/include/deal.II/lac/trilinos_tpetra_vector.templates.h +++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h @@ -47,6 +47,7 @@ namespace LinearAlgebra template Vector::Vector() : Subscriptor() + , compressed(true) , has_ghost(false) , vector(Utilities::Trilinos::internal::make_rcp( Utilities::Trilinos::internal::make_rcp( @@ -60,12 +61,10 @@ namespace LinearAlgebra template Vector::Vector(const Vector &V) : Subscriptor() + , compressed(V.compressed) , has_ghost(V.has_ghost) - , vector(Utilities::Trilinos::internal::make_rcp( - V.trilinos_vector(), - Teuchos::Copy)) - , locally_owned_map(V.locally_owned_map) - , locally_relevant_map(V.locally_relevant_map) + , vector(V.vector) + , nonlocal_vector(V.nonlocal_vector) {} @@ -73,20 +72,10 @@ namespace LinearAlgebra template Vector::Vector(const Teuchos::RCP V) : Subscriptor() - , has_ghost(false) + , compressed(true) + , has_ghost(V->getMap()->isOneToOne() == false) , vector(V) - { - if (vector->getMap()->isOneToOne()) - { - locally_owned_map = - Teuchos::rcp_const_cast(vector->getMap()); - } - else - { - locally_relevant_map = - Teuchos::rcp_const_cast(vector->getMap()); - } - } + {} @@ -94,12 +83,54 @@ namespace LinearAlgebra Vector::Vector(const IndexSet ¶llel_partitioner, const MPI_Comm communicator) : Subscriptor() + , compressed(true) , has_ghost(false) - , locally_owned_map( - parallel_partitioner.make_tpetra_map_rcp(communicator, false)) + , vector(Utilities::Trilinos::internal::make_rcp( + parallel_partitioner.make_tpetra_map_rcp(communicator, true))) + {} + + + + template + Vector::Vector(const IndexSet &locally_owned_entries, + const IndexSet &ghost_entries, + const MPI_Comm communicator, + const bool vector_writable) + : Subscriptor() { - vector = Utilities::Trilinos::internal::make_rcp( - locally_owned_map); + if (!vector_writable) + { + IndexSet parallel_partitioner = locally_owned_entries; + parallel_partitioner.add_indices(ghost_entries); + + vector = Utilities::Trilinos::internal::make_rcp( + parallel_partitioner.make_tpetra_map_rcp(communicator, true)); + + compressed = true; + } + else + { + Teuchos::RCP map = + locally_owned_entries.make_tpetra_map_rcp(communicator, false); + vector = Utilities::Trilinos::internal::make_rcp(map); + + IndexSet nonlocal_entries(ghost_entries); + nonlocal_entries.subtract_set(locally_owned_entries); + nonlocal_vector = Utilities::Trilinos::internal::make_rcp( + nonlocal_entries.make_tpetra_map_rcp(communicator, true)); + + compressed = false; + } + + has_ghost = (vector->getMap()->isOneToOne() == false); + +# ifdef DEBUG + MPI_Comm comm = Utilities::Trilinos::teuchos_comm_to_mpi_comm( + vector->getMap()->getComm()); + const size_type n_elements_global = + Utilities::MPI::sum(vector->getLocalLength(), comm); + Assert(has_ghost || n_elements_global == size(), ExcInternalError()); +# endif } @@ -108,28 +139,15 @@ namespace LinearAlgebra void Vector::reinit(const IndexSet ¶llel_partitioner, const MPI_Comm communicator, - const bool omit_zeroing_entries) + const bool /*omit_zeroing_entries*/) { - // release memory before reallocation vector.reset(); - locally_owned_map.reset(); - locally_relevant_map.reset(); - source_stored_elements.clear(); - - // Here we assume, that the index set is non-overlapping. - // If the index set is overlapping the function - // make_tpetra_map_rcp() will throw an assert. - has_ghost = false; - locally_owned_map = - parallel_partitioner.make_tpetra_map_rcp(communicator, false); - - if (!vector->getMap()->isSameAs(*locally_owned_map)) - vector = Utilities::Trilinos::internal::make_rcp( - locally_owned_map); - if (!omit_zeroing_entries) - { - vector->putScalar(0.); - } + nonlocal_vector.reset(); + + compressed = true; + has_ghost = false; + vector = Utilities::Trilinos::internal::make_rcp( + parallel_partitioner.make_tpetra_map_rcp(communicator, true)); } @@ -138,25 +156,47 @@ namespace LinearAlgebra void Vector::reinit(const IndexSet &locally_owned_entries, const IndexSet &ghost_entries, - const MPI_Comm communicator) + const MPI_Comm communicator, + const bool vector_writable) { // release memory before reallocation - vector.reset(); - locally_owned_map.reset(); - locally_relevant_map.reset(); + nonlocal_vector.reset(); + + if (!vector_writable) + { + vector.reset(); + + IndexSet parallel_partitioner = locally_owned_entries; + parallel_partitioner.add_indices(ghost_entries); - locally_owned_map = - locally_owned_entries.make_tpetra_map_rcp(communicator, false); + vector = Utilities::Trilinos::internal::make_rcp( + parallel_partitioner.make_tpetra_map_rcp(communicator, true)); - IndexSet parallel_partitioner = locally_owned_entries; - parallel_partitioner.add_indices(ghost_entries); - locally_relevant_map = - parallel_partitioner.make_tpetra_map_rcp(communicator, true); + compressed = true; + } + else + { + Teuchos::RCP map = + locally_owned_entries.make_tpetra_map_rcp(communicator, false); + + if (!vector->getMap()->isSameAs(*map)) + { + vector.reset(); + vector = Utilities::Trilinos::internal::make_rcp(map); + } + else + vector->putScalar(0); - vector = Utilities::Trilinos::internal::make_rcp( - locally_relevant_map); + IndexSet nonlocal_entries(ghost_entries); + nonlocal_entries.subtract_set(locally_owned_entries); - has_ghost = false; + nonlocal_vector = Utilities::Trilinos::internal::make_rcp( + nonlocal_entries.make_tpetra_map_rcp(communicator, true)); + + compressed = false; + } + + has_ghost = (vector->getMap()->isOneToOne() == false); } @@ -271,9 +311,8 @@ namespace LinearAlgebra V.vector->getMap()); Tpetra::deep_copy(*vector, *V.vector); + compressed = V.compressed; has_ghost = V.has_ghost; - locally_owned_map = V.locally_owned_map; - locally_relevant_map = V.locally_relevant_map; source_stored_elements = V.source_stored_elements; tpetra_comm_pattern = V.tpetra_comm_pattern; } @@ -287,10 +326,18 @@ namespace LinearAlgebra Vector & Vector::operator=(const Number s) { - Assert(s == Number(0.), + (void)s; + Assert(s == Number(0.0), ExcMessage("Only 0 can be assigned to a vector.")); - vector->putScalar(s); + // As checked above, we are only allowed to use d==0.0, so pass + // a constant zero (instead of a run-time value 'd' that *happens* to + // have a zero value) to the underlying class in hopes that the compiler + // can optimize this somehow. + vector->putScalar(/*s=*/0.0); + + if (nonlocal_vector.is_null()) + nonlocal_vector->putScalar(/*s=*/0.0); return *this; } @@ -474,6 +521,7 @@ namespace LinearAlgebra ExcDimensionMismatch(this->size(), V.size())); Assert(vector->getMap()->isSameAs(*V.trilinos_vector().getMap()), ExcDifferentParallelPartitioning()); + Assert(!has_ghost_elements(), ExcGhostsPresent()); return vector->dot(V.trilinos_vector()); } @@ -493,7 +541,7 @@ namespace LinearAlgebra // If the element is not present on the current processor, we can't // continue. This is the main difference to the el() function. - if (local_index == -1) + if (local_index == Teuchos::OrdinalTraits::invalid()) { # if DEAL_II_TRILINOS_VERSION_GTE(14, 0, 0) Assert( @@ -525,6 +573,10 @@ namespace LinearAlgebra { AssertIsFinite(a); + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); + # if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0) auto vector_2d = vector->template getLocalView( Tpetra::Access::ReadWrite); @@ -555,9 +607,14 @@ namespace LinearAlgebra Vector::add(const Number a, const Vector &V) { AssertIsFinite(a); + Assert(vector->getMap()->isSameAs(*(V.trilinos_vector().getMap())), ExcDifferentParallelPartitioning()); + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); + vector->update(a, V.trilinos_vector(), 1.); } @@ -570,12 +627,17 @@ namespace LinearAlgebra const Number b, const Vector &W) { + AssertIsFinite(a); + AssertIsFinite(b); + Assert(vector->getMap()->isSameAs(*(V.trilinos_vector().getMap())), ExcDifferentParallelPartitioning()); Assert(vector->getMap()->isSameAs(*(W.trilinos_vector().getMap())), ExcDifferentParallelPartitioning()); - AssertIsFinite(a); - AssertIsFinite(b); + + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); vector->update(a, V.trilinos_vector(), b, W.trilinos_vector(), 1.); } @@ -588,6 +650,13 @@ namespace LinearAlgebra const Number a, const Vector &V) { + AssertIsFinite(s); + AssertIsFinite(a); + + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); + *this *= s; Vector tmp(V); @@ -605,6 +674,10 @@ namespace LinearAlgebra *(scaling_factors.trilinos_vector().getMap())), ExcDifferentParallelPartitioning()); + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); + vector->elementWiseMultiply(1., *scaling_factors.vector, *vector, 0.); } @@ -614,6 +687,12 @@ namespace LinearAlgebra void Vector::equ(const Number a, const Vector &V) { + AssertIsFinite(a); + + // if we have ghost values, do not allow + // writing to this vector at all. + Assert(!has_ghost_elements(), ExcGhostsPresent()); + // If we don't have the same map, copy. if (vector->getMap()->isSameAs(*V.trilinos_vector().getMap()) == false) this->sadd(0., a, V); @@ -661,6 +740,8 @@ namespace LinearAlgebra Number Vector::mean_value() const { + Assert(!has_ghost_elements(), ExcGhostsPresent()); + return vector->meanValue(); } @@ -670,6 +751,8 @@ namespace LinearAlgebra typename Vector::real_type Vector::l1_norm() const { + Assert(!has_ghost_elements(), ExcGhostsPresent()); + return vector->norm1(); } @@ -679,6 +762,8 @@ namespace LinearAlgebra typename Vector::real_type Vector::l2_norm() const { + Assert(!has_ghost_elements(), ExcGhostsPresent()); + return vector->norm2(); } @@ -688,6 +773,8 @@ namespace LinearAlgebra typename Vector::real_type Vector::linfty_norm() const { + Assert(!has_ghost_elements(), ExcGhostsPresent()); + return vector->normInf(); } @@ -699,6 +786,8 @@ namespace LinearAlgebra const Vector &V, const Vector &W) { + AssertIsFinite(a); + this->add(a, V); return *this * W; @@ -769,13 +858,14 @@ namespace LinearAlgebra Vector::compress(const VectorOperation::values operation) { Assert(has_ghost == false, - ExcMessage("Calling compress() is only useful if a vector " - "has been written into, but this is a vector with ghost " - "elements and consequently is read-only. It does " - "not make sense to call compress() for such " - "vectors.")); - - if (!vector->getMap()->isOneToOne()) + ExcMessage( + "Calling compress() is only useful if a vector " + "has been written into, but this is a vector with ghost " + "elements and consequently is read-only. It does " + "not make sense to call compress() for such " + "vectors.")); + + if (!compressed) { Tpetra::CombineMode tpetra_operation = Tpetra::ZERO; if (operation == VectorOperation::insert) @@ -785,21 +875,11 @@ namespace LinearAlgebra else Assert(false, ExcNotImplemented()); - if (locally_relevant_map.is_null()) - locally_relevant_map = - Teuchos::rcp_const_cast(vector->getMap()); - - Assert(!locally_relevant_map.is_null(), ExcMissingIndexSet()); - Assert(!locally_owned_map.is_null(), ExcMissingIndexSet()); - - VectorType dummy = VectorType(locally_owned_map.getConst()); - Teuchos::RCP exporter = - Tpetra::createExport(locally_relevant_map.getConst(), - locally_owned_map.getConst()); - dummy.doExport(*vector, *exporter, tpetra_operation); + Tpetra::createExport(nonlocal_vector->getMap(), vector->getMap()); + vector->doExport(*nonlocal_vector, *exporter, tpetra_operation); - *vector = dummy; + compressed = true; } }