]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix some bugs in PETSc-deal.II vector copy code.
authorDavid Wells <wellsd2@rpi.edu>
Sat, 6 May 2017 22:53:03 +0000 (18:53 -0400)
committerDavid Wells <wellsd2@rpi.edu>
Mon, 8 May 2017 19:42:30 +0000 (15:42 -0400)
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
tests/petsc/copy_parallel_vector.cc [new file with mode: 0644]
tests/petsc/copy_parallel_vector.with_mpi=true.mpirun=4.output [new file with mode: 0644]

index 73ce84f6e1b2aef0e982da128bd95272dcaf713c..f4db5bc8c81061aedcc632493067a4315f3f28e1 100644 (file)
@@ -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<Vec, std::function<void(Vec *)>> 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 (file)
index 0000000..2be88ed
--- /dev/null
@@ -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 <deal.II/base/index_set.h>
+#include <deal.II/base/mpi.h>
+
+#include <deal.II/lac/petsc_parallel_vector.h>
+#include <deal.II/lac/vector.h>
+
+
+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<double> 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 (file)
index 0000000..cbe78fd
--- /dev/null
@@ -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
+

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.