From 254e961d8beec8438090c906a363e24e950edb48 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Wed, 22 May 2024 14:48:37 -0700 Subject: [PATCH 1/2] Allow creating subvectors only supplying local indices Also add some sanity checking about the parallel type of the subvector --- include/numerics/numeric_vector.h | 3 +- include/numerics/petsc_vector.h | 3 +- include/numerics/trilinos_epetra_vector.h | 3 -- src/numerics/petsc_vector.C | 51 +++++++++++++++++------ src/numerics/trilinos_epetra_vector.C | 9 ---- 5 files changed, 42 insertions(+), 27 deletions(-) diff --git a/include/numerics/numeric_vector.h b/include/numerics/numeric_vector.h index 12f19d7fdf2..c89efa28887 100644 --- a/include/numerics/numeric_vector.h +++ b/include/numerics/numeric_vector.h @@ -723,7 +723,8 @@ class NumericVector : public ReferenceCountedObject>, * PetscVectors. */ virtual void create_subvector(NumericVector & , - const std::vector &) const + const std::vector &, + bool = true) const { libmesh_not_implemented(); } diff --git a/include/numerics/petsc_vector.h b/include/numerics/petsc_vector.h index e81596fe276..3b01e3d02d4 100644 --- a/include/numerics/petsc_vector.h +++ b/include/numerics/petsc_vector.h @@ -327,7 +327,8 @@ class PetscVector final : public NumericVector virtual void print_matlab(const std::string & name = "") const override; virtual void create_subvector(NumericVector & subvector, - const std::vector & rows) const override; + const std::vector & rows, + bool supplying_global_rows = true) const override; virtual void swap (NumericVector & v) override; diff --git a/include/numerics/trilinos_epetra_vector.h b/include/numerics/trilinos_epetra_vector.h index 1a0f74ba253..e3af168f5c9 100644 --- a/include/numerics/trilinos_epetra_vector.h +++ b/include/numerics/trilinos_epetra_vector.h @@ -253,9 +253,6 @@ class EpetraVector final : public NumericVector virtual void pointwise_mult (const NumericVector & vec1, const NumericVector & vec2) override; - virtual void create_subvector (NumericVector & subvector, - const std::vector & rows) const override; - virtual void swap (NumericVector & v) override; virtual std::size_t max_allowed_id() const override; diff --git a/src/numerics/petsc_vector.C b/src/numerics/petsc_vector.C index 752d94e744c..d24c31ca8bf 100644 --- a/src/numerics/petsc_vector.C +++ b/src/numerics/petsc_vector.C @@ -1187,10 +1187,15 @@ void PetscVector::print_matlab (const std::string & name) const template void PetscVector::create_subvector(NumericVector & subvector, - const std::vector & rows) const + const std::vector & rows, + const bool supplying_global_rows) const { parallel_object_only(); + libmesh_error_msg_if( + subvector.type() == GHOSTED, + "We do not support scattering parallel information to ghosts for subvectors"); + this->_restore_array(); // PETSc data structures @@ -1204,21 +1209,32 @@ void PetscVector::create_subvector(NumericVector & subvector, // If not, we use the appropriate PETSc routines to initialize it. if (!petsc_subvector->initialized()) { - // Initialize the petsc_subvector to have enough space to hold - // the entries which will be scattered into it. Note: such an - // init() function (where we let PETSc decide the number of local - // entries) is not currently offered by the PetscVector - // class. Should we differentiate here between sequential and - // parallel vector creation based on this->n_processors() ? - ierr = VecCreateMPI(this->comm().get(), - PETSC_DECIDE, // n_local - cast_int(rows.size()), // n_global - &(petsc_subvector->_vec)); + libmesh_assert(petsc_subvector->_type == AUTOMATIC || petsc_subvector->_type == PARALLEL); + + if (supplying_global_rows) + // Initialize the petsc_subvector to have enough space to hold + // the entries which will be scattered into it. Note: such an + // init() function (where we let PETSc decide the number of local + // entries) is not currently offered by the PetscVector + // class. Should we differentiate here between sequential and + // parallel vector creation based on this->n_processors() ? + ierr = VecCreateMPI(this->comm().get(), + PETSC_DECIDE, // n_local + cast_int(rows.size()), // n_global + &(petsc_subvector->_vec)); + else + ierr = VecCreateMPI(this->comm().get(), + cast_int(rows.size()), + PETSC_DETERMINE, + &(petsc_subvector->_vec)); LIBMESH_CHKERR(ierr); ierr = VecSetFromOptions (petsc_subvector->_vec); LIBMESH_CHKERR(ierr); + // We created a parallel vector + petsc_subvector->_type = PARALLEL; + // Mark the subvector as initialized petsc_subvector->_is_initialized = true; } @@ -1227,9 +1243,16 @@ void PetscVector::create_subvector(NumericVector & subvector, petsc_subvector->_restore_array(); } - // Use iota to fill an array with entries [0,1,2,3,4,...rows.size()] std::vector idx(rows.size()); - std::iota (idx.begin(), idx.end(), 0); + if (supplying_global_rows) + std::iota (idx.begin(), idx.end(), 0); + else + { + PetscInt start; + ierr = VecGetOwnershipRange(petsc_subvector->_vec, &start, nullptr); + LIBMESH_CHKERR(ierr); + std::iota (idx.begin(), idx.end(), start); + } // Construct index sets WrappedPetsc parent_is; @@ -1258,6 +1281,8 @@ void PetscVector::create_subvector(NumericVector & subvector, // Actually perform the scatter VecScatterBeginEnd(this->comm(), scatter, this->_vec, petsc_subvector->_vec, INSERT_VALUES, SCATTER_FORWARD); + + petsc_subvector->_is_closed = true; } diff --git a/src/numerics/trilinos_epetra_vector.C b/src/numerics/trilinos_epetra_vector.C index a2e797ac1f9..187ae881cc9 100644 --- a/src/numerics/trilinos_epetra_vector.C +++ b/src/numerics/trilinos_epetra_vector.C @@ -576,15 +576,6 @@ void EpetraVector::localize_to_one (std::vector & v_local, } - -template -void EpetraVector::create_subvector(NumericVector & /* subvector */, - const std::vector & /* rows */) const -{ - libmesh_not_implemented(); -} - - /********************************************************************* * The following were copied (and slightly modified) from * Epetra_FEVector.h in order to allow us to use a standard From 6bdf108aa9e37a27f74f7ffebbdb4f2b3a9bca29 Mon Sep 17 00:00:00 2001 From: Alex Lindsay Date: Thu, 23 May 2024 09:21:42 -0700 Subject: [PATCH 2/2] Add doxygen for boolean param --- include/numerics/numeric_vector.h | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/include/numerics/numeric_vector.h b/include/numerics/numeric_vector.h index c89efa28887..0ac3de4cec7 100644 --- a/include/numerics/numeric_vector.h +++ b/include/numerics/numeric_vector.h @@ -717,10 +717,13 @@ class NumericVector : public ReferenceCountedObject>, } /** - * Fills in \p subvector from this vector using the indices in \p - * rows. Similar to the \p create_submatrix() routine for the - * SparseMatrix class, it is currently only implemented for - * PetscVectors. + * Fills in \p subvector from this vector using the indices in \p rows. + * Similar to the \p create_submatrix() routine for the SparseMatrix class, it + * is currently only implemented for PetscVectors. The boolean parameter + * communicates whether the supplied vector of rows corresponds to all the + * rows that should be used in the subvector's index set, e.g. whether the + * rows correspond to the global collective. If the rows supplied are only the + * local indices, then the boolean parameter should be set to false */ virtual void create_subvector(NumericVector & , const std::vector &,