diff --git a/cpp/dolfinx/fem/petsc.cpp b/cpp/dolfinx/fem/petsc.cpp
index 305a20e9c2f..323c4ac80ea 100644
--- a/cpp/dolfinx/fem/petsc.cpp
+++ b/cpp/dolfinx/fem/petsc.cpp
@@ -54,7 +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.");
+    la::petsc::error(ierr, __FILE__, "VecCreateGhost");
 
   return x;
 }
@@ -76,8 +76,10 @@ 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)
+    la::petsc::error(ierr, __FILE__, "VecCreateNest");
   return y;
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/la/petsc.cpp b/cpp/dolfinx/la/petsc.cpp
index cb23cfafc4e..eb41a6a0317 100644
--- a/cpp/dolfinx/la/petsc.cpp
+++ b/cpp/dolfinx/la/petsc.cpp
@@ -57,11 +57,15 @@ 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]);
-    PetscScalar* data;
-    VecGetArray(v[i], &data);
+    PetscErrorCode ierr
+        = VecCreateMPI(comm, x[i].size(), PETSC_DETERMINE, &v[i]);
+    CHECK_ERROR("VecCreateMPI");
+    PetscScalar* data = nullptr;
+    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;
@@ -140,7 +144,9 @@ std::vector<IS> la::petsc::create_index_sets(
     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);
+    PetscErrorCode ierr
+        = ISCreateStride(PETSC_COMM_SELF, bs * size, offset, 1, &_is);
+    CHECK_ERROR("ISCreateStride");
     is.push_back(_is);
     offset += bs * size;
   }
@@ -160,11 +166,14 @@ std::vector<std::vector<PetscScalar>> la::petsc::get_local_vectors(
 
   // Unwrap PETSc vector
   Vec x_local;
-  VecGhostGetLocalForm(x, &x_local);
+  PetscErrorCode 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,8 +194,10 @@ 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;
 }
@@ -205,11 +216,14 @@ void la::petsc::scatter_local_vectors(
     offset_owned += map.first.get().size_local() * map.second;
 
   Vec x_local;
-  VecGhostGetLocalForm(x, &x_local);
+  PetscErrorCode 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 +242,10 @@ 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,
@@ -245,8 +261,12 @@ Mat la::petsc::create_matrix(MPI_Comm comm, const SparsityPattern& sp,
   std::array maps = {sp.index_map(0), sp.index_map(1)};
   const std::array bs = {sp.block_size(0), sp.block_size(1)};
 
-  if (type)
-    MatSetType(A, type->c_str());
+  if (type and !type->empty())
+  {
+    ierr = MatSetType(A, type->c_str());
+    if (ierr != 0)
+      petsc::error(ierr, __FILE__, "MatSetType");
+  }
 
   // Get global and local dimensions
   const std::int64_t M = bs[0] * maps[0]->size_global();
@@ -294,7 +314,7 @@ Mat la::petsc::create_matrix(MPI_Comm comm, const SparsityPattern& sp,
   ierr = MatXAIJSetPreallocation(A, _bs, _nnz_diag.data(), _nnz_offdiag.data(),
                                  nullptr, nullptr);
   if (ierr != 0)
-    petsc::error(ierr, __FILE__, "MatXIJSetPreallocation");
+    petsc::error(ierr, __FILE__, "MatXAIJSetPreallocation");
 
   // Set block sizes
   ierr = MatSetBlockSizes(A, bs[0], bs[1]);
@@ -405,7 +425,10 @@ petsc::Vector::Vector(Vec x, bool inc_ref_count) : _x(x)
 {
   assert(x);
   if (inc_ref_count)
-    PetscObjectReference((PetscObject)_x);
+  {
+    PetscErrorCode ierr = PetscObjectReference((PetscObject)_x);
+    CHECK_ERROR("PetscObjectReference");
+  }
 }
 //-----------------------------------------------------------------------------
 petsc::Vector::Vector(Vector&& v) noexcept : _x(std::exchange(v._x, nullptr)) {}
@@ -413,7 +436,7 @@ petsc::Vector::Vector(Vector&& v) noexcept : _x(std::exchange(v._x, nullptr)) {}
 petsc::Vector::~Vector()
 {
   if (_x)
-    VecDestroy(&_x);
+    (void)VecDestroy(&_x);
 }
 //-----------------------------------------------------------------------------
 petsc::Vector& petsc::Vector::operator=(Vector&& v) noexcept
@@ -425,10 +448,13 @@ petsc::Vector& petsc::Vector::operator=(Vector&& v) noexcept
 petsc::Vector petsc::Vector::copy() const
 {
   Vec _y;
-  VecDuplicate(_x, &_y);
-  VecCopy(_x, _y);
+  PetscErrorCode ierr = VecDuplicate(_x, &_y);
+  CHECK_ERROR("VecDuplicate");
+  ierr = VecCopy(_x, _y);
+  CHECK_ERROR("VecCopy");
   Vector y(_y, true);
-  VecDestroy(&_y);
+  ierr = VecDestroy(&_y);
+  CHECK_ERROR("VecDestroy");
   return y;
 }
 //-----------------------------------------------------------------------------
@@ -499,7 +525,10 @@ petsc::Operator::Operator(Mat A, bool inc_ref_count) : _matA(A)
 {
   assert(A);
   if (inc_ref_count)
-    PetscObjectReference((PetscObject)_matA);
+  {
+    PetscErrorCode ierr = PetscObjectReference((PetscObject)_matA);
+    CHECK_ERROR("PetscObjectReference");
+  }
 }
 //-----------------------------------------------------------------------------
 petsc::Operator::Operator(Operator&& A) noexcept
@@ -512,7 +541,7 @@ petsc::Operator::~Operator()
   // Decrease reference count (PETSc will destroy object once reference
   // counts reached zero)
   if (_matA)
-    MatDestroy(&_matA);
+    (void)MatDestroy(&_matA);
 }
 //-----------------------------------------------------------------------------
 petsc::Operator& petsc::Operator::operator=(Operator&& A) noexcept
@@ -527,7 +556,7 @@ std::array<std::int64_t, 2> petsc::Operator::size() const
   PetscInt m(0), n(0);
   PetscErrorCode ierr = MatGetSize(_matA, &m, &n);
   if (ierr != 0)
-    petsc::error(ierr, __FILE__, "MetGetSize");
+    petsc::error(ierr, __FILE__, "MatGetSize");
   return {{m, n}};
 }
 //-----------------------------------------------------------------------------
@@ -620,21 +649,24 @@ void petsc::Matrix::apply(AssemblyType type)
 void petsc::Matrix::set_options_prefix(const std::string& options_prefix)
 {
   assert(_matA);
-  MatSetOptionsPrefix(_matA, options_prefix.c_str());
+  PetscErrorCode ierr = MatSetOptionsPrefix(_matA, options_prefix.c_str());
+  CHECK_ERROR("MatSetOptionsPrefix");
 }
 //-----------------------------------------------------------------------------
 std::string petsc::Matrix::get_options_prefix() const
 {
   assert(_matA);
   const char* prefix = nullptr;
-  MatGetOptionsPrefix(_matA, &prefix);
+  PetscErrorCode ierr = MatGetOptionsPrefix(_matA, &prefix);
+  CHECK_ERROR("MatGetOptionsPrefix");
   return std::string(prefix);
 }
 //-----------------------------------------------------------------------------
 void petsc::Matrix::set_from_options()
 {
   assert(_matA);
-  MatSetFromOptions(_matA);
+  PetscErrorCode ierr = MatSetFromOptions(_matA);
+  CHECK_ERROR("MatSetFromOptions");
 }
 //-----------------------------------------------------------------------------
 //-----------------------------------------------------------------------------
@@ -666,7 +698,7 @@ petsc::KrylovSolver::KrylovSolver(KrylovSolver&& solver) noexcept
 petsc::KrylovSolver::~KrylovSolver()
 {
   if (_ksp)
-    KSPDestroy(&_ksp);
+    (void)KSPDestroy(&_ksp);
 }
 //-----------------------------------------------------------------------------
 petsc::KrylovSolver&
@@ -695,11 +727,11 @@ int petsc::KrylovSolver::solve(Vec x, const Vec b, bool transpose) const
 
   // Get PETSc operators
   Mat _A, _P;
-  KSPGetOperators(_ksp, &_A, &_P);
+  PetscErrorCode ierr = KSPGetOperators(_ksp, &_A, &_P);
+  if (ierr != 0)
+    petsc::error(ierr, __FILE__, "KSPGetOperators");
   assert(_A);
 
-  PetscErrorCode ierr;
-
   // Solve linear system
   spdlog::info("PETSc Krylov solver starting to solve system.");
 
diff --git a/cpp/dolfinx/la/slepc.cpp b/cpp/dolfinx/la/slepc.cpp
index a942d5ccb18..5ee14191a2f 100644
--- a/cpp/dolfinx/la/slepc.cpp
+++ b/cpp/dolfinx/la/slepc.cpp
@@ -17,8 +17,22 @@
 using namespace dolfinx;
 using namespace dolfinx::la;
 
+namespace
+{
+//-----------------------------------------------------------------------------
+void check_petsc_error(PetscErrorCode ierr, const char* petsc_function)
+{
+  if (ierr != 0)
+    petsc::error(ierr, __FILE__, petsc_function);
+}
+//-----------------------------------------------------------------------------
+} // namespace
+
 //-----------------------------------------------------------------------------
-SLEPcEigenSolver::SLEPcEigenSolver(MPI_Comm comm) { EPSCreate(comm, &_eps); }
+SLEPcEigenSolver::SLEPcEigenSolver(MPI_Comm comm)
+{
+  check_petsc_error(EPSCreate(comm, &_eps), "EPSCreate");
+}
 //-----------------------------------------------------------------------------
 SLEPcEigenSolver::SLEPcEigenSolver(EPS eps, bool inc_ref_count) : _eps(eps)
 {
@@ -29,8 +43,7 @@ SLEPcEigenSolver::SLEPcEigenSolver(EPS eps, bool inc_ref_count) : _eps(eps)
   if (inc_ref_count)
   {
     ierr = PetscObjectReference((PetscObject)_eps);
-    if (ierr != 0)
-      petsc::error(ierr, __FILE__, "PetscObjectReference");
+    check_petsc_error(ierr, "PetscObjectReference");
   }
 }
 //-----------------------------------------------------------------------------
@@ -43,7 +56,7 @@ SLEPcEigenSolver::SLEPcEigenSolver(SLEPcEigenSolver&& solver) noexcept
 SLEPcEigenSolver::~SLEPcEigenSolver()
 {
   if (_eps)
-    EPSDestroy(&_eps);
+    (void)EPSDestroy(&_eps);
 }
 //-----------------------------------------------------------------------------
 SLEPcEigenSolver&
@@ -56,7 +69,7 @@ SLEPcEigenSolver::operator=(SLEPcEigenSolver&& solver) noexcept
 void SLEPcEigenSolver::set_operators(const Mat A, const Mat B)
 {
   assert(_eps);
-  EPSSetOperators(_eps, A, B);
+  check_petsc_error(EPSSetOperators(_eps, A, B), "EPSSetOperators");
 }
 //-----------------------------------------------------------------------------
 void SLEPcEigenSolver::solve()
@@ -64,10 +77,10 @@ void SLEPcEigenSolver::solve()
   // Get operators
   Mat A, B;
   assert(_eps);
-  EPSGetOperators(_eps, &A, &B);
+  check_petsc_error(EPSGetOperators(_eps, &A, &B), "EPSGetOperators");
 
   PetscInt m(0), n(0);
-  MatGetSize(A, &m, &n);
+  check_petsc_error(MatGetSize(A, &m, &n), "MatGetSize");
   solve(m);
 }
 //-----------------------------------------------------------------------------
@@ -77,35 +90,38 @@ void SLEPcEigenSolver::solve(std::int64_t n)
   // Get operators
   Mat A, B;
   assert(_eps);
-  EPSGetOperators(_eps, &A, &B);
+  check_petsc_error(EPSGetOperators(_eps, &A, &B), "EPSGetOperators");
 
   PetscInt _m(0), _n(0);
-  MatGetSize(A, &_m, &_n);
+  check_petsc_error(MatGetSize(A, &_m, &_n), "MatGetSize");
   assert(n <= _n);
 #endif
 
   // Set number of eigenpairs to compute
   assert(_eps);
-  EPSSetDimensions(_eps, n, PETSC_DECIDE, PETSC_DECIDE);
+  check_petsc_error(EPSSetDimensions(_eps, n, PETSC_DECIDE, PETSC_DECIDE),
+                    "EPSSetDimensions");
 
   // Set any options from the PETSc database
-  EPSSetFromOptions(_eps);
+  check_petsc_error(EPSSetFromOptions(_eps), "EPSSetFromOptions");
 
   // Solve eigenvalue problem
-  EPSSolve(_eps);
+  check_petsc_error(EPSSolve(_eps), "EPSSolve");
 
   // Check for convergence
   EPSConvergedReason reason;
-  EPSGetConvergedReason(_eps, &reason);
+  check_petsc_error(EPSGetConvergedReason(_eps, &reason),
+                    "EPSGetConvergedReason");
   if (reason < 0)
     spdlog::warn("Eigenvalue solver did not converge");
 
   // Report solver status
   PetscInt num_iterations = 0;
-  EPSGetIterationNumber(_eps, &num_iterations);
+  check_petsc_error(EPSGetIterationNumber(_eps, &num_iterations),
+                    "EPSGetIterationNumber");
 
   EPSType eps_type = nullptr;
-  EPSGetType(_eps, &eps_type);
+  check_petsc_error(EPSGetType(_eps, &eps_type), "EPSGetType");
   spdlog::info("Eigenvalue solver ({}) converged in {} iterations.", eps_type,
                num_iterations);
 }
@@ -116,17 +132,20 @@ std::complex<PetscReal> SLEPcEigenSolver::get_eigenvalue(int i) const
 
   // Get number of computed values
   PetscInt num_computed_eigenvalues;
-  EPSGetConverged(_eps, &num_computed_eigenvalues);
+  check_petsc_error(EPSGetConverged(_eps, &num_computed_eigenvalues),
+                    "EPSGetConverged");
 
   if (i < num_computed_eigenvalues)
   {
 #ifdef PETSC_USE_COMPLEX
     PetscScalar l;
-    EPSGetEigenvalue(_eps, i, &l, nullptr);
+    check_petsc_error(EPSGetEigenvalue(_eps, i, &l, nullptr),
+                      "EPSGetEigenvalue");
     return l;
 #else
     PetscScalar lr, li;
-    EPSGetEigenvalue(_eps, i, &lr, &li);
+    check_petsc_error(EPSGetEigenvalue(_eps, i, &lr, &li),
+                      "EPSGetEigenvalue");
     return std::complex<PetscReal>(lr, li);
 #endif
   }
@@ -145,9 +164,11 @@ void SLEPcEigenSolver::get_eigenpair(PetscScalar& lr, PetscScalar& lc, Vec r,
 
   // Get number of computed eigenvectors/values
   PetscInt num_computed_eigenvalues;
-  EPSGetConverged(_eps, &num_computed_eigenvalues);
+  check_petsc_error(EPSGetConverged(_eps, &num_computed_eigenvalues),
+                    "EPSGetConverged");
   if (ii < num_computed_eigenvalues)
-    EPSGetEigenpair(_eps, ii, &lr, &lc, r, c);
+    check_petsc_error(EPSGetEigenpair(_eps, ii, &lr, &lc, r, c),
+                      "EPSGetEigenpair");
   else
   {
     throw std::runtime_error("Requested eigenpair (" + std::to_string(i)
@@ -159,7 +180,7 @@ std::int64_t SLEPcEigenSolver::get_number_converged() const
 {
   PetscInt num_conv;
   assert(_eps);
-  EPSGetConverged(_eps, &num_conv);
+  check_petsc_error(EPSGetConverged(_eps, &num_conv), "EPSGetConverged");
   return num_conv;
 }
 //-----------------------------------------------------------------------------
@@ -167,8 +188,7 @@ void SLEPcEigenSolver::set_options_prefix(const std::string& options_prefix)
 {
   assert(_eps);
   PetscErrorCode ierr = EPSSetOptionsPrefix(_eps, options_prefix.c_str());
-  if (ierr != 0)
-    petsc::error(ierr, __FILE__, "EPSSetOptionsPrefix");
+  check_petsc_error(ierr, "EPSSetOptionsPrefix");
 }
 //-----------------------------------------------------------------------------
 std::string SLEPcEigenSolver::get_options_prefix() const
@@ -176,8 +196,7 @@ std::string SLEPcEigenSolver::get_options_prefix() const
   assert(_eps);
   const char* prefix = nullptr;
   PetscErrorCode ierr = EPSGetOptionsPrefix(_eps, &prefix);
-  if (ierr != 0)
-    petsc::error(ierr, __FILE__, "EPSGetOptionsPrefix");
+  check_petsc_error(ierr, "EPSGetOptionsPrefix");
   return std::string(prefix);
 }
 //-----------------------------------------------------------------------------
@@ -185,15 +204,15 @@ void SLEPcEigenSolver::set_from_options() const
 {
   assert(_eps);
   PetscErrorCode ierr = EPSSetFromOptions(_eps);
-  if (ierr != 0)
-    petsc::error(ierr, __FILE__, "EPSSetFromOptions");
+  check_petsc_error(ierr, "EPSSetFromOptions");
 }
 //-----------------------------------------------------------------------------
 int SLEPcEigenSolver::get_iteration_number() const
 {
   assert(_eps);
   PetscInt num_iter;
-  EPSGetIterationNumber(_eps, &num_iter);
+  check_petsc_error(EPSGetIterationNumber(_eps, &num_iter),
+                    "EPSGetIterationNumber");
   return num_iter;
 }
 //-----------------------------------------------------------------------------
@@ -203,7 +222,8 @@ MPI_Comm SLEPcEigenSolver::comm() const
 {
   assert(_eps);
   MPI_Comm mpi_comm = MPI_COMM_NULL;
-  PetscObjectGetComm((PetscObject)_eps, &mpi_comm);
+  check_petsc_error(PetscObjectGetComm((PetscObject)_eps, &mpi_comm),
+                    "PetscObjectGetComm");
   return mpi_comm;
 }
 //-----------------------------------------------------------------------------
diff --git a/cpp/dolfinx/nls/NewtonSolver.cpp b/cpp/dolfinx/nls/NewtonSolver.cpp
index 1a2984abc9e..80f84751d4c 100644
--- a/cpp/dolfinx/nls/NewtonSolver.cpp
+++ b/cpp/dolfinx/nls/NewtonSolver.cpp
@@ -19,6 +19,12 @@ using namespace dolfinx;
 namespace
 {
 //-----------------------------------------------------------------------------
+void check_petsc_error(PetscErrorCode ierr, const std::string& petsc_function)
+{
+  if (ierr != 0)
+    la::petsc::error(ierr, __FILE__, petsc_function);
+}
+//-----------------------------------------------------------------------------
 
 /// Convergence test
 /// @param solver The Newton solver
@@ -29,7 +35,7 @@ std::pair<double, bool> converged(const nls::petsc::NewtonSolver& solver,
                                   const Vec r)
 {
   PetscReal residual = 0;
-  VecNorm(r, NORM_2, &residual);
+  check_petsc_error(VecNorm(r, NORM_2, &residual), "VecNorm");
 
   // Relative residual
   const double relative_residual = residual / solver.residual0();
@@ -62,7 +68,7 @@ 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);
+  check_petsc_error(VecAXPY(x, -solver.relaxation_parameter, dx), "VecAXPY");
 }
 //-----------------------------------------------------------------------------
 } // namespace
@@ -85,13 +91,13 @@ nls::petsc::NewtonSolver::NewtonSolver(MPI_Comm comm)
 nls::petsc::NewtonSolver::~NewtonSolver()
 {
   if (_b)
-    VecDestroy(&_b);
+    (void)VecDestroy(&_b);
   if (_dx)
-    VecDestroy(&_dx);
+    (void)VecDestroy(&_dx);
   if (_matJ)
-    MatDestroy(&_matJ);
+    (void)MatDestroy(&_matJ);
   if (_matP)
-    MatDestroy(&_matP);
+    (void)MatDestroy(&_matP);
 }
 //-----------------------------------------------------------------------------
 void nls::petsc::NewtonSolver::setF(std::function<void(const Vec, Vec)> F,
@@ -99,7 +105,8 @@ void nls::petsc::NewtonSolver::setF(std::function<void(const Vec, Vec)> F,
 {
   _fnF = std::move(F);
   _b = b;
-  PetscObjectReference((PetscObject)_b);
+  check_petsc_error(PetscObjectReference((PetscObject)_b),
+                    "PetscObjectReference");
 }
 //-----------------------------------------------------------------------------
 void nls::petsc::NewtonSolver::setJ(std::function<void(const Vec, Mat)> J,
@@ -107,7 +114,8 @@ void nls::petsc::NewtonSolver::setJ(std::function<void(const Vec, Mat)> J,
 {
   _fnJ = std::move(J);
   _matJ = Jmat;
-  PetscObjectReference((PetscObject)_matJ);
+  check_petsc_error(PetscObjectReference((PetscObject)_matJ),
+                    "PetscObjectReference");
 }
 //-----------------------------------------------------------------------------
 void nls::petsc::NewtonSolver::setP(std::function<void(const Vec, Mat)> P,
@@ -115,7 +123,8 @@ void nls::petsc::NewtonSolver::setP(std::function<void(const Vec, Mat)> P,
 {
   _fnP = std::move(P);
   _matP = Pmat;
-  PetscObjectReference((PetscObject)_matP);
+  check_petsc_error(PetscObjectReference((PetscObject)_matP),
+                    "PetscObjectReference");
 }
 //-----------------------------------------------------------------------------
 const la::petsc::KrylovSolver&
@@ -195,7 +204,7 @@ std::pair<int, bool> nls::petsc::NewtonSolver::solve(Vec x)
     _solver.set_operators(_matJ, _matJ);
 
   if (!_dx)
-    MatCreateVecs(_matJ, &_dx, nullptr);
+    check_petsc_error(MatCreateVecs(_matJ, &_dx, nullptr), "MatCreateVecs");
 
   // Start iterations
   while (!newton_converged and _iteration < max_it)
@@ -224,7 +233,7 @@ std::pair<int, bool> nls::petsc::NewtonSolver::solve(Vec x)
     if (_iteration == 1)
     {
       PetscReal _r = 0;
-      VecNorm(_dx, NORM_2, &_r);
+      check_petsc_error(VecNorm(_dx, NORM_2, &_r), "VecNorm");
       _residual0 = _r;
     }
 
diff --git a/python/dolfinx/fem/assemble.py b/python/dolfinx/fem/assemble.py
index 4a1b89e017d..9fc2dca9ae5 100644
--- a/python/dolfinx/fem/assemble.py
+++ b/python/dolfinx/fem/assemble.py
@@ -26,6 +26,25 @@
 from dolfinx.fem.function import FunctionSpace
 
 
+def _check_bcs_for_different_trial_test_spaces(a: Form, bcs: Sequence[DirichletBC]) -> None:
+    """Reject DirichletBCs for unsupported Petrov-Galerkin assembly."""
+    if not bcs or len(a.function_spaces) != 2:
+        return
+
+    test_space, trial_space = a.function_spaces
+    if test_space is trial_space:
+        return
+    if test_space.mesh is not trial_space.mesh:
+        return
+
+    for bc in bcs:
+        if test_space.contains(bc.function_space) or trial_space.contains(bc.function_space):
+            raise RuntimeError(
+                "Dirichlet boundary conditions for bilinear forms with different test "
+                "and trial spaces are not supported by this matrix assembly path."
+            )
+
+
 def pack_constants(
     form: Form | Sequence[Form],
 ) -> npt.NDArray | Sequence[npt.NDArray]:
@@ -288,6 +307,7 @@ def assemble_matrix(
         accumulated.
     """
     bcs = [] if bcs is None else bcs
+    _check_bcs_for_different_trial_test_spaces(a, bcs)
     A: la.MatrixCSR = create_matrix(a, block_mode)
     _assemble_matrix_csr(A, a, bcs, diag, constants, coeffs)
     return A
@@ -325,15 +345,17 @@ def _assemble_matrix_csr(
         The returned matrix is not finalised, i.e. ghost values are not
         accumulated.
     """
-    bcs = [] if bcs is None else [bc._cpp_object for bc in bcs]
+    bcs = [] if bcs is None else bcs
+    _check_bcs_for_different_trial_test_spaces(a, bcs)
+    _bcs = [bc._cpp_object for bc in bcs]
     constants = pack_constants(a) if constants is None else constants  # type: ignore[assignment]
     coeffs = pack_coefficients(a) if coeffs is None else coeffs  # type: ignore[assignment]
-    _cpp.fem.assemble_matrix(A._cpp_object, a._cpp_object, constants, coeffs, bcs)
+    _cpp.fem.assemble_matrix(A._cpp_object, a._cpp_object, constants, coeffs, _bcs)
 
     # If matrix is a 'diagonal'block, set diagonal entry for constrained
     # dofs
     if a.function_spaces[0] is a.function_spaces[1]:
-        _cpp.fem.insert_diagonal(A._cpp_object, a.function_spaces[0], bcs, diag)
+        _cpp.fem.insert_diagonal(A._cpp_object, a.function_spaces[0], _bcs, diag)
     return A
 
 
diff --git a/python/dolfinx/fem/petsc.py b/python/dolfinx/fem/petsc.py
index e5dee8c5cf1..21e4dc60c8f 100644
--- a/python/dolfinx/fem/petsc.py
+++ b/python/dolfinx/fem/petsc.py
@@ -49,7 +49,10 @@
 from dolfinx.cpp.fem.petsc import discrete_gradient as _discrete_gradient
 from dolfinx.cpp.fem.petsc import interpolation_matrix as _interpolation_matrix
 from dolfinx.fem import IntegralType, pack_coefficients, pack_constants
-from dolfinx.fem.assemble import _assemble_vector_array
+from dolfinx.fem.assemble import (
+    _assemble_vector_array,
+    _check_bcs_for_different_trial_test_spaces,
+)
 from dolfinx.fem.assemble import apply_lifting as _apply_lifting
 from dolfinx.fem.bcs import DirichletBC
 from dolfinx.fem.bcs import bcs_by_block as _bcs_by_block
@@ -395,8 +398,12 @@ def assemble_matrix(
     Returns:
         Matrix representing the bilinear form.
     """  # noqa: D301
+    bcs = [] if bcs is None else bcs
+    if not isinstance(a, Sequence):
+        _check_bcs_for_different_trial_test_spaces(a, bcs)
+
     A = create_matrix(a, kind)
-    assemble_matrix(A, a, bcs, diag, constants, coeffs)  # type: ignore[arg-type]
+    assemble_matrix(A, a, bcs, diag, constants, coeffs, _check_bcs=False)  # type: ignore[arg-type]
     return A
 
 
@@ -412,6 +419,8 @@ def _(
         | Sequence[Sequence[dict[tuple[IntegralType, int], npt.NDArray]]]
         | None
     ) = None,
+    *,
+    _check_bcs: bool = True,
 ) -> PETSc.Mat:  # type: ignore[name-defined]
     """Assemble bilinear form into a matrix.
 
@@ -422,6 +431,10 @@ def _(
     The returned matrix is not finalised, i.e. ghost values are not
     accumulated.
     """
+    bcs = [] if bcs is None else bcs
+    if _check_bcs and not isinstance(a, Sequence):
+        _check_bcs_for_different_trial_test_spaces(a, bcs)
+
     if A.getType() == PETSc.Mat.Type.NEST:  # type: ignore[attr-defined]
         if not isinstance(a, Sequence):
             raise ValueError("Must provide a sequence of forms when assembling a nest matrix")
@@ -431,7 +444,7 @@ def _(
             for j, (a_block, const, coeff) in enumerate(zip(a_row, const_row, coeff_row)):
                 if a_block is not None:
                     Asub = A.getNestSubMatrix(i, j)
-                    assemble_matrix(Asub, a_block, bcs, diag, const, coeff)  # type: ignore[arg-type]
+                    assemble_matrix(Asub, a_block, bcs, diag, const, coeff, _check_bcs=False)  # type: ignore[arg-type]
                 elif i == j:
                     for bc in bcs:
                         row_forms = [row_form for row_form in a_row if row_form is not None]
@@ -462,7 +475,7 @@ def _(
             [(Vsub.dofmaps[0].index_map, Vsub.dofmaps[0].index_map_bs) for Vsub in V[1]]  # type: ignore[union-attr]
         )
 
-        _bcs = [bc._cpp_object for bc in bcs] if bcs is not None else []
+        _bcs = [bc._cpp_object for bc in bcs]
         for i, a_row in enumerate(a):
             for j, a_sub in enumerate(a_row):
                 if a_sub is not None:
@@ -501,7 +514,7 @@ def _(
     else:  # Non-blocked
         constants = pack_constants(a) if constants is None else constants  # type: ignore[assignment]
         coeffs = pack_coefficients(a) if coeffs is None else coeffs  # type: ignore[assignment]
-        _bcs = [bc._cpp_object for bc in bcs] if bcs is not None else []
+        _bcs = [bc._cpp_object for bc in bcs]
         _cpp.fem.petsc.assemble_matrix(A, a._cpp_object, constants, coeffs, _bcs)
         if a.function_spaces[0] is a.function_spaces[1]:
             A.assemblyBegin(PETSc.Mat.AssemblyType.FLUSH)  # type: ignore[attr-defined]
diff --git a/python/test/unit/fem/test_assembler.py b/python/test/unit/fem/test_assembler.py
index 55d77b88a51..2f2afa2b03f 100644
--- a/python/test/unit/fem/test_assembler.py
+++ b/python/test/unit/fem/test_assembler.py
@@ -178,6 +178,26 @@ def nest_matrix_norm(A):
     return math.sqrt(norm)
 
 
+def _different_trial_test_space_form_and_bc():
+    msh = create_unit_square(MPI.COMM_WORLD, 4, 4)
+    U = functionspace(msh, ("Lagrange", 1))
+    V = functionspace(msh, ("Lagrange", 2))
+    u = ufl.TrialFunction(U)
+    v = ufl.TestFunction(V)
+    a = form(inner(u, v) * dx)
+
+    facets = locate_entities_boundary(msh, msh.topology.dim - 1, lambda x: np.isclose(x[0], 0.0))
+    dofs = locate_dofs_topological(V, msh.topology.dim - 1, facets)
+    bc = dirichletbc(default_scalar_type(0), dofs, V)
+    return a, bc
+
+
+def test_assembly_rejects_bcs_for_different_trial_test_spaces():
+    a, bc = _different_trial_test_space_form_and_bc()
+    with pytest.raises(RuntimeError, match="different test and trial spaces"):
+        assemble_matrix(a, bcs=[bc])
+
+
 @pytest.mark.petsc4py
 def test_vector_single_space_as_block():
     from dolfinx.fem.petsc import create_vector as petsc_create_vector
@@ -193,6 +213,14 @@ def test_vector_single_space_as_block():
 class TestPETScAssemblers:
     """Test PETSc-based assemblers for matrices and vectors."""
 
+    def test_assembly_rejects_bcs_for_different_trial_test_spaces(self):
+        """Test that PETSc assembly rejects same-domain Petrov-Galerkin BCs."""
+        from dolfinx.fem.petsc import assemble_matrix as petsc_assemble_matrix
+
+        a, bc = _different_trial_test_space_form_and_bc()
+        with pytest.raises(RuntimeError, match="different test and trial spaces"):
+            petsc_assemble_matrix(a, bcs=[bc])
+
     @pytest.mark.parametrize("mode", [GhostMode.none, GhostMode.shared_facet])
     def test_basic_assembly_petsc_matrixcsr(self, mode):
         """Test basic assembly of PETSc Mat and compare with MatrixCSR assembly."""
diff --git a/python/test/unit/io/test_adios2.py b/python/test/unit/io/test_adios2.py
index aeeb526a0c6..70a49a6c597 100644
--- a/python/test/unit/io/test_adios2.py
+++ b/python/test/unit/io/test_adios2.py
@@ -251,9 +251,9 @@ def test_vtx_reuse_mesh(self, tempdir, dim, simplex, reuse):
 
     def test_dg_0_data(self, tempdir):
         """Test that we can mix DG-0 and other Lagrange functions."""
+        adios2 = import_adios2()
         from dolfinx.io import VTXWriter
 
-        adios2 = import_adios2()
         mesh = generate_mesh(2, False)
         v = Function(functionspace(mesh, ("Lagrange", 2, (2,))))
         filename = Path(tempdir, "v.bp")