diff --git a/src/analysis/history.cpp b/src/analysis/history.cpp
index d88e3d65..a55f8d1d 100644
--- a/src/analysis/history.cpp
+++ b/src/analysis/history.cpp
@@ -26,40 +26,37 @@ Real ReduceMassAccretionRate(MeshData<Real> *md) {
   const auto jb = md->GetBoundsJ(IndexDomain::interior);
   const auto kb = md->GetBoundsK(IndexDomain::interior);
 
+  namespace p = fluid_prim;
+  using parthenon::MakePackDescriptor;
   Mesh *pmesh = md->GetMeshPointer();
+  auto &resolved_pkgs = pmesh->resolved_packages;
+  static auto desc = MakePackDescriptor<p::density, p::velocity>(resolved_pkgs.get());
+  auto v = desc.GetPack(md);
+  const int nblocks = v.GetNBlocks();
+
   auto &pars = pmesh->packages.Get("geometry")->AllParams();
   const Real xh = pars.Get<Real>("xh");
 
-  namespace p = fluid_prim;
-  const std::vector<std::string> vars({p::density::name(), p::velocity::name()});
-
-  PackIndexMap imap;
-  auto pack = md->PackVariables(vars, imap);
-
-  const int prho = imap[p::density::name()].first;
-  const int pvel_lo = imap[p::velocity::name()].first;
-  const int pvel_hi = imap[p::velocity::name()].second;
-
   auto geom = Geometry::GetCoordinateSystem(md);
 
   Real result = 0.0;
   parthenon::par_reduce(
       parthenon::LoopPatternMDRange(), "Phoebus History for Mass Accretion Rate",
-      DevExecSpace(), 0, pack.GetDim(5) - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+      DevExecSpace(), 0, nblocks - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
       KOKKOS_LAMBDA(const int b, const int k, const int j, const int i, Real &lresult) {
-        const auto &coords = pack.GetCoords(b);
+        const auto &coords = v.GetCoordinates(b);
         if (coords.Xf<1>(i) <= xh && xh < coords.Xf<1>(i + 1)) {
           const Real dx1 = coords.CellWidthFA(X1DIR, k, j, i);
           const Real dx2 = coords.CellWidthFA(X2DIR, k, j, i);
           const Real dx3 = coords.CellWidthFA(X3DIR, k, j, i);
 
           // interp to make sure we're getting the horizon correct
-          auto m = (CalcMassFlux(pack, geom, prho, pvel_lo, pvel_hi, b, k, j, i + 1) -
-                    CalcMassFlux(pack, geom, prho, pvel_lo, pvel_hi, b, k, j, i - 1)) /
+          auto m = (CalcMassFlux(md, v, geom, b, k, j, i + 1) -
+                    CalcMassFlux(md, v, geom, b, k, j, i - 1)) /
                    (2.0 * dx1);
-          auto flux = (CalcMassFlux(pack, geom, prho, pvel_lo, pvel_hi, b, k, j, i) +
-                       (xh - coords.Xc<1>(i)) * m) *
-                      dx2 * dx3;
+          auto flux =
+              (CalcMassFlux(md, v, geom, b, k, j, i) + (xh - coords.Xc<1>(i)) * m) * dx2 *
+              dx3;
 
           lresult += flux;
         } else {
diff --git a/src/analysis/history_utils.hpp b/src/analysis/history_utils.hpp
index 5eba068b..0625afc9 100644
--- a/src/analysis/history_utils.hpp
+++ b/src/analysis/history_utils.hpp
@@ -22,6 +22,7 @@
 #include "phoebus_utils/relativity_utils.hpp"
 #include "phoebus_utils/robust.hpp"
 #include "phoebus_utils/variables.hpp"
+#include <interface/sparse_pack.hpp>
 #include <kokkos_abstraction.hpp>
 #include <parthenon/package.hpp>
 #include <utils/error_checking.hpp>
@@ -32,25 +33,29 @@ using namespace Geometry;
 namespace History {
 
 template <typename Pack, typename Geometry>
-KOKKOS_INLINE_FUNCTION Real CalcMassFlux(Pack &pack, Geometry &geom, const int prho,
-                                         const int pvel_lo, const int pvel_hi,
+KOKKOS_INLINE_FUNCTION Real CalcMassFlux(MeshData<Real> *md, Pack &pack, Geometry &geom,
                                          const int b, const int k, const int j,
                                          const int i) {
+  namespace p = fluid_prim;
+  using parthenon::MakePackDescriptor;
+  Mesh *pmesh = md->GetMeshPointer();
+  auto &resolved_pkgs = pmesh->resolved_packages;
+  static auto desc = MakePackDescriptor<p::density, p::velocity>(resolved_pkgs.get());
 
   Real gdet = geom.DetGamma(CellLocation::Cent, b, k, j, i);
   Real lapse = geom.Lapse(CellLocation::Cent, b, k, j, i);
   Real shift[3];
   geom.ContravariantShift(CellLocation::Cent, b, k, j, i, shift);
 
-  const Real vel[] = {pack(b, pvel_lo, k, j, i), pack(b, pvel_lo + 1, k, j, i),
-                      pack(b, pvel_hi, k, j, i)};
+  const Real vel[] = {pack(b, p::velocity(), k, j, i), pack(b, p::velocity(1), k, j, i),
+                      pack(b, p::velocity(2), k, j, i)};
 
   Real gcov4[4][4];
   geom.SpacetimeMetric(CellLocation::Cent, b, k, j, i, gcov4);
   const Real W = phoebus::GetLorentzFactor(vel, gcov4);
   const Real ucon = vel[0] - shift[0] * W / lapse;
 
-  return -lapse * gdet * pack(b, prho, k, j, i) * ucon;
+  return -lapse * gdet * pack(b, p::density(), k, j, i) * ucon;
 }
 
 template <typename Pack, typename Geometry>