Work around a bug in PETSc whereby PETSc does not zero-initialize ghost elements...

git-svn-id: https://svn.dealii.org/trunk@33142 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 635384ed356d8c3f027bc12f1f1ea14a371ae90a..5b1d5b9809aaa254b29c263742713c63f23800d2 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -148,6 +148,13 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> Fixed: PETSc up to at least version 3.5 has a bug where it does
+  not zero-initialize the ghost elements of a newly created ghosted
+  parallel vector. This is now worked around inside deal.II.
+  <br>
+  (Wolfgang Bangerth, Michal Wichrowski, 2014/07/12)
+  </li>
+
   <li> Improved: The Trilinos direct solver, TrilinosWrappers::SolverDirect,
   now takes a string to select among the available Amesos solvers. Moreover,
   the solver now also supports deal.II's distributed vectors.

diff --git a/deal.II/source/lac/petsc_parallel_vector.cc b/deal.II/source/lac/petsc_parallel_vector.cc
index ce866b663dc79fdc3ca8de275630ec4b40963ed8..749e4644105af680b0bf4b7e84accdc2eb13d4ee 100644 (file)
--- a/deal.II/source/lac/petsc_parallel_vector.cc
+++ b/deal.II/source/lac/petsc_parallel_vector.cc
@@ -1,7 +1,7 @@
 // ---------------------------------------------------------------------
 // $Id$
 //
-// Copyright (C) 2004 - 2013 by the deal.II authors
+// Copyright (C) 2004 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -332,29 +332,42 @@ namespace PETScWrappers
               ExcDimensionMismatch (size(), n));
 
 #if DEBUG
-      // test ghost allocation in debug mode
-      PetscInt begin, end;
+      {
+       // test ghost allocation in debug mode
+       PetscInt begin, end;
 
-      ierr = VecGetOwnershipRange (vector, &begin, &end);
+       ierr = VecGetOwnershipRange (vector, &begin, &end);
 
-      Assert(local_size==(size_type)(end-begin), ExcInternalError());
+       Assert(local_size==(size_type)(end-begin), ExcInternalError());
 
-      Vec l;
-      ierr = VecGhostGetLocalForm(vector, &l);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
-
-      PetscInt lsize;
-      ierr = VecGetSize(l, &lsize);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
+       Vec l;
+       ierr = VecGhostGetLocalForm(vector, &l);
+       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
-      ierr = VecGhostRestoreLocalForm(vector, &l);
-      AssertThrow (ierr == 0, ExcPETScError(ierr));
+       PetscInt lsize;
+       ierr = VecGetSize(l, &lsize);
+       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
-      Assert( lsize==end-begin+(PetscInt)ghost_indices.n_elements() ,ExcInternalError());
+       ierr = VecGhostRestoreLocalForm(vector, &l);
+       AssertThrow (ierr == 0, ExcPETScError(ierr));
 
+       Assert (lsize==end-begin+(PetscInt)ghost_indices.n_elements(),
+               ExcInternalError());
+      }
 #endif
 
 
+      // in PETSc versions up to 3.5, VecCreateGhost zeroed out the locally
+      // owned vector elements but forgot about the ghost elements. we need to
+      // do this ourselves
+      //
+      // see https://code.google.com/p/dealii/issues/detail?id=233
+#if DEAL_II_PETSC_VERSION_LT(3,6,0)
+      PETScWrappers::MPI::Vector zero;
+      zero.reinit (communicator, this->size(), local_size);
+      *this = zero;
+#endif
+
     }