From bec7fe052bd2a9e8858c24363ba6777fb41aab37 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Wed, 19 Apr 2023 13:50:44 -0600 Subject: [PATCH] Address the pesky 'MPI_Comm as void*' problem with the SUNDIALS interfaces. --- include/deal.II/sundials/n_vector.templates.h | 145 ++++++++++++------ 1 file changed, 99 insertions(+), 46 deletions(-) diff --git a/include/deal.II/sundials/n_vector.templates.h b/include/deal.II/sundials/n_vector.templates.h index a8cb07b451..535c2f41ac 100644 --- a/include/deal.II/sundials/n_vector.templates.h +++ b/include/deal.II/sundials/n_vector.templates.h @@ -92,6 +92,21 @@ namespace SUNDIALS const VectorType * get() const; + /** + * Return a reference to a copy of the communicator the vector uses. + * This function exists because the N_Vector + * interface requires a function that returns a `void*` pointing + * to the communicator object -- so somewhere, we need to have an + * address to point to. The issue is that our vectors typically + * return a *copy* of the communicator, rather than a reference to + * the communicator they use, and so there is only a temporary + * object and no address we can point to. To work around this + * requirement, this class stores a copy of the communicator, + * and this function here returns a reference to this copy. + */ + const MPI_Comm & + get_mpi_communicator() const; + private: using PointerType = std::unique_ptr>; @@ -107,6 +122,18 @@ namespace SUNDIALS */ PointerType vector; + /** + * A copy of the communicator the vector uses, initialized in the + * constructor of this class. We store this because the N_Vector + * interface requires a function that returns a `void*` pointing + * to the communicator object -- so somewhere, we need to have an + * address to point to. The issue is that our vectors typically + * return a *copy* of the communicator, rather than a reference to + * the communicator they use, and so there is only a temporary + * object and no address we can point to. + */ + MPI_Comm comm; + /** * Flag storing whether the stored pointer is to be treated as const. If * the pointer passed in the constructor was indeed const, it is cast away @@ -261,13 +288,7 @@ namespace SUNDIALS void add_constant(N_Vector x, realtype b, N_Vector z); - template ::value, int> = 0> - const MPI_Comm & - get_communicator(N_Vector v); - - template ::value, int> = 0> + template const MPI_Comm & get_communicator(N_Vector v); @@ -275,16 +296,9 @@ namespace SUNDIALS * Sundials likes a void* but we want to use the above functions * internally with a safe type. */ - template ::value, int> = 0> - inline void * - get_communicator_as_void_ptr(N_Vector v); - - template ::value, int> = 0> + template inline void * get_communicator_as_void_ptr(N_Vector v); - } // namespace NVectorOperations } // namespace internal } // namespace SUNDIALS @@ -297,18 +311,65 @@ namespace SUNDIALS { namespace internal { + namespace + { + template ::value, int> = 0> + MPI_Comm + get_mpi_communicator_from_vector(const VectorType &) + { + return MPI_COMM_SELF; + } + + + + template ::value && + !IsBlockVector::value, + int> = 0> + MPI_Comm + get_mpi_communicator_from_vector(const VectorType &v) + { +# ifndef DEAL_II_WITH_MPI + (void)v; + return MPI_COMM_SELF; +# else + return v.get_mpi_communicator(); +# endif + } + + + + template ::value && + IsBlockVector::value, + int> = 0> + MPI_Comm + get_mpi_communicator_from_vector(const VectorType &v) + { +# ifndef DEAL_II_WITH_MPI + (void)v; + return MPI_COMM_SELF; +# else + return v.block(0).get_mpi_communicator(); +# endif + } + } // namespace + + template NVectorContent::NVectorContent() : vector(typename VectorMemory::Pointer(mem)) + , comm(get_mpi_communicator_from_vector(*vector)) , is_const(false) {} - template NVectorContent::NVectorContent(VectorType *vector) : vector(vector, [](VectorType *) { /* not owning memory -> don't free*/ }) + , comm(get_mpi_communicator_from_vector(*vector)) , is_const(false) {} @@ -318,6 +379,7 @@ namespace SUNDIALS NVectorContent::NVectorContent(const VectorType *vector) : vector(const_cast(vector), [](VectorType *) { /* not owning memory -> don't free*/ }) + , comm(get_mpi_communicator_from_vector(*vector)) , is_const(true) {} @@ -349,6 +411,14 @@ namespace SUNDIALS + template + const MPI_Comm & + NVectorContent::get_mpi_communicator() const + { + return comm; + } + + # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0) template NVectorView @@ -528,41 +598,20 @@ namespace SUNDIALS - template ::value, int>> - const MPI_Comm & - get_communicator(N_Vector v) - { - return unwrap_nvector_const(v) - ->block(0) - .get_mpi_communicator(); - } - - - - template ::value, int>> + template const MPI_Comm & get_communicator(N_Vector v) { - return unwrap_nvector_const(v)->get_mpi_communicator(); - } - - - - template ::value, int>> - void *get_communicator_as_void_ptr(N_Vector) - { - // required by SUNDIALS: MPI-unaware vectors should return the nullptr - // as comm - return nullptr; + Assert(v != nullptr, ExcInternalError()); + Assert(v->content != nullptr, ExcInternalError()); + auto *pContent = + reinterpret_cast *>(v->content); + return pContent->get_mpi_communicator(); } - template ::value, int>> + template void * get_communicator_as_void_ptr(N_Vector v) { @@ -570,8 +619,12 @@ namespace SUNDIALS (void)v; return nullptr; # else - // We need to cast away const here, as SUNDIALS demands a pure `void *`. - return &(const_cast(get_communicator(v))); + if (is_serial_vector::value == false) + // We need to cast away const here, as SUNDIALS demands a pure + // `void*`. + return &(const_cast(get_communicator(v))); + else + return nullptr; # endif } -- 2.39.5