diff --git a/kharma/emhd/emhd.cpp b/kharma/emhd/emhd.cpp
index 34153ce0..3150a198 100644
--- a/kharma/emhd/emhd.cpp
+++ b/kharma/emhd/emhd.cpp
@@ -122,15 +122,16 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
     // Only enable limits internally if we're actually doing EMHD
     params.Add("enable_emhd_limits", enable_emhd_limits);
 
-    Metadata::AddUserFlag("EMHD");
+
+    Metadata::AddUserFlag("EMHDVar");
 
     // General options for primitive and conserved scalar variables in ImEx driver
     // EMHD is supported only with imex driver and implicit evolution,
     // synchronizing primitive variables
     Metadata m_con  = Metadata({Metadata::Real, Metadata::Cell, Metadata::Independent, Metadata::GetUserFlag("Implicit"),
-                                Metadata::WithFluxes, Metadata::Conserved, Metadata::GetUserFlag("EMHD")});
+                                Metadata::WithFluxes, Metadata::Conserved, Metadata::GetUserFlag("EMHDVar")});
     Metadata m_prim = Metadata({Metadata::Real, Metadata::Cell, Metadata::Derived, Metadata::GetUserFlag("Implicit"),
-                                Metadata::Restart, Metadata::FillGhost, Metadata::GetUserFlag("Primitive"), Metadata::GetUserFlag("EMHD")});
+                                Metadata::Restart, Metadata::FillGhost, Metadata::GetUserFlag("Primitive"), Metadata::GetUserFlag("EMHDVar")});
 
     // Heat conduction
     if (conduction) {
@@ -184,7 +185,7 @@ std::shared_ptr<KHARMAPackage> Initialize(ParameterInput *pin, std::shared_ptr<P
 //     auto pmb = rc->GetBlockPointer();
 
 //     PackIndexMap prims_map, cons_map;
-//     auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHD"), Metadata::Conserved}, cons_map);
+//     auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHDVar"), Metadata::Conserved}, cons_map);
 //     auto P = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
 //     const VarMap m_p(prims_map, false), m_u(cons_map, true);
 
@@ -216,7 +217,7 @@ void BlockPtoU(MeshBlockData<Real> *rc, IndexDomain domain, bool coarse)
     auto pmb = rc->GetBlockPointer();
 
     PackIndexMap prims_map, cons_map;
-    auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHD"), Metadata::Conserved}, cons_map);
+    auto U_E = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("EMHDVar"), Metadata::Conserved}, cons_map);
     auto P = rc->PackVariables(std::vector<MetadataFlag>{Metadata::GetUserFlag("Primitive")}, prims_map);
     const VarMap m_p(prims_map, false), m_u(cons_map, true);
 
diff --git a/kharma/grmhd/grmhd.hpp b/kharma/grmhd/grmhd.hpp
index 51736a2f..53878118 100644
--- a/kharma/grmhd/grmhd.hpp
+++ b/kharma/grmhd/grmhd.hpp
@@ -74,7 +74,7 @@ void FillOutput(MeshBlock *pmb, ParameterInput *pin);
 
 /**
  * Diagnostics performed after each step.
- * Currently finds any negative flags or 0/NaN values in ctop
+ * Currently just looks for negative density/internal energy
  */
 TaskStatus PostStepDiagnostics(const SimTime& tm, MeshData<Real> *rc);
 }
diff --git a/kharma/prob/post_initialize.cpp b/kharma/prob/post_initialize.cpp
index 6b338019..5944bbe1 100644
--- a/kharma/prob/post_initialize.cpp
+++ b/kharma/prob/post_initialize.cpp
@@ -92,8 +92,6 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptr<MeshData<Rea
     const bool use_b_cd = pmesh->packages.AllPackages().count("B_CD");
     const int verbose = pmesh->packages.Get("Globals")->Param<int>("verbose");
 
-    fprintf(stderr, "0.5");
-
     Flag("SeedBField");
     // Seed the magnetic field on each block
     for (auto &pmb : pmesh->block_list) {
@@ -108,8 +106,6 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptr<MeshData<Rea
     }
     EndFlag();
 
-    fprintf(stderr, "0.9");
-
     // Then, if we're in a torus problem or we explicitly ask for it,
     // normalize the magnetic field according to the density
     auto prob = pin->GetString("parthenon/job", "problem_id");
@@ -127,11 +123,8 @@ void KHARMA::SeedAndNormalizeB(ParameterInput *pin, std::shared_ptr<MeshData<Rea
         // Calculate current beta_min value
         Real bsq_max, p_max, beta_min;
         if (beta_calc_legacy) {
-            fprintf(stderr, "1");
             bsq_max = MPIReduce_once(MaxBsq(md.get()), MPI_MAX);
-            fprintf(stderr, "2");
             p_max = MPIReduce_once(MaxPressure(md.get()), MPI_MAX);
-            fprintf(stderr, "3");
             beta_min = p_max / (0.5 * bsq_max);
         } else {
             beta_min = MPIReduce_once(MinBeta(md.get()), MPI_MIN);
@@ -194,27 +187,22 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
     // If your problem requires custom boundary conditions, these should be implemented
     // with the problem and assigned to the relevant functions in the "Boundaries" package.
 
-    fprintf(stderr, "0.0");
     auto &md = pmesh->mesh_data.Get();
 
     auto& pkgs = pmesh->packages.AllPackages();
 
-    fprintf(stderr, "0.1");
     // Magnetic field operations
     if (pin->GetString("b_field", "solver") != "none") {
         // If we need to seed a field based on the problem's fluid initialization...
         if (pin->GetOrAddString("b_field", "type", "none") != "none" && !is_restart) {
             // B field init is not stencil-1, needs boundaries sync'd.
             // FreezeDirichlet ensures any Dirichlet conditions aren't overwritten by zeros
-            fprintf(stderr, "0.2");
             KBoundaries::FreezeDirichlet(md);
             KHARMADriver::SyncAllBounds(md);
 
-            fprintf(stderr, "0.3");
             // Then init B field on each block...
             KHARMA::SeedAndNormalizeB(pin, md);
         }
-        fprintf(stderr, "4");
 
         // Regardless, if evolving a field we should print max(divB)
         // divB is not stencil-1 and we may not have run the above.
@@ -222,8 +210,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
         KBoundaries::FreezeDirichlet(md);
         KHARMADriver::SyncAllBounds(md);
 
-        fprintf(stderr, "5");
-
         if (pkgs.count("B_FluxCT")) {
             B_FluxCT::PrintGlobalMaxDivB(md.get());
         } else if (pkgs.count("B_CT")) {
@@ -233,8 +219,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
         }
     }
 
-    fprintf(stderr, "6");
-
     // Add any hotspots.
     // Note any other modifications made when restarting should be made around here
     if (pin->GetOrAddBoolean("blob", "add_blob", false)) {
@@ -245,8 +229,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
         }
     }
 
-    fprintf(stderr, "7");
-
     // Any extra cleanup & init especially when restarting
     if (is_restart) {
         // Parthenon restores all parameters (global vars) when restarting,
@@ -254,8 +236,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
         KHARMA::ResetGlobals(pin, pmesh);
     }
 
-    fprintf(stderr, "8");
-
     // Clean the B field if we've introduced a divergence somewhere
     // We call this function any time the package is loaded:
     // if we decided to load it in kharma.cpp, we need to clean.
@@ -270,8 +250,6 @@ void KHARMA::PostInitialize(ParameterInput *pin, Mesh *pmesh, bool is_restart)
         B_Cleanup::CleanupDivergence(md);
     }
 
-    fprintf(stderr, "9");
-
     // Finally, synchronize boundary values.
     // Freeze any Dirichlet physical boundaries as they are now, after cleanup/sync/etc.
     KBoundaries::FreezeDirichlet(md);
diff --git a/kharma/reductions/reductions_impl.hpp b/kharma/reductions/reductions_impl.hpp
index 0575cd88..df42b1b8 100644
--- a/kharma/reductions/reductions_impl.hpp
+++ b/kharma/reductions/reductions_impl.hpp
@@ -73,10 +73,10 @@ void Reductions::Start(MeshData<Real> *md, int channel, T val, MPI_Op op)
     auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams();
     auto *reduce_pool = pars.GetMutable<std::vector<Reduce<T>>>(pool_name);
     while (reduce_pool->size() <= channel) reduce_pool->push_back(Reduce<T>());
-    auto& vector_int_reduce = (*reduce_pool)[channel];
+    auto& reduce = (*reduce_pool)[channel];
     // Fill with flags
-    vector_int_reduce.val = val;
-    vector_int_reduce.StartReduce(0, op);
+    reduce.val = val;
+    reduce.StartReduce(0, op);
 }
 template<typename T>
 void Reductions::StartToAll(MeshData<Real> *md, int channel, T val, MPI_Op op)
@@ -86,10 +86,10 @@ void Reductions::StartToAll(MeshData<Real> *md, int channel, T val, MPI_Op op)
     auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams();
     auto *allreduce_pool = pars.GetMutable<std::vector<AllReduce<T>>>(pool_name);
     while (allreduce_pool->size() <= channel) allreduce_pool->push_back(AllReduce<T>());
-    auto& vector_int_reduce = (*allreduce_pool)[channel];
+    auto& reduce = (*allreduce_pool)[channel];
     // Fill with flags
-    vector_int_reduce.val = val;
-    vector_int_reduce.StartReduce(op);
+    reduce.val = val;
+    reduce.StartReduce(op);
 }
 
 // MPI reduction checks
@@ -111,7 +111,7 @@ T Reductions::CheckOnAll(MeshData<Real> *md, int channel)
     // Get the relevant reducer and result
     const std::string pool_name = GetPoolName<T, true>();
     auto& pars = md->GetMeshPointer()->packages.Get("Reductions")->AllParams();
-    auto *reduce_pool = pars.GetMutable<std::vector<Reduce<T>>>(pool_name);
+    auto *reduce_pool = pars.GetMutable<std::vector<AllReduce<T>>>(pool_name);
     auto& reducer = (*reduce_pool)[channel];
 
     while (reducer.CheckReduce() == TaskStatus::incomplete);