From 340e86bfbdc532af7c8428d75babe170032b9194 Mon Sep 17 00:00:00 2001 From: Sivapriya Venkateswarar Date: Mon, 20 Apr 2026 11:14:21 +0400 Subject: [PATCH] Add PETSc error handling across modules --- cpp/dolfinx/fem/petsc.cpp | 11 ++++-- cpp/dolfinx/la/petsc.cpp | 64 ++++++++++++++++++++++++-------- cpp/dolfinx/la/slepc.cpp | 31 ++++++++++++---- cpp/dolfinx/nls/NewtonSolver.cpp | 19 +++++++--- 4 files changed, 94 insertions(+), 31 deletions(-) diff --git a/cpp/dolfinx/fem/petsc.cpp b/cpp/dolfinx/fem/petsc.cpp index 305a20e9c2f..5cd1aa10091 100644 --- a/cpp/dolfinx/fem/petsc.cpp +++ b/cpp/dolfinx/fem/petsc.cpp @@ -54,8 +54,7 @@ Vec fem::petsc::create_vector_block( = VecCreateGhost(maps[0].first.get().comm(), local_size, PETSC_DETERMINE, _ghosts.size(), _ghosts.data(), &x); if (ierr != 0) - throw std::runtime_error("Call to PETSc VecCreateGhost failed."); - + petsc::error(ierr, __FILE__, "VecCreateGhost"); return x; } //----------------------------------------------------------------------------- @@ -76,8 +75,12 @@ Vec fem::petsc::create_vector_nest( // Create nested (VecNest) vector Vec y; - VecCreateNest(vecs.front()->comm(), petsc_vecs.size(), nullptr, - petsc_vecs.data(), &y); + PetscErrorCode ierr + = VecCreateNest(vecs.front()->comm(), petsc_vecs.size(), nullptr, + petsc_vecs.data(), &y); + if (ierr != 0) + petsc::error(ierr, __FILE__, "VecCreateNest"); + return y; } //----------------------------------------------------------------------------- diff --git a/cpp/dolfinx/la/petsc.cpp b/cpp/dolfinx/la/petsc.cpp index cb23cfafc4e..a4b6b9b7d71 100644 --- a/cpp/dolfinx/la/petsc.cpp +++ b/cpp/dolfinx/la/petsc.cpp @@ -46,7 +46,7 @@ void la::petsc::error(PetscErrorCode error_code, const std::string& filename, desc); throw std::runtime_error("Failed to successfully call PETSc function '" + petsc_function + "'. PETSc error code is: " - + std ::to_string(error_code) + ", " + + std::to_string(error_code) + ", " + std::string(desc)); } //----------------------------------------------------------------------------- @@ -57,11 +57,19 @@ la::petsc::create_vectors(MPI_Comm comm, std::vector v(x.size()); for (std::size_t i = 0; i < v.size(); ++i) { - VecCreateMPI(comm, x[i].size(), PETSC_DETERMINE, &v[i]); + PetscErrorCode ierr; + + ierr = VecCreateMPI(comm, x[i].size(), PETSC_DETERMINE, &v[i]); + CHECK_ERROR("VecCreateMPI"); + PetscScalar* data; - VecGetArray(v[i], &data); + ierr = VecGetArray(v[i], &data); + CHECK_ERROR("VecGetArray"); + std::ranges::copy(x[i], data); - VecRestoreArray(v[i], &data); + + ierr = VecRestoreArray(v[i], &data); + CHECK_ERROR("VecRestoreArray"); } return v; @@ -136,17 +144,20 @@ std::vector la::petsc::create_index_sets( std::int64_t offset = 0; for (auto& map : maps) { + PetscErrorCode ierr; int bs = map.second; std::int32_t size = map.first.get().size_local() + map.first.get().num_ghosts(); IS _is; - ISCreateStride(PETSC_COMM_SELF, bs * size, offset, 1, &_is); + ierr = ISCreateStride(PETSC_COMM_SELF, bs * size, offset, 1, &_is); + CHECK_ERROR("ISCreateStride"); is.push_back(_is); offset += bs * size; } return is; } + //----------------------------------------------------------------------------- std::vector> la::petsc::get_local_vectors( const Vec x, @@ -158,13 +169,21 @@ std::vector> la::petsc::get_local_vectors( for (auto& map : maps) offset_owned += map.first.get().size_local() * map.second; + PetscErrorCode ierr; + // Unwrap PETSc vector Vec x_local; - VecGhostGetLocalForm(x, &x_local); + ierr = VecGhostGetLocalForm(x, &x_local); + CHECK_ERROR("VecGhostGetLocalForm"); + PetscInt n = 0; - VecGetSize(x_local, &n); + ierr = VecGetSize(x_local, &n); + CHECK_ERROR("VecGetSize"); + const PetscScalar* array = nullptr; - VecGetArrayRead(x_local, &array); + ierr = VecGetArrayRead(x_local, &array); + CHECK_ERROR("VecGetArrayRead"); + std::span _x(array, n); // Copy PETSc Vec data in to local vectors @@ -185,11 +204,15 @@ std::vector> la::petsc::get_local_vectors( offset_ghost += size_ghost; } - VecRestoreArrayRead(x_local, &array); - VecGhostRestoreLocalForm(x, &x_local); + ierr = VecRestoreArrayRead(x_local, &array); + CHECK_ERROR("VecRestoreArrayRead"); + + ierr = VecGhostRestoreLocalForm(x, &x_local); + CHECK_ERROR("VecGhostRestoreLocalForm"); return x_b; } + //----------------------------------------------------------------------------- void la::petsc::scatter_local_vectors( Vec x, const std::vector>& x_b, @@ -204,12 +227,20 @@ void la::petsc::scatter_local_vectors( for (auto& map : maps) offset_owned += map.first.get().size_local() * map.second; + PetscErrorCode ierr; + Vec x_local; - VecGhostGetLocalForm(x, &x_local); + ierr = VecGhostGetLocalForm(x, &x_local); + CHECK_ERROR("VecGhostGetLocalForm"); + PetscInt n = 0; - VecGetSize(x_local, &n); + ierr = VecGetSize(x_local, &n); + CHECK_ERROR("VecGetSize"); + PetscScalar* array = nullptr; - VecGetArray(x_local, &array); + ierr = VecGetArray(x_local, &array); + CHECK_ERROR("VecGetArray"); + std::span _x(array, n); // Copy local vectors into PETSc Vec @@ -228,8 +259,11 @@ void la::petsc::scatter_local_vectors( offset_ghost += size_ghost; } - VecRestoreArray(x_local, &array); - VecGhostRestoreLocalForm(x, &x_local); + ierr = VecRestoreArray(x_local, &array); + CHECK_ERROR("VecRestoreArray"); + + ierr = VecGhostRestoreLocalForm(x, &x_local); + CHECK_ERROR("VecGhostRestoreLocalForm"); } //----------------------------------------------------------------------------- Mat la::petsc::create_matrix(MPI_Comm comm, const SparsityPattern& sp, diff --git a/cpp/dolfinx/la/slepc.cpp b/cpp/dolfinx/la/slepc.cpp index a942d5ccb18..f72457d703a 100644 --- a/cpp/dolfinx/la/slepc.cpp +++ b/cpp/dolfinx/la/slepc.cpp @@ -18,7 +18,12 @@ using namespace dolfinx; using namespace dolfinx::la; //----------------------------------------------------------------------------- -SLEPcEigenSolver::SLEPcEigenSolver(MPI_Comm comm) { EPSCreate(comm, &_eps); } +SLEPcEigenSolver::SLEPcEigenSolver(MPI_Comm comm) +{ + PetscErrorCode ierr = EPSCreate(comm, &_eps); + if (ierr != 0) + petsc::error(ierr, __FILE__, "EPSCreate"); +} //----------------------------------------------------------------------------- SLEPcEigenSolver::SLEPcEigenSolver(EPS eps, bool inc_ref_count) : _eps(eps) { @@ -56,7 +61,9 @@ SLEPcEigenSolver::operator=(SLEPcEigenSolver&& solver) noexcept void SLEPcEigenSolver::set_operators(const Mat A, const Mat B) { assert(_eps); - EPSSetOperators(_eps, A, B); + PetscErrorCode ierr = EPSSetOperators(_eps, A, B); + if (ierr != 0) + petsc::error(ierr, __FILE__, "EPSSetOperators"); } //----------------------------------------------------------------------------- void SLEPcEigenSolver::solve() @@ -64,10 +71,15 @@ void SLEPcEigenSolver::solve() // Get operators Mat A, B; assert(_eps); - EPSGetOperators(_eps, &A, &B); + PetscErrorCode ierr = EPSGetOperators(_eps, &A, &B); + if (ierr != 0) + petsc::error(ierr, __FILE__, "EPSGetOperators"); PetscInt m(0), n(0); - MatGetSize(A, &m, &n); + ierr = MatGetSize(A, &m, &n); + if (ierr != 0) + petsc::error(ierr, __FILE__, "MatGetSize"); + solve(m); } //----------------------------------------------------------------------------- @@ -77,10 +89,15 @@ void SLEPcEigenSolver::solve(std::int64_t n) // Get operators Mat A, B; assert(_eps); - EPSGetOperators(_eps, &A, &B); + PetscErrorCode ierr = EPSGetOperators(_eps, &A, &B); + if (ierr != 0) + petsc::error(ierr, __FILE__, "EPSGetOperators"); PetscInt _m(0), _n(0); - MatGetSize(A, &_m, &_n); + ierr = MatGetSize(A, &_m, &_n); + if (ierr != 0) + petsc::error(ierr, __FILE__, "MatGetSize"); + assert(n <= _n); #endif @@ -208,4 +225,4 @@ MPI_Comm SLEPcEigenSolver::comm() const } //----------------------------------------------------------------------------- -#endif +#endif \ No newline at end of file diff --git a/cpp/dolfinx/nls/NewtonSolver.cpp b/cpp/dolfinx/nls/NewtonSolver.cpp index 1a2984abc9e..b97b7c09422 100644 --- a/cpp/dolfinx/nls/NewtonSolver.cpp +++ b/cpp/dolfinx/nls/NewtonSolver.cpp @@ -29,7 +29,9 @@ std::pair converged(const nls::petsc::NewtonSolver& solver, const Vec r) { PetscReal residual = 0; - VecNorm(r, NORM_2, &residual); + PetscErrorCode ierr = VecNorm(r, NORM_2, &residual); + if (ierr != 0) + petsc::error(ierr, __FILE__, "VecNorm"); // Relative residual const double relative_residual = residual / solver.residual0(); @@ -62,7 +64,9 @@ std::pair converged(const nls::petsc::NewtonSolver& solver, void update_solution(const nls::petsc::NewtonSolver& solver, const Vec dx, Vec x) { - VecAXPY(x, -solver.relaxation_parameter, dx); + PetscErrorCode ierr = VecAXPY(x, -solver.relaxation_parameter, dx); + if (ierr != 0) + petsc::error(ierr, __FILE__, "VecAXPY"); } //----------------------------------------------------------------------------- } // namespace @@ -195,8 +199,11 @@ std::pair nls::petsc::NewtonSolver::solve(Vec x) _solver.set_operators(_matJ, _matJ); if (!_dx) - MatCreateVecs(_matJ, &_dx, nullptr); - + { + PetscErrorCode ierr = MatCreateVecs(_matJ, &_dx, nullptr); + if (ierr != 0) + petsc::error(ierr, __FILE__, "MatCreateVecs"); + } // Start iterations while (!newton_converged and _iteration < max_it) { @@ -224,7 +231,9 @@ std::pair nls::petsc::NewtonSolver::solve(Vec x) if (_iteration == 1) { PetscReal _r = 0; - VecNorm(_dx, NORM_2, &_r); + PetscErrorCode ierr = VecNorm(_dx, NORM_2, &_r); + if (ierr != 0) + petsc::error(ierr, __FILE__, "VecNorm"); _residual0 = _r; }