diff --git a/src/problem/Formulation.cpp b/src/problem/Formulation.cpp index 5be2edffc24bc2dabb0399befa5910b3d2375bad..c22ddca56f8af75693f83bfb01d1fb0da481358d 100644 --- a/src/problem/Formulation.cpp +++ b/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) diff --git a/src/problem/Formulation.h b/src/problem/Formulation.h index ad8acf945ef5f2804b8f93180826eb3291d0d34e..ded0c343e7a5a59c2aa49c17f99fd9707cac1002 100644 --- a/src/problem/Formulation.h +++ b/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);