Building block RHS

src/problem/Formulation.cpp

@@ -1040,8 +1040,7 @@ namespace gmshddm
       PetscInt globSize;
 
       // Create Petsc distributed vector - let Petsc compute the global vector size
-      VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, &RHS);
-      VecCopy(_rhs.getPetsc(), RHS);
+      _buildPetscRHS(&RHS);
       VecGetSize(RHS, &globSize);
 
       VecDuplicate(RHS, &X); // current iterate
@@ -1624,6 +1623,58 @@ namespace gmshddm
       }
     }
 
+    template< class T_Scalar >
+    void Formulation< T_Scalar >::_buildPetscRHS(Vec *RHS)
+    {
+
+      auto locSize = _rhs.size();
+      VecCreateMPI(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, RHS);
+      VecCopy(_rhs.getPetsc(), *RHS);
+    }
+
+    template< class T_Scalar >
+    void Formulation< T_Scalar >::_buildPetscAllRHS(Mat *RHS)
+    {
+      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();
+
+      MatCreateDense(MPI_COMM_WORLD, locSize, PETSC_DETERMINE, PETSC_DETERMINE, nRHS, NULL, RHS);
+      MatSetUp(*RHS);
+
+      std::vector< PetscInt > rowIndices(locSize);
+      PetscInt startRow, endRow;
+      MatGetOwnershipRange(*RHS, &startRow, &endRow);
+
+      for(PetscInt i = 0; i < locSize; i++) {
+        rowIndices[i] = startRow + i;
+      }
+
+      // For each RHS, copy the vector into the column of the block
+      for(PetscInt j = 0; j < nRHS; j++) {
+        Vec bj = _multi_rhs[j].getPetsc();
+
+        PetscScalar *colValues;
+        VecGetArray(bj, &colValues);
+        MatSetValues(*RHS, locSize, rowIndices.data(), 1, &j, colValues, INSERT_VALUES);
+        VecRestoreArray(bj, &colValues);
+      }
+
+      MatAssemblyBegin(*RHS, MAT_FINAL_ASSEMBLY);
+      MatAssemblyEnd(*RHS, MAT_FINAL_ASSEMBLY);
+
+
+      time.tock();
+      if(mpi::getMPIRank() == 0) {
+        gmshfem::msg::info << "Done assembling block of RHS in " << time << "s" << gmshfem::msg::endl;
+      }
+    }
+
 
     template< class T_Scalar >
     int MatVectProductImpl(Mat A, Vec X, Vec Y, bool IA)

src/problem/Formulation.h

@@ -89,6 +89,9 @@ namespace gmshddm
       void _activateBilinear (const unsigned long long idom);
       void _deactivateBilinear (const unsigned long long idom);
 
+      void _buildPetscRHS(Vec* RHS);
+      void _buildPetscAllRHS(Mat* RHS);
+
      public:
       Formulation(const std::string &name, const std::vector< std::vector< unsigned int > > &topology);
       Formulation(const std::string &name, const std::vector< gmshfem::problem::Formulation< T_Scalar > * > &volume, const std::vector< std::unordered_map< unsigned int, gmshfem::problem::Formulation< T_Scalar > * > > &surface);