From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sat, 12 Jul 2014 12:57:51 +0000 (+0000)
Subject: Work around a bug in PETSc whereby PETSc does not zero-initialize ghost elements... 
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=b496a19fdde8945c1a6f4d0fb5f3c80dac9ef6b5;p=dealii-svn.git

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

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 635384ed35..5b1d5b9809 100644
--- 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 ce866b663d..749e464410 100644
--- 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
+
     }