From b20b75c89a37b784f95d47676d94d23b5c43d207 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 6 May 2017 18:53:03 -0400 Subject: [PATCH] Fix some bugs in PETSc-deal.II vector copy code. This new, correct implementation is based on the original, working implementation present before commit 3aa64f983d4 in the (now removed) PETSc serial vector class. --- include/deal.II/lac/vector.templates.h | 44 +++++---- tests/petsc/copy_parallel_vector.cc | 59 +++++++++++ ...allel_vector.with_mpi=true.mpirun=4.output | 99 +++++++++++++++++++ 3 files changed, 182 insertions(+), 20 deletions(-) create mode 100644 tests/petsc/copy_parallel_vector.cc create mode 100644 tests/petsc/copy_parallel_vector.with_mpi=true.mpirun=4.output diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 73ce84f6e1..f4db5bc8c8 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -104,40 +104,44 @@ namespace internal { // Create a sequential PETSc vector and then copy over the entries into // the deal.II vector. - // - // Wrap the sequential vector with a custom deleter. - std::unique_ptr> sequential_vector - (new Vec(), [](Vec *ptr) - { -#if DEAL_II_PETSC_VERSION_LT(3,2,0) - const PetscErrorCode ierr = VecDestroy (*ptr); -#else - const PetscErrorCode ierr = VecDestroy (ptr); -#endif - (void)ierr; - AssertNothrow (ierr == 0, ExcPETScError(ierr)); - }); + Vec sequential_vector; + VecScatter scatter_context; + + PetscErrorCode ierr = VecScatterCreateToAll(v, &scatter_context, &sequential_vector); + AssertThrow (ierr == 0, ExcPETScError(ierr)); - PetscErrorCode ierr = VecCreateSeq (v.get_mpi_communicator(), - v.size(), - sequential_vector.get()); + ierr = VecScatterBegin(scatter_context, v, sequential_vector, INSERT_VALUES, + SCATTER_FORWARD); AssertThrow (ierr == 0, ExcPETScError(ierr)); - ierr = VecCopy (v, *sequential_vector); + ierr = VecScatterEnd(scatter_context, v, sequential_vector, INSERT_VALUES, + SCATTER_FORWARD); AssertThrow (ierr == 0, ExcPETScError(ierr)); PetscScalar *start_ptr; - ierr = VecGetArray(*sequential_vector, &start_ptr); + ierr = VecGetArray(sequential_vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - const auto v_size = v.size(); + const PETScWrappers::VectorBase::size_type v_size = v.size(); if (out.size() != v_size) out.reinit (v_size, true); internal::VectorOperations::copy (start_ptr, start_ptr + out.size(), out.begin()); - ierr = VecRestoreArray (*sequential_vector, &start_ptr); + ierr = VecRestoreArray (sequential_vector, &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); + +#if DEAL_II_PETSC_VERSION_LT(3,2,0) + ierr = VecScatterDestroy(scatter_context); + AssertNothrow (ierr == 0, ExcPETScError(ierr)); + ierr = VecDestroy (sequential_vector); + AssertNothrow (ierr == 0, ExcPETScError(ierr)); +#else + ierr = VecScatterDestroy(&scatter_context); + AssertNothrow (ierr == 0, ExcPETScError(ierr)); + ierr = VecDestroy (&sequential_vector); + AssertNothrow (ierr == 0, ExcPETScError(ierr)); +#endif } } diff --git a/tests/petsc/copy_parallel_vector.cc b/tests/petsc/copy_parallel_vector.cc new file mode 100644 index 0000000000..2be88ede28 --- /dev/null +++ b/tests/petsc/copy_parallel_vector.cc @@ -0,0 +1,59 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Check that we can correctly copy a distributed PETSc parallel vector onto a +// single deal.II vector. + +#include "../tests.h" + +#include +#include + +#include +#include + + +int main(int argc, char **argv) +{ + using namespace dealii; + + Utilities::MPI::MPI_InitFinalize mpi_context(argc, argv, 1); + MPILogInitAll mpi_log; + + const auto rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD); + const auto n_mpi_processes = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD); + + constexpr types::global_dof_index dofs_per_process = 5; + IndexSet owned_indices(dofs_per_process*n_mpi_processes); + owned_indices.add_range(dofs_per_process*rank, dofs_per_process*(rank + 1)); + owned_indices.compress(); + + deallog << "rank is: " << rank << std::endl; + owned_indices.print(deallog); + + PETScWrappers::MPI::Vector petsc_vector(owned_indices, MPI_COMM_WORLD); + for (const auto dof : owned_indices) + { + petsc_vector(dof) = double(dof); + } + + Vector deal_vector(petsc_vector); + deallog << rank << ": deal vector size: " << deal_vector.size() << std::endl; + + for (const auto value : deal_vector) + { + deallog << value << std::endl; + } +} diff --git a/tests/petsc/copy_parallel_vector.with_mpi=true.mpirun=4.output b/tests/petsc/copy_parallel_vector.with_mpi=true.mpirun=4.output new file mode 100644 index 0000000000..cbe78fdd8c --- /dev/null +++ b/tests/petsc/copy_parallel_vector.with_mpi=true.mpirun=4.output @@ -0,0 +1,99 @@ + +DEAL:0::rank is: 0 +DEAL:0::{[0,4]} +DEAL:0::0: deal vector size: 20 +DEAL:0::0.00000 +DEAL:0::1.00000 +DEAL:0::2.00000 +DEAL:0::3.00000 +DEAL:0::4.00000 +DEAL:0::5.00000 +DEAL:0::6.00000 +DEAL:0::7.00000 +DEAL:0::8.00000 +DEAL:0::9.00000 +DEAL:0::10.0000 +DEAL:0::11.0000 +DEAL:0::12.0000 +DEAL:0::13.0000 +DEAL:0::14.0000 +DEAL:0::15.0000 +DEAL:0::16.0000 +DEAL:0::17.0000 +DEAL:0::18.0000 +DEAL:0::19.0000 + +DEAL:1::rank is: 1 +DEAL:1::{[5,9]} +DEAL:1::1: deal vector size: 20 +DEAL:1::0.00000 +DEAL:1::1.00000 +DEAL:1::2.00000 +DEAL:1::3.00000 +DEAL:1::4.00000 +DEAL:1::5.00000 +DEAL:1::6.00000 +DEAL:1::7.00000 +DEAL:1::8.00000 +DEAL:1::9.00000 +DEAL:1::10.0000 +DEAL:1::11.0000 +DEAL:1::12.0000 +DEAL:1::13.0000 +DEAL:1::14.0000 +DEAL:1::15.0000 +DEAL:1::16.0000 +DEAL:1::17.0000 +DEAL:1::18.0000 +DEAL:1::19.0000 + + +DEAL:2::rank is: 2 +DEAL:2::{[10,14]} +DEAL:2::2: deal vector size: 20 +DEAL:2::0.00000 +DEAL:2::1.00000 +DEAL:2::2.00000 +DEAL:2::3.00000 +DEAL:2::4.00000 +DEAL:2::5.00000 +DEAL:2::6.00000 +DEAL:2::7.00000 +DEAL:2::8.00000 +DEAL:2::9.00000 +DEAL:2::10.0000 +DEAL:2::11.0000 +DEAL:2::12.0000 +DEAL:2::13.0000 +DEAL:2::14.0000 +DEAL:2::15.0000 +DEAL:2::16.0000 +DEAL:2::17.0000 +DEAL:2::18.0000 +DEAL:2::19.0000 + + +DEAL:3::rank is: 3 +DEAL:3::{[15,19]} +DEAL:3::3: deal vector size: 20 +DEAL:3::0.00000 +DEAL:3::1.00000 +DEAL:3::2.00000 +DEAL:3::3.00000 +DEAL:3::4.00000 +DEAL:3::5.00000 +DEAL:3::6.00000 +DEAL:3::7.00000 +DEAL:3::8.00000 +DEAL:3::9.00000 +DEAL:3::10.0000 +DEAL:3::11.0000 +DEAL:3::12.0000 +DEAL:3::13.0000 +DEAL:3::14.0000 +DEAL:3::15.0000 +DEAL:3::16.0000 +DEAL:3::17.0000 +DEAL:3::18.0000 +DEAL:3::19.0000 + -- 2.39.5