From f0e5d9571000f745ead2c476eeed75a91f957c78 Mon Sep 17 00:00:00 2001
From: Boris Martin <boris.martin@uliege.be>
Date: Mon, 21 Aug 2023 19:01:36 +0200
Subject: [PATCH] added getGlobalRHS()
---
src/problem/Formulation.cpp | 43 ++++++++++++++++++++++++++++++++++---
src/problem/Formulation.h | 1 +
2 files changed, 41 insertions(+), 3 deletions(-)
diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp
index c22ddca5..bccd1b18 100644
--- a/src/problem/Formulation.cpp
+++ b/src/problem/Formulation.cpp
@@ -1432,6 +1432,41 @@ namespace gmshddm
return _rhs;
}
+ template< class T_Scalar >
+ gmshfem::algebra::Vector< T_Scalar > Formulation< T_Scalar >::getGlobalRHS()
+ {
+ Vec GlobalRHS;
+ VecScatter ctx;
+ Vec DistributedRHS;
+ _buildPetscRHS(&DistributedRHS);
+
+ // Allocate locally the full vector
+ PetscInt n, N;
+ VecGetLocalSize(DistributedRHS, &n);
+ VecGetSize(DistributedRHS, &N);
+ //VecCreateSeq(PETSC_COMM_SELF, N, &GlobalRHS);
+
+ // Scatter: RHS (distributed) -> GlobalRHS (full copy on each process)
+ VecScatterCreateToAll(DistributedRHS, &ctx, &GlobalRHS);
+ VecScatterBegin(ctx, DistributedRHS, GlobalRHS, INSERT_VALUES, SCATTER_FORWARD);
+ VecScatterEnd(ctx, DistributedRHS, GlobalRHS, INSERT_VALUES, SCATTER_FORWARD);
+
+ // At this point, GlobalRHS on each process contains the entire vector
+ const PetscScalar *petscArray;
+ gmshfem::algebra::Vector<T_Scalar> result(N);
+ PetscInt petscArraySize;
+ VecGetArrayRead(GlobalRHS, &petscArray);
+ for (unsigned i = 0; i < N; ++i) {
+ result[i] = petscArray[i];
+ }
+ VecRestoreArrayRead(GlobalRHS, &petscArray);
+
+ // Clean up
+ VecScatterDestroy(&ctx);
+ VecDestroy(&GlobalRHS);
+ return result;
+ }
+
template< class T_Scalar >
unsigned int Formulation< T_Scalar >::numberOfIterations() const
{
@@ -1637,13 +1672,15 @@ namespace gmshddm
{
gmshfem::common::Timer time;
time.tick();
- if(mpi::getMPIRank() == 0) {
- gmshfem::msg::info << "Starting to build a block of RHS." << gmshfem::msg::endl;
- }
auto locSize = _rhs.size();
auto nRHS = _activeShots.size();
+ if(mpi::getMPIRank() == 0) {
+ gmshfem::msg::info << "Starting to build a block of " << nRHS << "RHS." << gmshfem::msg::endl;
+ }
+
+
MatCreateDense(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, PETSC_DETERMINE, nRHS, NULL, RHS);
MatSetUp(*RHS);
diff --git a/src/problem/Formulation.h b/src/problem/Formulation.h
index ded0c343..66a49b92 100644
--- a/src/problem/Formulation.h
+++ b/src/problem/Formulation.h
@@ -138,6 +138,7 @@ namespace gmshddm
void setSolutionIntoInterfaceFields(const std::vector< T_Scalar > &solution);
const gmshfem::algebra::Vector< T_Scalar > &getRHS() const;
+ gmshfem::algebra::Vector< T_Scalar > getGlobalRHS(); //Built on the fly, so not a reference. Modifies internally the Petsc representation of the local RHS.
unsigned int numberOfIterations() const;
--
GitLab