diff --git a/example/sparse_advection/parthenon_app_inputs.cpp b/example/sparse_advection/parthenon_app_inputs.cpp
index c34a2041eb99..7c4740b8b5a9 100644
--- a/example/sparse_advection/parthenon_app_inputs.cpp
+++ b/example/sparse_advection/parthenon_app_inputs.cpp
@@ -85,6 +85,43 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
           });
     }
   }
+
+  if (pkg->Param<bool>("restart_test")) {
+    bool any_nonzero = false;
+
+    for (int k = kb.s; k <= kb.e; k++) {
+      for (int j = jb.s; j <= jb.e; j++) {
+        for (int i = ib.s; i <= ib.e; i++) {
+          auto x = coords.x1v(i);
+          auto y = coords.x2v(j);
+          auto z = coords.x3v(k);
+          auto r2 = x * x + y * y + z * z;
+          if (r2 < 0.5 * size) {
+            any_nonzero = true;
+          }
+        }
+      }
+    }
+
+    if (any_nonzero) {
+      data->AllocSparseID("z_shape_shift", 1);
+      data->AllocSparseID("z_shape_shift", 3);
+      data->AllocSparseID("z_shape_shift", 4);
+
+      auto v = data->PackVariables(
+          std::vector<std::string>{"z_dense_A", "z_dense_B", "z_shape_shift"});
+
+      pmb->par_for(
+          "SparseAdvection::ProblemGenerator", 0, v.GetDim(4) - 1, kb.s, kb.e, jb.s, jb.e,
+          ib.s, ib.e, KOKKOS_LAMBDA(const int n, const int k, const int j, const int i) {
+            auto x = coords.x1v(i);
+            auto y = coords.x2v(j);
+            auto z = coords.x3v(k);
+            auto r2 = x * x + y * y + z * z;
+            v(n, k, j, i) = (r2 < 0.5 * size ? 1.0 : 0.0);
+          });
+    }
+  }
 }
 
 //========================================================================================
diff --git a/example/sparse_advection/parthinput.sparse_advection b/example/sparse_advection/parthinput.sparse_advection
index c07f3ebf0024..441b9501bda7 100644
--- a/example/sparse_advection/parthinput.sparse_advection
+++ b/example/sparse_advection/parthinput.sparse_advection
@@ -58,6 +58,8 @@ integrator = rk2
 ncycle_out_mesh = -10000
 
 <SparseAdvection>
+restart_test = true
+
 cfl = 0.45
 speed = 1.5
 
@@ -71,7 +73,7 @@ dt = 0.5
 <parthenon/output0>
 file_type = hdf5
 dt = 0.1
-variables = sparse
+variables = sparse, z_dense_A, z_dense_B, z_shape_shift
 
 <parthenon/output3>
 file_type = hst
diff --git a/example/sparse_advection/sparse_advection_package.cpp b/example/sparse_advection/sparse_advection_package.cpp
index c10a4674a459..6f7f72b698ac 100644
--- a/example/sparse_advection/sparse_advection_package.cpp
+++ b/example/sparse_advection/sparse_advection_package.cpp
@@ -22,6 +22,7 @@
 #include <parthenon/package.hpp>
 
 #include "defs.hpp"
+#include "interface/metadata.hpp"
 #include "interface/sparse_pool.hpp"
 #include "kokkos_abstraction.hpp"
 #include "reconstruct/dc_inline.hpp"
@@ -42,6 +43,9 @@ using parthenon::UserHistoryOperation;
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   auto pkg = std::make_shared<StateDescriptor>("sparse_advection_package");
 
+  bool restart_test = pin->GetOrAddBoolean("SparseAdvection", "restart_test", false);
+  pkg->AddParam("restart_test", restart_test);
+
   Real cfl = pin->GetOrAddReal("SparseAdvection", "cfl", 0.45);
   pkg->AddParam("cfl", cfl);
   Real refine_tol = pin->GetOrAddReal("SparseAdvection", "refine_tol", 0.3);
@@ -64,15 +68,37 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   pkg->AddParam("vz", RealArr_t{0.0, 0.0, 0.0, 0.0});
 
   // add sparse field
-  Metadata m({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes,
-              Metadata::FillGhost, Metadata::Sparse},
-             std::vector<int>({1}));
-  parthenon::SparsePool pool("sparse", m);
+  {
+    Metadata m({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes,
+                Metadata::FillGhost, Metadata::Sparse},
+               std::vector<int>({1}));
+    SparsePool pool("sparse", m);
+
+    for (int sid = 0; sid < NUM_FIELDS; ++sid) {
+      pool.Add(sid);
+    }
+    pkg->AddSparsePool(pool);
+  }
+
+  // add fields for restart test ("Z" prefix so they are after sparse in alphabetical
+  // list, helps with reusing velocity vectors)
+  if (restart_test) {
+    Metadata m_dense({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes,
+                      Metadata::FillGhost});
+    pkg->AddField("z_dense_A", m_dense);
+    pkg->AddField("z_dense_B", m_dense);
+
+    Metadata m_sparse({Metadata::Cell, Metadata::Independent, Metadata::WithFluxes,
+                       Metadata::FillGhost, Metadata::Sparse});
 
-  for (int sid = 0; sid < NUM_FIELDS; ++sid) {
-    pool.Add(sid);
+    SparsePool pool("z_shape_shift", m_sparse);
+    pool.Add(1, std::vector<int>{1}, std::vector<std::string>{"scalar"});
+    pool.Add(3, std::vector<int>{3}, Metadata::Vector,
+             std::vector<std::string>{"vec_x", "vec_y", "vec_z"});
+    pool.Add(4, std::vector<int>{4}, Metadata::Vector);
+
+    pkg->AddSparsePool(pool);
   }
-  pkg->AddSparsePool(pool);
 
   pkg->CheckRefinementBlock = CheckRefinement;
   pkg->EstimateTimestepBlock = EstimateTimestepBlock;
@@ -182,22 +208,19 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
       KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int k, const int j) {
         parthenon::ScratchPad2D<Real> ql(member.team_scratch(scratch_level), nvar, nx1);
         parthenon::ScratchPad2D<Real> qr(member.team_scratch(scratch_level), nvar, nx1);
+
         // get reconstructed state on faces
         parthenon::DonorCellX1(member, k, j, ib.s - 1, ib.e + 1, v, ql, qr);
+
         // Sync all threads in the team so that scratch memory is consistent
         member.team_barrier();
 
         for (int n = 0; n < nvar; n++) {
           if (!v.IsAllocated(n)) continue;
-          if (vx[n] > 0.0) {
-            par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
-              v.flux(X1DIR, n, k, j, i) = ql(n, i) * vx[n];
-            });
-          } else {
-            par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
-              v.flux(X1DIR, n, k, j, i) = qr(n, i) * vx[n];
-            });
-          }
+          const auto this_v = vx[n % NUM_FIELDS];
+          par_for_inner(member, ib.s, ib.e + 1, [&](const int i) {
+            v.flux(X1DIR, n, k, j, i) = (this_v > 0.0 ? ql(n, i) : qr(n, i)) * this_v;
+          });
         }
       });
 
@@ -206,31 +229,28 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
     pmb->par_for_outer(
         "x2 flux", 3 * scratch_size_in_bytes, scratch_level, kb.s, kb.e, jb.s, jb.e + 1,
         KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int k, const int j) {
-          // the overall algorithm/use of scratch pad here is clear inefficient and kept
+          // the overall algorithm/use of scratch pad here is clearly inefficient and kept
           // just for demonstrating purposes. The key point is that we cannot reuse
           // reconstructed arrays for different `j` with `j` being part of the outer
           // loop given that this loop can be handled by multiple threads simultaneously.
-
           parthenon::ScratchPad2D<Real> ql(member.team_scratch(scratch_level), nvar, nx1);
           parthenon::ScratchPad2D<Real> qr(member.team_scratch(scratch_level), nvar, nx1);
           parthenon::ScratchPad2D<Real> q_unused(member.team_scratch(scratch_level), nvar,
                                                  nx1);
+
           // get reconstructed state on faces
           parthenon::DonorCellX2(member, k, j - 1, ib.s, ib.e, v, ql, q_unused);
           parthenon::DonorCellX2(member, k, j, ib.s, ib.e, v, q_unused, qr);
+
           // Sync all threads in the team so that scratch memory is consistent
           member.team_barrier();
+
           for (int n = 0; n < nvar; n++) {
             if (!v.IsAllocated(n)) continue;
-            if (vy[n] > 0.0) {
-              par_for_inner(member, ib.s, ib.e, [&](const int i) {
-                v.flux(X2DIR, n, k, j, i) = ql(n, i) * vy[n];
-              });
-            } else {
-              par_for_inner(member, ib.s, ib.e, [&](const int i) {
-                v.flux(X2DIR, n, k, j, i) = qr(n, i) * vy[n];
-              });
-            }
+            const auto this_v = vy[n % NUM_FIELDS];
+            par_for_inner(member, ib.s, ib.e, [&](const int i) {
+              v.flux(X2DIR, n, k, j, i) = (this_v > 0.0 ? ql(n, i) : qr(n, i)) * this_v;
+            });
           }
         });
   }
diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py
index 4d31313dc4cd..987cb5c31ead 100644
--- a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py
+++ b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py
@@ -172,11 +172,11 @@ def compare_attribute_group(f0, f1, name):
 
             print("\nFirst file:")
             for k in diffs:
-                print("%20s: " % k, group0[k])
+                print("%20s: " % k, group0[k] if k in group0 else "<does not exist>")
 
             print("\nSecond file:")
             for k in diffs:
-                print("%20s: " % k, group1[k])
+                print("%20s: " % k, group1[k] if k in group1 else "<does not exist>")
         else:
             print("  %20s: no diffs" % name)
 
diff --git a/src/interface/meshblock_data.hpp b/src/interface/meshblock_data.hpp
index b428c371b41c..76a5994ae2d1 100644
--- a/src/interface/meshblock_data.hpp
+++ b/src/interface/meshblock_data.hpp
@@ -142,10 +142,8 @@ class MeshBlockData {
 
   std::shared_ptr<CellVariable<T>> GetCellVarPtr(const std::string &label) const {
     auto it = varMap_.find(label);
-    if (it == varMap_.end()) {
-      PARTHENON_THROW(std::string("\n") + std::string(label) +
-                      std::string(" array not found in Get()\n"));
-    }
+    PARTHENON_REQUIRE_THROWS(it != varMap_.end(),
+                             "Couldn't find variable '" + label + "'");
     return it->second;
   }
 
diff --git a/src/interface/metadata.cpp b/src/interface/metadata.cpp
index 5bd9ebe3d348..bab299b61fbc 100644
--- a/src/interface/metadata.cpp
+++ b/src/interface/metadata.cpp
@@ -120,7 +120,6 @@ std::vector<MetadataFlag> Metadata::Flags() const {
 std::array<int, 6> Metadata::GetArrayDims(std::weak_ptr<MeshBlock> wpmb,
                                           bool coarse) const {
   std::array<int, 6> arrDims;
-
   const auto &shape = shape_;
   const int N = shape.size();
 
diff --git a/src/interface/metadata.hpp b/src/interface/metadata.hpp
index d2a92e90c5c0..d5683ec06ed0 100644
--- a/src/interface/metadata.hpp
+++ b/src/interface/metadata.hpp
@@ -225,6 +225,16 @@ class Metadata {
       PARTHENON_REQUIRE_THROWS(
           shape_.size() <= 3,
           "Variables tied to mesh entities can only have a shape of rank <= 3");
+
+      int num_comp = 1;
+      for (auto s : shape) {
+        num_comp *= s;
+      }
+
+      PARTHENON_REQUIRE_THROWS(component_labels.size() == 0 ||
+                                   (component_labels.size() == num_comp),
+                               "Must provide either 0 component labels or the same "
+                               "number as the number of components");
     }
   }
 
diff --git a/src/interface/update.cpp b/src/interface/update.cpp
index 953a80cef7ad..c2cdd810b7bf 100644
--- a/src/interface/update.cpp
+++ b/src/interface/update.cpp
@@ -20,6 +20,7 @@
 #include "globals.hpp"
 #include "interface/meshblock_data.hpp"
 #include "interface/metadata.hpp"
+#include "interface/variable_pack.hpp"
 #include "mesh/mesh.hpp"
 #include "mesh/meshblock.hpp"
 
@@ -147,12 +148,9 @@ TaskStatus SparseDeallocCheck(MeshData<Real> *md) {
   const IndexRange jb = md->GetBoundsJ(IndexDomain::entire);
   const IndexRange kb = md->GetBoundsK(IndexDomain::entire);
 
+  PackIndexMap map;
   const std::vector<MetadataFlag> flags({Metadata::Sparse});
-  const auto &pack = md->PackVariables(flags);
-
-  // get list of variables in same order as they appear in the pack, since variable
-  // metadata is the same on all blocks, we can just use the first block
-  const auto &var_list = md->GetBlockData(0)->GetVariablesByFlag(flags, true);
+  const auto &pack = md->PackVariables(flags, map);
 
   const int num_blocks = pack.GetDim(5);
   const int num_vars = pack.GetDim(4);
@@ -205,10 +203,19 @@ TaskStatus SparseDeallocCheck(MeshData<Real> *md) {
   auto is_zero_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), is_zero);
 
   for (int b = 0; b < num_blocks; ++b) {
-    for (int v = 0; v < num_vars; ++v) {
-      auto &counter = md->GetBlockData(b)->Get(var_list.labels()[v]).dealloc_count;
-      if (is_zero_h(b, v)) {
-        // this var is zero, increment dealloc counter
+    for (auto var_itr : map.Map()) {
+      const auto label = var_itr.first;
+      auto &counter = md->GetBlockData(b)->Get(label).dealloc_count;
+      bool all_zero = true;
+      for (int v = var_itr.second.first; v <= var_itr.second.second; ++v) {
+        if (!is_zero_h(b, v)) {
+          all_zero = false;
+          break;
+        }
+      }
+
+      if (all_zero) {
+        // all components of this var are zero, increment dealloc counter
         counter += 1;
       } else {
         counter = 0;
@@ -217,9 +224,9 @@ TaskStatus SparseDeallocCheck(MeshData<Real> *md) {
       if (counter > Globals::sparse_config.deallocation_count) {
         // this variable has been flagged for deallocation deallocation_count times in a
         // row, now deallocate it
-        // printf("Deallocating var %s on block %i\n", var_list.labels()[v].c_str(),
+        // printf("Deallocating var %s on block %i\n", label.c_str(),
         //        md->GetBlockData(b)->GetBlockPointer()->gid);
-        md->GetBlockData(b)->GetBlockPointer()->DeallocateSparse(var_list.labels()[v]);
+        md->GetBlockData(b)->GetBlockPointer()->DeallocateSparse(label);
       }
     }
   }
diff --git a/src/interface/variable.hpp b/src/interface/variable.hpp
index 90a00f50c136..dfcaae7b269c 100644
--- a/src/interface/variable.hpp
+++ b/src/interface/variable.hpp
@@ -72,6 +72,9 @@ class CellVariable {
     }
 
     if (IsSet(Metadata::FillGhost)) {
+      PARTHENON_REQUIRE_THROWS(
+          GetDim(4) == NumComponents(),
+          "CellCenteredBoundaryVariable currently only supports rank-1 variables");
       vbvar = std::make_shared<CellCenteredBoundaryVariable>(wpmb.lock(), IsSparse(),
                                                              label(), GetDim(4));
     }
diff --git a/src/interface/variable_pack.hpp b/src/interface/variable_pack.hpp
index 57f433ae8b77..1f1ea632ced2 100644
--- a/src/interface/variable_pack.hpp
+++ b/src/interface/variable_pack.hpp
@@ -93,6 +93,8 @@ class PackIndexMap {
  public:
   PackIndexMap() = default;
 
+  const auto &Map() const { return map_; }
+
   const auto &get(const std::string &base_name, int sparse_id = InvalidSparseID) const {
     const auto &key = MakeVarLabel(base_name, sparse_id);
     auto itr = map_.find(key);
diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp
index 0aeab9a078fd..ea515dee385d 100644
--- a/src/outputs/parthenon_hdf5.cpp
+++ b/src/outputs/parthenon_hdf5.cpp
@@ -107,7 +107,7 @@ struct VarInfo {
   }
 
   explicit VarInfo(const std::shared_ptr<CellVariable<Real>> &var)
-      : VarInfo(var->label(), var->metadata().getComponentLabels(), var->GetDim(4),
+      : VarInfo(var->label(), var->metadata().getComponentLabels(), var->NumComponents(),
                 var->IsSparse(), var->IsSet(Metadata::Vector)) {}
 };
 
@@ -137,23 +137,25 @@ static void writeXdmfArrayRef(std::ofstream &fid, const std::string &prefix,
 }
 
 static void writeXdmfSlabVariableRef(std::ofstream &fid, const std::string &name,
+                                     const std::vector<std::string> &component_labels,
                                      std::string &hdfFile, int iblock, const int &vlen,
                                      int &ndims, hsize_t *dims,
                                      const std::string &dims321, bool isVector) {
   // writes a slab reference to file
-
   std::vector<std::string> names;
   int nentries = 1;
-  int vector_size = 1;
   if (vlen == 1 || isVector) {
+    // we only make one entry, because either vlen == 1, or we write this as a vector
     names.push_back(name);
   } else {
     nentries = vlen;
     for (int i = 0; i < vlen; i++) {
-      names.push_back(name + "_" + std::to_string(i));
+      names.push_back(
+          name + "_" +
+          (component_labels.empty() ? std::to_string(i) : component_labels[i]));
     }
   }
-  if (isVector) vector_size = vlen;
+  const int vector_size = isVector ? vlen : 1;
 
   const std::string prefix = "      ";
   for (int i = 0; i < nentries; i++) {
@@ -261,8 +263,8 @@ void genXDMF(std::string hdfFile, Mesh *pm, SimTime *tm, int nx1, int nx2, int n
     for (const auto &vinfo : var_list) {
       const int vlen = vinfo.vlen;
       dims[4] = vlen;
-      writeXdmfSlabVariableRef(xdmf, vinfo.label, hdfFile, ib, vlen, ndims, dims, dims321,
-                               vinfo.is_vector);
+      writeXdmfSlabVariableRef(xdmf, vinfo.label, vinfo.component_labels, hdfFile, ib,
+                               vlen, ndims, dims, dims321, vinfo.is_vector);
     }
     xdmf << "      </Grid>" << std::endl;
   }
@@ -325,7 +327,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
   filename.append(restart_ ? ".rhdf" : ".phdf");
 
   // After file has been opened with the current number, already advance output
-  // parameters so that for restarts the file is not immediatly overwritten again.
+  // parameters so that for restarts the file is not immediately overwritten again.
   output_params.file_number++;
   output_params.next_time += output_params.dt;
   pin->SetInteger(output_params.block_name, "file_number", output_params.file_number);
diff --git a/src/outputs/parthenon_hdf5.hpp b/src/outputs/parthenon_hdf5.hpp
index d546ad363e9a..3d78ad246a26 100644
--- a/src/outputs/parthenon_hdf5.hpp
+++ b/src/outputs/parthenon_hdf5.hpp
@@ -186,13 +186,15 @@ std::vector<T> HDF5ReadAttributeVec(hid_t location, const std::string &name) {
   auto type = getHDF5Type(res.data());
 
   // check if attribute exists
-  PARTHENON_HDF5_CHECK(H5Aexists(location, name.c_str()));
+  auto status = PARTHENON_HDF5_CHECK(H5Aexists(location, name.c_str()));
+  PARTHENON_REQUIRE_THROWS(status > 0, "Attribute '" + name + "' does not exist");
 
   const H5A attr = H5A::FromHIDCheck(H5Aopen(location, name.c_str(), H5P_DEFAULT));
 
   // check data type
   const H5T hdf5_type = H5T::FromHIDCheck(H5Aget_type(attr));
-  PARTHENON_HDF5_CHECK(H5Tequal(type, hdf5_type));
+  status = PARTHENON_HDF5_CHECK(H5Tequal(type, hdf5_type));
+  PARTHENON_REQUIRE_THROWS(status > 0, "Type mismatch for attribute " + name);
 
   // Allocate array of correct size
   const H5S dataspace = H5S::FromHIDCheck(H5Aget_space(attr));
diff --git a/src/outputs/restart.cpp b/src/outputs/restart.cpp
index a7bcfbd6b662..1fd948ee798c 100644
--- a/src/outputs/restart.cpp
+++ b/src/outputs/restart.cpp
@@ -25,6 +25,7 @@
 #include "mesh/mesh.hpp"
 #include "mesh/meshblock.hpp"
 #include "outputs/outputs.hpp"
+#include "outputs/parthenon_hdf5.hpp"
 #include "outputs/restart.hpp"
 #include "utils/error_checking.hpp"
 
@@ -40,12 +41,48 @@ RestartReader::RestartReader(const char *filename) : filename_(filename) {
       << "Executable not configured for HDF5 outputs, but HDF5 file format "
       << "is required for restarts" << std::endl;
   PARTHENON_FAIL(msg);
-#else
+#else  // HDF5 enabled
   // Open the HDF file in read only mode
   fh_ = H5F::FromHIDCheck(H5Fopen(filename, H5F_ACC_RDONLY, H5P_DEFAULT));
 
   hasGhost = GetAttr<int>("Info", "IncludesGhost");
-#endif
+#endif // ENABLE_HDF5
+}
+
+RestartReader::SparseInfo RestartReader::GetSparseInfo() const {
+#ifndef ENABLE_HDF5
+  PARTHENON_FAIL("Restart functionality is not available because HDF5 is disabled");
+#else  // HDF5 enabled
+  SparseInfo info;
+
+  // check if SparseInfo exists, if not, return the default-constructed SparseInfo
+  // instance
+  auto status = PARTHENON_HDF5_CHECK(H5Oexists_by_name(fh_, "SparseInfo", H5P_DEFAULT));
+  if (status > 0) {
+    // SparseInfo exists, read its contents
+    auto hdl = OpenDataset<bool>("SparseInfo");
+    PARTHENON_REQUIRE_THROWS(hdl.rank == 2, "SparseInfo expected to have rank 2");
+
+    info.labels = HDF5ReadAttributeVec<std::string>(hdl.dataset, "SparseFields");
+    info.num_sparse = static_cast<int>(info.labels.size());
+    PARTHENON_REQUIRE_THROWS(info.num_sparse == static_cast<int>(hdl.dims[1]),
+                             "Mismatch in number of sparse fields");
+
+    // Note: We cannot use ReadData, because std::vector<bool> doesn't have a data()
+    // member
+    info.allocated.reset(new hbool_t[hdl.count]);
+    info.num_blocks = static_cast<int>(hdl.dims[0]);
+
+    const H5S memspace =
+        H5S::FromHIDCheck(H5Screate_simple(hdl.rank, hdl.dims.data(), NULL));
+
+    // Read data from file
+    PARTHENON_HDF5_CHECK(H5Dread(hdl.dataset, hdl.type, memspace, hdl.dataspace,
+                                 H5P_DEFAULT, static_cast<void *>(info.allocated.get())));
+  }
+
+  return info;
+#endif // ENABLE_HDF5
 }
 
 } // namespace parthenon
diff --git a/src/outputs/restart.hpp b/src/outputs/restart.hpp
index 11367c2f8bcb..1b65ad32d014 100644
--- a/src/outputs/restart.hpp
+++ b/src/outputs/restart.hpp
@@ -38,97 +38,145 @@ class RestartReader {
  public:
   explicit RestartReader(const char *theFile);
 
-  // Gets data for all blocks on current rank.
-  // Assumes blocks are contiguous
-  // fills internal data for given pointer
-  // returns NBlocks on success, -1 on failure
+  struct SparseInfo {
+    // labels of sparse fields (full label, i.e. base name and sparse id)
+    std::vector<std::string> labels;
+
+    // allocation status of sparse fields (2D array outer dimension: block, inner
+    // dimension: sparse field)
+    // can't use std::vector here because std::vector<hbool_t> is the same as
+    // std::vector<bool> and it doesn't have .data() member
+    std::unique_ptr<hbool_t[]> allocated;
+
+    int num_blocks = 0;
+    int num_sparse = 0;
+
+    bool IsAllocated(int block, int sparse_field_idx) const {
+      PARTHENON_REQUIRE_THROWS(allocated != nullptr,
+                               "Tried to get allocation status but no data present");
+      PARTHENON_REQUIRE_THROWS((block >= 0) && (block < num_blocks),
+                               "Invalid block index in SparseInfo::IsAllocated");
+      PARTHENON_REQUIRE_THROWS((sparse_field_idx >= 0) && (sparse_field_idx < num_sparse),
+                               "Invalid sparse field index in SparseInfo::IsAllocated");
+
+      return allocated[block * num_sparse + sparse_field_idx];
+    }
+  };
+
+  SparseInfo GetSparseInfo() const;
+
+ private:
+  struct DatasetHandle {
+    hid_t type;
+    H5D dataset;
+    H5S dataspace;
+    int rank;
+    hsize_t count;
+    std::vector<hsize_t> dims;
+  };
+
+  // internal convenience function to open a dataset, perform some checks, and get
+  // dimensions
   template <typename T>
-  int ReadBlocks(const char *name, IndexRange range, std::vector<T> &dataVec,
-                 const std::vector<size_t> &bsize, size_t vlen = 1) {
+  DatasetHandle OpenDataset(const std::string &name) const {
 #ifndef ENABLE_HDF5
     PARTHENON_FAIL("Restart functionality is not available because HDF5 is disabled");
-#else
-    try {
-      // dataVec is assumed to be of the correct size
-      T *data = dataVec.data();
-
-      // compute block size, probaby could cache this.
-      hid_t const theHdfType = getHDF5Type(data);
-
-      H5D const dataset = H5D::FromHIDCheck(H5Dopen2(fh_, name, H5P_DEFAULT));
-      H5S const dataspace = H5S::FromHIDCheck(H5Dget_space(dataset));
-
-      /** Define hyperslab in dataset **/
-      hsize_t offset[5] = {static_cast<hsize_t>(range.s), 0, 0, 0, 0};
-      hsize_t count[5] = {static_cast<hsize_t>(range.e - range.s + 1), bsize[2], bsize[1],
-                          bsize[0], vlen};
-      PARTHENON_HDF5_CHECK(
-          H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
-
-      /** Define memory dataspace **/
-      H5S const memspace = H5S::FromHIDCheck(H5Screate_simple(5, count, NULL));
-      PARTHENON_HDF5_CHECK(
-          H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
-
-      // Read data from file
-      PARTHENON_HDF5_CHECK(
-          H5Dread(dataset, H5T_NATIVE_DOUBLE, memspace, dataspace, H5P_DEFAULT, data));
-
-      return static_cast<int>(count[0]);
-    } catch (const std::exception &e) {
-      std::cout << e.what();
-      return -1;
+#else  // HDF5 enabled
+    DatasetHandle handle;
+
+    // make sure dataset exists
+    auto status = PARTHENON_HDF5_CHECK(H5Oexists_by_name(fh_, name.c_str(), H5P_DEFAULT));
+    PARTHENON_REQUIRE_THROWS(
+        status > 0, "Dataset '" + name + "' does not exist in HDF5 file " + filename_);
+
+    // open dataset
+    handle.dataset = H5D::FromHIDCheck(H5Dopen2(fh_, name.c_str(), H5P_DEFAULT));
+    handle.dataspace = H5S::FromHIDCheck(H5Dget_space(handle.dataset));
+
+    // get the HDF5 type from the template parameter and make sure it matches the dataset
+    // type
+    T *typepointer = nullptr;
+    handle.type = getHDF5Type(typepointer);
+    const H5T dset_type = H5T::FromHIDCheck(H5Dget_type(handle.dataset));
+    status = PARTHENON_HDF5_CHECK(H5Tequal(handle.type, dset_type));
+    PARTHENON_REQUIRE_THROWS(status > 0, "Type mismatch for dataset " + name);
+
+    // get rank and dims
+    const H5S filespace = H5S::FromHIDCheck(H5Dget_space(handle.dataset));
+    handle.rank = PARTHENON_HDF5_CHECK(H5Sget_simple_extent_ndims(filespace));
+
+    handle.dims.resize(handle.rank);
+    PARTHENON_HDF5_CHECK(H5Sget_simple_extent_dims(filespace, handle.dims.data(), NULL));
+    handle.count = 1;
+    for (int idir = 0; idir < handle.rank; idir++) {
+      handle.count = handle.count * handle.dims[idir];
     }
-#endif
+
+    return handle;
+#endif // ENABLE_HDF5
   }
 
-  // Reads an array dataset from file as a 1D vector.
-  // Returns number of items read in count if provided
+ public:
+  // Gets data for all blocks on current rank.
+  // Assumes blocks are contiguous
+  // fills internal data for given pointer
   template <typename T>
-  std::vector<T> ReadDataset(const char *name, size_t *count = nullptr) {
-    // Returns entire 1D array.
-    // status, never checked.  We should...
+  void ReadBlocks(const std::string &name, IndexRange range, std::vector<T> &dataVec,
+                  const std::vector<size_t> &bsize, size_t vlen = 1) const {
 #ifndef ENABLE_HDF5
     PARTHENON_FAIL("Restart functionality is not available because HDF5 is disabled");
-#else
-    T *typepointer = nullptr;
-    hid_t const theHdfType = getHDF5Type(typepointer);
+#else  // HDF5 enabled
+    auto hdl = OpenDataset<T>(name);
 
-    H5D const dataset = H5D::FromHIDCheck(H5Dopen2(fh_, name, H5P_DEFAULT));
-    H5S const dataspace = H5S::FromHIDCheck(H5Dget_space(dataset));
+    PARTHENON_REQUIRE_THROWS(hdl.rank == 5, "Expected data set of rank 5, but dataset " +
+                                                name + " has rank " +
+                                                std::to_string(hdl.rank));
 
-    // Allocate array of correct size
-    H5S const filespace = H5S::FromHIDCheck(H5Dget_space(dataset));
+    PARTHENON_REQUIRE_THROWS(dataVec.size() >= hdl.count,
+                             "Buffer too small for dataset " + name);
 
-    int rank = PARTHENON_HDF5_CHECK(H5Sget_simple_extent_ndims(filespace));
+    /** Select hyperslab in dataset **/
+    hsize_t offset[5] = {static_cast<hsize_t>(range.s), 0, 0, 0, 0};
+    hsize_t count[5] = {static_cast<hsize_t>(range.e - range.s + 1), bsize[2], bsize[1],
+                        bsize[0], vlen};
+    PARTHENON_HDF5_CHECK(
+        H5Sselect_hyperslab(hdl.dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
 
-    std::vector<hsize_t> dims(rank);
-    PARTHENON_HDF5_CHECK(H5Sget_simple_extent_dims(filespace, dims.data(), NULL));
-    hsize_t isize = 1;
-    for (int idir = 0; idir < rank; idir++) {
-      isize = isize * dims[idir];
-    }
-    if (count != nullptr) {
-      *count = isize;
-    }
+    const H5S memspace = H5S::FromHIDCheck(H5Screate_simple(5, count, NULL));
+    PARTHENON_HDF5_CHECK(
+        H5Sselect_hyperslab(hdl.dataspace, H5S_SELECT_SET, offset, NULL, count, NULL));
+
+    // Read data from file
+    PARTHENON_HDF5_CHECK(H5Dread(hdl.dataset, hdl.type, memspace, hdl.dataspace,
+                                 H5P_DEFAULT, dataVec.data()));
+#endif // ENABLE_HDF5
+  }
 
-    std::vector<T> data(isize);
-    /** Define memory dataspace **/
-    H5S const memspace = H5S::FromHIDCheck(H5Screate_simple(rank, dims.data(), NULL));
+  // Reads an array dataset from file as a 1D vector.
+  template <typename T>
+  std::vector<T> ReadDataset(const std::string &name) const {
+#ifndef ENABLE_HDF5
+    PARTHENON_FAIL("Restart functionality is not available because HDF5 is disabled");
+#else  // HDF5 enabled
+    auto hdl = OpenDataset<T>(name);
+
+    std::vector<T> data(hdl.count);
+    const H5S memspace =
+        H5S::FromHIDCheck(H5Screate_simple(hdl.rank, hdl.dims.data(), NULL));
 
     // Read data from file
-    PARTHENON_HDF5_CHECK(H5Dread(dataset, theHdfType, memspace, dataspace, H5P_DEFAULT,
-                                 static_cast<void *>(data.data())));
+    PARTHENON_HDF5_CHECK(H5Dread(hdl.dataset, hdl.type, memspace, hdl.dataspace,
+                                 H5P_DEFAULT, static_cast<void *>(data.data())));
 
     return data;
-#endif
+#endif // ENABLE_HDF5
   }
 
   template <typename T>
-  std::vector<T> GetAttrVec(const std::string &location, const std::string &name) {
+  std::vector<T> GetAttrVec(const std::string &location, const std::string &name) const {
 #ifndef ENABLE_HDF5
     PARTHENON_FAIL("Restart functionality is not available because HDF5 is disabled");
-#else
+#else  // HDF5 enabled
     // check if the location exists in the file
     PARTHENON_HDF5_CHECK(H5Oexists_by_name(fh_, location.c_str(), H5P_DEFAULT));
 
@@ -136,11 +184,11 @@ class RestartReader {
     const H5O obj = H5O::FromHIDCheck(H5Oopen(fh_, location.c_str(), H5P_DEFAULT));
 
     return HDF5ReadAttributeVec<T>(obj, name);
-#endif
+#endif // ENABLE_HDF5
   }
 
   template <typename T>
-  T GetAttr(const std::string &location, const std::string &name) {
+  T GetAttr(const std::string &location, const std::string &name) const {
     // Note: We don't need a template specialization for std::string, since that case will
     // be handled by HDF5ReadAttributeVec
     auto res = GetAttrVec<T>(location, name);
@@ -166,7 +214,7 @@ class RestartReader {
   // Currently all restarts are HDF5 files
   // when that changes, this will be revisited
   H5F fh_;
-#endif
+#endif // ENABLE_HDF5
 };
 
 } // namespace parthenon
diff --git a/src/parthenon_manager.cpp b/src/parthenon_manager.cpp
index f2bbc841a26a..25340fc6feb9 100644
--- a/src/parthenon_manager.cpp
+++ b/src/parthenon_manager.cpp
@@ -13,7 +13,10 @@
 
 #include "parthenon_manager.hpp"
 
+#include <algorithm>
+#include <exception>
 #include <string>
+#include <unordered_map>
 #include <utility>
 #include <vector>
 
@@ -25,6 +28,7 @@
 #include "mesh/domain.hpp"
 #include "mesh/meshblock.hpp"
 #include "refinement/refinement.hpp"
+#include "utils/error_checking.hpp"
 
 namespace parthenon {
 
@@ -220,69 +224,87 @@ void ParthenonManager::RestartPackages(Mesh &rm, RestartReader &resfile) {
 
   size_t nCells = bsize[0] * bsize[1] * bsize[2];
 
-  // Get list of variables, assumed same for all blocks
-
-  // TODO(JL) this won't work when reading true sparse variables, will be updated in PR
-  // #383
+  // Get list of variables, they are the same for all blocks (since all blocks have the
+  // same variable metadata)
   const auto indep_restart_vars =
       mb.meshblock_data.Get()
           ->GetVariablesByFlag(
               {parthenon::Metadata::Independent, parthenon::Metadata::Restart}, false)
           .vars();
 
+  const auto sparse_info = resfile.GetSparseInfo();
+  // create map of sparse field labels to index in the SparseInfo table
+  std::unordered_map<std::string, int> sparse_idxs;
+  for (int i = 0; i < sparse_info.num_sparse; ++i) {
+    sparse_idxs.insert({sparse_info.labels[i], i});
+  }
+
   // Allocate space based on largest vector
-  size_t vlen = 1;
+  int max_vlen = 1;
+  int num_sparse = 0;
   for (auto &v : indep_restart_vars) {
-    if (v->GetDim(4) > vlen) {
-      vlen = v->GetDim(4);
+    // check that variable is in the list of sparse fields if and only if it is sparse
+    if (v->IsSparse()) {
+      ++num_sparse;
+      PARTHENON_REQUIRE_THROWS(sparse_idxs.count(v->label()) == 1,
+                               "Sparse field " + v->label() +
+                                   " is not marked as sparse in restart file");
+    } else {
+      PARTHENON_REQUIRE_THROWS(sparse_idxs.count(v->label()) == 0,
+                               "Dense field " + v->label() +
+                                   " is marked as sparse in restart file");
     }
+    max_vlen = std::max(max_vlen, v->GetDim(4));
   }
-  std::vector<Real> tmp(static_cast<size_t>(nb) * nCells * vlen);
-  for (auto &v : indep_restart_vars) {
-    const size_t v4 = v->GetDim(4);
-    const std::string vName = v->label();
-
-    if (Globals::my_rank == 0) std::cout << "Var:" << vName << ":" << v4 << std::endl;
-    // Read relevant data from the hdf file
-    int stat = resfile.ReadBlocks(vName.c_str(), myBlocks, tmp, bsize, v4);
-    if (stat < 0) {
-      std::cout << " WARNING: Variable " << v->label() << " Not found in restart file";
+
+  // make sure we have all sparse variables that are in the restart file
+  PARTHENON_REQUIRE_THROWS(
+      num_sparse == sparse_info.num_sparse,
+      "Mismatch between sparse fields in simulation and restart file");
+
+  std::vector<Real> tmp(static_cast<size_t>(nb) * nCells * max_vlen);
+  for (auto &v_info : indep_restart_vars) {
+    const auto vlen = v_info->GetDim(4);
+    const auto &label = v_info->label();
+
+    if (Globals::my_rank == 0) std::cout << "Var:" << label << ":" << vlen << std::endl;
+    // Read relevant data from the hdf file, this works for dense and sparse variables
+    try {
+      resfile.ReadBlocks(label, myBlocks, tmp, bsize, vlen);
+    } catch (std::exception &ex) {
+      std::cout << " WARNING: Failed to read variable " << label
+                << " from restart file. Error message: " << ex.what() << std::endl;
       continue;
     }
 
     size_t index = 0;
     for (auto &pmb : rm.block_list) {
-      bool found = false;
-      const auto this_indep_restart_vars =
-          pmb->meshblock_data.Get()
-              ->GetVariablesByFlag(
-                  {parthenon::Metadata::Independent, parthenon::Metadata::Restart}, false)
-              .vars();
-      for (auto &v : this_indep_restart_vars) {
-        if (vName.compare(v->label()) == 0) {
-          auto v_h = v->data.GetHostMirror();
-
-          // Note index l transposed to interior
-          for (int k = out_kb.s; k <= out_kb.e; ++k) {
-            for (int j = out_jb.s; j <= out_jb.e; ++j) {
-              for (int i = out_ib.s; i <= out_ib.e; ++i) {
-                for (int l = 0; l < v_h.GetDim(4); ++l) {
-                  v_h(l, k, j, i) = tmp[index++];
-                }
-              }
+      if (v_info->IsSparse()) {
+        // check if the sparse variable is allocated on this block
+        if (sparse_info.IsAllocated(pmb->gid, sparse_idxs.at(label))) {
+          pmb->AllocateSparse(label);
+        } else {
+          // nothing to read for this block, advance reading index
+          index += nCells * vlen;
+          continue;
+        }
+      }
+
+      auto v = pmb->meshblock_data.Get()->GetCellVarPtr(label);
+      auto v_h = v->data.GetHostMirror();
+
+      // Note index l transposed to interior
+      for (int k = out_kb.s; k <= out_kb.e; ++k) {
+        for (int j = out_jb.s; j <= out_jb.e; ++j) {
+          for (int i = out_ib.s; i <= out_ib.e; ++i) {
+            for (int l = 0; l < vlen; ++l) {
+              v_h(l, k, j, i) = tmp[index++];
             }
           }
-
-          v->data.DeepCopy(v_h);
-          found = true;
-          break;
         }
       }
-      if (!found) {
-        std::stringstream msg;
-        msg << "### ERROR: Unable to find variable " << vName << std::endl;
-        PARTHENON_FAIL(msg);
-      }
+
+      v->data.DeepCopy(v_h);
     }
   }
 }
diff --git a/tst/regression/CMakeLists.txt b/tst/regression/CMakeLists.txt
index cd998757040c..fd8ca82f812c 100644
--- a/tst/regression/CMakeLists.txt
+++ b/tst/regression/CMakeLists.txt
@@ -64,13 +64,13 @@ list(APPEND EXTRA_TEST_LABELS "")
 
 if (ENABLE_HDF5)
 
-  # h5py is needed for restart and hdff5 test
+  # h5py is needed for restart and hdf5 test
   list(APPEND REQUIRED_PYTHON_MODULES h5py)
 
   # Restart
   list(APPEND TEST_DIRS restart)
   list(APPEND TEST_PROCS ${NUM_MPI_PROC_TESTING})
-  list(APPEND TEST_ARGS "--driver ${PROJECT_BINARY_DIR}/example/advection/advection-example \
+  list(APPEND TEST_ARGS "--driver ${PROJECT_BINARY_DIR}/example/sparse_advection/sparse_advection-example \
     --driver_input ${CMAKE_CURRENT_SOURCE_DIR}/test_suites/restart/parthinput.restart \
     --num_steps 2")
   list(APPEND EXTRA_TEST_LABELS "")
diff --git a/tst/regression/test_suites/advection_outflow/advection_outflow.py b/tst/regression/test_suites/advection_outflow/advection_outflow.py
index ce74cff35b4a..1e6aa34292e1 100644
--- a/tst/regression/test_suites/advection_outflow/advection_outflow.py
+++ b/tst/regression/test_suites/advection_outflow/advection_outflow.py
@@ -25,12 +25,10 @@
 
 class TestCase(utils.test_case.TestCaseAbs):
     def Prepare(self, parameters, step):
-
         parameters.coverage_status = "both"
         return parameters
 
     def Analyse(self, parameters):
-
         sys.path.insert(
             1,
             parameters.parthenon_path
diff --git a/tst/regression/test_suites/restart/parthinput.restart b/tst/regression/test_suites/restart/parthinput.restart
index d037ec9cb8fc..bf39cc03ecff 100644
--- a/tst/regression/test_suites/restart/parthinput.restart
+++ b/tst/regression/test_suites/restart/parthinput.restart
@@ -19,20 +19,20 @@ refinement = adaptive
 numlevel = 3
 
 nx1 = 64
-x1min = -0.5
-x1max = 0.5
+x1min = -1
+x1max = 1
 ix1_bc = periodic
 ox1_bc = periodic
 
 nx2 = 64
-x2min = -0.5
-x2max = 0.5
+x2min = -1
+x2max = 1
 ix2_bc = periodic
 ox2_bc = periodic
 
 nx3 = 1
-x3min = -0.5
-x3max = 0.5
+x3min = -1
+x3max = 1
 ix3_bc = periodic
 ox3_bc = periodic
 
@@ -45,7 +45,9 @@ nx3 = 1
 tlim = 0.5
 integrator = rk2
 
-<Advection>
+<SparseAdvection>
+restart_test = true
+
 cfl = 0.45
 vx = 1.0
 vy = 1.0
@@ -53,7 +55,8 @@ vz = 1.0
 profile = hard_sphere
 
 refine_tol = 0.3    # control the package specific refinement tagging function
-derefine_tol = 0.03
+derefine_tol = 0.01 # with a larger value there will be bigger diffs due to
+                    # de-refinement counters being reset to 0 at restart
 compute_error = false
 
 <parthenon/output0>
diff --git a/tst/regression/test_suites/restart/restart.py b/tst/regression/test_suites/restart/restart.py
index 1b2874bbdeaf..8b357799e67d 100644
--- a/tst/regression/test_suites/restart/restart.py
+++ b/tst/regression/test_suites/restart/restart.py
@@ -16,11 +16,7 @@
 # ========================================================================================
 
 # Modules
-import h5py
-import math
-import numpy as np
 import sys
-import os
 import utils.test_case
 
 
@@ -30,7 +26,6 @@
 
 class TestCase(utils.test_case.TestCaseAbs):
     def Prepare(self, parameters, step):
-
         # enable coverage testing on pass where restart
         # files are both read and written
         parameters.coverage_status = "both"
@@ -47,16 +42,23 @@ def Prepare(self, parameters, step):
         return parameters
 
     def Analyse(self, parameters):
-        # spotcheck one variable
-        goldFile = "gold.out0.00002.rhdf"
-        silverFile = "silver.out0.00002.rhdf"
-        varName = "advected"
-        with h5py.File(goldFile, "r") as gold, h5py.File(silverFile, "r") as silver:
-            goldData = np.zeros(gold[varName].shape, dtype=np.float64)
-            gold[varName].read_direct(goldData)
-            silverData = np.zeros(silver[varName].shape, dtype=np.float64)
-            silver[varName].read_direct(silverData)
-
-        maxdiff = np.abs(goldData - silverData).max()
-        print("Variable: %s, diff=%g, N=%d" % (varName, maxdiff, len(goldData)))
-        return maxdiff == 0.0
+        sys.path.insert(
+            1,
+            parameters.parthenon_path
+            + "/scripts/python/packages/parthenon_tools/parthenon_tools",
+        )
+
+        try:
+            from phdf_diff import compare
+        except ModuleNotFoundError:
+            print("Couldn't find module to compare Parthenon hdf5 files.")
+            return False
+
+        delta = compare(
+            [
+                "gold.out0.00002.rhdf",
+                "silver.out0.00002.rhdf",
+            ]
+        )
+
+        return delta == 0