diff --git a/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx b/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx
index 881a0b7fc0..b3554a05c6 100644
--- a/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx
+++ b/src/invert/laplacexy/impls/petsc/laplacexy-petsc.cxx
@@ -90,15 +90,15 @@ LaplaceXYpetsc::LaplaceXYpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
   const int localN = localSize();
 
   // Create Vectors
-  VecCreate(comm, &xs);
-  VecSetSizes(xs, localN, PETSC_DETERMINE);
-  VecSetFromOptions(xs);
-  VecDuplicate(xs, &bs);
+  BOUT_DO_PETSC(VecCreate(comm, &xs));
+  BOUT_DO_PETSC(VecSetSizes(xs, localN, PETSC_DETERMINE));
+  BOUT_DO_PETSC(VecSetFromOptions(xs));
+  BOUT_DO_PETSC(VecDuplicate(xs, &bs));
 
   // Set size of Matrix on each processor to localN x localN
-  MatCreate(comm, &MatA);
-  MatSetSizes(MatA, localN, localN, PETSC_DETERMINE, PETSC_DETERMINE);
-  MatSetFromOptions(MatA);
+  BOUT_DO_PETSC(MatCreate(comm, &MatA));
+  BOUT_DO_PETSC(MatSetSizes(MatA, localN, localN, PETSC_DETERMINE, PETSC_DETERMINE));
+  BOUT_DO_PETSC(MatSetFromOptions(MatA));
 
   //////////////////////////////////////////////////
   // Specify local indices. This creates a mapping
@@ -187,16 +187,16 @@ LaplaceXYpetsc::LaplaceXYpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
     setPreallocationFiniteDifference(d_nnz, o_nnz);
   }
   // Pre-allocate
-  MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz);
-  MatSetUp(MatA);
+  BOUT_DO_PETSC(MatMPIAIJSetPreallocation(MatA, 0, d_nnz, 0, o_nnz));
+  BOUT_DO_PETSC(MatSetUp(MatA));
 
-  PetscFree(d_nnz);
-  PetscFree(o_nnz);
+  BOUT_DO_PETSC(PetscFree(d_nnz));
+  BOUT_DO_PETSC(PetscFree(o_nnz));
 
   // Determine which row/columns of the matrix are locally owned
   int Istart = 0;
   int Iend = 0;
-  MatGetOwnershipRange(MatA, &Istart, &Iend);
+  BOUT_DO_PETSC(MatGetOwnershipRange(MatA, &Istart, &Iend));
 
   // Convert indexXY from local index to global index
   indexXY += Istart;
@@ -209,14 +209,14 @@ LaplaceXYpetsc::LaplaceXYpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
   // Set up KSP
 
   // Declare KSP Context
-  KSPCreate(comm, &ksp);
+  BOUT_DO_PETSC(KSPCreate(comm, &ksp));
 
   // Configure Linear Solver
 
   const bool direct = (*opt)["direct"].doc("Use a direct LU solver").withDefault(false);
 
   if (direct) {
-    KSPGetPC(ksp, &pc);
+    BOUT_DO_PETSC(KSPGetPC(ksp, &pc));
     PCSetType(pc, PCLU);
 #if PETSC_VERSION_GE(3, 9, 0)
     PCFactorSetMatSolverType(pc, "mumps");
@@ -248,12 +248,12 @@ LaplaceXYpetsc::LaplaceXYpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
     const std::string pctype =
         (*opt)["pctype"].doc("Preconditioner type").withDefault("none");
 
-    KSPSetType(ksp, ksptype.c_str());
-    KSPSetTolerances(ksp, rtol, atol, dtol, maxits);
+    BOUT_DO_PETSC(KSPSetType(ksp, ksptype.c_str()));
+    BOUT_DO_PETSC(KSPSetTolerances(ksp, rtol, atol, dtol, maxits));
 
-    KSPSetInitialGuessNonzero(ksp, static_cast<PetscBool>(true));
+    BOUT_DO_PETSC(KSPSetInitialGuessNonzero(ksp, static_cast<PetscBool>(true)));
 
-    KSPGetPC(ksp, &pc);
+    BOUT_DO_PETSC(KSPGetPC(ksp, &pc));
     PCSetType(pc, pctype.c_str());
 
     if (pctype == "shell") {
@@ -264,9 +264,9 @@ LaplaceXYpetsc::LaplaceXYpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
       const bool rightprec =
           (*opt)["rightprec"].doc("Use right preconditioning?").withDefault(true);
       if (rightprec) {
-        KSPSetPCSide(ksp, PC_RIGHT); // Right preconditioning
+        BOUT_DO_PETSC(KSPSetPCSide(ksp, PC_RIGHT)); // Right preconditioning
       } else {
-        KSPSetPCSide(ksp, PC_LEFT); // Left preconditioning
+        BOUT_DO_PETSC(KSPSetPCSide(ksp, PC_LEFT)); // Left preconditioning
       }
     }
   }
@@ -578,10 +578,10 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
         int row = globalIndex(localmesh->xstart - 1, y);
         PetscScalar val = 0.5;
-        MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
         int col = globalIndex(localmesh->xstart, y);
-        MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
         // Preconditioner
         bcoef(y - localmesh->ystart, 0) = 0.5;
@@ -594,11 +594,11 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
         int row = globalIndex(localmesh->xstart - 1, y);
         PetscScalar val = 1.0;
-        MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
         int col = globalIndex(localmesh->xstart, y);
         val = -1.0;
-        MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
         // Preconditioner
         bcoef(y - localmesh->ystart, 0) = 1.0;
@@ -612,10 +612,10 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
     for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
       int row = globalIndex(localmesh->xend + 1, y);
       PetscScalar val = 0.5;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       int col = globalIndex(localmesh->xend, y);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
       // Preconditioner
       acoef(y - localmesh->ystart, localmesh->xend + 1 - xstart) = 0.5;
@@ -632,10 +632,10 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       }
       int row = globalIndex(it.ind, localmesh->ystart - 1);
       PetscScalar val = 0.5;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       int col = globalIndex(it.ind, localmesh->ystart);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
     }
 
     for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) {
@@ -645,10 +645,10 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       }
       int row = globalIndex(it.ind, localmesh->yend + 1);
       PetscScalar val = 0.5;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       int col = globalIndex(it.ind, localmesh->yend);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
     }
   } else if (y_bndry == "neumann") {
     // Neumann on Y boundaries
@@ -659,12 +659,12 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       }
       int row = globalIndex(it.ind, localmesh->ystart - 1);
       PetscScalar val = 1.0;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       val = -1.0;
       int col = globalIndex(it.ind, localmesh->ystart);
 
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
     }
 
     for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) {
@@ -674,11 +674,11 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       }
       int row = globalIndex(it.ind, localmesh->yend + 1);
       PetscScalar val = 1.0;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       val = -1.0;
       int col = globalIndex(it.ind, localmesh->yend);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
     }
   } else if (y_bndry == "free_o3") {
     // 'free_o3' extrapolating boundary condition on Y boundaries
@@ -689,19 +689,19 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       }
       int row = globalIndex(it.ind, localmesh->ystart - 1);
       PetscScalar val = 1.0;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       val = -3.0;
       int col = globalIndex(it.ind, localmesh->ystart);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
       val = 3.0;
       col = globalIndex(it.ind, localmesh->ystart + 1);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
       val = -1.0;
       col = globalIndex(it.ind, localmesh->ystart + 2);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
     }
 
     for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) {
@@ -711,19 +711,19 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
       }
       int row = globalIndex(it.ind, localmesh->yend + 1);
       PetscScalar val = 1.0;
-      MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
       val = -3.0;
       int col = globalIndex(it.ind, localmesh->yend);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
       val = 3.0;
       col = globalIndex(it.ind, localmesh->yend - 1);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
       val = -1.0;
       col = globalIndex(it.ind, localmesh->yend - 2);
-      MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
     }
   } else {
     throw BoutException("Unsupported option for y_bndry");
@@ -741,11 +741,11 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // f(xs-1,ys-1) = f(xs-1,ys)
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -1.0;
           int col = globalIndex(localmesh->xstart - 1, localmesh->ystart);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") {
           // 'free_o3' extrapolating boundary condition on Y boundaries
           // f(xs-1,ys-1) = 3*f(xs-1,ys) - 3*f(xs-1,ys+1) + f(xs-1,ys+2)
@@ -754,19 +754,19 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // what value to pass for the corner
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xstart - 1, localmesh->ystart - 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -3.0;
           int col = globalIndex(localmesh->xstart - 1, localmesh->ystart);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = 3.0;
           col = globalIndex(localmesh->xstart - 1, localmesh->ystart + 1);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = -1.0;
           col = globalIndex(localmesh->xstart - 1, localmesh->ystart + 2);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else {
           throw BoutException("Unsupported option for y_bndry");
         }
@@ -777,11 +777,11 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // f(xe+1,ys-1) = f(xe+1,ys)
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xend + 1, localmesh->ystart - 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -1.0;
           int col = globalIndex(localmesh->xend + 1, localmesh->ystart);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") {
           // 'free_o3' extrapolating boundary condition on Y boundaries
           // f(xe+1,ys-1) = 3*f(xe+1,ys) - 3*f(xe+1,ys+1) + f(xe+1,ys+2)
@@ -790,19 +790,19 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // what value to pass for the corner
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xend + 1, localmesh->ystart - 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -3.0;
           int col = globalIndex(localmesh->xend + 1, localmesh->ystart);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = 3.0;
           col = globalIndex(localmesh->xend + 1, localmesh->ystart + 1);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = -1.0;
           col = globalIndex(localmesh->xend + 1, localmesh->ystart + 2);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else {
           throw BoutException("Unsupported option for y_bndry");
         }
@@ -815,11 +815,11 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // f(xs-1,ys-1) = f(xs-1,ys)
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xstart - 1, localmesh->yend + 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -1.0;
           int col = globalIndex(localmesh->xstart - 1, localmesh->yend);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") {
           // 'free_o3' extrapolating boundary condition on Y boundaries
           // f(xs-1,ys-1) = 3*f(xs-1,ys) - 3*f(xs-1,ys+1) + f(xs-1,ys+2)
@@ -828,19 +828,19 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // what value to pass for the corner
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xstart - 1, localmesh->yend + 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -3.0;
           int col = globalIndex(localmesh->xstart - 1, localmesh->yend);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = 3.0;
           col = globalIndex(localmesh->xstart - 1, localmesh->yend - 1);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = -1.0;
           col = globalIndex(localmesh->xstart - 1, localmesh->yend - 2);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else {
           throw BoutException("Unsupported option for y_bndry");
         }
@@ -851,11 +851,11 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // f(xe+1,ys-1) = f(xe+1,ys)
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xend + 1, localmesh->yend + 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -1.0;
           int col = globalIndex(localmesh->xend + 1, localmesh->yend);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else if (y_bndry == "free_o3" or y_bndry == "dirichlet") {
           // 'free_o3' extrapolating boundary condition on Y boundaries
           // f(xe+1,ys-1) = 3*f(xe+1,ys) - 3*f(xe+1,ys+1) + f(xe+1,ys+2)
@@ -864,19 +864,19 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
           // what value to pass for the corner
           PetscScalar val = 1.0;
           int row = globalIndex(localmesh->xend + 1, localmesh->yend + 1);
-          MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           val = -3.0;
           int col = globalIndex(localmesh->xend + 1, localmesh->yend);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = 3.0;
           col = globalIndex(localmesh->xend + 1, localmesh->yend - 1);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           val = -1.0;
           col = globalIndex(localmesh->xend + 1, localmesh->yend - 2);
-          MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
         } else {
           throw BoutException("Unsupported option for y_bndry");
         }
@@ -885,14 +885,14 @@ void LaplaceXYpetsc::setCoefs(const Field2D& A, const Field2D& B) {
   }
 
   // Assemble Matrix
-  MatAssemblyBegin(MatA, MAT_FINAL_ASSEMBLY);
-  MatAssemblyEnd(MatA, MAT_FINAL_ASSEMBLY);
+  BOUT_DO_PETSC(MatAssemblyBegin(MatA, MAT_FINAL_ASSEMBLY));
+  BOUT_DO_PETSC(MatAssemblyEnd(MatA, MAT_FINAL_ASSEMBLY));
 
   // Set the operator
 #if PETSC_VERSION_GE(3, 5, 0)
-  KSPSetOperators(ksp, MatA, MatA);
+  BOUT_DO_PETSC(KSPSetOperators(ksp, MatA, MatA));
 #else
-  KSPSetOperators(ksp, MatA, MatA, DIFFERENT_NONZERO_PATTERN);
+  BOUT_DO_PETSC(KSPSetOperators(ksp, MatA, MatA, DIFFERENT_NONZERO_PATTERN));
 #endif
 
   // Set coefficients for preconditioner
@@ -981,24 +981,24 @@ void LaplaceXYpetsc::setMatrixElementsFiniteVolume(const Field2D& A, const Field
       int row = globalIndex(x, y);
 
       // Set the centre (diagonal)
-      MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES));
 
       // X + 1
       int col = globalIndex(x + 1, y);
-      MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES));
 
       // X - 1
       col = globalIndex(x - 1, y);
-      MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES));
 
       if (include_y_derivs) {
         // Y + 1
         col = globalIndex(x, y + 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES));
 
         // Y - 1
         col = globalIndex(x, y - 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES));
       }
     }
   }
@@ -1113,40 +1113,40 @@ void LaplaceXYpetsc::setMatrixElementsFiniteDifference(const Field2D& A,
       int row = globalIndex(x, y);
 
       // Set the centre (diagonal)
-      MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &row, &c, INSERT_VALUES));
 
       // X + 1
       int col = globalIndex(x + 1, y);
-      MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xp, INSERT_VALUES));
 
       // X - 1
       col = globalIndex(x - 1, y);
-      MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES);
+      BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xm, INSERT_VALUES));
 
       if (include_y_derivs) {
         // Y + 1
         col = globalIndex(x, y + 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &yp, INSERT_VALUES));
 
         // Y - 1
         col = globalIndex(x, y - 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &ym, INSERT_VALUES));
 
         // X + 1, Y + 1
         col = globalIndex(x + 1, y + 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &xpyp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xpyp, INSERT_VALUES));
 
         // X + 1, Y - 1
         col = globalIndex(x + 1, y - 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &xpym, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xpym, INSERT_VALUES));
 
         // X - 1, Y + 1
         col = globalIndex(x - 1, y + 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &xmyp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xmyp, INSERT_VALUES));
 
         // X - 1, Y - 1
         col = globalIndex(x - 1, y - 1);
-        MatSetValues(MatA, 1, &row, 1, &col, &xmym, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(MatA, 1, &row, 1, &col, &xmym, INSERT_VALUES));
       }
     }
   }
@@ -1158,12 +1158,12 @@ LaplaceXYpetsc::~LaplaceXYpetsc() {
 
   if (!is_finalised) {
     // PetscFinalize may already have destroyed this object
-    KSPDestroy(&ksp);
+    BOUT_DO_PETSC(KSPDestroy(&ksp));
   }
 
-  VecDestroy(&xs);
-  VecDestroy(&bs);
-  MatDestroy(&MatA);
+  BOUT_DO_PETSC(VecDestroy(&xs));
+  BOUT_DO_PETSC(VecDestroy(&bs));
+  BOUT_DO_PETSC(MatDestroy(&MatA));
 }
 
 Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
@@ -1181,10 +1181,10 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
       int ind = globalIndex(x, y);
 
       PetscScalar val = x0(x, y);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       val = rhs(x, y);
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
   }
 
@@ -1195,18 +1195,18 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
   }
 
   // Assemble RHS Vector
-  VecAssemblyBegin(bs);
-  VecAssemblyEnd(bs);
+  BOUT_DO_PETSC(VecAssemblyBegin(bs));
+  BOUT_DO_PETSC(VecAssemblyEnd(bs));
 
   // Assemble Trial Solution Vector
-  VecAssemblyBegin(xs);
-  VecAssemblyEnd(xs);
+  BOUT_DO_PETSC(VecAssemblyBegin(xs));
+  BOUT_DO_PETSC(VecAssemblyEnd(xs));
 
   // Solve the system
-  KSPSolve(ksp, bs, xs);
+  BOUT_DO_PETSC(KSPSolve(ksp, bs, xs));
 
   KSPConvergedReason reason;
-  KSPGetConvergedReason(ksp, &reason);
+  BOUT_DO_PETSC(KSPGetConvergedReason(ksp, &reason));
 
   if (reason <= 0) {
     throw BoutException("LaplaceXYpetsc failed to converge. Reason {} ({:d})",
@@ -1218,7 +1218,7 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
     n_calls++;
 
     int iterations = 0;
-    KSPGetIterationNumber(ksp, &iterations);
+    BOUT_DO_PETSC(KSPGetIterationNumber(ksp, &iterations));
 
     average_iterations = BoutReal(n_calls - 1) / BoutReal(n_calls) * average_iterations
                          + BoutReal(iterations) / BoutReal(n_calls);
@@ -1236,7 +1236,7 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
       int ind = globalIndex(x, y);
 
       PetscScalar val;
-      VecGetValues(xs, 1, &ind, &val);
+      BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
       result(x, y) = val;
     }
   }
@@ -1246,7 +1246,7 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
     for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
       int ind = globalIndex(localmesh->xstart - 1, y);
       PetscScalar val;
-      VecGetValues(xs, 1, &ind, &val);
+      BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
       for (int x = localmesh->xstart - 1; x >= 0; x--) {
         result(x, y) = val;
       }
@@ -1258,7 +1258,7 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
     for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
       int ind = globalIndex(localmesh->xend + 1, y);
       PetscScalar val;
-      VecGetValues(xs, 1, &ind, &val);
+      BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
       for (int x = localmesh->xend + 1; x < localmesh->LocalNx; x++) {
         result(x, y) = val;
       }
@@ -1278,7 +1278,7 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
     }
     int ind = globalIndex(it.ind, localmesh->ystart - 1);
     PetscScalar val;
-    VecGetValues(xs, 1, &ind, &val);
+    BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
     for (int y = localmesh->ystart - 1; y >= 0; y--) {
       result(it.ind, y) = val;
     }
@@ -1297,7 +1297,7 @@ Field2D LaplaceXYpetsc::solve(const Field2D& rhs, const Field2D& x0) {
     }
     int ind = globalIndex(it.ind, localmesh->yend + 1);
     PetscScalar val;
-    VecGetValues(xs, 1, &ind, &val);
+    BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
     for (int y = localmesh->yend + 1; y < localmesh->LocalNy; y++) {
       result(it.ind, y) = val;
     }
@@ -1315,10 +1315,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
         int ind = globalIndex(localmesh->xstart - 1, y);
 
         PetscScalar val = x0(localmesh->xstart - 1, y);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         val = 0.5 * (x0(localmesh->xstart - 1, y) + x0(localmesh->xstart, y));
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     } else {
       // Inner X boundary (Neumann)
@@ -1326,10 +1326,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
         int ind = globalIndex(localmesh->xstart - 1, y);
 
         PetscScalar val = x0(localmesh->xstart - 1, y);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         val = 0.0; // x0(localmesh->xstart-1,y) - x0(localmesh->xstart,y);
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
   }
@@ -1340,10 +1340,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
       int ind = globalIndex(localmesh->xend + 1, y);
 
       PetscScalar val = x0(localmesh->xend + 1, y);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       val = 0.5 * (x0(localmesh->xend, y) + x0(localmesh->xend + 1, y));
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
   }
 
@@ -1356,10 +1356,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
       int ind = globalIndex(it.ind, localmesh->ystart - 1);
 
       PetscScalar val = x0(it.ind, localmesh->ystart - 1);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       val = 0.5 * (x0(it.ind, localmesh->ystart - 1) + x0(it.ind, localmesh->ystart));
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
 
     for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) {
@@ -1370,10 +1370,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
       int ind = globalIndex(it.ind, localmesh->yend + 1);
 
       PetscScalar val = x0(it.ind, localmesh->yend + 1);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       val = 0.5 * (x0(it.ind, localmesh->yend + 1) + x0(it.ind, localmesh->yend));
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
   } else if (y_bndry == "neumann" or y_bndry == "free_o3") {
     // Y boundaries Neumann
@@ -1385,10 +1385,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
       int ind = globalIndex(it.ind, localmesh->ystart - 1);
 
       PetscScalar val = x0(it.ind, localmesh->ystart - 1);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       val = 0.0;
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
 
     for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) {
@@ -1399,10 +1399,10 @@ void LaplaceXYpetsc::solveFiniteVolume(const Field2D& x0) {
       int ind = globalIndex(it.ind, localmesh->yend + 1);
 
       PetscScalar val = x0(it.ind, localmesh->yend + 1);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       val = 0.0;
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
   } else {
     throw BoutException("Unsupported option for y_bndry");
@@ -1422,11 +1422,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         // For the boundary value of the initial guess, use the value that would be set by
         // applying the boundary condition to the initial guess
         PetscScalar val = 2. * x0(localmesh->xstart - 1, y) - x0(localmesh->xstart, y);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Pass the value from boundary cell of x0 as the boundary condition to the rhs
         val = x0(localmesh->xstart - 1, y);
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     } else {
       // Inner X boundary (Neumann)
@@ -1436,10 +1436,10 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         // For the boundary value of the initial guess, use the value that would be set by
         // applying the boundary condition to the initial guess
         PetscScalar val = x0(localmesh->xstart, y);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
   }
@@ -1452,11 +1452,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       // For the boundary value of the initial guess, use the value that would be set by
       // applying the boundary condition to the initial guess
       PetscScalar val = 2. * x0(localmesh->xend + 1, y) - x0(localmesh->xend, y);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Pass the value from boundary cell of x0 as the boundary condition to the rhs
       val = x0(localmesh->xend + 1, y);
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
   }
 
@@ -1472,11 +1472,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       // applying the boundary condition to the initial guess
       PetscScalar val =
           2. * x0(it.ind, localmesh->ystart - 1) - x0(it.ind, localmesh->ystart);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Pass the value from boundary cell of x0 as the boundary condition to the rhs
       val = x0(it.ind, localmesh->ystart - 1);
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
 
     for (RangeIterator it = localmesh->iterateBndryUpperY(); !it.isDone(); it++) {
@@ -1490,11 +1490,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       // applying the boundary condition to the initial guess
       PetscScalar val =
           2. * x0(it.ind, localmesh->yend + 1) - x0(it.ind, localmesh->yend);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Pass the value from boundary cell of x0 as the boundary condition to the rhs
       val = x0(it.ind, localmesh->yend + 1);
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
 
     // Use free_o3 for the corner boundary cells
@@ -1507,11 +1507,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->ystart)
                           - 3. * x0(localmesh->xstart - 1, localmesh->ystart + 1)
                           + x0(localmesh->xstart - 1, localmesh->ystart + 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
       if (localmesh->lastX()) {
         int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1);
@@ -1521,11 +1521,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->ystart)
                           - 3. * x0(localmesh->xend + 1, localmesh->ystart + 1)
                           + x0(localmesh->xend + 1, localmesh->ystart + 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
     if (localmesh->hasBndryUpperY()) {
@@ -1537,11 +1537,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->yend)
                           - 3. * x0(localmesh->xstart - 1, localmesh->yend - 1)
                           + x0(localmesh->xstart - 1, localmesh->yend - 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
       if (localmesh->lastX()) {
         int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1);
@@ -1551,11 +1551,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->yend)
                           - 3. * x0(localmesh->xend + 1, localmesh->yend - 1)
                           + x0(localmesh->xend + 1, localmesh->yend - 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
   } else if (y_bndry == "neumann") {
@@ -1570,11 +1570,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       // For the boundary value of the initial guess, use the value that would be set by
       // applying the boundary condition to the initial guess
       PetscScalar val = x0(it.ind, localmesh->ystart);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Set the value for the rhs of the boundary condition to zero
       val = 0.0;
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
     if (localmesh->hasBndryLowerY()) {
       if (localmesh->firstX()) {
@@ -1583,11 +1583,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         // For the boundary value of the initial guess, use the value that would be set by
         // applying the boundary condition to the initial guess
         PetscScalar val = x0(localmesh->xstart - 1, localmesh->ystart);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
       if (localmesh->lastX()) {
         int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1);
@@ -1595,11 +1595,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         // For the boundary value of the initial guess, use the value that would be set by
         // applying the boundary condition to the initial guess
         PetscScalar val = x0(localmesh->xend + 1, localmesh->ystart);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
 
@@ -1613,11 +1613,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       // For the boundary value of the initial guess, use the value that would be set by
       // applying the boundary condition to the initial guess
       PetscScalar val = x0(it.ind, localmesh->yend);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Set the value for the rhs of the boundary condition to zero
       val = 0.0;
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
     if (localmesh->hasBndryUpperY()) {
       if (localmesh->firstX()) {
@@ -1626,11 +1626,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         // For the boundary value of the initial guess, use the value that would be set by
         // applying the boundary condition to the initial guess
         PetscScalar val = x0(localmesh->xstart - 1, localmesh->yend);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
       if (localmesh->lastX()) {
         int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1);
@@ -1638,11 +1638,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         // For the boundary value of the initial guess, use the value that would be set by
         // applying the boundary condition to the initial guess
         PetscScalar val = x0(localmesh->xend + 1, localmesh->yend);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
   } else if (y_bndry == "free_o3") {
@@ -1659,11 +1659,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       PetscScalar val = 3. * x0(it.ind, localmesh->ystart)
                         - 3. * x0(it.ind, localmesh->ystart + 1)
                         + x0(it.ind, localmesh->ystart + 2);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Set the value for the rhs of the boundary condition to zero
       val = 0.0;
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
     if (localmesh->hasBndryLowerY()) {
       if (localmesh->firstX()) {
@@ -1674,11 +1674,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->ystart)
                           - 3. * x0(localmesh->xstart - 1, localmesh->ystart + 1)
                           + x0(localmesh->xstart - 1, localmesh->ystart + 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
       if (localmesh->lastX()) {
         int ind = globalIndex(localmesh->xend + 1, localmesh->ystart - 1);
@@ -1688,11 +1688,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->ystart)
                           - 3. * x0(localmesh->xend + 1, localmesh->ystart + 1)
                           + x0(localmesh->xend + 1, localmesh->ystart + 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
 
@@ -1708,11 +1708,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
       PetscScalar val = 3. * x0(it.ind, localmesh->yend)
                         - 3. * x0(it.ind, localmesh->yend - 1)
                         + x0(it.ind, localmesh->yend - 2);
-      VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
       // Set the value for the rhs of the boundary condition to zero
       val = 0.0;
-      VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
     }
     if (localmesh->hasBndryUpperY()) {
       if (localmesh->firstX()) {
@@ -1723,11 +1723,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xstart - 1, localmesh->yend)
                           - 3. * x0(localmesh->xstart - 1, localmesh->yend - 1)
                           + x0(localmesh->xstart - 1, localmesh->yend - 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
       if (localmesh->lastX()) {
         int ind = globalIndex(localmesh->xend + 1, localmesh->yend + 1);
@@ -1737,11 +1737,11 @@ void LaplaceXYpetsc::solveFiniteDifference(const Field2D& x0) {
         PetscScalar val = 3. * x0(localmesh->xend + 1, localmesh->yend)
                           - 3. * x0(localmesh->xend + 1, localmesh->yend - 1)
                           + x0(localmesh->xend + 1, localmesh->yend - 2);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         // Set the value for the rhs of the boundary condition to zero
         val = 0.0;
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
       }
     }
   } else {
@@ -1764,8 +1764,8 @@ int LaplaceXYpetsc::precon(Vec input, Vec result) {
     }
     const int ind = globalIndex(itdwn.ind, localmesh->ystart - 1);
     PetscScalar val;
-    VecGetValues(input, 1, &ind, &val);
-    VecSetValues(result, 1, &ind, &val, INSERT_VALUES);
+    BOUT_DO_PETSC(VecGetValues(input, 1, &ind, &val));
+    BOUT_DO_PETSC(VecSetValues(result, 1, &ind, &val, INSERT_VALUES));
   }
   for (auto itup = localmesh->iterateBndryUpperY(); !itup.isDone(); itup++) {
     // Should not go into corner cells, LaplaceXYpetsc stencil does not include them
@@ -1774,8 +1774,8 @@ int LaplaceXYpetsc::precon(Vec input, Vec result) {
     }
     const int ind = globalIndex(itup.ind, localmesh->yend + 1);
     PetscScalar val;
-    VecGetValues(input, 1, &ind, &val);
-    VecSetValues(result, 1, &ind, &val, INSERT_VALUES);
+    BOUT_DO_PETSC(VecGetValues(input, 1, &ind, &val));
+    BOUT_DO_PETSC(VecSetValues(result, 1, &ind, &val, INSERT_VALUES));
   }
 
   // Load vector x into bvals array
@@ -1783,7 +1783,7 @@ int LaplaceXYpetsc::precon(Vec input, Vec result) {
     for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
       const int ind = globalIndex(x, y);
       PetscScalar val;
-      VecGetValues(input, 1, &ind, &val);
+      BOUT_DO_PETSC(VecGetValues(input, 1, &ind, &val));
       bvals(y - localmesh->ystart, x - xstart) = val;
     }
   }
@@ -1796,11 +1796,11 @@ int LaplaceXYpetsc::precon(Vec input, Vec result) {
     for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
       const int ind = globalIndex(x, y);
       PetscScalar val = xvals(y - localmesh->ystart, x - xstart);
-      VecSetValues(result, 1, &ind, &val, INSERT_VALUES);
+      BOUT_DO_PETSC(VecSetValues(result, 1, &ind, &val, INSERT_VALUES));
     }
   }
-  VecAssemblyBegin(result);
-  VecAssemblyEnd(result);
+  BOUT_DO_PETSC(VecAssemblyBegin(result));
+  BOUT_DO_PETSC(VecAssemblyEnd(result));
   return 0;
 }
 
diff --git a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx
index 3e55262b6d..93e2004720 100644
--- a/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx
+++ b/src/invert/laplacexz/impls/petsc/laplacexz-petsc.cxx
@@ -9,6 +9,7 @@
  */
 
 #include "bout/build_defines.hxx"
+#include "bout/petsclib.hxx"
 
 #include "laplacexz-petsc.hxx"
 
@@ -167,10 +168,10 @@ LaplaceXZpetsc::LaplaceXZpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
   }
 
   // Create Vectors
-  VecCreate(comm, &xs);
-  VecSetSizes(xs, localN, PETSC_DETERMINE);
-  VecSetFromOptions(xs);
-  VecDuplicate(xs, &bs);
+  BOUT_DO_PETSC(VecCreate(comm, &xs));
+  BOUT_DO_PETSC(VecSetSizes(xs, localN, PETSC_DETERMINE));
+  BOUT_DO_PETSC(VecSetFromOptions(xs));
+  BOUT_DO_PETSC(VecDuplicate(xs, &bs));
 
   for (int y = localmesh->ystart; y <= localmesh->yend; y++) {
     YSlice data;
@@ -178,9 +179,10 @@ LaplaceXZpetsc::LaplaceXZpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
     data.yindex = y;
 
     // Set size of Matrix on each processor to localN x localN
-    MatCreate(comm, &data.MatA);
-    MatSetSizes(data.MatA, localN, localN, PETSC_DETERMINE, PETSC_DETERMINE);
-    MatSetFromOptions(data.MatA);
+    BOUT_DO_PETSC(MatCreate(comm, &data.MatA));
+    BOUT_DO_PETSC(
+        MatSetSizes(data.MatA, localN, localN, PETSC_DETERMINE, PETSC_DETERMINE));
+    BOUT_DO_PETSC(MatSetFromOptions(data.MatA));
 
     //////////////////////////////////////////////////
     // Pre-allocate space for matrix elements
@@ -223,28 +225,28 @@ LaplaceXZpetsc::LaplaceXZpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
       }
     }
 
-    MatMPIAIJSetPreallocation(data.MatA, 0, d_nnz, 0, o_nnz);
-    MatSetUp(data.MatA);
-    PetscFree(d_nnz);
-    PetscFree(o_nnz);
+    BOUT_DO_PETSC(MatMPIAIJSetPreallocation(data.MatA, 0, d_nnz, 0, o_nnz));
+    BOUT_DO_PETSC(MatSetUp(data.MatA));
+    BOUT_DO_PETSC(PetscFree(d_nnz));
+    BOUT_DO_PETSC(PetscFree(o_nnz));
 
     //////////////////////////////////////////////////
     // Declare KSP Context
-    KSPCreate(comm, &data.ksp);
+    BOUT_DO_PETSC(KSPCreate(comm, &data.ksp));
 
     // Set KSP type
-    KSPSetType(data.ksp, ksptype.c_str());
+    BOUT_DO_PETSC(KSPSetType(data.ksp, ksptype.c_str()));
     // Set KSP tolerances
-    KSPSetTolerances(data.ksp, rtol, atol, dtol, maxits);
+    BOUT_DO_PETSC(KSPSetTolerances(data.ksp, rtol, atol, dtol, maxits));
 
     // Set KSP preconditioner
     PC pc;
-    KSPGetPC(data.ksp, &pc);
-    PCSetType(pc, pctype.c_str());
+    BOUT_DO_PETSC(KSPGetPC(data.ksp, &pc));
+    BOUT_DO_PETSC(PCSetType(pc, pctype.c_str()));
 #if PETSC_VERSION_GE(3, 9, 0)
-    PCFactorSetMatSolverType(pc, factor_package.c_str());
+    BOUT_DO_PETSC(PCFactorSetMatSolverType(pc, factor_package.c_str()));
 #else
-    PCFactorSetMatSolverPackage(pc, factor_package.c_str());
+    BOUT_DO_PETSC(PCFactorSetMatSolverPackage(pc, factor_package.c_str()));
 #endif
 
     lib.setOptionsFromInputFile(data.ksp);
@@ -257,22 +259,22 @@ LaplaceXZpetsc::LaplaceXZpetsc(Mesh* m, Options* opt, const CELL_LOC loc)
 LaplaceXZpetsc::~LaplaceXZpetsc() {
 
   PetscBool petsc_is_finalised;
-  PetscFinalized(&petsc_is_finalised);
+  BOUT_DO_PETSC(PetscFinalized(&petsc_is_finalised));
 
   for (auto& it : slice) {
-    MatDestroy(&it.MatA);
+    BOUT_DO_PETSC(MatDestroy(&it.MatA));
     if (coefs_set) {
-      MatDestroy(&it.MatP);
+      BOUT_DO_PETSC(MatDestroy(&it.MatP));
     }
 
     if (!petsc_is_finalised) {
       // PetscFinalize may already have destroyed this object
-      KSPDestroy(&it.ksp);
+      BOUT_DO_PETSC(KSPDestroy(&it.ksp));
     }
   }
 
-  VecDestroy(&bs);
-  VecDestroy(&xs);
+  BOUT_DO_PETSC(VecDestroy(&bs));
+  BOUT_DO_PETSC(VecDestroy(&xs));
 }
 
 void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
@@ -316,7 +318,7 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
     int y = it.yindex;
 
     int Istart, Iend;
-    MatGetOwnershipRange(it.MatA, &Istart, &Iend);
+    BOUT_DO_PETSC(MatGetOwnershipRange(it.MatA, &Istart, &Iend));
 
     ////////////////////////////////////////////////
     // Inner X boundary (see note about BC in LaplaceXZ constructor)
@@ -329,11 +331,11 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
          */
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row + (localmesh->LocalNz); // +1 in X
           val = -1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -341,11 +343,11 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         // Setting BC from x0
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row + (localmesh->LocalNz); // +1 in X
           val = 0.0;
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -353,11 +355,11 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         // Setting BC from b
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row + (localmesh->LocalNz); // +1 in X
           val = 0.0;
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -368,10 +370,10 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
          */
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 0.5;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row + (localmesh->LocalNz); // +1 in X
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -597,15 +599,15 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         //
 
         // Set the centre (diagonal)
-        MatSetValues(it.MatA, 1, &row, 1, &row, &c, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &c, INSERT_VALUES));
 
         // X + 1
         int col = row + (localmesh->LocalNz);
-        MatSetValues(it.MatA, 1, &row, 1, &col, &xp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &xp, INSERT_VALUES));
 
         // X - 1
         col = row - (localmesh->LocalNz);
-        MatSetValues(it.MatA, 1, &row, 1, &col, &xm, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &xm, INSERT_VALUES));
 
         // Z + 1
         col = row + 1;
@@ -613,30 +615,30 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
           col -= localmesh->LocalNz; // Wrap around
         }
 
-        MatSetValues(it.MatA, 1, &row, 1, &col, &zp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &zp, INSERT_VALUES));
 
         // X + 1, Z + 1
         const int xpzp_col = col + localmesh->LocalNz;
-        MatSetValues(it.MatA, 1, &row, 1, &xpzp_col, &xpzp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &xpzp_col, &xpzp, INSERT_VALUES));
 
         // X - 1, Z + 1
         const int xmzp_col = col - localmesh->LocalNz;
-        MatSetValues(it.MatA, 1, &row, 1, &xmzp_col, &xmzp, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &xmzp_col, &xmzp, INSERT_VALUES));
 
         // Z - 1
         col = row - 1;
         if (z == 0) {
           col += localmesh->LocalNz; // Wrap around
         }
-        MatSetValues(it.MatA, 1, &row, 1, &col, &zm, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &zm, INSERT_VALUES));
 
         // X + 1, Z - 1
         const int xpzm_col = col + localmesh->LocalNz;
-        MatSetValues(it.MatA, 1, &row, 1, &xpzm_col, &xpzm, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &xpzm_col, &xpzm, INSERT_VALUES));
 
         // X - 1, Z - 1
         const int xmzm_col = col - localmesh->LocalNz;
-        MatSetValues(it.MatA, 1, &row, 1, &xmzm_col, &xmzm, INSERT_VALUES);
+        BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &xmzm_col, &xmzm, INSERT_VALUES));
 
         row++;
       }
@@ -649,11 +651,11 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         // Neumann 0
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row - (localmesh->LocalNz); // -1 in X
           val = -1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -661,11 +663,11 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         // Setting BC from x0
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row - (localmesh->LocalNz); // -1 in X
           val = 0.0;
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -673,11 +675,11 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         // Setting BC from b
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           PetscScalar val = 1.0;
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row - (localmesh->LocalNz); // -1 in X
           val = 0.0;
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -686,10 +688,10 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
         PetscScalar val = 0.5;
 
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
-          MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &row, &val, INSERT_VALUES));
 
           int col = row - (localmesh->LocalNz); // -1 in X
-          MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(MatSetValues(it.MatA, 1, &row, 1, &col, &val, INSERT_VALUES));
 
           row++;
         }
@@ -699,8 +701,8 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
     ASSERT1(row == Iend); // Check that row is currently on the last row
 
     // Assemble Matrix
-    MatAssemblyBegin(it.MatA, MAT_FINAL_ASSEMBLY);
-    MatAssemblyEnd(it.MatA, MAT_FINAL_ASSEMBLY);
+    BOUT_DO_PETSC(MatAssemblyBegin(it.MatA, MAT_FINAL_ASSEMBLY));
+    BOUT_DO_PETSC(MatAssemblyEnd(it.MatA, MAT_FINAL_ASSEMBLY));
   }
 
   // Increase reuse count
@@ -714,9 +716,9 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
       // Copy matrix into preconditioner
       if (coefs_set) {
         // Preconditioner already set
-        MatDestroy(&it.MatP);
+        BOUT_DO_PETSC(MatDestroy(&it.MatP));
       }
-      MatConvert(it.MatA, MATSAME, MAT_INITIAL_MATRIX, &it.MatP);
+      BOUT_DO_PETSC(MatConvert(it.MatA, MATSAME, MAT_INITIAL_MATRIX, &it.MatP));
     }
 
     // Set operators
@@ -724,9 +726,9 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
 
       // Note: This is a hack to force update of the preconditioner matrix
 #if PETSC_VERSION_GE(3, 5, 0)
-      KSPSetOperators(it.ksp, it.MatA, it.MatP);
+      BOUT_DO_PETSC(KSPSetOperators(it.ksp, it.MatA, it.MatP));
 #else
-      KSPSetOperators(it.ksp, it.MatA, it.MatP, SAME_NONZERO_PATTERN);
+      BOUT_DO_PETSC(KSPSetOperators(it.ksp, it.MatA, it.MatP, SAME_NONZERO_PATTERN));
 #endif
     }
   } else {
@@ -734,9 +736,9 @@ void LaplaceXZpetsc::setCoefs(const Field3D& Ain, const Field3D& Bin) {
       /// Reuse the preconditioner, even if the operator changes
 
 #if PETSC_VERSION_GE(3, 5, 0)
-      KSPSetReusePreconditioner(it.ksp, PETSC_TRUE);
+      BOUT_DO_PETSC(KSPSetReusePreconditioner(it.ksp, PETSC_TRUE));
 #else
-      KSPSetOperators(it.ksp, it.MatA, it.MatP, SAME_PRECONDITIONER);
+      BOUT_DO_PETSC(KSPSetOperators(it.ksp, it.MatA, it.MatP, SAME_PRECONDITIONER));
 #endif
     }
   }
@@ -780,12 +782,12 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
     int y = it.yindex;
 
     /// Specify non-zero starting guess for solution (from input x0)
-    KSPSetInitialGuessNonzero(it.ksp, PETSC_TRUE);
+    BOUT_DO_PETSC(KSPSetInitialGuessNonzero(it.ksp, PETSC_TRUE));
 
     //////////////////////////
     // Load initial guess x0 into xs and rhs into bs
     int Istart, Iend;
-    MatGetOwnershipRange(it.MatA, &Istart, &Iend);
+    BOUT_DO_PETSC(MatGetOwnershipRange(it.MatA, &Istart, &Iend));
 
     // Starting index
     int ind = Istart;
@@ -797,11 +799,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xstart - 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = 0.0;
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
           ind++;
         }
       } else if (inner_boundary_flags & INVERT_SET) {
@@ -809,11 +811,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xstart - 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = x0(localmesh->xstart, y, z);
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
           ind++;
         }
       } else if (inner_boundary_flags & INVERT_RHS) {
@@ -821,11 +823,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xstart - 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = b(localmesh->xstart, y, z);
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
           ind++;
         }
       } else {
@@ -833,11 +835,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xstart - 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = x0(localmesh->xstart - 1, y, z) - x0(localmesh->xstart, y, z);
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
           ind++;
         }
       }
@@ -847,10 +849,10 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
     for (int x = localmesh->xstart; x <= localmesh->xend; x++) {
       for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
         PetscScalar val = x0(x, y, z);
-        VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
         val = b(x, y, z);
-        VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+        BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
         ind++;
       }
     }
@@ -862,11 +864,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xend + 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = 0.0;
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
 
           ind++;
         }
@@ -875,11 +877,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xend + 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = x0(localmesh->xend + 1, y, z);
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
 
           ind++;
         }
@@ -888,11 +890,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xend + 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = b(localmesh->xend + 1, y, z);
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
 
           ind++;
         }
@@ -901,11 +903,11 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
         for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
           // Setting the initial guess x0
           PetscScalar val = x0(localmesh->xend + 1, y, z);
-          VecSetValues(xs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(xs, 1, &ind, &val, INSERT_VALUES));
 
           // Setting the solution b
           val = 0.5 * (x0(localmesh->xend, y, z) + x0(localmesh->xend + 1, y, z));
-          VecSetValues(bs, 1, &ind, &val, INSERT_VALUES);
+          BOUT_DO_PETSC(VecSetValues(bs, 1, &ind, &val, INSERT_VALUES));
 
           ind++;
         }
@@ -915,21 +917,21 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
     ASSERT1(ind == Iend); // Reached end of range
 
     // Assemble RHS Vector
-    VecAssemblyBegin(bs);
-    VecAssemblyEnd(bs);
+    BOUT_DO_PETSC(VecAssemblyBegin(bs));
+    BOUT_DO_PETSC(VecAssemblyEnd(bs));
 
     // Assemble Trial Solution Vector
-    VecAssemblyBegin(xs);
-    VecAssemblyEnd(xs);
+    BOUT_DO_PETSC(VecAssemblyBegin(xs));
+    BOUT_DO_PETSC(VecAssemblyEnd(xs));
 
     //////////////////////////
     // Solve the system
 
-    KSPSolve(it.ksp, bs, xs);
+    BOUT_DO_PETSC(KSPSolve(it.ksp, bs, xs));
 
     // Check if the solve converged
     KSPConvergedReason reason;
-    KSPGetConvergedReason(it.ksp, &reason);
+    BOUT_DO_PETSC(KSPGetConvergedReason(it.ksp, &reason));
 
     if (reason <= 0) {
       throw BoutException("LaplaceXZ failed to converge. Reason {} ({:d})",
@@ -944,7 +946,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
     if (localmesh->firstX()) {
       for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
         PetscScalar val;
-        VecGetValues(xs, 1, &ind, &val);
+        BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
         for (int x = localmesh->xstart - 1; x >= 0; --x) {
           result(x, y, z) = val;
         }
@@ -955,7 +957,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
     for (int x = localmesh->xstart; x <= localmesh->xend; x++) {
       for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
         PetscScalar val;
-        VecGetValues(xs, 1, &ind, &val);
+        BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
         result(x, y, z) = val;
         ind++;
       }
@@ -965,7 +967,7 @@ Field3D LaplaceXZpetsc::solve(const Field3D& bin, const Field3D& x0in) {
     if (localmesh->lastX()) {
       for (int z = localmesh->zstart; z <= localmesh->zend; z++) {
         PetscScalar val;
-        VecGetValues(xs, 1, &ind, &val);
+        BOUT_DO_PETSC(VecGetValues(xs, 1, &ind, &val));
         for (int x = localmesh->xend + 1; x < localmesh->LocalNx; ++x) {
           result(x, y, z) = val;
         }
diff --git a/src/solver/impls/imex-bdf2/imex-bdf2.cxx b/src/solver/impls/imex-bdf2/imex-bdf2.cxx
index 26721da95b..4bf235feea 100644
--- a/src/solver/impls/imex-bdf2/imex-bdf2.cxx
+++ b/src/solver/impls/imex-bdf2/imex-bdf2.cxx
@@ -12,6 +12,7 @@
 #include <bout/mesh.hxx>
 #include <bout/options.hxx>
 #include <bout/petsc_interface.hxx>
+#include <bout/petsclib.hxx>
 #include <bout/solver.hxx>
 #include <bout/utils.hxx>
 
@@ -492,15 +493,15 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
       }
 
       // Pre-allocate
-      MatMPIAIJSetPreallocation(Jmf, 0, d_nnz, 0, o_nnz);
-      MatSetUp(Jmf);
-      MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
-      PetscFree(d_nnz);
-      PetscFree(o_nnz);
+      BOUT_DO_PETSC(MatMPIAIJSetPreallocation(Jmf, 0, d_nnz, 0, o_nnz));
+      BOUT_DO_PETSC(MatSetUp(Jmf));
+      BOUT_DO_PETSC(MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
+      BOUT_DO_PETSC(PetscFree(d_nnz));
+      BOUT_DO_PETSC(PetscFree(o_nnz));
 
       // Determine which row/columns of the matrix are locally owned
       int Istart, Iend;
-      MatGetOwnershipRange(Jmf, &Istart, &Iend);
+      BOUT_DO_PETSC(MatGetOwnershipRange(Jmf, &Istart, &Iend));
 
       // Convert local into global indices
       index += Istart;
@@ -547,7 +548,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
                 PetscInt col = ind2 + j;
 
                 // output.write("SETTING 1: {:d}, {:d}\n", row, col);
-                MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES);
+                BOUT_DO_PETSC(MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES));
               }
             }
           }
@@ -567,7 +568,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
               for (int j = 0; j < n2d; j++) {
                 PetscInt col = ind0 + j;
                 // output.write("SETTING 2: {:d}, {:d}\n", row, col);
-                MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES);
+                BOUT_DO_PETSC(MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES));
               }
 
               // 5 point star pattern
@@ -593,7 +594,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
                 for (int j = 0; j < n3d; j++) {
                   PetscInt col = ind2 + j;
                   // output.write("SETTING 3: {:d}, {:d}\n", row, col);
-                  MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES);
+                  BOUT_DO_PETSC(MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES));
                 }
               }
 
@@ -610,7 +611,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
                 for (int j = 0; j < n3d; j++) {
                   PetscInt col = ind2 + j;
                   // output.write("SETTING 4: {:d}, {:d}\n", row, col);
-                  MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES);
+                  BOUT_DO_PETSC(MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES));
                 }
 
                 int zm = (z - 1 + nz) % nz;
@@ -621,7 +622,7 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
                 for (int j = 0; j < n3d; j++) {
                   PetscInt col = ind2 + j;
                   // output.write("SETTING 5: {:d}, {:d}\n", row, col);
-                  MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES);
+                  BOUT_DO_PETSC(MatSetValues(Jmf, 1, &row, 1, &col, &val, INSERT_VALUES));
                 }
               }
             }
@@ -631,31 +632,32 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
       // Finished marking non-zero entries
 
       // Assemble Matrix
-      MatAssemblyBegin(Jmf, MAT_FINAL_ASSEMBLY);
-      MatAssemblyEnd(Jmf, MAT_FINAL_ASSEMBLY);
+      BOUT_DO_PETSC(MatAssemblyBegin(Jmf, MAT_FINAL_ASSEMBLY));
+      BOUT_DO_PETSC(MatAssemblyEnd(Jmf, MAT_FINAL_ASSEMBLY));
 
       ISColoring iscoloring;
       MatColoring coloring; // This new in PETSc 3.5
-      MatColoringCreate(Jmf, &coloring);
-      MatColoringSetType(coloring, MATCOLORINGSL);
-      MatColoringSetFromOptions(coloring);
+      BOUT_DO_PETSC(MatColoringCreate(Jmf, &coloring));
+      BOUT_DO_PETSC(MatColoringSetType(coloring, MATCOLORINGSL));
+      BOUT_DO_PETSC(MatColoringSetFromOptions(coloring));
       // Calculate index sets
-      MatColoringApply(coloring, &iscoloring);
-      MatColoringDestroy(&coloring);
+      BOUT_DO_PETSC(MatColoringApply(coloring, &iscoloring));
+      BOUT_DO_PETSC(MatColoringDestroy(&coloring));
 
       // Create data structure for SNESComputeJacobianDefaultColor
-      MatFDColoringCreate(Jmf, iscoloring, &fdcoloring);
+      BOUT_DO_PETSC(MatFDColoringCreate(Jmf, iscoloring, &fdcoloring));
       // Set the function to difference
-      MatFDColoringSetFunction(fdcoloring,
-                               bout::cast_MatFDColoringFn(FormFunctionForColoring), this);
-      MatFDColoringSetFromOptions(fdcoloring);
-      MatFDColoringSetUp(Jmf, iscoloring, fdcoloring);
-      ISColoringDestroy(&iscoloring);
+      BOUT_DO_PETSC(MatFDColoringSetFunction(
+          fdcoloring, bout::cast_MatFDColoringFn(FormFunctionForColoring), this));
+      BOUT_DO_PETSC(MatFDColoringSetFromOptions(fdcoloring));
+      BOUT_DO_PETSC(MatFDColoringSetUp(Jmf, iscoloring, fdcoloring));
+      BOUT_DO_PETSC(ISColoringDestroy(&iscoloring));
 
-      SNESSetJacobian(*snesIn, Jmf, Jmf, SNESComputeJacobianDefaultColor, fdcoloring);
+      BOUT_DO_PETSC(SNESSetJacobian(*snesIn, Jmf, Jmf, SNESComputeJacobianDefaultColor,
+                                    fdcoloring));
 
       // Re-use Jacobian
-      SNESSetLagJacobian(*snesIn, lag_jacobian);
+      BOUT_DO_PETSC(SNESSetLagJacobian(*snesIn, lag_jacobian));
 
       // MatView(Jmf, PETSC_VIEWER_DRAW_WORLD);
       // MatView(Jmf, PETSC_VIEWER_STDOUT_WORLD);
@@ -663,40 +665,41 @@ void IMEXBDF2::constructSNES(SNES* snesIn) {
       // Brute force calculation
       // NOTE: Slow!
 
-      MatCreateAIJ(
+      BOUT_DO_PETSC(MatCreateAIJ(
           BoutComm::get(), nlocal, nlocal,  // Local sizes
           PETSC_DETERMINE, PETSC_DETERMINE, // Global sizes
           3, // Number of nonzero entries in diagonal portion of local submatrix
           nullptr,
           0, // Number of nonzeros per row in off-diagonal portion of local submatrix
-          nullptr, &Jmf);
+          nullptr, &Jmf));
 
-      SNESSetJacobian(*snesIn, Jmf, Jmf, SNESComputeJacobianDefault, this);
+      BOUT_DO_PETSC(SNESSetJacobian(*snesIn, Jmf, Jmf, SNESComputeJacobianDefault, this));
 
-      MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+      BOUT_DO_PETSC(MatSetOption(Jmf, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
     }
   }
 
   /////////////////////////////////////////////////////
   // Set tolerances
-  SNESSetTolerances(*snesIn, atol, rtol, PETSC_DEFAULT, max_nonlinear_it, PETSC_DEFAULT);
+  BOUT_DO_PETSC(SNESSetTolerances(*snesIn, atol, rtol, PETSC_DEFAULT, max_nonlinear_it,
+                                  PETSC_DEFAULT));
 
   /////////////////////////////////////////////////////
   // Preconditioner
 
   // Get KSP context from SNES
   KSP ksp;
-  SNESGetKSP(*snesIn, &ksp);
+  BOUT_DO_PETSC(SNESGetKSP(*snesIn, &ksp));
 
   if (kspsetinitialguessnonzero) {
-    KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
+    BOUT_DO_PETSC(KSPSetInitialGuessNonzero(ksp, PETSC_TRUE));
   }
 
-  KSPSetTolerances(ksp,
-                   PETSC_DEFAULT, // rtol
-                   PETSC_DEFAULT, // abstol
-                   PETSC_DEFAULT, // dtol (divergence tolerance)
-                   maxl);         // Maximum number of iterations
+  BOUT_DO_PETSC(KSPSetTolerances(ksp,
+                                 PETSC_DEFAULT, // rtol
+                                 PETSC_DEFAULT, // abstol
+                                 PETSC_DEFAULT, // dtol (divergence tolerance)
+                                 maxl));        // Maximum number of iterations
 
   // Get PC context from KSP
   PC pc;