diff --git a/src/system/DistributedContext.cpp b/src/system/DistributedContext.cpp
index 99c37e11acc650d6866bfc0741340fc8fbc3a1a5..988a8f3aad29cede3aa51bb4b83d8ebe0103a6eb 100644
--- a/src/system/DistributedContext.cpp
+++ b/src/system/DistributedContext.cpp
@@ -61,6 +61,16 @@ namespace gmshfem::system
       return _distributedExtendedVector;
       #endif
     }
+    PetscErrorCode sumDistributedRHS(Vec localRHSdistributed, Vec* out) const {
+      PetscFunctionBeginUser;
+      Vec res;
+      PetscCall(VecDuplicate(_distributedVector, &res));
+      VecView(localRHSdistributed, PETSC_VIEWER_STDOUT_WORLD);
+      PetscCall(VecScatterBegin(_scatter, localRHSdistributed, res, ADD_VALUES, SCATTER_REVERSE));
+      PetscCall(VecScatterEnd(_scatter, localRHSdistributed, res, ADD_VALUES, SCATTER_REVERSE));
+      *out = res;
+      PetscFunctionReturn(PETSC_SUCCESS);
+    };
     ~DistributedContextImpl() {
       #ifdef HAVE_PETSC
       VecDestroy(&_distributedExtendedVector);
@@ -125,6 +135,16 @@ namespace gmshfem::system
 
     PetscFunctionReturn(PETSC_SUCCESS);
   }
+
+  template< class T_Scalar >
+  PetscErrorCode DistributedContext< T_Scalar >::sumDistributedRHS(Vec localRHSdistributed, Vec *out) const
+  {
+    PetscFunctionBeginUser;
+    if (!_impl)
+      throw common::Exception("Uninitialized DistributedContextImpl");
+    PetscCall(_impl->sumDistributedRHS(localRHSdistributed, out));
+    PetscFunctionReturn(PETSC_SUCCESS);
+  }
 #endif
 
   template <class T_Scalar>
diff --git a/src/system/DistributedContext.h b/src/system/DistributedContext.h
index b037e1630864cde287030e5ee888822adf8a8ade..eff68c99f6c975f4ee08311a4da25b045e8768ab 100644
--- a/src/system/DistributedContext.h
+++ b/src/system/DistributedContext.h
@@ -45,9 +45,12 @@ class DistributedContext {
     const Indices& localToGlobal() const;
     const Indices& localIDofOwned() const;
     const Indices& localIDofNonOwned() const;
+    // Vec is a distbuted vector. We scatter to an extended vector (duplicated overlaps) then set the local vals
     void readScatteredData(std::vector<T_Scalar>& values, Vec sol) const;
     #ifdef HAVE_PETSC
     PetscErrorCode getSubdomainsIS(IS* is) const;
+    // Map an extended vector into a global vector with sums. Output becomes a standalone vector
+    PetscErrorCode sumDistributedRHS(Vec localRHSdistributed, Vec* out) const;
     #endif
     DistributedContext(Indices&& localToGlobal, Indices&& owned, Indices&& nonOwned);
     ~DistributedContext();
diff --git a/src/system/VectorFactory.h b/src/system/VectorFactory.h
index e27d5013b4a8440ba78d75180b7c174b0e5fce70..77bc2dbe0060b2d14730004a6f83b480ce214180 100644
--- a/src/system/VectorFactory.h
+++ b/src/system/VectorFactory.h
@@ -59,6 +59,7 @@ namespace gmshfem::system
     Vec getPetsc() const;
     #ifdef HAVE_MPI
     Vec getPetscDistributed(const DistributedContext<T_Scalar>& distributedContext) const;
+    PetscErrorCode getPetscDistributedSummedRHS(const DistributedContext<T_Scalar>& distributedContext, Vec* out) const;
     // Return PetscScalar* but not declared yet
     void* generateRawLocalData(const DistributedContext<T_Scalar>& distributedContext) const;