diff --git a/framework/src/utils/PetscDMMoose.C b/framework/src/utils/PetscDMMoose.C
index 3f71cf0a3115..0a835bc445ae 100644
--- a/framework/src/utils/PetscDMMoose.C
+++ b/framework/src/utils/PetscDMMoose.C
@@ -140,7 +140,7 @@ DMMooseValidityCheck(DM dm)
                      "Got DM of type %s, not of type %s",
                      ((PetscObject)dm)->type_name,
                      DMMOOSE);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -159,7 +159,7 @@ DMMooseGetContacts(DM dm,
     contact_names.push_back(it.second);
     displaced.push_back((*dmm->_contact_displaced)[it.second]);
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -178,7 +178,7 @@ DMMooseGetUnContacts(DM dm,
     uncontact_names.push_back(it.second);
     displaced.push_back((*dmm->_uncontact_displaced)[it.second]);
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -192,7 +192,7 @@ DMMooseGetSides(DM dm, std::vector<std::string> & side_names)
   DM_Moose * dmm = (DM_Moose *)dm->data;
   for (const auto & it : *(dmm->_side_ids))
     side_names.push_back(it.first);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -206,7 +206,7 @@ DMMooseGetUnSides(DM dm, std::vector<std::string> & side_names)
   DM_Moose * dmm = (DM_Moose *)dm->data;
   for (const auto & it : *(dmm->_unside_ids))
     side_names.push_back(it.first);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -220,7 +220,7 @@ DMMooseGetBlocks(DM dm, std::vector<std::string> & block_names)
   DM_Moose * dmm = (DM_Moose *)dm->data;
   for (const auto & it : *(dmm->_block_ids))
     block_names.push_back(it.first);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -234,7 +234,7 @@ DMMooseGetVariables(DM dm, std::vector<std::string> & var_names)
   DM_Moose * dmm = (DM_Moose *)(dm->data);
   for (const auto & it : *(dmm->_var_ids))
     var_names.push_back(it.first);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -251,7 +251,7 @@ DMMooseSetNonlinearSystem(DM dm, NonlinearSystemBase & nl)
             "Cannot reset the NonlinearSystem after DM has been set up.");
   DM_Moose * dmm = (DM_Moose *)(dm->data);
   dmm->_nl = &nl;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -268,7 +268,7 @@ DMMooseSetName(DM dm, const std::string & dm_name)
             "Cannot reset the MOOSE DM name after DM has been set up.");
   DM_Moose * dmm = (DM_Moose *)(dm->data);
   *dmm->_name = dm_name;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -286,7 +286,7 @@ DMMooseSetParentDM(DM dm, DM_Moose * parent)
 
   DM_Moose * dmm = (DM_Moose *)(dm->data);
   dmm->_parent = parent;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -317,7 +317,7 @@ DMMooseSetVariables(DM dm, const std::set<std::string> & vars)
   }
 
   dmm->_vars = new std::set<std::string>(std::move(processed_vars));
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -334,7 +334,7 @@ DMMooseSetBlocks(DM dm, const std::set<std::string> & blocks)
   if (dmm->_blocks)
     delete dmm->_blocks;
   dmm->_blocks = new std::set<std::string>(blocks);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -351,7 +351,7 @@ DMMooseSetSides(DM dm, const std::set<std::string> & sides)
   if (dmm->_sides)
     delete dmm->_sides;
   dmm->_sides = new std::set<std::string>(sides);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -368,7 +368,7 @@ DMMooseSetUnSides(DM dm, const std::set<std::string> & unsides)
   if (dmm->_unsides)
     delete dmm->_unsides;
   dmm->_unsides = new std::set<std::string>(unsides);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -385,7 +385,7 @@ DMMooseSetUnSideByVar(DM dm, const std::set<std::string> & unside_by_var)
   if (dmm->_unside_by_var)
     delete dmm->_unside_by_var;
   dmm->_unside_by_var = new std::set<std::string>(unside_by_var);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -417,7 +417,7 @@ DMMooseSetContacts(DM dm,
     dmm->_contacts->insert(contacts[i]);
     dmm->_contact_displaced->insert(std::make_pair(contacts[i], displaced[i]));
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -450,7 +450,7 @@ DMMooseSetUnContacts(DM dm,
     dmm->_uncontacts->insert(uncontacts[i]);
     dmm->_uncontact_displaced->insert(std::make_pair(uncontacts[i], displaced[i]));
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -463,7 +463,7 @@ DMMooseGetNonlinearSystem(DM dm, NonlinearSystemBase *& nl)
   CHKERRQ(ierr);
   DM_Moose * dmm = (DM_Moose *)(dm->data);
   nl = dmm->_nl;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -504,7 +504,7 @@ DMMooseSetSplitNames(DM dm, const std::vector<std::string> & split_names)
     (*dmm->_splits)[name] = info;
     dmm->_splitlocs->insert(std::make_pair(name, i));
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -527,7 +527,7 @@ DMMooseGetSplitNames(DM dm, std::vector<std::string> & split_names)
       unsigned int sloc = lit.second;
       split_names[sloc] = sname;
     }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -538,7 +538,7 @@ DMMooseGetEmbedding_Private(DM dm, IS * embedding)
 
   PetscFunctionBegin;
   if (!embedding)
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   if (!dmm->_embedding)
   {
     // The rules interpreting the coexistence of blocks (un)sides/(un)contacts are these
@@ -870,7 +870,7 @@ DMMooseGetEmbedding_Private(DM dm, IS * embedding)
   CHKERRQ(ierr);
   *embedding = dmm->_embedding;
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -886,7 +886,7 @@ DMCreateFieldDecomposition_Moose(
 
   /* Only called after DMSetUp(). */
   if (!dmm->_splitlocs)
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   *len = dmm->_splitlocs->size();
   if (namelist)
   {
@@ -1014,7 +1014,7 @@ DMCreateFieldDecomposition_Moose(
               split_size_sum,
               dmm->_nl->nonlinearSolver()->system().get_system_matrix().local_m());
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1030,7 +1030,7 @@ DMCreateDomainDecomposition_Moose(
     *outerislist = LIBMESH_PETSC_NULLPTR; /* FIX: allow mesh-based overlap. */
   ierr = DMCreateFieldDecomposition_Moose(dm, len, namelist, innerislist, dmlist);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1102,7 +1102,7 @@ DMMooseFunction(DM dm, Vec x, Vec r)
     mooseError(err.str());
   }
   R.close();
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1114,7 +1114,7 @@ SNESFunction_DMMoose(SNES, Vec x, Vec r, void * ctx)
   PetscFunctionBegin;
   ierr = DMMooseFunction(dm, x, r);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1196,7 +1196,7 @@ DMMooseJacobian(DM dm, Vec x, Mat jac, Mat pc)
   }
   the_pc.close();
   Jac.close();
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1208,7 +1208,7 @@ SNESJacobian_DMMoose(SNES, Vec x, Mat jac, Mat pc, void * ctx)
   PetscFunctionBegin;
   ierr = DMMooseJacobian(dm, x, jac, pc);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1235,7 +1235,7 @@ DMVariableBounds_Moose(DM dm, Vec xl, Vec xu)
   else
     SETERRQ(
         ((PetscObject)dm)->comm, PETSC_ERR_ARG_WRONG, "No bounds calculation in this Moose object");
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1285,7 +1285,7 @@ DMCreateGlobalVector_Moose(DM dm, Vec * x)
   ierr = VecSetDM(*x, dm);
   CHKERRQ(ierr);
 #endif
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1342,7 +1342,7 @@ DMCreateMatrix_Moose(DM dm, Mat * A)
    * settings made here. */
   ierr = MatSetUp(*A);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1475,7 +1475,7 @@ DMView_Moose(DM dm, PetscViewer viewer)
   else
     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Non-ASCII viewers are not supported");
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1498,7 +1498,7 @@ DMMooseGetMeshBlocks_Private(DM dm, std::set<subdomain_id_type> & blocks)
     blocks.insert(elem->subdomain_id());
   // Some subdomains may only live on other processors
   mesh.comm().set_union(blocks);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1757,7 +1757,7 @@ DMSetUp_Moose_Pre(DM dm)
   }
   ierr = PetscObjectSetName((PetscObject)dm, name.c_str());
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1771,7 +1771,7 @@ DMMooseReset(DM dm)
   ierr = PetscObjectTypeCompare((PetscObject)dm, DMMOOSE, &ismoose);
   CHKERRQ(ierr);
   if (!ismoose)
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   if (!dmm->_nl)
     SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONGSTATE, "No Moose system set for DM_Moose");
   ierr = ISDestroy(&dmm->_embedding);
@@ -1788,7 +1788,7 @@ DMMooseReset(DM dm)
     }
   }
   dm->setupcalled = PETSC_FALSE;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -1851,7 +1851,7 @@ DMSetUp_Moose(DM dm)
       ierr = DMSetVariableBounds(dm, DMVariableBounds_Moose);
     CHKERRQ(ierr);
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 #if !PETSC_VERSION_LESS_THAN(3, 18, 0)
@@ -2236,7 +2236,7 @@ DMSetFromOptions_Moose(PetscOptions * /*options*/, DM dm) // >= 3.6.0
   ierr = DMSetUp_Moose_Pre(dm);
   CHKERRQ(ierr); /* Need some preliminary set up because, strangely enough, DMView() is called in
                     DMSetFromOptions(). */
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 static PetscErrorCode
@@ -2291,7 +2291,7 @@ DMDestroy_Moose(DM dm)
   CHKERRQ(ierr);
   ierr = PetscFree(dm->data);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -2308,7 +2308,7 @@ DMCreateMoose(MPI_Comm comm, NonlinearSystemBase & nl, const std::string & dm_na
   CHKERRQ(ierr);
   ierr = DMMooseSetName(*dm, dm_name);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 EXTERN_C_BEGIN
@@ -2370,7 +2370,7 @@ DMCreate_Moose(DM dm)
   dm->ops->view = DMView_Moose;
   dm->ops->setfromoptions = DMSetFromOptions_Moose;
   dm->ops->setup = DMSetUp_Moose;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 EXTERN_C_END
 
@@ -2425,7 +2425,7 @@ SNESUpdateDMMoose(SNES snes, PetscInt iteration)
     ierr = PCDestroy(&pc);
     CHKERRQ(ierr);
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -2441,5 +2441,5 @@ DMMooseRegisterAll()
     CHKERRQ(ierr);
     DMMooseRegisterAllCalled = PETSC_TRUE;
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
diff --git a/framework/src/utils/PetscSupport.C b/framework/src/utils/PetscSupport.C
index 4bbbc632787d..3d4a315dc1d9 100644
--- a/framework/src/utils/PetscSupport.C
+++ b/framework/src/utils/PetscSupport.C
@@ -268,6 +268,7 @@ petscSetOptions(const PetscOptions & po, const SolverParams & solver_params)
 PetscErrorCode
 petscSetupOutput(CommandLine * cmd_line)
 {
+  PetscFunctionBegin;
   char code[10] = {45, 45, 109, 111, 111, 115, 101};
   const std::vector<std::string> argv = cmd_line->getArguments();
   for (const auto & arg : argv)
@@ -278,7 +279,7 @@ petscSetupOutput(CommandLine * cmd_line)
       break;
     }
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -290,6 +291,7 @@ petscNonlinearConverged(SNES snes,
                         SNESConvergedReason * reason,
                         void * ctx)
 {
+  PetscFunctionBegin;
   FEProblemBase & problem = *static_cast<FEProblemBase *>(ctx);
 
   // Let's be nice and always check PETSc error codes.
@@ -346,7 +348,7 @@ petscNonlinearConverged(SNES snes,
   if (domainerror)
   {
     *reason = SNES_DIVERGED_FUNCTION_DOMAIN;
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
 
   // Error message that will be set by the FEProblemBase.
@@ -418,7 +420,7 @@ petscNonlinearConverged(SNES snes,
       break;
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PCSide
diff --git a/framework/src/utils/SlepcSupport.C b/framework/src/utils/SlepcSupport.C
index 741c13afb375..9b198ca61366 100644
--- a/framework/src/utils/SlepcSupport.C
+++ b/framework/src/utils/SlepcSupport.C
@@ -520,13 +520,13 @@ mooseEPSFormMatrices(EigenProblem & eigen_problem, EPS eps, Vec x, void * ctx)
   PetscFunctionBegin;
 
   if (eigen_problem.constantMatrices() && eigen_problem.wereMatricesFormed())
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
 
   if (eigen_problem.onLinearSolver())
     // We reach here during linear iteration when solve type is PJFNKMO.
     // We will use the matrices assembled at the beginning of this Newton
     // iteration for the following residual evaluation.
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
 
   NonlinearEigenSystem & eigen_nl = eigen_problem.getCurrentNonlinearEigenSystem();
   SNES snes = eigen_nl.getSNES();
@@ -552,7 +552,7 @@ mooseEPSFormMatrices(EigenProblem & eigen_problem, EPS eps, Vec x, void * ctx)
   moosePetscSNESFormMatricesTags(
       snes, x, mats, ctx, {eigen_nl.nonEigenMatrixTag(), eigen_nl.eigenMatrixTag()});
   eigen_problem.wereMatricesFormed(true);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 void
@@ -644,7 +644,7 @@ mooseSlepcEigenFormFunctionMFFD(void * ctx, Vec x, Vec r)
 
   eigen_problem->onLinearSolver(false);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -683,7 +683,7 @@ mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
       ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
       CHKERRQ(ierr);
     }
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
 
   ierr = PetscObjectTypeCompare((PetscObject)pc, MATSHELL, &pisshell);
@@ -700,7 +700,7 @@ mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
     ierr = MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY);
     CHKERRQ(ierr);
 
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
 
   // Jacobian and precond matrix are the same
@@ -709,7 +709,7 @@ mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
     if (!pisshell)
       moosePetscSNESFormMatrixTag(snes, x, pc, ctx, eigen_nl.precondMatrixTag());
 
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
   else
   {
@@ -718,7 +718,7 @@ mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
       std::vector<Mat> mats = {jac, pc};
       moosePetscSNESFormMatricesTags(
           snes, x, mats, ctx, {eigen_nl.nonEigenMatrixTag(), eigen_nl.precondMatrixTag()});
-      return PETSC_SUCCESS;
+      PetscFunctionReturn(PETSC_SUCCESS);
     }
     if (!pisshell) // We need to form only precond matrix
     {
@@ -727,7 +727,7 @@ mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
       CHKERRQ(ierr);
       ierr = MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
       CHKERRQ(ierr);
-      return PETSC_SUCCESS;
+      PetscFunctionReturn(PETSC_SUCCESS);
     }
     if (!jisshell && !jismffd) // We need to form only Jacobian matrix
     {
@@ -736,10 +736,10 @@ mooseSlepcEigenFormJacobianA(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
       CHKERRQ(ierr);
       ierr = MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY);
       CHKERRQ(ierr);
-      return PETSC_SUCCESS;
+      PetscFunctionReturn(PETSC_SUCCESS);
     }
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -774,7 +774,7 @@ mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
     ierr = MatAssemblyEnd(pc, MAT_FINAL_ASSEMBLY);
     CHKERRQ(ierr);
 
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
 
   if (jac != pc && (!jshell && !jshell))
@@ -790,7 +790,7 @@ mooseSlepcEigenFormJacobianB(SNES snes, Vec x, Mat jac, Mat pc, void * ctx)
     CHKERRQ(ierr);
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 void
@@ -850,12 +850,12 @@ mooseSlepcEigenFormFunctionA(SNES snes, Vec x, Vec r, void * ctx)
     ierr = MatMult(A, x, r);
     CHKERRQ(ierr);
 
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
 
   moosePetscSNESFormFunction(snes, x, r, ctx, eigen_nl.nonEigenVectorTag());
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -898,7 +898,7 @@ mooseSlepcEigenFormFunctionB(SNES snes, Vec x, Vec r, void * ctx)
     CHKERRQ(ierr);
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -942,7 +942,7 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
       CHKERRQ(ierr);
     }
 
-    return PETSC_SUCCESS;
+    PetscFunctionReturn(PETSC_SUCCESS);
   }
 
   PetscVector<Number> X_global(x, sys.comm()), AX(Ax, sys.comm()), BX(Bx, sys.comm());
@@ -978,7 +978,7 @@ mooseSlepcEigenFormFunctionAB(SNES /*snes*/, Vec x, Vec Ax, Vec Bx, void * ctx)
     CHKERRQ(ierr);
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -987,7 +987,7 @@ mooseSlepcEigenFormNorm(SNES /*snes*/, Vec /*Bx*/, PetscReal * norm, void * ctx)
   PetscFunctionBegin;
   auto * const eigen_problem = static_cast<EigenProblem *>(ctx);
   *norm = eigen_problem->formNorm();
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 void
@@ -1098,7 +1098,7 @@ mooseMatMult_Eigen(Mat mat, Vec x, Vec r)
     CHKERRQ(ierr);
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1117,7 +1117,7 @@ mooseMatMult_NonEigen(Mat mat, Vec x, Vec r)
 
   mooseMatMult(*eigen_problem, x, r, eigen_nl.nonEigenVectorTag());
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 void
@@ -1141,7 +1141,7 @@ registerPCToPETSc()
   ierr = PCRegister("moosepc", PCCreate_MoosePC);
   CHKERRQ(ierr);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PETSC_EXTERN PetscErrorCode
@@ -1154,7 +1154,7 @@ PCCreate_MoosePC(PC pc)
   pc->ops->setup = PCSetUp_MoosePC;
   pc->ops->apply = PCApply_MoosePC;
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1163,7 +1163,7 @@ PCDestroy_MoosePC(PC /*pc*/)
   PetscFunctionBegin;
   /* We do not need to do anything right now, but later we may have some data we need to free here
    */
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1180,7 +1180,7 @@ PCView_MoosePC(PC /*pc*/, PetscViewer viewer)
     ierr = PetscViewerASCIIPrintf(viewer, "  %s\n", "moosepc");
     CHKERRQ(ierr);
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1191,6 +1191,7 @@ PCApply_MoosePC(PC pc, Vec x, Vec y)
   PetscContainer container;
   PetscErrorCode ierr;
 
+  PetscFunctionBegin;
   ierr = PCGetOperators(pc, &Amat, &Pmat);
   CHKERRQ(ierr);
   ierr = PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container);
@@ -1217,7 +1218,7 @@ PCApply_MoosePC(PC pc, Vec x, Vec y)
 
   preconditioner->apply(x_vec, y_vec);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1228,6 +1229,7 @@ PCSetUp_MoosePC(PC pc)
   Mat Amat, Pmat;
   PetscContainer container;
 
+  PetscFunctionBegin;
   ierr = PCGetOperators(pc, &Amat, &Pmat);
   CHKERRQ(ierr);
   ierr = PetscObjectQuery((PetscObject)Pmat, "formFunctionCtx", (PetscObject *)&container);
@@ -1253,7 +1255,7 @@ PCSetUp_MoosePC(PC pc)
 
   preconditioner->setup();
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1268,6 +1270,7 @@ mooseSlepcStoppingTest(EPS eps,
   PetscErrorCode ierr;
   EigenProblem * eigen_problem = static_cast<EigenProblem *>(ctx);
 
+  PetscFunctionBegin;
   ierr = EPSStoppingBasic(eps, its, max_it, nconv, nev, reason, NULL);
   LIBMESH_CHKERR(ierr);
 
@@ -1282,7 +1285,7 @@ mooseSlepcStoppingTest(EPS eps,
     *reason = EPS_CONVERGED_USER;
     eps->nconv = 1;
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1291,6 +1294,7 @@ mooseSlepcEPSGetSNES(EPS eps, SNES * snes)
   PetscErrorCode ierr;
   PetscBool same, nonlinear;
 
+  PetscFunctionBegin;
   ierr = PetscObjectTypeCompare((PetscObject)eps, EPSPOWER, &same);
   LIBMESH_CHKERR(ierr);
 
@@ -1306,7 +1310,7 @@ mooseSlepcEPSGetSNES(EPS eps, SNES * snes)
   ierr = EPSPowerGetSNES(eps, snes);
   LIBMESH_CHKERR(ierr);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1316,6 +1320,7 @@ mooseSlepcEPSSNESSetUpOptionPrefix(EPS eps)
   SNES snes;
   const char * prefix = nullptr;
 
+  PetscFunctionBegin;
   ierr = mooseSlepcEPSGetSNES(eps, &snes);
   LIBMESH_CHKERR(ierr);
   // There is an extra "eps_power" in snes that users do not like it.
@@ -1327,7 +1332,7 @@ mooseSlepcEPSSNESSetUpOptionPrefix(EPS eps)
   ierr = SNESSetOptionsPrefix(snes, prefix);
   LIBMESH_CHKERR(ierr);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1338,6 +1343,7 @@ mooseSlepcEPSSNESSetCustomizePC(EPS eps)
   KSP ksp;
   PC pc;
 
+  PetscFunctionBegin;
   // Get SNES from EPS
   ierr = mooseSlepcEPSGetSNES(eps, &snes);
   LIBMESH_CHKERR(ierr);
@@ -1350,7 +1356,7 @@ mooseSlepcEPSSNESSetCustomizePC(EPS eps)
   // Set PC type
   ierr = PCSetType(pc, "moosepc");
   LIBMESH_CHKERR(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1360,6 +1366,7 @@ mooseSlepcEPSSNESKSPSetPCSide(FEProblemBase & problem, EPS eps)
   SNES snes;
   KSP ksp;
 
+  PetscFunctionBegin;
   // Get SNES from EPS
   ierr = mooseSlepcEPSGetSNES(eps, &snes);
   LIBMESH_CHKERR(ierr);
@@ -1370,7 +1377,7 @@ mooseSlepcEPSSNESKSPSetPCSide(FEProblemBase & problem, EPS eps)
   Moose::PetscSupport::petscSetDefaultPCSide(problem, ksp);
 
   Moose::PetscSupport::petscSetDefaultKSPNormType(problem, ksp);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -1387,6 +1394,7 @@ mooseSlepcEPSMonitor(EPS eps,
   auto ierr = (PetscErrorCode)0;
   PetscScalar eigenr, eigeni;
 
+  PetscFunctionBegin;
   EigenProblem * eigen_problem = static_cast<EigenProblem *>(mctx);
   auto & console = eigen_problem->console();
 
@@ -1405,7 +1413,7 @@ mooseSlepcEPSMonitor(EPS eps,
   console << " Iteration " << its << std::setprecision(10) << std::fixed
           << (inverse ? " k-eigenvalue = " : " eigenvalue = ") << eigenvalue << std::endl;
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 } // namespace SlepcSupport
diff --git a/modules/external_petsc_solver/src/petscsolver/PETScDiffusionFDM.C b/modules/external_petsc_solver/src/petscsolver/PETScDiffusionFDM.C
index 00b7b67c9501..d843033e2d72 100644
--- a/modules/external_petsc_solver/src/petscsolver/PETScDiffusionFDM.C
+++ b/modules/external_petsc_solver/src/petscsolver/PETScDiffusionFDM.C
@@ -83,7 +83,7 @@ PETScExternalSolverCreate(MPI_Comm comm, TS * ts)
   CHKERRQ(ierr);
   ierr = TSSetFromOptions(*ts);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -94,7 +94,7 @@ PETScExternalSolverDestroy(TS ts)
   PetscFunctionBeginUser;
   ierr = TSDestroy(&ts);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 /*
@@ -185,7 +185,7 @@ externalPETScDiffusionFDMSolve(
   CHKERRQ(ierr);
   *converged = reason > 0 ? PETSC_TRUE : PETSC_FALSE;
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 /* --------------------------------------------------------------------- */
@@ -335,7 +335,7 @@ FormIFunction(TS ts, PetscReal /*t*/, Vec U, Vec Udot, Vec F, void * /*ctx*/)
   CHKERRQ(ierr);
   ierr = PetscLogFlops(11.0 * ym * xm);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 /* --------------------------------------------------------------------- */
@@ -462,7 +462,7 @@ FormIJacobian(
     CHKERRQ(ierr);
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 /* ------------------------------------------------------------------- */
@@ -524,5 +524,5 @@ FormInitialSolution(TS ts, Vec U, void * /*ptr*/)
   /* Restore vectors */
   ierr = DMDAVecRestoreArray(da, U, &u);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
diff --git a/modules/fsi/test/tests/2d-finite-strain-steady/tests b/modules/fsi/test/tests/2d-finite-strain-steady/tests
index 09fc023fc123..273ccdec19a5 100644
--- a/modules/fsi/test/tests/2d-finite-strain-steady/tests
+++ b/modules/fsi/test/tests/2d-finite-strain-steady/tests
@@ -8,7 +8,7 @@
     csvdiff = thermal-me.csv
     method = '!dbg'
     valgrind = 'none'
-    rel_err = 1e-4
+    rel_err = 2e-4
     # skip test if test is being run out-of-tree. Issue Ref: #25279
     installation_type = in_tree
   []
diff --git a/modules/optimization/src/executioners/OptimizeSolve.C b/modules/optimization/src/executioners/OptimizeSolve.C
index 0aa8d0602f40..d6e7068acef8 100644
--- a/modules/optimization/src/executioners/OptimizeSolve.C
+++ b/modules/optimization/src/executioners/OptimizeSolve.C
@@ -87,6 +87,7 @@ OptimizeSolve::taoSolve()
   // Petsc error code to be checked after each petsc call
   auto ierr = (PetscErrorCode)0;
 
+  PetscFunctionBegin;
   // Initialize tao object
   ierr = TaoCreate(_my_comm.get(), &_tao);
   CHKERRQ(ierr);
@@ -339,31 +340,34 @@ OptimizeSolve::monitor(Tao tao, void * ctx)
   PetscInt its;
   PetscReal f, gnorm, cnorm, xdiff;
 
+  PetscFunctionBegin;
   auto ierr = TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, &reason);
   CHKERRQ(ierr);
 
   auto * solver = static_cast<OptimizeSolve *>(ctx);
   solver->setTaoSolutionStatus((double)f, (int)its, (double)gnorm, (double)cnorm, (double)xdiff);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::objectiveFunctionWrapper(Tao /*tao*/, Vec x, Real * objective, void * ctx)
 {
+  PetscFunctionBegin;
   auto * solver = static_cast<OptimizeSolve *>(ctx);
 
   libMesh::PetscVector<Number> param(x, solver->_my_comm);
   *solver->_parameters = param;
 
   (*objective) = solver->objectiveFunction();
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::objectiveAndGradientFunctionWrapper(
     Tao /*tao*/, Vec x, Real * objective, Vec gradient, void * ctx)
 {
+  PetscFunctionBegin;
   auto * solver = static_cast<OptimizeSolve *>(ctx);
 
   libMesh::PetscVector<Number> param(x, solver->_my_comm);
@@ -374,12 +378,13 @@ OptimizeSolve::objectiveAndGradientFunctionWrapper(
   libMesh::PetscVector<Number> grad(gradient, solver->_my_comm);
 
   solver->gradientFunction(grad);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::hessianFunctionWrapper(Tao /*tao*/, Vec x, Mat /*hessian*/, Mat /*pc*/, void * ctx)
 {
+  PetscFunctionBegin;
   // Define Hessian-vector multiplication routine
   auto * solver = static_cast<OptimizeSolve *>(ctx);
   libMesh::PetscVector<Number> param(x, solver->_my_comm);
@@ -387,13 +392,15 @@ OptimizeSolve::hessianFunctionWrapper(Tao /*tao*/, Vec x, Mat /*hessian*/, Mat /
   PetscErrorCode ierr = MatShellSetOperation(
       solver->_hessian, MATOP_MULT, (void (*)(void))OptimizeSolve::applyHessianWrapper);
   CHKERRQ(ierr);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::applyHessianWrapper(Mat H, Vec s, Vec Hs)
 {
   void * ctx;
+
+  PetscFunctionBegin;
   auto ierr = MatShellGetContext(H, &ctx);
   CHKERRQ(ierr);
 
@@ -406,6 +413,7 @@ OptimizeSolve::applyHessianWrapper(Mat H, Vec s, Vec Hs)
 PetscErrorCode
 OptimizeSolve::variableBoundsWrapper(Tao tao, Vec /*xl*/, Vec /*xu*/, void * ctx)
 {
+  PetscFunctionBegin;
   auto * solver = static_cast<OptimizeSolve *>(ctx);
 
   PetscErrorCode ierr = solver->variableBounds(tao);
@@ -456,6 +464,7 @@ OptimizeSolve::gradientFunction(libMesh::PetscVector<Number> & gradient)
 PetscErrorCode
 OptimizeSolve::applyHessian(libMesh::PetscVector<Number> & s, libMesh::PetscVector<Number> & Hs)
 {
+  PetscFunctionBegin;
   TIME_SECTION("applyHessian", 2, "Hessian forward/adjoint solve");
   // What happens for material inversion when the Hessian
   // is dependent on the parameters? Deal with it later???
@@ -486,12 +495,13 @@ OptimizeSolve::applyHessian(libMesh::PetscVector<Number> & s, libMesh::PetscVect
 
   _obj_function->computeGradient(Hs);
   _hess_iterate++;
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::variableBounds(Tao tao)
 {
+  PetscFunctionBegin;
   unsigned int sz = _obj_function->getNumParams();
 
   libMesh::PetscVector<Number> xl(_my_comm, sz);
@@ -512,19 +522,21 @@ OptimizeSolve::variableBounds(Tao tao)
 PetscErrorCode
 OptimizeSolve::equalityFunctionWrapper(Tao /*tao*/, Vec /*x*/, Vec ce, void * ctx)
 {
+  PetscFunctionBegin;
   // grab the solver
   auto * solver = static_cast<OptimizeSolve *>(ctx);
   libMesh::PetscVector<Number> eq_con(ce, solver->_my_comm);
   // use the OptimizationReporterBase class to actually compute equality constraints
   OptimizationReporterBase * obj_func = solver->getObjFunction();
   obj_func->computeEqualityConstraints(eq_con);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::equalityGradientFunctionWrapper(
     Tao /*tao*/, Vec /*x*/, Mat gradient_e, Mat /*gradient_epre*/, void * ctx)
 {
+  PetscFunctionBegin;
   // grab the solver
   auto * solver = static_cast<OptimizeSolve *>(ctx);
   libMesh::PetscMatrix<Number> grad_eq(gradient_e, solver->_my_comm);
@@ -532,25 +544,27 @@ OptimizeSolve::equalityGradientFunctionWrapper(
   // constraints gradient
   OptimizationReporterBase * obj_func = solver->getObjFunction();
   obj_func->computeEqualityGradient(grad_eq);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::inequalityFunctionWrapper(Tao /*tao*/, Vec /*x*/, Vec ci, void * ctx)
 {
+  PetscFunctionBegin;
   // grab the solver
   auto * solver = static_cast<OptimizeSolve *>(ctx);
   libMesh::PetscVector<Number> ineq_con(ci, solver->_my_comm);
   // use the OptimizationReporterBase class to actually compute equality constraints
   OptimizationReporterBase * obj_func = solver->getObjFunction();
   obj_func->computeInequalityConstraints(ineq_con);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::inequalityGradientFunctionWrapper(
     Tao /*tao*/, Vec /*x*/, Mat gradient_i, Mat /*gradient_ipre*/, void * ctx)
 {
+  PetscFunctionBegin;
   // grab the solver
   auto * solver = static_cast<OptimizeSolve *>(ctx);
   libMesh::PetscMatrix<Number> grad_ineq(gradient_i, solver->_my_comm);
@@ -558,7 +572,7 @@ OptimizeSolve::inequalityGradientFunctionWrapper(
   // constraints gradient
   OptimizationReporterBase * obj_func = solver->getObjFunction();
   obj_func->computeInequalityGradient(grad_ineq);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
@@ -566,6 +580,7 @@ OptimizeSolve::taoALCreate()
 {
   auto ierr = (PetscErrorCode)0;
 
+  PetscFunctionBegin;
   if (_obj_function->getNumEqCons())
   {
     // Create equality vector
@@ -633,13 +648,15 @@ OptimizeSolve::taoALCreate()
         _tao, _gradient_i, _gradient_i, inequalityGradientFunctionWrapper, this);
     CHKERRQ(ierr);
   }
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 OptimizeSolve::taoALDestroy()
 {
   auto ierr = (PetscErrorCode)0;
+
+  PetscFunctionBegin;
   if (_obj_function->getNumEqCons())
   {
     ierr = VecDestroy(&_ce);
@@ -656,5 +673,5 @@ OptimizeSolve::taoALDestroy()
     CHKERRQ(ierr);
   }
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
diff --git a/modules/stochastic_tools/src/utils/GaussianProcessHandler.C b/modules/stochastic_tools/src/utils/GaussianProcessHandler.C
index 52d5c89eb2db..fdd40a50e298 100644
--- a/modules/stochastic_tools/src/utils/GaussianProcessHandler.C
+++ b/modules/stochastic_tools/src/utils/GaussianProcessHandler.C
@@ -158,6 +158,7 @@ GaussianProcessHandler::tuneHyperParamsTAO(const RealEigenMatrix & training_para
   PetscErrorCode ierr;
   Tao tao;
 
+  PetscFunctionBegin;
   _training_params = &training_params;
   _training_data = &training_data;
 
@@ -213,16 +214,17 @@ GaussianProcessHandler::tuneHyperParamsTAO(const RealEigenMatrix & training_para
   ierr = TaoDestroy(&tao);
   CHKERRQ(ierr);
 
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 PetscErrorCode
 GaussianProcessHandler::formInitialGuessTAO(Vec theta_vec)
 {
+  PetscFunctionBegin;
   libMesh::PetscVector<Number> theta(theta_vec, _tao_comm);
   _covariance_function->buildHyperParamMap(_hyperparam_map, _hyperparam_vec_map);
   mapToPetscVec(_tuning_data, _hyperparam_map, _hyperparam_vec_map, theta);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 void
@@ -243,9 +245,10 @@ PetscErrorCode
 GaussianProcessHandler::formFunctionGradientWrapper(
     Tao tao, Vec theta_vec, PetscReal * f, Vec grad_vec, void * ptr)
 {
+  PetscFunctionBegin;
   GaussianProcessHandler * GP_ptr = (GaussianProcessHandler *)ptr;
   GP_ptr->formFunctionGradient(tao, theta_vec, f, grad_vec);
-  return PETSC_SUCCESS;
+  PetscFunctionReturn(PETSC_SUCCESS);
 }
 
 void