From: Stefano Zampini Date: Sat, 28 Jan 2023 08:47:25 +0000 (+0300) Subject: Fix PETSc ghost vectors X-Git-Tag: v9.5.0-rc1~369^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=bee31456661b8984f8ce20df26cf9973ed0c6b3f;p=dealii.git Fix PETSc ghost vectors --- diff --git a/include/deal.II/lac/petsc_vector_base.h b/include/deal.II/lac/petsc_vector_base.h index 4a1cf5bea1..a5ccbd0d46 100644 --- a/include/deal.II/lac/petsc_vector_base.h +++ b/include/deal.II/lac/petsc_vector_base.h @@ -426,10 +426,13 @@ namespace PETScWrappers has_ghost_elements() const; /** - * This function only exists for compatibility with the @p - * LinearAlgebra::distributed::Vector class and does nothing: this class - * implements ghost value updates in a different way that is a better fit - * with the underlying PETSc vector object. + * Return the IndexSet of ghost elements. + */ + const IndexSet & + ghost_elements() const; + + /** + * Update ghosted elements. */ void update_ghost_values() const; @@ -847,6 +850,12 @@ namespace PETScWrappers const size_type * indices, const PetscScalar *values, const bool add_values); + + /** + * Determine ghost indices from the internal PETSc Vec + */ + void + determine_ghost_indices(); }; @@ -1116,10 +1125,26 @@ namespace PETScWrappers } + inline const IndexSet & + VectorBase::ghost_elements() const + { + return ghost_indices; + } + inline void VectorBase::update_ghost_values() const - {} + { + if (ghosted) + { + PetscErrorCode ierr; + + ierr = VecGhostUpdateBegin(vector, INSERT_VALUES, SCATTER_FORWARD); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGhostUpdateEnd(vector, INSERT_VALUES, SCATTER_FORWARD); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } + } diff --git a/source/lac/petsc_vector_base.cc b/source/lac/petsc_vector_base.cc index 1160f17aeb..2d28ec1d9d 100644 --- a/source/lac/petsc_vector_base.cc +++ b/source/lac/petsc_vector_base.cc @@ -149,11 +149,11 @@ namespace PETScWrappers , ghosted(false) , last_action(VectorOperation::unknown) { - /* TODO GHOSTED */ const PetscErrorCode ierr = PetscObjectReference(reinterpret_cast(vector)); AssertNothrow(ierr == 0, ExcPETScError(ierr)); (void)ierr; + this->determine_ghost_indices(); } @@ -170,7 +170,6 @@ namespace PETScWrappers void VectorBase::reinit(Vec v) { - /* TODO GHOSTED */ AssertThrow(last_action == VectorOperation::unknown, ExcMessage("Cannot assign a new Vec")); PetscErrorCode ierr = @@ -179,8 +178,75 @@ namespace PETScWrappers ierr = VecDestroy(&vector); AssertThrow(ierr == 0, ExcPETScError(ierr)); vector = v; + this->determine_ghost_indices(); } + void + VectorBase::determine_ghost_indices() + { + // Reset ghost data + ghosted = false; + ghost_indices.clear(); + + // There's no API to infer ghost indices from a PETSc Vec + PetscErrorCode ierr; + Vec ghosted_vec; + ierr = VecGhostGetLocalForm(vector, &ghosted_vec); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + if (ghosted_vec && ghosted_vec != vector) + { + Vec tvector; + PetscScalar *array; + PetscInt st, en, N, ln; + + ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + + ierr = VecGetSize(vector, &N); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGetOwnershipRange(vector, &st, &en); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecDuplicate(vector, &tvector); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGetArray(tvector, &array); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + for (PetscInt i = 0; i < en - st; i++) + array[i] = st + i; + ierr = VecRestoreArray(tvector, &array); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGhostUpdateBegin(tvector, INSERT_VALUES, SCATTER_FORWARD); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGhostUpdateEnd(tvector, INSERT_VALUES, SCATTER_FORWARD); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGhostGetLocalForm(tvector, &ghosted_vec); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGetLocalSize(ghosted_vec, &ln); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGetArrayRead(ghosted_vec, (const PetscScalar **)&array); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + + // Populate ghosted and ghost_indices + ghosted = true; + ghost_indices.set_size(N); + for (PetscInt i = en - st; i < ln; i++) + ghost_indices.add_index(static_cast(array[i])); + ghost_indices.compress(); + + ierr = VecRestoreArrayRead(ghosted_vec, (const PetscScalar **)&array); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecGhostRestoreLocalForm(tvector, &ghosted_vec); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ierr = VecDestroy(&tvector); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } + else + { + ierr = VecGhostRestoreLocalForm(vector, &ghosted_vec); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + } + } + + void VectorBase::clear() { @@ -443,6 +509,9 @@ namespace PETScWrappers // indicate that we're back to a // pristine state last_action = VectorOperation::unknown; + + // update ghost values if needed + update_ghost_values(); } @@ -710,6 +779,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecScale(vector, a); AssertThrow(ierr == 0, ExcPETScError(ierr)); + // update ghost values if needed + update_ghost_values(); + return *this; } @@ -727,6 +799,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecScale(vector, factor); AssertThrow(ierr == 0, ExcPETScError(ierr)); + // update ghost values if needed + update_ghost_values(); + return *this; } @@ -739,6 +814,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecAXPY(vector, 1, v); AssertThrow(ierr == 0, ExcPETScError(ierr)); + // update ghost values if needed + update_ghost_values(); + return *this; } @@ -751,6 +829,8 @@ namespace PETScWrappers const PetscErrorCode ierr = VecAXPY(vector, -1, v); AssertThrow(ierr == 0, ExcPETScError(ierr)); + // update ghost values if needed + update_ghost_values(); return *this; } @@ -764,6 +844,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecShift(vector, s); AssertThrow(ierr == 0, ExcPETScError(ierr)); + + // update ghost values if needed + update_ghost_values(); } @@ -776,6 +859,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecAXPY(vector, a, v); AssertThrow(ierr == 0, ExcPETScError(ierr)); + + // update ghost values if needed + update_ghost_values(); } @@ -795,6 +881,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecMAXPY(vector, 2, weights, addends); AssertThrow(ierr == 0, ExcPETScError(ierr)); + + // update ghost values if needed + update_ghost_values(); } @@ -807,6 +896,9 @@ namespace PETScWrappers const PetscErrorCode ierr = VecAYPX(vector, s, v); AssertThrow(ierr == 0, ExcPETScError(ierr)); + + // update ghost values if needed + update_ghost_values(); } @@ -821,7 +913,7 @@ namespace PETScWrappers AssertIsFinite(a); // there is nothing like a AXPAY - // operation in Petsc, so do it in two + // operation in PETSc, so do it in two // steps *this *= s; add(a, v); @@ -835,6 +927,9 @@ namespace PETScWrappers Assert(!has_ghost_elements(), ExcGhostsPresent()); const PetscErrorCode ierr = VecPointwiseMult(vector, factors, vector); AssertThrow(ierr == 0, ExcPETScError(ierr)); + + // update ghost values if needed + update_ghost_values(); } diff --git a/tests/petsc/copy_to_dealvec_block.cc b/tests/petsc/copy_to_dealvec_block.cc index cb6132c559..2c4803c80b 100644 --- a/tests/petsc/copy_to_dealvec_block.cc +++ b/tests/petsc/copy_to_dealvec_block.cc @@ -65,6 +65,17 @@ test() vb_one *= 2.0; v_one = vb_one; + // Create a copy of v_one using the internal Vec + PETScWrappers::MPI::Vector v_one_from_vec(v_one.petsc_vector()); + Assert(v_one.size() == v_one_from_vec.size(), ExcInternalError()); + Assert(v_one.locally_owned_size() == v_one_from_vec.locally_owned_size(), + ExcInternalError()); + Assert(v_one.has_ghost_elements() == v_one_from_vec.has_ghost_elements(), + ExcInternalError()); + Assert(v_one.ghost_elements() == v_one_from_vec.ghost_elements(), + ExcInternalError()); + + PETScWrappers::MPI::BlockVector vb, v; vb.reinit(2); v.reinit(2);