Merge pull request #3857 from lindsayad/subvector-from-local-indices

Allow creating subvectors only supplying local indices
libMesh · May 24, 2024 · 6cd7a31 · 6cd7a31
2 parents e7bc93d + 6bdf108
commit 6cd7a31
Show file tree

Hide file tree

Showing 5 changed files with 49 additions and 31 deletions.
diff --git a/include/numerics/numeric_vector.h b/include/numerics/numeric_vector.h
@@ -717,13 +717,17 @@ class NumericVector : public ReferenceCountedObject<NumericVector<T>>,
   }
 
   /**
-   * 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<T> & ,
-                                const std::vector<numeric_index_type> &) const
+                                const std::vector<numeric_index_type> &,
+                                bool = true) const
   {
     libmesh_not_implemented();
   }

diff --git a/include/numerics/petsc_vector.h b/include/numerics/petsc_vector.h
@@ -327,7 +327,8 @@ class PetscVector final : public NumericVector<T>
   virtual void print_matlab(const std::string & name = "") const override;
 
   virtual void create_subvector(NumericVector<T> & subvector,
-                                const std::vector<numeric_index_type> & rows) const override;
+                                const std::vector<numeric_index_type> & rows,
+                                bool supplying_global_rows = true) const override;
 
   virtual void swap (NumericVector<T> & v) override;
 

diff --git a/include/numerics/trilinos_epetra_vector.h b/include/numerics/trilinos_epetra_vector.h
@@ -253,9 +253,6 @@ class EpetraVector final : public NumericVector<T>
   virtual void pointwise_mult (const NumericVector<T> & vec1,
                                const NumericVector<T> & vec2) override;
 
-  virtual void create_subvector (NumericVector<T> & subvector,
-                                 const std::vector<numeric_index_type> & rows) const override;
-
   virtual void swap (NumericVector<T> & 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
@@ -1187,10 +1187,15 @@ void PetscVector<T>::print_matlab (const std::string & name) const
 
 template <typename T>
 void PetscVector<T>::create_subvector(NumericVector<T> & subvector,
-                                      const std::vector<numeric_index_type> & rows) const
+                                      const std::vector<numeric_index_type> & 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<T>::create_subvector(NumericVector<T> & 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<PetscInt>(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<PetscInt>(rows.size()), // n_global
+                            &(petsc_subvector->_vec));
+      else
+        ierr = VecCreateMPI(this->comm().get(),
+                            cast_int<PetscInt>(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<T>::create_subvector(NumericVector<T> & subvector,
       petsc_subvector->_restore_array();
     }
 
-  // Use iota to fill an array with entries [0,1,2,3,4,...rows.size()]
   std::vector<PetscInt> 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<IS> parent_is;
@@ -1258,6 +1281,8 @@ void PetscVector<T>::create_subvector(NumericVector<T> & 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
@@ -576,15 +576,6 @@ void EpetraVector<T>::localize_to_one (std::vector<T> &  v_local,
 }
 
 
-
-template <typename T>
-void EpetraVector<T>::create_subvector(NumericVector<T> & /* subvector */,
-                                       const std::vector<numeric_index_type> & /* 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