diff --git a/include/numerics/numeric_vector.h b/include/numerics/numeric_vector.h index fb90bc33706..0eea10e1912 100644 --- a/include/numerics/numeric_vector.h +++ b/include/numerics/numeric_vector.h @@ -139,7 +139,8 @@ class NumericVector : public ReferenceCountedObject>, */ static std::unique_ptr> build(const Parallel::Communicator & comm, - const SolverPackage solver_package = libMesh::default_solver_package()); + SolverPackage solver_package = libMesh::default_solver_package(), + ParallelType parallel_type = AUTOMATIC); /** * \returns \p true if the vector has been initialized, @@ -154,8 +155,25 @@ class NumericVector : public ReferenceCountedObject>, /** * \returns The type (SERIAL, PARALLEL, GHOSTED) of the vector. + * + * \deprecated because it is dangerous to change the ParallelType + * of an already-initialized NumericVector. See NumericVector::set_type() + * for a safer, non-deprecated setter. */ +#ifdef LIBMESH_ENABLE_DEPRECATED ParallelType & type() { return _type; } +#endif + + /** + * Allow the user to change the ParallelType of the NumericVector + * under some circumstances. If the NumericVector has not been + * initialized yet, then it is generally safe to change the + * ParallelType. otherwise, if the NumericVector has already been + * initialized with a specific type, we cannot change it without + * doing some extra copying/reinitialization work, and we currently + * throw an error if this is requested. + */ + void set_type(ParallelType t); /** * \returns \p true if the vector is closed and ready for diff --git a/src/numerics/numeric_vector.C b/src/numerics/numeric_vector.C index 4fd45749ad1..8c95827aa4e 100644 --- a/src/numerics/numeric_vector.C +++ b/src/numerics/numeric_vector.C @@ -47,7 +47,9 @@ namespace libMesh // Full specialization for Real datatypes template std::unique_ptr> -NumericVector::build(const Parallel::Communicator & comm, const SolverPackage solver_package) +NumericVector::build(const Parallel::Communicator & comm, + SolverPackage solver_package, + ParallelType parallel_type) { // Build the appropriate vector switch (solver_package) @@ -55,31 +57,62 @@ NumericVector::build(const Parallel::Communicator & comm, const SolverPackage #ifdef LIBMESH_HAVE_LASPACK case LASPACK_SOLVERS: - return std::make_unique>(comm, AUTOMATIC); + return std::make_unique>(comm, parallel_type); #endif #ifdef LIBMESH_HAVE_PETSC case PETSC_SOLVERS: - return std::make_unique>(comm, AUTOMATIC); + return std::make_unique>(comm, parallel_type); #endif #ifdef LIBMESH_TRILINOS_HAVE_EPETRA case TRILINOS_SOLVERS: - return std::make_unique>(comm, AUTOMATIC); + return std::make_unique>(comm, parallel_type); #endif #ifdef LIBMESH_HAVE_EIGEN case EIGEN_SOLVERS: - return std::make_unique>(comm, AUTOMATIC); + return std::make_unique>(comm, parallel_type); #endif default: - return std::make_unique>(comm, AUTOMATIC); + return std::make_unique>(comm, parallel_type); } } +template +void NumericVector::set_type(ParallelType t) +{ + // Check for no-op + if (_type == t) + return; + + // If the NumericVector is not yet initialized, then it is generally + // safe to change the ParallelType, with minor restrictions. + if (!this->initialized()) + { + // If ghosted vectors are not enabled and the user requested a + // GHOSTED vector, fall back on SERIAL. +#ifndef LIBMESH_ENABLE_GHOSTED + if (t == GHOSTED) + { + _type = SERIAL; + return; + } +#endif + + _type = t; + return; + } + + // If we made it here, then the NumericVector was already + // initialized and we don't currently allow the ParallelType to be + // changed, although this could potentially be added later. + libmesh_not_implemented(); +} + template void NumericVector::insert (const T * v, const std::vector & dof_indices) diff --git a/src/systems/system.C b/src/systems/system.C index 5c4fa61021a..dabcdbf06b2 100644 --- a/src/systems/system.C +++ b/src/systems/system.C @@ -792,7 +792,9 @@ NumericVector & System::add_vector (std::string_view vec_name, vec.swap(*new_vec); } else - vec.type() = type; + // The PARALLEL vec is not yet initialized, so we can + // just "upgrade" it to GHOSTED. + vec.set_type(type); } } @@ -800,11 +802,19 @@ NumericVector & System::add_vector (std::string_view vec_name, return vec; } - // Otherwise build the vector - auto pr = _vectors.emplace(vec_name, NumericVector::build(this->comm())); + // Otherwise, build the vector. The following emplace() is + // guaranteed to succeed because, if we made it here, we don't + // already have a vector named "vec_name". We pass the user's + // requested ParallelType directly to NumericVector::build() so + // that, even if the vector is not initialized now, it will get the + // right type when it is initialized later. + auto pr = + _vectors.emplace(vec_name, + NumericVector::build(this->comm(), + libMesh::default_solver_package(), + type)); auto buf = pr.first->second.get(); _vector_projections.emplace(vec_name, projections); - buf->type() = type; // Vectors are primal by default _vector_is_adjoint.emplace(vec_name, -1);