Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 7 additions & 4 deletions cpp/dolfinx/fem/petsc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
//-----------------------------------------------------------------------------
Expand All @@ -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;
}
//-----------------------------------------------------------------------------
Expand Down
64 changes: 49 additions & 15 deletions cpp/dolfinx/la/petsc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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));
}
//-----------------------------------------------------------------------------
Expand All @@ -57,11 +57,19 @@ la::petsc::create_vectors(MPI_Comm comm,
std::vector<Vec> 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;
Expand Down Expand Up @@ -136,17 +144,20 @@ std::vector<IS> 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<std::vector<PetscScalar>> la::petsc::get_local_vectors(
const Vec x,
Expand All @@ -158,13 +169,21 @@ std::vector<std::vector<PetscScalar>> 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
Expand All @@ -185,11 +204,15 @@ std::vector<std::vector<PetscScalar>> 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<std::span<const PetscScalar>>& x_b,
Expand All @@ -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
Expand All @@ -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,
Expand Down
31 changes: 24 additions & 7 deletions cpp/dolfinx/la/slepc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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)
{
Expand Down Expand Up @@ -56,18 +61,25 @@ 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()
{
// 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);
}
//-----------------------------------------------------------------------------
Expand All @@ -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

Expand Down Expand Up @@ -208,4 +225,4 @@ MPI_Comm SLEPcEigenSolver::comm() const
}
//-----------------------------------------------------------------------------

#endif
#endif
19 changes: 14 additions & 5 deletions cpp/dolfinx/nls/NewtonSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,9 @@ std::pair<double, bool> 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();
Expand Down Expand Up @@ -62,7 +64,9 @@ std::pair<double, bool> 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
Expand Down Expand Up @@ -195,8 +199,11 @@ std::pair<int, bool> 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)
{
Expand Down Expand Up @@ -224,7 +231,9 @@ std::pair<int, bool> 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;
}

Expand Down