Correctly handle ghosted PETSc vectors and make ghosted vectors read-only

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

diff --git a/deal.II/include/deal.II/base/exceptions.h b/deal.II/include/deal.II/base/exceptions.h
index 5ef1d87a262b8234129580fd9cb9609e09b378be..fb1d11f7b1567a8289c6947905a3a4f0917fb65a 100644 (file)
--- a/deal.II/include/deal.II/base/exceptions.h
+++ b/deal.II/include/deal.II/base/exceptions.h
@@ -838,6 +838,11 @@ namespace StandardExceptions
    */
   DeclException0 (ExcNeedsNetCDF);
 
+  /**
+   * Parallel vectors with ghost elements are read-only vectors.
+   */
+  DeclException0 (ExcGhostsPresent);
+
   /**
    * A configuration option disabled this feature. In order to use it,
    * you must reconfigure and recompile the libraries.

diff --git a/deal.II/include/deal.II/lac/petsc_parallel_vector.h b/deal.II/include/deal.II/lac/petsc_parallel_vector.h
index e1fab26b3dc248962ff03baaa43b641310e6ec20..52a1399c5a76864a7b3825d03a111af479177fa9 100644 (file)
--- a/deal.II/include/deal.II/lac/petsc_parallel_vector.h
+++ b/deal.II/include/deal.II/lac/petsc_parallel_vector.h
@@ -229,7 +229,7 @@ namespace PETScWrappers
 
 
       /**
-       * Constructs a new parallel PETSc
+       * Constructs a new parallel ghosted PETSc
        * vector from an Indexset. Note that
        * @p local must be contiguous and
        * the global size of the vector is
@@ -245,11 +245,21 @@ namespace PETScWrappers
        * way, the ghost parameter can equal
        * the set of locally relevant
        * degrees of freedom, see step-32.
+       *
+       * @Note: This always creates a ghosted
+       * vector.
        */
       explicit Vector (const MPI_Comm     &communicator,
                        const IndexSet   &local,
-                       const IndexSet &ghost = IndexSet(0));
+                       const IndexSet &ghost);
 
+      /**
+       * Constructs a new parallel PETSc
+       * vector from an Indexset. This creates a non
+       * ghosted vector.
+       */
+      explicit Vector (const MPI_Comm     &communicator,
+                       const IndexSet   &local);
 
       /**
        * Copy the given vector. Resize the

diff --git a/deal.II/include/deal.II/lac/petsc_vector_base.h b/deal.II/include/deal.II/lac/petsc_vector_base.h
index 0588ad8c6337c72a3a28d0ee50310b76ec827db0..07882b876a42bafbb8e861592f5c6d3e5ec68b9c 100644 (file)
--- a/deal.II/include/deal.II/lac/petsc_vector_base.h
+++ b/deal.II/include/deal.II/lac/petsc_vector_base.h
@@ -957,6 +957,8 @@ namespace PETScWrappers
               ExcWrongMode (VectorOperation::insert,
                             vector.last_action));
 
+      Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
+
 #ifdef PETSC_USE_64BIT_INDICES
       const PetscInt petsc_i = index;
 #else
@@ -984,6 +986,8 @@ namespace PETScWrappers
               ExcWrongMode (VectorOperation::add,
                             vector.last_action));
 
+      Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
+
       vector.last_action = VectorOperation::add;
 
       // we have to do above actions in any
@@ -1022,6 +1026,8 @@ namespace PETScWrappers
               ExcWrongMode (VectorOperation::add,
                             vector.last_action));
 
+      Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
+
       vector.last_action = VectorOperation::add;
 
       // we have to do above actions in any
@@ -1061,6 +1067,8 @@ namespace PETScWrappers
               ExcWrongMode (VectorOperation::insert,
                             vector.last_action));
 
+      Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
+
       vector.last_action = VectorOperation::insert;
 
       // we have to do above actions in any
@@ -1101,6 +1109,8 @@ namespace PETScWrappers
               ExcWrongMode (VectorOperation::insert,
                             vector.last_action));
 
+      Assert (!vector.has_ghost_elements(), ExcGhostsPresent());
+
       vector.last_action = VectorOperation::insert;
 
       // we have to do above actions in any

diff --git a/deal.II/source/lac/petsc_parallel_vector.cc b/deal.II/source/lac/petsc_parallel_vector.cc
index e2302f37b8bda60728e6ec90de445361a283c844..fdadfa6ae4957b4291c26505e51d28a3b08fde69 100644 (file)
--- a/deal.II/source/lac/petsc_parallel_vector.cc
+++ b/deal.II/source/lac/petsc_parallel_vector.cc
@@ -85,7 +85,14 @@ namespace PETScWrappers
       Vector::create_vector(local.size(), local.n_elements(), ghost_set);
     }
 
-
+    Vector::Vector (const MPI_Comm     &communicator,
+                     const IndexSet   &local)
+    :
+          communicator (communicator)
+    {
+      Assert(local.is_contiguous(), ExcNotImplemented());
+      Vector::create_vector(local.size(), local.n_elements());
+    }
 
     void
     Vector::reinit (const MPI_Comm    &comm,

diff --git a/deal.II/source/lac/petsc_vector_base.cc b/deal.II/source/lac/petsc_vector_base.cc
index 4e3622e0e1049eca4a9eb662cc21842d603e0387..1102071c14532e25a98ef6a34934e0f5fd4580b4 100644 (file)
--- a/deal.II/source/lac/petsc_vector_base.cc
+++ b/deal.II/source/lac/petsc_vector_base.cc
@@ -567,7 +567,7 @@ namespace PETScWrappers
   VectorBase::normalize () const
   {
     real_type d;
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr = VecNormalize (vector, &d);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -604,6 +604,8 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::abs ()
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
+
     const int ierr = VecAbs (vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -615,6 +617,8 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::conjugate ()
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
+
     const int ierr = VecConjugate (vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -626,6 +630,8 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::mult ()
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
+
     const int ierr = VecPointwiseMult (vector,vector,vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -636,6 +642,7 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::mult (const VectorBase &v)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr = VecPointwiseMult (vector,vector,v);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -647,6 +654,7 @@ namespace PETScWrappers
   VectorBase::mult (const VectorBase &u,
                     const VectorBase &v)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr = VecPointwiseMult (vector,u,v);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -741,7 +749,7 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator *= (const PetscScalar a)
   {
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
 
     const int ierr = VecScale (vector, a);
@@ -755,7 +763,7 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator /= (const PetscScalar a)
   {
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
 
     const PetscScalar factor = 1./a;
@@ -772,6 +780,7 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator += (const VectorBase &v)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr = VecAXPY (vector, 1, v);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -783,6 +792,7 @@ namespace PETScWrappers
   VectorBase &
   VectorBase::operator -= (const VectorBase &v)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr = VecAXPY (vector, -1, v);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
 
@@ -794,6 +804,7 @@ namespace PETScWrappers
   void
   VectorBase::add (const PetscScalar s)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(s), ExcNumberNotFinite());
 
     const int ierr = VecShift (vector, s);
@@ -814,6 +825,7 @@ namespace PETScWrappers
   VectorBase::add (const PetscScalar a,
                    const VectorBase     &v)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
 
     const int ierr = VecAXPY (vector, a, v);
@@ -828,6 +840,7 @@ namespace PETScWrappers
                    const PetscScalar b,
                    const VectorBase &w)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
     Assert (numbers::is_finite(b), ExcNumberNotFinite());
 
@@ -844,6 +857,7 @@ namespace PETScWrappers
   VectorBase::sadd (const PetscScalar s,
                     const VectorBase &v)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(s), ExcNumberNotFinite());
 
     const int ierr = VecAYPX (vector, s, v);
@@ -857,7 +871,7 @@ namespace PETScWrappers
                     const PetscScalar a,
                     const VectorBase     &v)
   {
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(s), ExcNumberNotFinite());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
 
@@ -877,6 +891,7 @@ namespace PETScWrappers
                     const PetscScalar b,
                     const VectorBase     &w)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(s), ExcNumberNotFinite());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
     Assert (numbers::is_finite(b), ExcNumberNotFinite());
@@ -903,7 +918,7 @@ namespace PETScWrappers
                     const PetscScalar c,
                     const VectorBase     &x)
   {
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(s), ExcNumberNotFinite());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
     Assert (numbers::is_finite(b), ExcNumberNotFinite());
@@ -925,6 +940,7 @@ namespace PETScWrappers
   void
   VectorBase::scale (const VectorBase &factors)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr
       = VecPointwiseMult (vector, factors, vector);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
@@ -936,7 +952,7 @@ namespace PETScWrappers
   VectorBase::equ (const PetscScalar a,
                    const VectorBase &v)
   {
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
 
     Assert (size() == v.size(),
@@ -959,7 +975,7 @@ namespace PETScWrappers
                    const PetscScalar b,
                    const VectorBase &w)
   {
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
     Assert (numbers::is_finite(b), ExcNumberNotFinite());
 
@@ -981,6 +997,7 @@ namespace PETScWrappers
   VectorBase::ratio (const VectorBase &a,
                      const VectorBase &b)
   {
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     const int ierr = VecPointwiseDivide (vector, a, b);
     AssertThrow (ierr == 0, ExcPETScError(ierr));
   }
@@ -1091,7 +1108,7 @@ namespace PETScWrappers
             (last_action == ::dealii::VectorOperation::unknown),
             internal::VectorReference::ExcWrongMode (action,
                                                      last_action));
-
+    Assert (!has_ghost_elements(), ExcGhostsPresent());
     // VecSetValues complains if we
     // come with an empty
     // vector. however, it is not a