From: heister Date: Sat, 16 Feb 2013 21:03:23 +0000 (+0000) Subject: Correctly handle ghosted PETSc vectors and make ghosted vectors read-only X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e17e71209965c6e86c7092bab5d0edb73dacf179;p=dealii-svn.git 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 5ef1d87a26..fb1d11f7b1 100644 --- 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 e1fab26b3d..52a1399c5a 100644 --- 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 0588ad8c63..07882b876a 100644 --- 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 e2302f37b8..fdadfa6ae4 100644 --- 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 4e3622e0e1..1102071c14 100644 --- 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