diff --git a/doc/sphinx/src/sparse_packs.rst b/doc/sphinx/src/sparse_packs.rst
index 026e437b60866..43a0f4470405e 100644
--- a/doc/sphinx/src/sparse_packs.rst
+++ b/doc/sphinx/src/sparse_packs.rst
@@ -85,6 +85,83 @@ snippet.
    care should be taken not to assume that scratch variables can persist between tasks, even when they are directly
    dependent on each other. Issues resolved by setting ``-DPARTHENON_DEBUG_SCRATCH=ON`` can be an indication of this issue.
 
+Slicing into ``SparsePack``\ s with ``SubPack``\ s
+--------------------------------------------------
+
+A `SubPack` provdies a view into a slice of a `SparsePack` along a given dimension(s).
+`SubPack`\ s are constructed with a block + `kji` index to give a slice into the 
+fields at the meshblock + cell index. When a `SubPack` is constructed with the
+`Axis` template parameters then the `SubPack` also providies a view into slices
+of the `SparsePack` along the provided axes offset from the provided `ijk` indices.
+
+.. code:: c++
+
+  const int ni = ib.e - ib.s + 1;
+  const int ic = ib.s + ni / 2;
+  par_for(
+      PARTHENON_AUTO_LABEL, 0, sparse_pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+      KOKKOS_LAMBDA(int b, int k, int j) {
+        int ltot = 0;
+        int lo = sparse_pack.GetLowerBound(b, v3());
+        int hi = sparse_pack.GetUpperBound(b, v3());
+        auto sub_pack = parthenon::SubPack<Axis::I>(sparse_pack, b, k, j, ic);
+
+        for (int i = ib.s - ni / 2; i <= ib.e - ni / 2; i++) {
+          for (int c = 0; c <= hi - lo; ++c) {
+            Real n = i + ic + 1e1 * j + 1e2 * k + 1e4 * c + 1e5 * v + 1e3 * b;
+            // indexes into sparse_pack(b, v3(c), k, j, ic + i)
+            if (n != sub_pack(v3(c), i)) ltot += 1;
+          }
+        }
+      });
+
+Using Virtual Fields
+--------------------
+
+Type-based ``SparsePack``\ s, when indexed with special ``virtual_variable_t`` types will return a value that is
+calculated from other variables in the pack from the ``evaluate`` function registered to the type. Virtual variable types
+need to export the other fields that their ``evaluate`` method depends on, allowing the virtual variable type, when included
+in the ``PackDescriptor``, to pack the real variables. 
+
+At the simplest a virtual variable's ``evaluate`` method has its signature take in a ``SparsePack<Ts...>`` as well as the ``i,j,k``
+indices and returns a ``Real`` as the value for the virtual field. 
+
+Additionally a virtual variable can declare a ``pack_type`` to refer to a 0D, 1D, 2D, or 3D ``SubPack``, which case the ``evaluate``
+signature only takes in the desired ``SubPack`` centered at teh ``i,j,k`` index.
+
+.. code:: c++
+
+   struct d1 : parthenon::variable_names::virtual_variable_t<v1, v3> {
+     KOKKOS_INLINE_FUNCTION d1() : x(0.0) {}
+
+     KOKKOS_INLINE_FUNCTION d1(const Real &xx) : x(xx) {}
+
+     template <typename Pack_t>
+     KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack, const int b, const int k,
+                                          const int j, const int i) const {
+       return pack(b, v1(), k, j, i) * pack(b, v3(1), k, j, i) * pack(b, v3(2), k, j, i) + x;
+     }
+
+     Real x;
+   };
+
+   // use a subpack to index into the sparse pack
+   struct d1_subpack : public parthenon::variable_names::virtual_variable_t<v1, v3> {
+     // declare the type of subpack we want to use for our evaluate method
+     using pack_type = parthenon::SubPack0D;
+
+     KOKKOS_INLINE_FUNCTION d1_subpack() : x(0.0) {}
+
+     KOKKOS_INLINE_FUNCTION d1_subpack(const Real &xx) : x(xx) {}
+
+     template <typename Pack_t>
+     KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack) const {
+       return pack(v1()) * pack(v3(1)) * pack(v3(2)) + x;
+     }
+
+     Real x;
+   };
+
 
 Building and Using a ``SparsePack``
 -----------------------------------
diff --git a/src/basic_types.hpp b/src/basic_types.hpp
index e55b6d1c6edc1..c039dcc235b7e 100644
--- a/src/basic_types.hpp
+++ b/src/basic_types.hpp
@@ -177,6 +177,9 @@ inline bool operator<(const GridIdentifier &lhs, const GridIdentifier &rhs) {
   return lhs.logical_level() < rhs.logical_level();
 }
 
+// Enumeration for accessing spatial axes of a meshblock
+enum class Axis { I = 0, J = 1, K = 2 };
+
 // Enumeration for accessing a field on different locations of the grid:
 // CC = cell center of (i, j, k)
 // F1 = x-normal face at (i - 1/2, j, k)
diff --git a/src/pack/pack_utils.hpp b/src/pack/pack_utils.hpp
index 0c1087c7c0a95..21cd2603a6b7d 100644
--- a/src/pack/pack_utils.hpp
+++ b/src/pack/pack_utils.hpp
@@ -22,7 +22,12 @@
 #include <utility>
 #include <vector>
 
+#include <Kokkos_Core.hpp>
+
 #include "basic_types.hpp"
+#include "utils/concepts_lite.hpp"
+#include "utils/error_checking.hpp"
+#include "utils/type_list.hpp"
 
 // SFINAE for block iter so that Sparse/SwarmPacks can work for MeshBlockData and MeshData
 namespace {
@@ -161,6 +166,65 @@ struct any_nonautoflux : public base_t<true> {
   }
 };
 using any = any_nonautoflux;
+
+// Concept to state that a typed-field is dependent on other
+// fields, and is not itself an actual indexable field.
+template <typename T>
+concept DependentVariable =
+    requires { typename T::independent_vars; }; // NOLINT(readability/braces)
+
+template <typename... Ts>
+struct virtual_variable_t {
+  using type = virtual_variable_t<Ts...>;
+  using independent_vars = TypeList<Ts...>;
+};
+
+// Concept to check that a type shares an ancestor with virtual_variable_t
+template <typename T>
+concept VirtualVariable = requires {
+  typename T::type;
+  requires is_specialization_of<typename T::type, virtual_variable_t>::value;
+}; // NOLINT(readability/braces)
+
+// Concept to check that a type wants a subpack to evaluate the virtual field
+template <typename T>
+concept VirtualSubPack = requires {
+  typename T::pack_type;
+  requires std::is_same_v<decltype(T::pack_type::Naxes), const int>;
+}; // NOLINT(readability/braces)
+
+namespace impl {
+
+struct AllIndependentVariables {
+  template <typename T>
+    requires(!DependentVariable<T>)
+  static auto get(T) {
+    return TypeList<T>();
+  }
+
+  template <template <typename...> typename TL, typename... Ts>
+  static auto get(TL<Ts...>) {
+    return AllIndependentVariables::get(Ts()...);
+  }
+
+  template <typename T>
+    requires(DependentVariable<T>)
+  static auto get(T) {
+    return AllIndependentVariables::get(typename T::independent_vars());
+  }
+
+  template <typename T, typename... Ts>
+  static auto get(T, Ts...) {
+    return union_type_lists_t<decltype(AllIndependentVariables::get(T())),
+                              decltype(AllIndependentVariables::get(Ts()...))>();
+  }
+};
+} // namespace impl
+
+// Get a TypeList of all the independent variables that make up the requested types.
+template <typename... Ts>
+using all_independent_variables_t = decltype(impl::AllIndependentVariables::get(Ts()...));
+
 } // namespace variable_names
 
 // Namespace in which to put swarm variable name types that are used for indexing into
diff --git a/src/pack/sparse_pack/make_pack_descriptor.hpp b/src/pack/sparse_pack/make_pack_descriptor.hpp
index e63e2af96e437..46ab039ea59b8 100644
--- a/src/pack/sparse_pack/make_pack_descriptor.hpp
+++ b/src/pack/sparse_pack/make_pack_descriptor.hpp
@@ -24,8 +24,10 @@
 #include "interface/metadata.hpp"
 #include "interface/state_descriptor.hpp"
 #include "mesh/mesh.hpp"
+#include "pack/pack_utils.hpp"
 #include "pack/sparse_pack/pack_descriptor.hpp"
 #include "pack/sparse_pack/sparse_pack.hpp"
+#include "utils/concepts_lite.hpp"
 #include "utils/type_list.hpp"
 
 namespace parthenon {
@@ -63,6 +65,7 @@ inline auto MakePackDescriptor(MT *pmd, const std::vector<std::string> &vars,
 }
 
 template <class... Ts, class MT>
+  requires(!variable_names::DependentVariable<Ts> && ...)
 inline auto MakePackDescriptor(MT *pmd, const std::vector<MetadataFlag> &flags = {},
                                const std::set<PDOpt> &options = {}) {
   const std::vector<std::string> vars{Ts::name()...};
@@ -72,6 +75,7 @@ inline auto MakePackDescriptor(MT *pmd, const std::vector<MetadataFlag> &flags =
 }
 
 template <class... Ts, class MT>
+  requires(!variable_names::DependentVariable<Ts> && ...)
 inline auto MakePackDescriptor(SparsePack<Ts...> pack, MT *pmd,
                                const std::vector<MetadataFlag> &flags = {},
                                const std::set<PDOpt> &options = {}) {
@@ -110,6 +114,14 @@ inline auto MakePackDescriptorFromTypeList(Args &&...args) {
   return MakePackDescriptorFromTypeList(TL(), std::forward<Args>(args)...);
 }
 
+template <class... Ts, class MT>
+  requires(variable_names::DependentVariable<Ts> || ...)
+inline auto MakePackDescriptor(MT *pmd, const std::vector<MetadataFlag> &flags = {},
+                               const std::set<PDOpt> &options = {}) {
+  return MakePackDescriptorFromTypeList(
+      variable_names::all_independent_variables_t<Ts...>(), pmd, flags, options);
+}
+
 struct PackDescriptorCacheBase {
   virtual ~PackDescriptorCacheBase() = default;
 };
diff --git a/src/pack/sparse_pack/sparse_pack.hpp b/src/pack/sparse_pack/sparse_pack.hpp
index 6dd5362c8ca6c..d61a2394f0929 100644
--- a/src/pack/sparse_pack/sparse_pack.hpp
+++ b/src/pack/sparse_pack/sparse_pack.hpp
@@ -21,6 +21,7 @@
 #include "pack/block_selector.hpp"
 #include "pack/pack_utils.hpp"
 #include "pack/sparse_pack/sparse_pack_base.hpp"
+#include "pack/subpack.hpp"
 #include "utils/concepts_lite.hpp"
 #include "utils/type_list.hpp"
 
@@ -340,6 +341,51 @@ class SparsePack : public SparsePackBase {
     return pack_(static_cast<int>(el) % 3, b, vidx)(k, j, i);
   }
 
+  template <typename Tin, typename... Args>
+    requires(variable_names::VirtualVariable<Tin> && !variable_names::VirtualSubPack<Tin>)
+  KOKKOS_INLINE_FUNCTION Real operator()(const int b, const Tin &t, const int k,
+                                         const int j, const int i, Args &&...args) const {
+    static_assert(
+        TypeList<Ts...>::IsIn(typename Tin::independent_vars()),
+        "SparsePack must pack all independent_vars needed for virtual variable");
+    return t.evaluate(*this, b, k, j, i, std::forward<Args>(args)...);
+  }
+
+  template <typename Tin, typename... Args>
+    requires(variable_names::VirtualVariable<Tin> && variable_names::VirtualSubPack<Tin>)
+  KOKKOS_INLINE_FUNCTION Real operator()(const int b, const Tin &t, const int k,
+                                         const int j, const int i, Args &&...args) const {
+    static_assert(
+        TypeList<Ts...>::IsIn(typename Tin::independent_vars()),
+        "SparsePack must pack all independent_vars needed for virtual variable");
+    using sub_pack_type = decltype(SubPack<typename Tin::pack_type>(*this, b, k, j, i));
+    return t.evaluate(SubPack<typename Tin::pack_type>(*this, b, k, j, i),
+                      std::forward<Args>(args)...);
+  }
+  template <typename Tin, typename... Args>
+    requires(variable_names::VirtualVariable<Tin> && variable_names::VirtualSubPack<Tin>)
+  KOKKOS_INLINE_FUNCTION Real operator()(const int b, const TE el, const Tin &t,
+                                         const int k, const int j, const int i,
+                                         Args &&...args) const {
+    static_assert(
+        TypeList<Ts...>::IsIn(typename Tin::independent_vars()),
+        "SparsePack must pack all independent_vars needed for virtual variable");
+    using sub_pack_type = decltype(SubPack<typename Tin::pack_type>(*this, b, k, j, i));
+    return t.evaluate(SubPack<typename Tin::pack_type>(*this, b, k, j, i), el,
+                      std::forward<Args>(args)...);
+  }
+
+  template <typename Tin, typename... Args>
+    requires(variable_names::VirtualVariable<Tin> && !variable_names::VirtualSubPack<Tin>)
+  KOKKOS_INLINE_FUNCTION Real operator()(const int b, const TE el, const Tin &t,
+                                         const int k, const int j, const int i,
+                                         Args &&...args) const {
+    static_assert(
+        TypeList<Ts...>::IsIn(typename Tin::independent_vars()),
+        "SparsePack must pack all independent_vars needed for virtual variable");
+    return t.evaluate(*this, b, el, k, j, i, std::forward<Args>(args)...);
+  }
+
   // flux() overloads
   template <class... Args>
   KOKKOS_INLINE_FUNCTION auto &flux(const int b, const TopologicalElement te,
diff --git a/src/pack/subpack.hpp b/src/pack/subpack.hpp
new file mode 100644
index 0000000000000..199183b2cb7f9
--- /dev/null
+++ b/src/pack/subpack.hpp
@@ -0,0 +1,224 @@
+//========================================================================================
+// (C) (or copyright) 2020-2025. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+#ifndef PACK_SUBPACK_HPP_
+#define PACK_SUBPACK_HPP_
+
+#include "basic_types.hpp"
+#include "utils/concepts_lite.hpp"
+#include "utils/type_list.hpp"
+
+namespace parthenon {
+
+namespace impl {
+
+template <typename PackType, Axis... axes>
+struct SubPack_impl {
+  KOKKOS_INLINE_FUNCTION SubPack_impl(PackType &pack, const int &b, const int &k,
+                                      const int &j, const int &i)
+      : pack_(pack), b_(b), k_(k), j_(j), i_(i) {}
+
+  template <typename Var_t>
+  KOKKOS_INLINE_FUNCTION Real &operator()(const Var_t &var) const {
+    return pack_(b_, var, k_, j_, i_);
+  }
+
+  template <typename Var_t>
+  KOKKOS_INLINE_FUNCTION Real &operator()(TopologicalElement te, const Var_t &var) const {
+    return pack_(b_, te, var, k_, j_, i_);
+  }
+
+  template <typename Var_t>
+  KOKKOS_INLINE_FUNCTION Real &flux(TopologicalElement te, const Var_t &var) const {
+    return pack_.flux(b_, te, var, k_, j_, i_);
+  }
+
+  template <typename V>
+  KOKKOS_INLINE_FUNCTION std::size_t GetSize(const V &var) const {
+    return pack_.GetSize(b_, var);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  const auto &GetCoordinates() const { return pack_.GetCoordinates(b_); }
+
+ private:
+  PackType &pack_;
+  const int b_, k_, j_, i_;
+};
+
+template <typename PackType, Axis... axes>
+struct StencilSubPack_impl {
+  KOKKOS_INLINE_FUNCTION StencilSubPack_impl(PackType &pack, const int &b, const int &k,
+                                             const int &j, const int &i)
+      : pack_(pack), b_(b), ijk_({i, j, k}) {}
+
+  template <typename Var_t, typename... Is>
+  KOKKOS_INLINE_FUNCTION Real &operator()(const Var_t &var, Is &&...idxs) const {
+    static_assert(sizeof...(Is) == sizeof...(axes),
+                  "number of indices passed to sub pack must match number of axes.");
+    Kokkos::Array<int, 3> ijk = ijk_;
+    ([&]() { ijk[static_cast<int>(axes)] += idxs; }(), ...);
+    return pack_(b_, var, ijk[2], ijk[1], ijk[0]);
+  }
+
+  template <typename Var_t, typename... Is>
+  KOKKOS_INLINE_FUNCTION Real &operator()(TopologicalElement te, const Var_t &var,
+                                          Is &&...idxs) const {
+    static_assert(sizeof...(Is) == sizeof...(axes),
+                  "number of indices passed to sub pack must match number of axes.");
+    Kokkos::Array<int, 3> ijk = ijk_;
+    ([&]() { ijk[static_cast<int>(axes)] += idxs; }(), ...);
+    return pack_(b_, te, var, ijk[2], ijk[1], ijk[0]);
+  }
+
+  template <typename Var_t, typename... Is>
+  KOKKOS_INLINE_FUNCTION Real &flux(TopologicalElement te, const Var_t &var,
+                                    Is &&...idxs) const {
+    static_assert(sizeof...(Is) == sizeof...(axes),
+                  "number of indices passed to sub pack must match number of axes.");
+    Kokkos::Array<int, 3> ijk = ijk_;
+    ([&]() { ijk[static_cast<int>(axes)] += idxs; }(), ...);
+    return pack_.flux(b_, te, var, ijk[2], ijk[1], ijk[0]);
+  }
+
+  template <typename V>
+  KOKKOS_INLINE_FUNCTION std::size_t GetSize(const V &var) const {
+    return pack_.GetSize(b_, var);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  const auto &GetCoordinates() const { return pack_.GetCoordinates(b_); }
+
+ private:
+  const PackType &pack_;
+  const Kokkos::Array<int, 3> ijk_;
+  const int b_;
+};
+
+template <typename Var_t, typename PackType, Axis... axes>
+struct VarStencilSubPack_impl {
+  KOKKOS_INLINE_FUNCTION VarStencilSubPack_impl(PackType &pack, const int &b,
+                                                const Var_t &var, const int &k,
+                                                const int &j, const int &i)
+      : pack_(pack), b_(b), var_(var), ijk_({i, j, k}) {}
+
+  template <typename... Is>
+  KOKKOS_INLINE_FUNCTION Real &operator()(Is &&...idxs) const {
+    static_assert(sizeof...(Is) == sizeof...(axes),
+                  "number of indices passed to sub pack must match number of axes.");
+    Kokkos::Array<int, 3> ijk = ijk_;
+    ([&]() { ijk[static_cast<int>(axes)] += idxs; }(), ...);
+    return pack_(b_, var_, ijk[2], ijk[1], ijk[0]);
+  }
+
+  template <typename... Is>
+  KOKKOS_INLINE_FUNCTION Real &operator()(TopologicalElement te, Is &&...idxs) const {
+    static_assert(sizeof...(Is) == sizeof...(axes),
+                  "number of indices passed to sub pack must match number of axes.");
+    Kokkos::Array<int, 3> ijk = ijk_;
+    ([&]() { ijk[static_cast<int>(axes)] += idxs; }(), ...);
+    return pack_(b_, te, var_, ijk[2], ijk[1], ijk[0]);
+  }
+
+  template <typename... Is>
+  KOKKOS_INLINE_FUNCTION Real &flux(TopologicalElement te, Is &&...idxs) const {
+    static_assert(sizeof...(Is) == sizeof...(axes),
+                  "number of indices passed to sub pack must match number of axes.");
+    Kokkos::Array<int, 3> ijk = ijk_;
+    ([&]() { ijk[static_cast<int>(axes)] += idxs; }(), ...);
+    return pack_.flux(b_, te, var_, ijk[2], ijk[1], ijk[0]);
+  }
+
+  KOKKOS_INLINE_FUNCTION
+  const auto &GetCoordinates() const { return pack_.GetCoordinates(b_); }
+
+ private:
+  const PackType &pack_;
+  const Kokkos::Array<int, 3> ijk_;
+  const Var_t var_;
+  const int b_;
+};
+} // namespace impl
+
+template <typename... Ts>
+class SparsePack;
+
+template <typename PackType>
+constexpr bool is_sparse_pack =
+    is_specialization_of<base_type<PackType>, SparsePack>::value;
+
+// Some convenience types to differentiate between the types of subpacks
+struct SubPack0D {
+  static constexpr int Naxes = 0;
+};
+
+template <Axis axis1>
+struct SubPack1D {
+  static constexpr int Naxes = 1;
+  static constexpr Axis axis = axis1;
+};
+
+template <Axis a1, Axis a2>
+struct SubPack2D {
+  static constexpr int Naxes = 2;
+  static constexpr Axis axis1 = a1;
+  static constexpr Axis axis2 = a2;
+};
+
+template <Axis a1, Axis a2, Axis a3>
+struct SubPack3D {
+  static constexpr int Naxes = 3;
+  static constexpr Axis axis1 = a1;
+  static constexpr Axis axis2 = a2;
+  static constexpr Axis axis3 = a3;
+};
+
+template <Axis axis, Axis... axes, typename Var_t, typename SparsePackType>
+  requires(is_sparse_pack<SparsePackType>)
+KOKKOS_INLINE_FUNCTION auto SubPack(SparsePackType &pack, const int &b, const Var_t &var,
+                                    const int &k, const int &j, const int &i) {
+  return VarStencilSubPack_impl<Var_t, SparsePackType, axis, axes...>(pack, b, var, k, j,
+                                                                      i);
+}
+
+template <Axis axis, Axis... axes, typename SparsePackType>
+  requires(is_sparse_pack<SparsePackType>)
+KOKKOS_INLINE_FUNCTION auto SubPack(SparsePackType &pack, const int &b, const int &k,
+                                    const int &j, const int &i) {
+  return StencilSubPack_impl<SparsePackType, axis, axes...>(pack, b, k, j, i);
+}
+
+template <typename SparsePackType>
+  requires(is_sparse_pack<SparsePackType>)
+KOKKOS_INLINE_FUNCTION auto SubPack(SparsePackType &pack, const int &b, const int &k,
+                                    const int &j, const int &i) {
+  return SubPack_impl<SparsePackType>(pack, b, k, j, i);
+}
+
+template <typename PackType, typename SparsePackType>
+  requires(is_sparse_pack<SparsePackType>)
+KOKKOS_INLINE_FUNCTION auto SubPack(SparsePackType &pack, const int &b, const int &k,
+                                    const int &j, const int &i) {
+  constexpr int Naxes = PackType::Naxes;
+  if constexpr (Naxes == 3) {
+    return SubPack<PackType::axis1, PackType::axis2, PackType::axis3>(pack, b, k, j, i);
+  } else if constexpr (Naxes == 2) {
+    return SubPack<PackType::axis1, PackType::axis2>(pack, b, k, j, i);
+  } else if constexpr (Naxes == 1) {
+    return SubPack<PackType::axis>(pack, b, k, j, i);
+  } else {
+    return SubPack(pack, b, k, j, i);
+  }
+}
+
+} // namespace parthenon
+#endif // PACK_SUBPACK_HPP_
diff --git a/src/utils/type_list.hpp b/src/utils/type_list.hpp
index 88e1d6b1ecedd..cead515f9a6d0 100644
--- a/src/utils/type_list.hpp
+++ b/src/utils/type_list.hpp
@@ -67,6 +67,11 @@ struct TypeList {
     }
   }
 
+  template <template <typename...> typename TL, typename... Ts>
+  static constexpr bool IsIn(TL<Ts...>) {
+    return (IsIn<Ts>() && ...);
+  }
+
   template <class F>
   static void IterateTypes(F func) {
     (func(Args()), ...);
@@ -96,6 +101,28 @@ auto ConcatenateTypeLists(TypeList<Args1...>, TypeList<Args2...>, Args...) {
   return ConcatenateTypeLists(TypeList<Args1..., Args2...>(), Args()...);
 }
 
+template <typename, typename, typename...>
+struct UnionTypeLists {};
+
+template <typename... Ts>
+struct UnionTypeLists<TypeList<Ts...>, TypeList<>> {
+  using type = TypeList<Ts...>;
+};
+
+template <typename... Ts, typename V, typename... Vs>
+struct UnionTypeLists<TypeList<Ts...>, TypeList<V, Vs...>> {
+  using TL = TypeList<Ts...>;
+  using type = typename UnionTypeLists<
+      std::conditional_t<TL::template IsIn<V>(), TL, TypeList<Ts..., V>>,
+      TypeList<Vs...>>::type;
+};
+
+template <typename T, typename U, typename V, typename... Args>
+struct UnionTypeLists<T, U, V, Args...> {
+  using type =
+      typename UnionTypeLists<typename UnionTypeLists<T, U>::type, V, Args...>::type;
+};
+
 template <class T, std::size_t I, class... Ts>
 static auto InsertTypeImpl(TypeList<Ts...>) {
   if constexpr (I == 0) {
@@ -115,6 +142,9 @@ template <class... TLs>
 using concatenate_type_lists_t =
     decltype(impl::ConcatenateTypeLists(std::declval<TLs>()...));
 
+template <typename... TLs>
+using union_type_lists_t = typename impl::UnionTypeLists<TLs...>::type;
+
 template <class T, class TL, std::size_t I = TL::n_types>
 using insert_type_list_t = decltype(impl::InsertTypeImpl<T, I>(TL()));
 
diff --git a/tst/unit/test_sparse_pack.cpp b/tst/unit/test_sparse_pack.cpp
index b2c5d3cb5a92a..09d75d6b5af2b 100644
--- a/tst/unit/test_sparse_pack.cpp
+++ b/tst/unit/test_sparse_pack.cpp
@@ -1,5 +1,5 @@
-//========================================================================================
-// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
+//=======================================================================================
+//(C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -27,11 +27,15 @@
 #include "kokkos_abstraction.hpp"
 #include "loop_abstraction/loop_abstraction.hpp"
 #include "mesh/meshblock.hpp"
+#include "pack/pack_utils.hpp"
 #include "pack/sparse_pack/make_pack_descriptor.hpp"
 #include "pack/sparse_pack/sparse_pack.hpp"
 
 // TODO(jcd): can't call the MeshBlock constructor without mesh_refinement.hpp???
 #include "mesh/mesh_refinement.hpp"
+#include "pack/subpack.hpp"
+#include "utils/error_checking.hpp"
+#include "utils/type_list.hpp"
 
 using parthenon::BlockList_t;
 using parthenon::DevExecSpace;
@@ -103,6 +107,102 @@ struct v7 : public parthenon::variable_names::base_t<false, ANYDIM, 3> {
   static constexpr bool is_sparse() { return false; }
 };
 
+struct d1 : parthenon::variable_names::virtual_variable_t<v1, v3> {
+  KOKKOS_INLINE_FUNCTION d1() : x(0.0) {}
+
+  KOKKOS_INLINE_FUNCTION d1(const Real &xx) : x(xx) {}
+
+  template <typename Pack_t>
+  KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack, const int b, const int k,
+                                       const int j, const int i) const {
+    return pack(b, v1(), k, j, i) * pack(b, v3(1), k, j, i) * pack(b, v3(2), k, j, i) + x;
+  }
+
+  Real x;
+};
+
+// use a subpack to index into the sparse pack
+struct d1_subpack : public parthenon::variable_names::virtual_variable_t<v1, v3> {
+  // declare the type of subpack we want to use for our evaluate method
+  using pack_type = parthenon::SubPack0D;
+
+  KOKKOS_INLINE_FUNCTION d1_subpack() {}
+
+  template <typename Pack_t>
+  KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack, const Real &x) const {
+    return pack(v1()) * pack(v3(1)) * pack(v3(2)) + x;
+  }
+};
+
+// use 1D subpack to get the gradient of v1()
+template <parthenon::Axis axis>
+struct gradient : public parthenon::variable_names::virtual_variable_t<v1> {
+  using pack_type = parthenon::SubPack1D<axis>;
+
+  template <typename Pack_t>
+  KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack) const {
+    return 0.5 * (pack(v1(), 1) - pack(v1(), -1));
+  }
+};
+
+// gradient normal to face from te
+struct gradient_te : public parthenon::variable_names::virtual_variable_t<v1> {
+  template <typename Pack_t>
+  KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack, const int b,
+                                       const parthenon::TopologicalElement el,
+                                       const int k, const int j, const int i) const {
+    using TE = parthenon::TopologicalElement;
+    switch (el) {
+    case (TE::F1):
+      return (pack(b, v1(), k, j, i) - pack(b, v1(), k, j, i - 1));
+    case (TE::F2):
+      return (pack(b, v1(), k, j, i) - pack(b, v1(), k, j - 1, i));
+    case (TE::F3):
+      return (pack(b, v1(), k, j, i) - pack(b, v1(), k - 1, j, i));
+    default:
+      return 0.;
+    }
+  }
+};
+
+// use 2d subpack to get the curl of a vector
+template <parthenon::Axis axis>
+struct curl : public parthenon::variable_names::virtual_variable_t<v3> {
+  static constexpr parthenon::Axis axis2 =
+      static_cast<parthenon::Axis>((static_cast<int>(axis) + 1) % 3);
+  static constexpr parthenon::Axis axis3 =
+      static_cast<parthenon::Axis>((static_cast<int>(axis) + 2) % 3);
+  using pack_type = parthenon::SubPack2D<axis2, axis3>;
+
+  template <typename Pack_t>
+  KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack) const {
+    constexpr int ax2 = static_cast<int>(axis2);
+    constexpr int ax3 = static_cast<int>(axis3);
+    return 0.5 * ((pack(v3(ax3), 1, 0) - pack(v3(ax3), -1, 0)) -
+                  (pack(v3(ax2), 0, 1) - pack(v3(ax2), 0, -1)));
+  }
+};
+
+// use 3d subpack to get the third derivative
+// d_1 d_2 d_3 v1
+struct der3 : public parthenon::variable_names::virtual_variable_t<v1> {
+  using pack_type =
+      parthenon::SubPack3D<parthenon::Axis::I, parthenon::Axis::J, parthenon::Axis::K>;
+
+  template <typename Pack_t>
+  KOKKOS_INLINE_FUNCTION Real evaluate(const Pack_t &pack) const {
+    Real d3 = 0.;
+    const parthenon::IndexRange pm{0, 1};
+    parthenon::seq_for(pm, pm, pm, [&](const int k, const int j, const int i) {
+      const int kk = 2 * k - 1;
+      const int jj = 2 * j - 1;
+      const int ii = 2 * i - 1;
+      d3 += kk * jj * ii * pack(v1(), kk, jj, ii);
+    });
+    return d3 * 0.125;
+  }
+};
+
 using namespace parthenon::loop_abstraction;
 struct PackViewSpec {
   int nblocks;
@@ -488,6 +588,125 @@ TEST_CASE("Test behavior of sparse packs", "[SparsePack]") {
         REQUIRE(nwrong == 0);
       }
 
+      THEN("A sub pack correctly loads this data and can be read from v3 on all "
+           "blocks") {
+        // Create a pack use type variables
+        auto desc =
+            parthenon::MakePackDescriptor<v5, v3>(pkg.get(), {Metadata::WithFluxes});
+        auto sparse_pack = desc.GetPack(&mesh_data);
+
+        const int v = 1; // v3 is the second variable in the loop above so v = 1 there
+        int nwrong = 0;
+        par_reduce(
+            loop_pattern_mdrange_tag, "check vector", DevExecSpace(), 0,
+            sparse_pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e, ib.s, ib.e,
+            KOKKOS_LAMBDA(int b, int k, int j, int i, int &ltot) {
+              int lo = sparse_pack.GetLowerBound(b, v3());
+              int hi = sparse_pack.GetUpperBound(b, v3());
+              for (int c = 0; c <= hi - lo; ++c) {
+                Real n = i + 1e1 * j + 1e2 * k + 1e4 * c + 1e5 * v + 1e3 * b;
+                auto sub_pack = parthenon::SubPack(sparse_pack, b, k, j, i);
+                if (n != sub_pack(v3(c))) ltot += 1;
+              }
+            },
+            nwrong);
+        REQUIRE(nwrong == 0);
+
+        using Axis = parthenon::Axis;
+        AND_THEN("1D Stencil subpacks can correctly access the data") {
+          const int ni = ib.e - ib.s + 1;
+          const int ic = ib.s + ni / 2;
+          par_reduce(
+              loop_pattern_mdrange_tag, "check vector", DevExecSpace(), 0,
+              sparse_pack.GetNBlocks() - 1, kb.s, kb.e, jb.s, jb.e,
+              KOKKOS_LAMBDA(int b, int k, int j, int &ltot) {
+                int lo = sparse_pack.GetLowerBound(b, v3());
+                int hi = sparse_pack.GetUpperBound(b, v3());
+                int vc = (hi - lo) / 2;
+                auto sub_pack_v =
+                    parthenon::SubPack<Axis::I>(sparse_pack, b, v3(vc), k, j, ic);
+                auto sub_pack = parthenon::SubPack<Axis::I>(sparse_pack, b, k, j, ic);
+
+                for (int i = ib.s - ni / 2; i <= ib.e - ni / 2; i++) {
+                  for (int c = 0; c <= hi - lo; ++c) {
+                    Real n = i + ic + 1e1 * j + 1e2 * k + 1e4 * c + 1e5 * v + 1e3 * b;
+                    if (n != sub_pack(v3(c), i)) ltot += 1;
+                    if (n != sub_pack_v(i) && c == vc) ltot += 1;
+                  }
+                }
+              },
+              nwrong);
+          REQUIRE(nwrong == 0);
+        }
+
+        AND_THEN("2D Stencil subpacks can correctly access the data") {
+          const int ni = ib.e - ib.s + 1;
+          const int ic = ib.s + ni / 2;
+          const int nj = jb.e - jb.s + 1;
+          const int jc = jb.s + nj / 2;
+          par_reduce(
+              loop_pattern_mdrange_tag, "check vector", DevExecSpace(), 0,
+              sparse_pack.GetNBlocks() - 1, kb.s, kb.e,
+              KOKKOS_LAMBDA(int b, int k, int &ltot) {
+                int lo = sparse_pack.GetLowerBound(b, v3());
+                int hi = sparse_pack.GetUpperBound(b, v3());
+                int vc = (hi - lo) / 2;
+                auto sub_pack_v = parthenon::SubPack<Axis::I, Axis::J>(sparse_pack, b,
+                                                                       v3(vc), k, jc, ic);
+                auto sub_pack =
+                    parthenon::SubPack<Axis::I, Axis::J>(sparse_pack, b, k, jc, ic);
+
+                for (int j = jb.s - nj / 2; j <= jb.e - nj / 2; j++) {
+                  for (int i = ib.s - ni / 2; i <= ib.e - ni / 2; i++) {
+                    for (int c = 0; c <= hi - lo; ++c) {
+                      Real n =
+                          i + ic + 1e1 * (j + jc) + 1e2 * k + 1e4 * c + 1e5 * v + 1e3 * b;
+                      if (n != sub_pack(v3(c), i, j)) ltot += 1;
+                      if (n != sub_pack_v(i, j) && c == vc) ltot += 1;
+                    }
+                  }
+                }
+              },
+              nwrong);
+          REQUIRE(nwrong == 0);
+        }
+
+        AND_THEN("3D Stencil subpacks can correctly access the data") {
+          const int ni = ib.e - ib.s + 1;
+          const int ic = ib.s + ni / 2;
+          const int nj = jb.e - jb.s + 1;
+          const int jc = jb.s + nj / 2;
+          const int nk = kb.e - kb.s + 1;
+          const int kc = kb.s + nk / 2;
+          parthenon::par_reduce(
+              "check vector", 0, sparse_pack.GetNBlocks() - 1,
+              KOKKOS_LAMBDA(int b, int &ltot) {
+                int lo = sparse_pack.GetLowerBound(b, v3());
+                int hi = sparse_pack.GetUpperBound(b, v3());
+                int vc = (hi - lo) / 2;
+                auto sub_pack_v = parthenon::SubPack<Axis::I, Axis::J, Axis::K>(
+                    sparse_pack, b, v3(vc), kc, jc, ic);
+                auto sub_pack = parthenon::SubPack<Axis::I, Axis::J, Axis::K>(
+                    sparse_pack, b, kc, jc, ic);
+
+                for (int k = kb.s - nk / 2; k <= kb.e - nk / 2; k++) {
+                  for (int j = jb.s - nj / 2; j <= jb.e - nj / 2; j++) {
+                    for (int i = ib.s - ni / 2; i <= ib.e - ni / 2; i++) {
+                      for (int c = 0; c <= hi - lo; ++c) {
+                        Real n = i + ic + 1e1 * (j + jc) + 1e2 * (k + kc) + 1e4 * c +
+                                 1e5 * v + 1e3 * b;
+                        if (n != sub_pack(v3(c), i, j, k)) ltot += 1;
+                        if (n != sub_pack_v(i, j, k) && c == vc) ltot += 1;
+                      }
+                    }
+                  }
+                }
+              },
+              nwrong);
+          REQUIRE(nwrong == 0);
+        }
+      }
+
       THEN("A bovi sparse pack view works through the loop abstraction on memory spans") {
         using namespace parthenon::loop_abstraction;
         using IS = IndexSpace<loop_tag::bovi, inner_tag::memory>;
@@ -698,6 +917,154 @@ TEST_CASE("Test behavior of sparse packs", "[SparsePack]") {
       }
     }
   }
+  GIVEN("A pair of fields on a mesh") {
+    const std::vector<int> scalar_shape{N, N, N};
+    Metadata m({Metadata::Independent}, scalar_shape);
+    const std::vector<int> vector_shape{N, N, N, 3};
+    Metadata m_vector({Metadata::Independent, Metadata::Vector}, vector_shape);
+
+    auto pkg = std::make_shared<StateDescriptor>("Test package");
+    pkg->AddField(v1::name(), m);
+    pkg->AddField(v3::name(), m_vector);
+    BlockList_t block_list = MakeBlockList(pkg, NBLOCKS, N, NDIM);
+
+    MeshData<Real> mesh_data("base");
+    mesh_data.Initialize(block_list, nullptr);
+    auto ib = block_list[0]->cellbounds.GetBoundsI(IndexDomain::entire);
+    auto jb = block_list[0]->cellbounds.GetBoundsJ(IndexDomain::entire);
+    auto kb = block_list[0]->cellbounds.GetBoundsK(IndexDomain::entire);
+
+    const parthenon::IndexRange ibi{ib.s + 1, ib.e - 1}, jbi{jb.s + 1, jb.e - 1},
+        kbi{kb.s + 1, kb.e - 1};
+
+    WHEN("We get a sparse pack for a virtual variable dependent on "
+         "our pair of fields, and initialize the dependent vars.") {
+      using Axis = parthenon::Axis;
+      auto desc = parthenon::MakePackDescriptor<d1, d1_subpack, gradient<Axis::I>,
+                                                curl<Axis::J>, der3>(pkg.get());
+      auto pack = desc.GetPack(&mesh_data);
+      par_for(
+          "initialize d1", 0, NBLOCKS - 1, kb, jb, ib,
+          KOKKOS_LAMBDA(int b, int k, int j, int i) {
+            Real n = b + k * j * i;
+            Real m = i * i + j * j + k * k + b * b;
+            pack(b, v1(), k, j, i) = n;
+            pack(b, v3(0), k, j, i) = 0.;
+            pack(b, v3(1), k, j, i) = m;
+            pack(b, v3(2), k, j, i) = m + n;
+          });
+
+      THEN("We can correctly evaluate the virtual field.") {
+        int nwrong = 0;
+        par_reduce(
+            "check virtual", 0, NBLOCKS - 1, kb, jb, ib,
+            KOKKOS_LAMBDA(int b, int k, int j, int i, int &ltot) {
+              const Real x = 3.8;
+              const Real answer = pack(b, v1(), k, j, i) * pack(b, v3(1), k, j, i) *
+                                      pack(b, v3(2), k, j, i) +
+                                  x;
+              if (pack(b, d1(x), k, j, i) != answer) {
+                ltot += 1;
+              }
+              if (pack(b, d1(x), k, j, i) != pack(b, d1_subpack(), k, j, i, x)) {
+                ltot += 1;
+              }
+            },
+            nwrong);
+        REQUIRE(nwrong == 0);
+      }
+      THEN("We can correctly use subpacks to get gradients.") {
+        par_for(
+            "initialize v1", 0, NBLOCKS - 1, kb, jb, ib,
+            KOKKOS_LAMBDA(int b, int k, int j, int i) {
+              pack(b, v1(), k, j, i) = 11.0 * i + 22.0 * j + 33.0 * k;
+            });
+        int nwrong = 0;
+        par_reduce(
+            "check virtual", 0, NBLOCKS - 1, kbi, jbi, ibi,
+            KOKKOS_LAMBDA(int b, int k, int j, int i, int &ltot) {
+              if (std::abs(pack(b, gradient<parthenon::Axis::I>(), k, j, i) - 11.0) >=
+                  1.e-12) {
+                ltot += 1;
+              }
+              if (std::abs(pack(b, gradient<parthenon::Axis::J>(), k, j, i) - 22.0) >=
+                  1.e-12) {
+                ltot += 1;
+              }
+              if (std::abs(pack(b, gradient<parthenon::Axis::K>(), k, j, i) - 33.0) >=
+                  1.e-12) {
+                ltot += 1;
+              }
+              using TE = parthenon::TopologicalElement;
+              if (std::abs(pack(b, gradient<parthenon::Axis::I>(), k, j, i) -
+                           pack(b, TE::F1, gradient_te(), k, j, i)) >= 1.e-12) {
+                ltot += 1;
+              }
+              if (std::abs(pack(b, gradient<parthenon::Axis::J>(), k, j, i) -
+                           pack(b, TE::F2, gradient_te(), k, j, i)) >= 1.e-12) {
+                ltot += 1;
+              }
+              if (std::abs(pack(b, gradient<parthenon::Axis::K>(), k, j, i) -
+                           pack(b, TE::F3, gradient_te(), k, j, i)) >= 1.e-12) {
+                ltot += 1;
+              }
+            },
+            nwrong);
+        REQUIRE(nwrong == 0);
+      }
+      THEN("We can use 2D subpacks") {
+        const Real axy = 12.;
+        const Real azx = 48.32;
+        const Real ayz = 3.14;
+        par_for(
+            "Initialize a simple vector with curl", 0, NBLOCKS - 1, kb, jb, ib,
+            KOKKOS_LAMBDA(int b, int k, int j, int i) {
+              pack(b, v3(0), k, j, i) = ayz * static_cast<Real>(j * k);
+              pack(b, v3(1), k, j, i) = azx * static_cast<Real>(k * i);
+              pack(b, v3(2), k, j, i) = axy * static_cast<Real>(i * j);
+            });
+
+        int nwrong = 0;
+        par_reduce(
+            "check virtual", 0, NBLOCKS - 1, kbi, jbi, ibi,
+            KOKKOS_LAMBDA(int b, int k, int j, int i, int &ltot) {
+              if (std::abs(pack(b, curl<parthenon::Axis::K>(), k, j, i) -
+                           (azx - ayz) * k) >= 1.e-12) {
+                ltot += 1;
+              }
+              if (std::abs(pack(b, curl<parthenon::Axis::J>(), k, j, i) -
+                           (ayz - axy) * j) >= 1.e-12) {
+                ltot += 1;
+              }
+              if (std::abs(pack(b, curl<parthenon::Axis::I>(), k, j, i) -
+                           (axy - azx) * i) >= 1.e-12) {
+                ltot += 1;
+              }
+            },
+            nwrong);
+        REQUIRE(nwrong == 0);
+      }
+      THEN("We can use 3D subpacks") {
+        const Real a = 12.345;
+        par_for(
+            "Initialize a simple vector with curl", 0, NBLOCKS - 1, kb, jb, ib,
+            KOKKOS_LAMBDA(int b, int k, int j, int i) {
+              pack(b, v1(), k, j, i) = a * static_cast<Real>(i * j * k);
+            });
+
+        int nwrong = 0;
+        par_reduce(
+            "check virtual", 0, NBLOCKS - 1, kbi, jbi, ibi,
+            KOKKOS_LAMBDA(int b, int k, int j, int i, int &ltot) {
+              if (std::abs(pack(b, der3(), k, j, i) - a) >= 1.e-12) {
+                ltot += 1;
+              }
+            },
+            nwrong);
+        REQUIRE(nwrong == 0);
+      }
+    }
+  }
 }
 
 TEST_CASE("Pack views preserve the loop abstraction contract through pack_view",