diff --git a/CHANGELOG.md b/CHANGELOG.md
index 8777c52338934..286d17381102e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -41,6 +41,7 @@
 
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 1398]](https://github.com/parthenon-hpc-lab/parthenon/pull/1398) Generalize sparse identifiers and add sparse control groups. Downstream sparse-state access may need to change from `sparse_idx` to `sparse_idx()`.
 - [[PR 1351]](https://github.com/parthenon-hpc-lab/parthenon/pull/1351) Bump Kokkos 5 & C++20
 - [[PR 1377]](https://github.com/parthenon-hpc-lab/parthenon/pull/1377) Extend Initialization Hierarchy
 - [[PR 1376]](https://github.com/parthenon-hpc-lab/parthenon/pull/1376) Refactor SwarmPacks
diff --git a/dev/plan_histories/1398.md b/dev/plan_histories/1398.md
new file mode 100644
index 0000000000000..40110e8bf30e2
--- /dev/null
+++ b/dev/plan_histories/1398.md
@@ -0,0 +1,53 @@
+<!-- This file was made in part with generative AI. -->
+
+## Goal
+
+Generalize sparse identifiers so Parthenon can represent more than a single integer
+index per sparse field, while keeping the existing scalar behavior and restart behavior
+unchanged for current codes. At the same time, introduce a cleaner control-group model
+for sparse fields so allocation and deallocation can operate on groups of controlling
+fields rather than only a single controller label.
+
+## Functional Plan
+
+- Introduce a dedicated `SparseID` value type and thread it through the sparse field
+  model instead of storing raw integers directly in the core interfaces.
+- Keep the sparse label format backward compatible for current scalar IDs so existing
+  sparse names continue to look the same.
+- Add support for multi-component sparse IDs in a way that remains device-visible and
+  inexpensive enough to be used in hot code paths.
+- Mark the sparse-ID comparison operators as device-available so pack selection and
+  other device code can reason about sparse ordering without special handling.
+- Keep the host-facing sparse APIs compatible with existing scalar call sites by
+  allowing integer inputs to be converted to `SparseID` at the API boundary.
+- Update sparse pools so they store and operate on typed sparse IDs internally rather
+  than treating sparse IDs as anonymous integers.
+- Add a sparse-ID projection rule for sparse pools so a pool can group fields by a
+  controller key derived from the sparse index space.
+- Extend `StateDescriptor` so it resolves sparse-pool control groups into the same
+  controller-reverse-map machinery used by ordinary fields.
+- Treat control groups as the source of truth for grouped sparse control, while keeping
+  the legacy single-controller accessor available as a compatibility shim.
+- Keep per-field sparse deallocation counters intact, but evaluate deallocation at the
+  group level by requiring every control member to be ready before the group is removed.
+- Make sparse deallocation ignore dense-only controller groups so dense fields do not
+  enter the sparse dealloc path.
+- Preserve restart exactness by keeping the sparse deallocation metadata and sparse
+  allocation behavior compatible with existing restart files.
+- Keep the existing sparse-advection example and restart flow working as the primary
+  exactness check while the sparse-control machinery changes underneath it.
+- Add focused unit coverage for the new sparse-ID type and for grouped sparse control.
+- Update sparse documentation so the split between sparse IDs, sparse pools, and
+  controlling groups is described clearly.
+
+## Scope Boundaries
+
+- Do not change the public downstream package API unless it is required for backward
+  compatibility with existing scalar sparse call sites.
+- Keep dense fields and sparse fields conceptually separate; only unify their internal
+  registration and control bookkeeping where that reduces duplication.
+- Do not change the numerical update kernels in the example problems unless a restart or
+  sparse-control bug forces a targeted fix.
+- Preserve the current sparse-advection restart gold output as the exactness oracle for
+  this PR.
+
diff --git a/example/advection/advection_package.cpp b/example/advection/advection_package.cpp
index 9f96e90166dfb..9c21868afee87 100644
--- a/example/advection/advection_package.cpp
+++ b/example/advection/advection_package.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <algorithm>
 #include <cmath>
 #include <iostream>
@@ -390,8 +392,10 @@ void PostFill(MeshBlockData<Real> *rc) {
     IndexRange kb = pmb->cellbounds.GetBoundsK(IndexDomain::entire);
 
     // check that we have the sparse indices we want
-    pmb->AllocSparseID("one_minus_sqrt_one_minus_advected_sq", 12);
-    pmb->AllocSparseID("one_minus_sqrt_one_minus_advected_sq", 37);
+    pmb->AllocSparseID("one_minus_sqrt_one_minus_advected_sq",
+                       parthenon::SparseID::Scalar(12));
+    pmb->AllocSparseID("one_minus_sqrt_one_minus_advected_sq",
+                       parthenon::SparseID::Scalar(37));
 
     // packing in principle unnecessary/convoluted here and just done for demonstration
     std::vector<std::string> vars(
@@ -402,7 +406,9 @@ void PostFill(MeshBlockData<Real> *rc) {
     const int in = imap.get("one_minus_advected_sq").first;
     // we can get sparse fields either by specifying base name and sparse id, or the full
     // name
-    const int out12 = imap.get("one_minus_sqrt_one_minus_advected_sq", 12).first;
+    const int out12 =
+        imap.get("one_minus_sqrt_one_minus_advected_sq", parthenon::SparseID::Scalar(12))
+            .first;
     const int out37 = imap.get("one_minus_sqrt_one_minus_advected_sq_37").first;
     const auto num_vars = rc->Get("advected").data.GetDim(4);
     pmb->par_for(
diff --git a/example/calculate_pi/calculate_pi.cpp b/example/calculate_pi/calculate_pi.cpp
index 3fe83b71343aa..74322891ff3be 100644
--- a/example/calculate_pi/calculate_pi.cpp
+++ b/example/calculate_pi/calculate_pi.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 // Standard Includes
 #include <iostream>
 #include <memory>
@@ -72,8 +74,8 @@ void SetInOrOut(MeshBlockData<Real> *rc) {
       return;
     }
 
-    pmb->AllocSparseID("in_or_out", 0);
-    v = rc->Get("in_or_out", 0).data;
+    pmb->AllocSparseID("in_or_out", parthenon::SparseID::Scalar(0));
+    v = rc->Get("in_or_out", parthenon::SparseID::Scalar(0)).data;
   } else {
     v = rc->Get("in_or_out").data;
   }
diff --git a/example/fine_advection/parthenon_app_inputs.cpp b/example/fine_advection/parthenon_app_inputs.cpp
index 8025970f35ca9..a30eb10845d91 100644
--- a/example/fine_advection/parthenon_app_inputs.cpp
+++ b/example/fine_advection/parthenon_app_inputs.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <memory>
 #include <sstream>
 #include <string>
@@ -65,7 +67,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   if (do_regular_advection) {
     const int sparse_size = pkg->Param<int>("sparse_size");
     for (int s = 0; s < sparse_size; ++s)
-      pmb->AllocSparseID(phi::name(), s);
+      pmb->AllocSparseID(phi::name(), SparseID::Scalar(s));
   }
   static auto desc = parthenon::MakePackDescriptor<phi, phi_fine, C, D>(data.get());
   auto pack = desc.GetPack(data.get());
diff --git a/example/sparse_advection/parthenon_app_inputs.cpp b/example/sparse_advection/parthenon_app_inputs.cpp
index f55e1dfc2aab4..13d5332695c1b 100644
--- a/example/sparse_advection/parthenon_app_inputs.cpp
+++ b/example/sparse_advection/parthenon_app_inputs.cpp
@@ -10,6 +10,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <algorithm>
 #include <cstdio>
 #include <limits>
@@ -89,15 +91,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
       VariablePack<Real> v;
 
       if (restart_test) {
-        pmb->AllocSparseID("shape_shift", 1);
-        pmb->AllocSparseID("shape_shift", 3);
-        pmb->AllocSparseID("shape_shift", 4);
+        pmb->AllocSparseID("shape_shift", parthenon::SparseID::Scalar(1));
+        pmb->AllocSparseID("shape_shift", parthenon::SparseID::Scalar(3));
+        pmb->AllocSparseID("shape_shift", parthenon::SparseID::Scalar(4));
 
         v = data->PackVariables(
             std::vector<std::string>{"dense_A", "dense_B", "shape_shift"});
       } else {
-        pmb->AllocSparseID("sparse", f);
-        v = data->PackVariables(std::vector<std::string>{MakeVarLabel("sparse", f)});
+        pmb->AllocSparseID("sparse", parthenon::SparseID::Scalar(f));
+        v = data->PackVariables(std::vector<std::string>{
+            parthenon::MakeVarLabel("sparse", parthenon::SparseID::Scalar(f))});
       }
 
       pmb->par_for(
@@ -127,7 +130,7 @@ void PostStepDiagnosticsInLoop(Mesh *mesh, ParameterInput *pin, const SimTime &t
   for (auto &pmb : mesh->block_list) {
     auto rc = pmb->meshblock_data.Get(); // get base container
     for (int i = 0; i < n; ++i) {
-      if (rc->IsAllocated("sparse", i)) {
+      if (rc->IsAllocated("sparse", parthenon::SparseID::Scalar(i))) {
         num_allocated[i] += 1;
       }
     }
diff --git a/example/sparse_advection/sparse_advection_package.cpp b/example/sparse_advection/sparse_advection_package.cpp
index 8351a0ffb4382..13927d492bbc3 100644
--- a/example/sparse_advection/sparse_advection_package.cpp
+++ b/example/sparse_advection/sparse_advection_package.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <algorithm>
 #include <cmath>
 #include <limits>
@@ -76,7 +78,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
     SparsePool pool("sparse", m);
 
     for (int sid = 0; sid < NUM_FIELDS; ++sid) {
-      pool.Add(sid);
+      pool.Add(parthenon::SparseID::Scalar(sid));
     }
     pkg->AddSparsePool(pool);
   }
@@ -93,10 +95,11 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
                        Metadata::FillGhost, Metadata::Sparse});
 
     SparsePool pool("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,
+    pool.Add(parthenon::SparseID::Scalar(1), std::vector<int>{1},
+             std::vector<std::string>{"scalar"});
+    pool.Add(parthenon::SparseID::Scalar(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);
+    pool.Add(parthenon::SparseID::Scalar(4), std::vector<int>{4}, Metadata::Vector);
 
     pkg->AddSparsePool(pool);
   }
@@ -216,7 +219,7 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
 
         for (int n = 0; n < nvar; n++) {
           if (!v.IsAllocated(n)) continue;
-          const auto this_v = vx[v(n).sparse_id % NUM_FIELDS];
+          const auto this_v = vx[v(n).sparse_id() % 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;
           });
@@ -246,7 +249,7 @@ TaskStatus CalculateFluxes(std::shared_ptr<MeshBlockData<Real>> &rc) {
 
           for (int n = 0; n < nvar; n++) {
             if (!v.IsAllocated(n)) continue;
-            const auto this_v = vy[v(n).sparse_id % NUM_FIELDS];
+            const auto this_v = vy[v(n).sparse_id() % 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/src/interface/meshblock_data.cpp b/src/interface/meshblock_data.cpp
index 2dc949dd4cf69..32fda50b671c2 100644
--- a/src/interface/meshblock_data.cpp
+++ b/src/interface/meshblock_data.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include "interface/meshblock_data.hpp"
 
 #include <algorithm>
@@ -44,7 +46,7 @@ namespace parthenon {
 /// @param sparse_id the sparse id of the variable
 template <typename T>
 void MeshBlockData<T>::AddField(const std::string &base_name, const Metadata &metadata,
-                                int sparse_id) {
+                                SparseID sparse_id) {
   auto pvar = std::make_shared<Variable<T>>(base_name, metadata, sparse_id, pmy_block);
   Add(pvar);
 
diff --git a/src/interface/meshblock_data.hpp b/src/interface/meshblock_data.hpp
index c23fa87fe5ab9..4fc9efa98d980 100644
--- a/src/interface/meshblock_data.hpp
+++ b/src/interface/meshblock_data.hpp
@@ -13,7 +13,7 @@
 #ifndef INTERFACE_MESHBLOCK_DATA_HPP_
 #define INTERFACE_MESHBLOCK_DATA_HPP_
 
-// This file was made in part with generative AI
+// This file was made in part with generative AI.
 
 #include <algorithm>
 #include <limits>
@@ -266,7 +266,8 @@ class MeshBlockData {
 
   const auto &GetUidMap() const { return varUidMap_; }
 
-  Variable<T> &Get(const std::string &base_name, int sparse_id = InvalidSparseID) const {
+  Variable<T> &Get(const std::string &base_name,
+                   SparseID sparse_id = InvalidSparseID) const {
     return *GetVarPtr(MakeVarLabel(base_name, sparse_id));
   }
   Variable<T> &Get(const Uid_t &uid) const { return *(varUidMap_.at(uid)); }
@@ -286,7 +287,8 @@ class MeshBlockData {
     return it->second->IsAllocated();
   }
 
-  inline bool IsAllocated(std::string const &base_name, int sparse_id) const noexcept {
+  inline bool IsAllocated(std::string const &base_name,
+                          SparseID sparse_id) const noexcept {
     return IsAllocated(MakeVarLabel(base_name, sparse_id));
   }
 
@@ -585,14 +587,14 @@ class MeshBlockData {
 
  private:
   void AddField(const std::string &base_name, const Metadata &metadata,
-                int sparse_id = InvalidSparseID);
+                SparseID sparse_id);
 
   void Add(std::shared_ptr<Variable<T>> var) noexcept;
 
   std::shared_ptr<Variable<T>> AllocateSparse(std::string const &label,
                                               bool flag_uninitialized = false);
   std::shared_ptr<Variable<T>> AllocSparseID(std::string const &base_name,
-                                             const int sparse_id) {
+                                             const SparseID sparse_id) {
     return AllocateSparse(MakeVarLabel(base_name, sparse_id));
   }
   void DeallocateSparse(std::string const &label);
diff --git a/src/interface/sparse_pool.cpp b/src/interface/sparse_pool.cpp
index 43814eb601663..f2138df866972 100644
--- a/src/interface/sparse_pool.cpp
+++ b/src/interface/sparse_pool.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <memory>
 #include <string>
 #include <vector>
@@ -23,7 +25,7 @@
 namespace parthenon {
 
 SparsePool::SparsePool(const std::string &base_name, const Metadata &metadata,
-                       const std::vector<int> &sparse_ids,
+                       const std::vector<SparseID> &sparse_ids,
                        const std::vector<std::vector<int>> &shapes,
                        const std::vector<MetadataFlag> &vector_tensor_flags,
                        const std::vector<std::vector<std::string>> &component_labels,
@@ -92,7 +94,7 @@ MakeSparseVarMetadataImpl(Metadata *in, const std::vector<int> &shape,
   return this_metadata;
 }
 
-const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector<int> &shape,
+const Metadata &SparsePool::AddImpl(SparseID sparse_id, const std::vector<int> &shape,
                                     const MetadataFlag *vector_tensor,
                                     const std::vector<std::string> &component_labels) {
   PARTHENON_REQUIRE_THROWS(sparse_id != InvalidSparseID,
@@ -107,22 +109,22 @@ const Metadata &SparsePool::AddImpl(int sparse_id, const std::vector<int> &shape
   }
 
   const auto ins = pool_.insert({sparse_id, *this_metadata});
-  PARTHENON_REQUIRE_THROWS(ins.second, "Tried to add sparse ID " +
-                                           std::to_string(sparse_id) +
-                                           " to sparse pool '" + base_name_ +
+  PARTHENON_REQUIRE_THROWS(ins.second, "Tried to add sparse field '" +
+                                           MakeVarLabel(base_name_, sparse_id) +
+                                           "' to sparse pool '" + base_name_ +
                                            "', but this sparse ID already exists");
 
   return ins.first->second;
 }
 
-const Metadata &SparsePool::Add(int sparse_id, const Metadata &md) {
+const Metadata &SparsePool::Add(SparseID sparse_id, const Metadata &md) {
   PARTHENON_REQUIRE_THROWS(sparse_id != InvalidSparseID,
                            "Tried to add InvalidSparseID to sparse pool " + base_name_);
 
   const auto ins = pool_.insert({sparse_id, md});
-  PARTHENON_REQUIRE_THROWS(ins.second, "Tried to add sparse ID " +
-                                           std::to_string(sparse_id) +
-                                           " to sparse pool '" + base_name_ +
+  PARTHENON_REQUIRE_THROWS(ins.second, "Tried to add sparse field '" +
+                                           MakeVarLabel(base_name_, sparse_id) +
+                                           "' to sparse pool '" + base_name_ +
                                            "', but this sparse ID already exists");
 
   return ins.first->second;
diff --git a/src/interface/sparse_pool.hpp b/src/interface/sparse_pool.hpp
index 1c3323da06733..93ea65d0be8f1 100644
--- a/src/interface/sparse_pool.hpp
+++ b/src/interface/sparse_pool.hpp
@@ -13,7 +13,10 @@
 #ifndef INTERFACE_SPARSE_POOL_HPP_
 #define INTERFACE_SPARSE_POOL_HPP_
 
+// This file was made in part with generative AI.
+
 #include <map>
+#include <set>
 #include <string>
 #include <type_traits>
 #include <utility>
@@ -29,6 +32,17 @@ class MetadataFlag;
 
 class SparsePool {
  public:
+  enum class ControlSparseIDMode { Identity, FirstComponent };
+
+  static std::vector<SparseID> ToSparseIDs(const std::vector<int> &sparse_ids) {
+    std::vector<SparseID> typed_ids;
+    typed_ids.reserve(sparse_ids.size());
+    for (const auto sparse_id : sparse_ids) {
+      typed_ids.emplace_back(SparseID::Scalar(sparse_id));
+    }
+    return typed_ids;
+  }
+
   // Create an empty sparse pool
   SparsePool(const std::string &base_name, const Metadata &metadata,
              const std::string &controller_base_name = "")
@@ -41,12 +55,13 @@ class SparsePool {
   // Create a copy of the sparse pool with a different name
   SparsePool(const std::string &new_base_name, const SparsePool &src)
       : base_name_(new_base_name), controller_base_name_(src.controller_base_name_),
-        shared_metadata_(src.shared_metadata()), pool_(src.pool()) {}
+        shared_metadata_(src.shared_metadata()),
+        control_sparse_id_mode_(src.control_sparse_id_mode_), pool_(src.pool()) {}
 
   // Create a sparse pool with given sparse ids, shapes, Vector/Tensor flags, and optional
   // component labels
   SparsePool(const std::string &base_name, const Metadata &metadata,
-             const std::vector<int> &sparse_ids,
+             const std::vector<SparseID> &sparse_ids,
              const std::vector<std::vector<int>> &shapes,
              const std::vector<MetadataFlag> &vector_tensor_flags,
              const std::vector<std::vector<std::string>> &component_labels = {},
@@ -55,7 +70,8 @@ class SparsePool {
   // Create a sparse pool with given sparse ids and controlling base name and optional
   // shapes and component labels
   SparsePool(const std::string &base_name, const Metadata &metadata,
-             const std::string &controller_base_name, const std::vector<int> &sparse_ids,
+             const std::string &controller_base_name,
+             const std::vector<SparseID> &sparse_ids,
              const std::vector<std::vector<int>> &shapes = {},
              const std::vector<std::vector<std::string>> &component_labels = {})
       : SparsePool(base_name, metadata, sparse_ids, shapes, {}, component_labels,
@@ -63,16 +79,36 @@ class SparsePool {
 
   // Create a sparse pool with given sparse ids and optional shapes and component labels
   SparsePool(const std::string &base_name, const Metadata &metadata,
-             const std::vector<int> &sparse_ids,
+             const std::vector<SparseID> &sparse_ids,
              const std::vector<std::vector<int>> &shapes = {},
              const std::vector<std::vector<std::string>> &component_labels = {})
       : SparsePool(base_name, metadata, sparse_ids, shapes, {}, component_labels, "") {}
 
-  // Create a sparse pool with given sparse ids and component labels
+  // Bridge constructor that accepts scalar sparse ids and converts them to SparseID.
   SparsePool(const std::string &base_name, const Metadata &metadata,
              const std::vector<int> &sparse_ids,
-             const std::vector<std::vector<std::string>> &component_labels)
-      : SparsePool(base_name, metadata, sparse_ids, {}, {}, component_labels, "") {}
+             const std::vector<std::vector<int>> &shapes,
+             const std::vector<MetadataFlag> &vector_tensor_flags,
+             const std::vector<std::vector<std::string>> &component_labels = {},
+             const std::string &controller_base_name = "")
+      : SparsePool(base_name, metadata, ToSparseIDs(sparse_ids), shapes,
+                   vector_tensor_flags, component_labels, controller_base_name) {}
+
+  // Bridge constructor that accepts scalar sparse ids and converts them to SparseID.
+  SparsePool(const std::string &base_name, const Metadata &metadata,
+             const std::string &controller_base_name, const std::vector<int> &sparse_ids,
+             const std::vector<std::vector<int>> &shapes = {},
+             const std::vector<std::vector<std::string>> &component_labels = {})
+      : SparsePool(base_name, metadata, controller_base_name, ToSparseIDs(sparse_ids),
+                   shapes, component_labels) {}
+
+  // Bridge constructor that accepts scalar sparse ids and converts them to SparseID.
+  SparsePool(const std::string &base_name, const Metadata &metadata,
+             const std::vector<int> &sparse_ids,
+             const std::vector<std::vector<int>> &shapes = {},
+             const std::vector<std::vector<std::string>> &component_labels = {})
+      : SparsePool(base_name, metadata, ToSparseIDs(sparse_ids), shapes,
+                   component_labels) {}
 
   // template on variable type
   template <typename T, typename... Args>
@@ -83,15 +119,35 @@ class SparsePool {
   const std::string &base_name() const { return base_name_; }
   const std::string &controller_base_name() const { return controller_base_name_; }
   const Metadata &shared_metadata() const { return shared_metadata_; }
-  const std::map<int, Metadata> &pool() const { return pool_; }
+  const std::map<SparseID, Metadata> &pool() const { return pool_; }
   auto size() const { return pool_.size(); }
+  ControlSparseIDMode control_sparse_id_mode() const { return control_sparse_id_mode_; }
+  void SetControlSparseIDMode(ControlSparseIDMode mode) {
+    control_sparse_id_mode_ = mode;
+  }
+  SparseID ControlSparseID(SparseID sparse_id) const {
+    if (control_sparse_id_mode_ == ControlSparseIDMode::FirstComponent) {
+      return SparseID::Scalar(sparse_id(0));
+    }
+    return sparse_id;
+  }
+  ControlGroup ControlGroupFor(SparseID sparse_id) const {
+    ControlGroup control_group;
+    const auto control_id = ControlSparseID(sparse_id);
+    for (const auto &pair : pool_) {
+      if (ControlSparseID(pair.first) == control_id) {
+        control_group.emplace(base_name_, pair.first);
+      }
+    }
+    return control_group;
+  }
 
   // Add a new sparse ID to the pool with optional arguments:
   // shape: use this shape if not {}, otherwise use shape from shared metadata (the
   // Vector/Tensor flag will be copied from shared metadata)
   // component_labels: use these component labels if not {}, otherwise use component
   // labels from shared metadata
-  const Metadata &Add(int sparse_id, const std::vector<int> &shape = {},
+  const Metadata &Add(SparseID sparse_id, const std::vector<int> &shape = {},
                       const std::vector<std::string> &component_labels = {}) {
     return AddImpl(sparse_id, shape, nullptr, component_labels);
   }
@@ -99,25 +155,26 @@ class SparsePool {
   // As above, but explicitly set Vector/Tensor metadata flag. Valid values for
   // vector_tensor are: None (unset both Vector and Tensor flag), Vector (set only Vector
   // flag), Tensor (set only Tensor flag)
-  const Metadata &Add(int sparse_id, const std::vector<int> &shape,
+  const Metadata &Add(SparseID sparse_id, const std::vector<int> &shape,
                       MetadataFlag vector_tensor,
                       const std::vector<std::string> &component_labels = {}) {
     return AddImpl(sparse_id, shape, &vector_tensor, component_labels);
   }
 
-  const Metadata &Add(int sparse_id, const std::vector<std::string> &component_labels) {
+  const Metadata &Add(SparseID sparse_id,
+                      const std::vector<std::string> &component_labels) {
     return AddImpl(sparse_id, {}, nullptr, component_labels);
   }
 
   // Let someone specify arbitrary metadata for this field
-  const Metadata &Add(int sparse_id, const Metadata &md);
+  const Metadata &Add(SparseID sparse_id, const Metadata &md);
 
  private:
   // TODO(JL) Once we have C++17 with std::optional, we can use
   // std::optional<MetadataFlag> instead of a pointer. We need to differentiate between
   // the getting a value form the user and not getting a value, but there is no good
   // default value
-  const Metadata &AddImpl(int sparse_id, const std::vector<int> &shape,
+  const Metadata &AddImpl(SparseID sparse_id, const std::vector<int> &shape,
                           const MetadataFlag *vector_tensor,
                           const std::vector<std::string> &component_labels);
 
@@ -125,10 +182,10 @@ class SparsePool {
   const std::string controller_base_name_;
 
   Metadata shared_metadata_;
+  ControlSparseIDMode control_sparse_id_mode_ = ControlSparseIDMode::Identity;
   // Metadata per sparse id
-  // JMM: note that this map SHOULD be ordered as sparse ids, being
-  // integers, have an implicit ordering.
-  std::map<int, Metadata> pool_;
+  // Sparse IDs provide an explicit ordering for deterministic iteration.
+  std::map<SparseID, Metadata> pool_;
 };
 
 } // namespace parthenon
diff --git a/src/interface/state_descriptor.cpp b/src/interface/state_descriptor.cpp
index d970887e50aea..4b82149db4266 100644
--- a/src/interface/state_descriptor.cpp
+++ b/src/interface/state_descriptor.cpp
@@ -11,6 +11,9 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
+#include <algorithm>
 #include <iomanip>
 #include <iostream>
 #include <memory>
@@ -167,8 +170,7 @@ class FieldProvider : public VariableProvider {
       const auto &src_pool = pkg->GetSparsePool(base_name);
       added = state_->AddSparsePool(new_name, src_pool);
     } else {
-      auto controller = packages_.Get(package)->GetFieldController(base_name);
-      added = state_->AddField(new_name, metadata, controller);
+      added = state_->AddField(new_name, metadata, pkg->GetFieldControlGroup(base_name));
     }
 
     PARTHENON_REQUIRE_THROWS(added, "Couldn't add private field '" + base_name +
@@ -182,8 +184,7 @@ class FieldProvider : public VariableProvider {
       const auto &pool = pkg->GetSparsePool(base_name);
       added = state_->AddSparsePool(pool);
     } else {
-      auto controller = packages_.Get(package)->GetFieldController(base_name);
-      added = state_->AddField(base_name, metadata, controller);
+      added = state_->AddField(base_name, metadata, pkg->GetFieldControlGroup(base_name));
     }
 
     PARTHENON_REQUIRE_THROWS(added, "Couldn't add provided field '" + base_name +
@@ -209,8 +210,8 @@ class FieldProvider : public VariableProvider {
       for (auto &pair : packages_.AllPackages()) {
         auto &package = pair.second;
         if (package->FieldPresent(base_name)) {
-          auto controller = package->GetFieldController(base_name);
-          added = state_->AddField(base_name, metadata, controller);
+          added = state_->AddField(base_name, metadata,
+                                   package->GetFieldControlGroup(base_name));
           break;
         }
       }
@@ -288,19 +289,37 @@ bool StateDescriptor::AddSwarmValue(const std::string &value_name,
 
 bool StateDescriptor::AddField(const std::string &field_name, const Metadata &m_in,
                                const std::string &controlling_field) {
+  if (controlling_field == "") {
+    return AddField(field_name, m_in, ControlGroup{VarID(field_name)});
+  }
+  return AddField(field_name, m_in, std::vector<std::string>{controlling_field});
+}
+
+bool StateDescriptor::AddField(const std::string &field_name, const Metadata &m_in,
+                               const std::vector<std::string> &controlling_fields) {
+  ControlGroup control_group;
+  for (const auto &controlling_field : controlling_fields) {
+    control_group.emplace(controlling_field);
+  }
+  if (control_group.empty()) control_group.emplace(field_name);
+  return AddField(field_name, m_in, control_group);
+}
+
+bool StateDescriptor::AddField(const std::string &field_name, const Metadata &m_in,
+                               const ControlGroup &controlling_fields) {
   Metadata m = m_in; // so we can modify it
   if (m.IsSet(Metadata::Sparse)) {
     PARTHENON_THROW(
         "Tried to add a sparse field with AddField, use AddSparsePool instead");
   }
   if (!m.IsSet(GetMetadataFlag())) m.Set(GetMetadataFlag());
-  VarID controller = VarID(controlling_field);
-  if (controlling_field == "") controller = VarID(field_name);
-  return AddFieldImpl_(VarID(field_name), m, controller);
+  ControlGroup control_group = controlling_fields;
+  if (control_group.empty()) control_group.emplace(VarID(field_name));
+  return AddFieldImpl_(VarID(field_name), m, control_group);
 }
 
 bool StateDescriptor::AddFieldImpl_(const VarID &vid, const Metadata &m_in,
-                                    const VarID &control_vid) {
+                                    const ControlGroup &control_group) {
   Metadata m = m_in; // Force const correctness
 
   const std::string &assoc = m.getAssociated();
@@ -313,13 +332,13 @@ bool StateDescriptor::AddFieldImpl_(const VarID &vid, const Metadata &m_in,
     if (m.IsSet(Metadata::WithFluxes) && m.GetFluxName() == "") {
       auto fId = VarID{internal_fluxname + internal_varname_seperator + vid.base_name,
                        vid.sparse_id};
-      AddFieldImpl_(fId, *(m.GetSPtrFluxMetadata()), control_vid);
+      AddFieldImpl_(fId, *(m.GetSPtrFluxMetadata()), control_group);
       m.SetFluxName(fId.label());
     }
     labelToVidMap_.insert({vid.label(), vid});
     metadataMap_.insert({vid, m});
     refinementFuncMaps_.Register(m, vid.label());
-    allocControllerReverseMap_.insert({vid, control_vid});
+    allocControllerReverseMap_.insert({vid, control_group});
     // Add this variable to the set of unique IDs at the
     // earliest possible time
     Variable<Real>::GetUniqueID(vid.label());
@@ -341,12 +360,12 @@ bool StateDescriptor::AddSparsePoolImpl_(const SparsePool &pool) {
   sparsePoolMap_.insert({pool.base_name(), pool});
   refinementFuncMaps_.Register(pool.shared_metadata(), pool.base_name());
 
-  std::string controller_base = pool.controller_base_name();
-  if (controller_base == "") controller_base = pool.base_name();
   // add all the sparse fields
   for (const auto itr : pool.pool()) {
-    if (!AddFieldImpl_(VarID(pool.base_name(), itr.first), itr.second,
-                       VarID(controller_base, itr.first))) {
+    const auto control_group = pool.controller_base_name().empty()
+                                   ? pool.ControlGroupFor(itr.first)
+                                   : ControlGroup{};
+    if (!AddFieldImpl_(VarID(pool.base_name(), itr.first), itr.second, control_group)) {
       // a field with this name already exists, this would leave the StateDescriptor in an
       // inconsistent state, so throw
       PARTHENON_THROW("Couldn't add sparse field " +
@@ -389,12 +408,21 @@ bool StateDescriptor::FlagsPresent(std::vector<MetadataFlag> const &flags,
 }
 
 std::string StateDescriptor::GetFieldController(const std::string &field_name) {
+  const auto &controller_group = GetFieldControlGroup(field_name);
+  PARTHENON_REQUIRE_THROWS(
+      controller_group.size() == 1,
+      "Asking for legacy controlling field for grouped control set " + field_name +
+          "; use GetFieldControlGroup instead");
+  return controller_group.begin()->label();
+}
+
+const ControlGroup &StateDescriptor::GetFieldControlGroup(const std::string &field_name) {
   VarID field_id(field_name);
   auto controller = allocControllerReverseMap_.find(field_id);
   PARTHENON_REQUIRE(controller != allocControllerReverseMap_.end(),
                     "Asking for controlling field that is not in this package (" +
                         field_name + ")");
-  return controller->second.label();
+  return controller->second;
 }
 
 bool StateDescriptor::SwarmValuePresent(const std::string &value_name,
@@ -420,7 +448,7 @@ std::vector<std::string> StateDescriptor::GetControlVariables() {
 
 // retrieve metadata for a specific field
 const Metadata &StateDescriptor::FieldMetadata(const std::string &base_name,
-                                               int sparse_id) const {
+                                               SparseID sparse_id) const {
   const auto itr = metadataMap_.find(VarID(base_name, sparse_id));
   PARTHENON_REQUIRE_THROWS(itr != metadataMap_.end(),
                            "FieldMetadata: Non-existent field: " +
@@ -465,30 +493,64 @@ std::ostream &operator<<(std::ostream &os, const StateDescriptor &sd) {
   return os;
 }
 
-// Take a map going from variable 1 -> variable that controls variable 1 and invert it
-// to give control variable -> list of variables controlled by control variable.
-// TODO(LFR): Currently, calling this repeatedly will just add a controlled variable
-// to the vector of a controlling variable repeatedly. I think this shouldn't cause any
-// issues since allocating or deallocating a variable multiple times in a row is the
-// same as allocating it or deallocating it once. That being said, it could be switched
-// to an unordered_set from a vector so that variable names can only show up once. Also,
-// it is not clear to me exactly what behavior this should have if invert control map is
-// called more than once (I think the normal use case would be for just a single call
-// during resolution of the combined state descriptor). It may be that we should be
-// calling allocControllerMap_.clear() before starting the for_each loop.
+// Take the reverse controller map, which stores the controlling group for each variable,
+// and invert it into the runtime fan-out map used by allocation and deallocation.
+// The forward map is keyed by individual controller labels, so grouped controllers end up
+// sharing the same controlled-variable list. We deduplicate controlled labels because a
+// variable may appear in more than one controller member's reverse entry, but repeated
+// allocation/deallocation of the same sparse field is harmless.
 void StateDescriptor::InvertControllerMap() {
-  std::for_each(allocControllerReverseMap_.begin(), allocControllerReverseMap_.end(),
-                [this](const auto &pair) {
-                  auto var = pair.first.label();
-                  auto cont = pair.second.label();
-                  auto iter = allocControllerMap_.find(cont);
-                  if (iter == allocControllerMap_.end()) {
-                    allocControllerMap_.emplace(
-                        std::make_pair(cont, std::vector<std::string>{var}));
-                  } else {
-                    iter->second.push_back(var);
-                  }
-                });
+  allocControllerMap_.clear();
+  allocControlGroups_.clear();
+  for (const auto &pair : allocControllerReverseMap_) {
+    const auto var = pair.first.label();
+    const auto &control_group = pair.second;
+
+    allocControlGroups_.insert(control_group);
+    for (const auto &cont : control_group) {
+      auto iter = allocControllerMap_.find(cont.label());
+      if (iter == allocControllerMap_.end()) {
+        allocControllerMap_.emplace(
+            std::make_pair(cont.label(), std::vector<std::string>{var}));
+      } else if (std::find(iter->second.begin(), iter->second.end(), var) ==
+                 iter->second.end()) {
+        iter->second.push_back(var);
+      }
+    }
+  }
+}
+
+// Sparse pools may project their sparse ID space into a smaller controller key space.
+// That projection is resolved here once all sparse pools are known, so every sparse
+// variable ends up with the concrete control group of real field labels that owns it.
+void StateDescriptor::ResolveSparseControllerGroups() {
+  for (const auto &pool_pair : sparsePoolMap_) {
+    const auto &pool = pool_pair.second;
+    const bool has_controller_pool = !pool.controller_base_name().empty();
+    std::string controller_base = pool.controller_base_name();
+    if (controller_base == "") controller_base = pool.base_name();
+
+    auto controller_pool = sparsePoolMap_.find(controller_base);
+    PARTHENON_REQUIRE_THROWS(controller_pool != sparsePoolMap_.end(),
+                             "Couldn't find controlling sparse pool '" + controller_base +
+                                 "' for sparse pool '" + pool.base_name() + "'");
+
+    for (const auto &entry : pool.pool()) {
+      const auto vid = VarID(pool.base_name(), entry.first);
+      const auto control_group =
+          controller_pool->second.ControlGroupFor(pool.ControlSparseID(entry.first));
+      auto [iter, inserted] = allocControllerReverseMap_.emplace(vid, control_group);
+      if (!inserted) {
+        if (has_controller_pool) {
+          iter->second = control_group;
+        } else {
+          PARTHENON_REQUIRE_THROWS(iter->second == control_group,
+                                   "Sparse field '" + vid.label() +
+                                       "' resolved to more than one control group");
+        }
+      }
+    }
+  }
 }
 
 // Takes all packages and combines them into a single state descriptor
@@ -561,6 +623,7 @@ StateDescriptor::CreateResolvedStateDescriptor(Packages_t &packages) {
   field_tracker.CheckOverridable(&field_provider);
   swarm_tracker.CheckOverridable(&swarm_provider);
 
+  state->ResolveSparseControllerGroups();
   state->InvertControllerMap();
 
   return state;
@@ -577,8 +640,9 @@ StateDescriptor::CreateResolvedStateDescriptor(Packages_t &packages) {
 std::vector<std::string>
 StateDescriptor::GetVariableNames(const std::vector<std::string> &requested_names,
                                   const Metadata::FlagCollection &flags,
-                                  const std::vector<int> &sparse_ids) {
-  std::unordered_set<int> sparse_ids_set(sparse_ids.begin(), sparse_ids.end());
+                                  const std::vector<SparseID> &sparse_ids) {
+  std::unordered_set<SparseID, SparseIDHasher> sparse_ids_set(sparse_ids.begin(),
+                                                              sparse_ids.end());
   std::unordered_set<std::string> names;
   std::vector<std::string> names_vec;
   // first add names that are present
@@ -636,20 +700,22 @@ StateDescriptor::GetVariableNames(const std::vector<std::string> &requested_name
 std::vector<std::string>
 StateDescriptor::GetVariableNames(const std::vector<std::string> &requested_names,
                                   const std::vector<int> &sparse_ids) {
-  return GetVariableNames(requested_names, Metadata::FlagCollection(), sparse_ids);
+  return GetVariableNames(requested_names, Metadata::FlagCollection(),
+                          SparsePool::ToSparseIDs(sparse_ids));
 }
 std::vector<std::string>
 StateDescriptor::GetVariableNames(const Metadata::FlagCollection &flags,
                                   const std::vector<int> &sparse_ids) {
-  return GetVariableNames({}, flags, sparse_ids);
+  return GetVariableNames({}, flags, SparsePool::ToSparseIDs(sparse_ids));
 }
 std::vector<std::string>
 StateDescriptor::GetVariableNames(const std::vector<std::string> &requested_names) {
-  return GetVariableNames(requested_names, Metadata::FlagCollection(), {});
+  return GetVariableNames(requested_names, Metadata::FlagCollection(),
+                          std::vector<SparseID>{});
 }
 std::vector<std::string>
 StateDescriptor::GetVariableNames(const Metadata::FlagCollection &flags) {
-  return GetVariableNames({}, flags, {});
+  return GetVariableNames({}, flags, std::vector<SparseID>{});
 }
 
 // Get the total length of this StateDescriptor's variables when packed
diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp
index ba24573f33893..0be18adf94dc8 100644
--- a/src/interface/state_descriptor.hpp
+++ b/src/interface/state_descriptor.hpp
@@ -20,6 +20,7 @@
 #include <limits>
 #include <map>
 #include <memory>
+#include <set>
 #include <sstream>
 #include <string>
 #include <unordered_map>
@@ -231,10 +232,22 @@ class StateDescriptor {
   // field addition / retrieval routines
   bool AddField(const std::string &field_name, const Metadata &m_in,
                 const std::string &controlling_field = "");
+  bool AddField(const std::string &field_name, const Metadata &m_in,
+                const std::vector<std::string> &controlling_fields);
+  bool AddField(const std::string &field_name, const Metadata &m_in,
+                const ControlGroup &controlling_fields);
   template <typename T>
   bool AddField(const Metadata &m, const std::string &controlling_field = "") {
     return AddField(T::name(), m, controlling_field);
   }
+  template <typename T>
+  bool AddField(const Metadata &m, const std::vector<std::string> &controlling_fields) {
+    return AddField(T::name(), m, controlling_fields);
+  }
+  template <typename T>
+  bool AddField(const Metadata &m, const ControlGroup &controlling_fields) {
+    return AddField(T::name(), m, controlling_fields);
+  }
 
   template <typename T>
   bool AddScratch() {
@@ -313,11 +326,16 @@ class StateDescriptor {
   }
   std::vector<std::string> GetVariableNames(const std::vector<std::string> &req_names,
                                             const Metadata::FlagCollection &flags,
-                                            const std::vector<int> &sparse_ids);
+                                            const std::vector<SparseID> &sparse_ids);
   std::vector<std::string> GetVariableNames(const std::vector<std::string> &req_names,
                                             const std::vector<int> &sparse_ids);
   std::vector<std::string> GetVariableNames(const Metadata::FlagCollection &flags,
                                             const std::vector<int> &sparse_ids);
+  std::vector<std::string> GetVariableNames(const std::vector<std::string> &req_names,
+                                            const Metadata::FlagCollection &flags,
+                                            const std::vector<int> &sparse_ids) {
+    return GetVariableNames(req_names, flags, SparsePool::ToSparseIDs(sparse_ids));
+  }
   std::vector<std::string> GetVariableNames(const std::vector<std::string> &req_names);
   std::vector<std::string> GetVariableNames(const Metadata::FlagCollection &flags);
 
@@ -344,7 +362,7 @@ class StateDescriptor {
   }
 
   bool FieldPresent(const std::string &base_name,
-                    int sparse_id = InvalidSparseID) const noexcept {
+                    SparseID sparse_id = InvalidSparseID) const noexcept {
     return metadataMap_.count(VarID(base_name, sparse_id)) > 0;
   }
   bool FieldPresent(const VarID &var_id) const noexcept {
@@ -359,14 +377,19 @@ class StateDescriptor {
   bool SwarmValuePresent(const std::string &value_name,
                          const std::string &swarm_name) const noexcept;
 
+  // Legacy single-controller accessor. Only valid when the control group has one member.
   std::string GetFieldController(const std::string &field_name);
+  // The full controlling group for a field. This is the preferred API for grouped sparse
+  // control and also reflects the underlying model for ordinary single-controller cases.
+  const ControlGroup &GetFieldControlGroup(const std::string &field_name);
   bool ControlVariablesSet() { return (allocControllerMap_.size() > 0); }
+  const std::set<ControlGroup> &GetControlGroups() const { return allocControlGroups_; }
   const std::vector<std::string> &GetControlledVariables(const std::string &field_name);
   std::vector<std::string> GetControlVariables();
 
   // retrieve metadata for a specific field
   const Metadata &FieldMetadata(const std::string &base_name,
-                                int sparse_id = InvalidSparseID) const;
+                                SparseID sparse_id = InvalidSparseID) const;
   // retrieve metadata for a specific swarm
   Metadata &SwarmMetadata(const std::string &swarm_name) noexcept {
     return swarmMetadataMap_[swarm_name];
@@ -503,10 +526,12 @@ class StateDescriptor {
  protected:
   // internal function to add dense/sparse fields. Private because outside classes must
   // use the public interface below
-  bool AddFieldImpl_(const VarID &vid, const Metadata &m, const VarID &control_vid);
+  bool AddFieldImpl_(const VarID &vid, const Metadata &m,
+                     const ControlGroup &control_group);
 
   // add a sparse pool
   bool AddSparsePoolImpl_(const SparsePool &pool);
+  void ResolveSparseControllerGroups();
 
   void InvertControllerMap();
 
@@ -516,8 +541,9 @@ class StateDescriptor {
   // for each variable label (full label for sparse variables) hold metadata
   std::unordered_map<std::string, VarID> labelToVidMap_;
   std::unordered_map<VarID, Metadata, VarIDHasher> metadataMap_;
-  std::unordered_map<VarID, VarID, VarIDHasher> allocControllerReverseMap_;
+  std::unordered_map<VarID, ControlGroup, VarIDHasher> allocControllerReverseMap_;
   std::unordered_map<std::string, std::vector<std::string>> allocControllerMap_;
+  std::set<ControlGroup> allocControlGroups_;
   const std::vector<std::string> nullControl_{};
 
   // for each sparse base name hold its sparse pool
diff --git a/src/interface/var_id.hpp b/src/interface/var_id.hpp
index cdb6aaa3bad44..58018e14aaad4 100644
--- a/src/interface/var_id.hpp
+++ b/src/interface/var_id.hpp
@@ -13,14 +13,85 @@
 #ifndef INTERFACE_VAR_ID_HPP_
 #define INTERFACE_VAR_ID_HPP_
 
+// This file was made in part with generative AI.
+
+#include <cstddef>
 #include <functional>
+#include <limits>
+#include <set>
 #include <string>
 
+#include <Kokkos_Core.hpp>
+
+#include "utils/error_checking.hpp"
+
 namespace parthenon {
 
-inline std::string MakeVarLabel(const std::string &base_name, int sparse_id) {
-  return base_name +
-         (sparse_id == InvalidSparseID ? "" : "_" + std::to_string(sparse_id));
+inline constexpr int InvalidSparseIDValue = std::numeric_limits<int>::min();
+
+struct SparseID {
+  int id_ = InvalidSparseIDValue;
+  int id2_ = InvalidSparseIDValue;
+
+  KOKKOS_INLINE_FUNCTION constexpr SparseID() = default;
+  KOKKOS_INLINE_FUNCTION constexpr SparseID(int id)
+      : id_(id), id2_(InvalidSparseIDValue) {}
+  KOKKOS_INLINE_FUNCTION constexpr SparseID(int id, int id2) : id_(id), id2_(id2) {}
+
+  KOKKOS_INLINE_FUNCTION static constexpr SparseID Scalar(const int id) {
+    return SparseID{id};
+  }
+  KOKKOS_INLINE_FUNCTION static constexpr SparseID Pair(const int id, const int id2) {
+    return SparseID{id, id2};
+  }
+
+  KOKKOS_INLINE_FUNCTION constexpr int operator()() const {
+    PARTHENON_DEBUG_REQUIRE(id2_ == InvalidSparseIDValue,
+                            "SparseID() only valid for scalar sparse IDs");
+    return id_;
+  }
+  KOKKOS_INLINE_FUNCTION constexpr int operator()(const int idx) const {
+    return (idx == 0) ? id_ : id2_;
+  }
+};
+
+inline constexpr SparseID InvalidSparseID{InvalidSparseIDValue, InvalidSparseIDValue};
+
+KOKKOS_INLINE_FUNCTION constexpr bool operator==(const SparseID lhs, const SparseID rhs) {
+  return lhs(0) == rhs(0) && lhs(1) == rhs(1);
+}
+
+KOKKOS_INLINE_FUNCTION constexpr bool operator!=(const SparseID lhs, const SparseID rhs) {
+  return !(lhs == rhs);
+}
+
+KOKKOS_INLINE_FUNCTION constexpr bool operator<(const SparseID lhs, const SparseID rhs) {
+  return (lhs(0) < rhs(0)) || (lhs(0) == rhs(0) && lhs(1) < rhs(1));
+}
+
+KOKKOS_INLINE_FUNCTION constexpr bool operator>(const SparseID lhs, const SparseID rhs) {
+  return rhs < lhs;
+}
+
+KOKKOS_INLINE_FUNCTION constexpr bool operator<=(const SparseID lhs, const SparseID rhs) {
+  return !(rhs < lhs);
+}
+
+KOKKOS_INLINE_FUNCTION constexpr bool operator>=(const SparseID lhs, const SparseID rhs) {
+  return !(lhs < rhs);
+}
+
+inline constexpr bool IsValidSparseID(const SparseID sparse_id) {
+  return sparse_id(0) != InvalidSparseIDValue;
+}
+
+inline std::string MakeVarLabel(const std::string &base_name, const SparseID sparse_id) {
+  if (sparse_id == InvalidSparseID) return base_name;
+  if (sparse_id(1) == InvalidSparseIDValue) {
+    return base_name + "_" + std::to_string(sparse_id(0));
+  }
+  return base_name + "_" + std::to_string(sparse_id(0)) + "_" +
+         std::to_string(sparse_id(1));
 }
 
 /// We uniquely identify a variable by its full label, i.e. base name plus sparse ID.
@@ -35,22 +106,31 @@ inline std::string MakeVarLabel(const std::string &base_name, int sparse_id) {
 /// prolongation/restriction operators.
 struct VarID {
   std::string base_name;
-  int sparse_id;
+  SparseID sparse_id;
 
-  explicit VarID(const std::string base_name, int sparse_id = InvalidSparseID)
+  explicit VarID(const std::string base_name, SparseID sparse_id = InvalidSparseID)
       : base_name(base_name), sparse_id(sparse_id) {}
 
   std::string label() const { return MakeVarLabel(base_name, sparse_id); }
 
   bool operator==(const VarID &other) const { return (label() == other.label()); }
+  bool operator<(const VarID &other) const { return label() < other.label(); }
 };
 
+using ControlGroup = std::set<VarID>;
+
 struct VarIDHasher {
   auto operator()(const VarID &vid) const {
     return std::hash<std::string>{}(vid.label());
   }
 };
 
+struct SparseIDHasher {
+  auto operator()(const SparseID sparse_id) const {
+    return std::hash<int>{}(sparse_id(0)) ^ (std::hash<int>{}(sparse_id(1)) << 1);
+  }
+};
+
 } // namespace parthenon
 
 #endif // INTERFACE_VAR_ID_HPP_
diff --git a/src/interface/variable.cpp b/src/interface/variable.cpp
index 65340a89dd79c..5a87c2dd91c8a 100644
--- a/src/interface/variable.cpp
+++ b/src/interface/variable.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include "interface/variable.hpp"
 
 #include <cstdio>
@@ -31,7 +33,7 @@ namespace parthenon {
 
 template <typename T>
 Variable<T>::Variable(const std::string &base_name, const Metadata &metadata,
-                      int sparse_id, std::weak_ptr<MeshBlock> wpmb)
+                      SparseID sparse_id, std::weak_ptr<MeshBlock> wpmb)
     : m_(metadata), base_name_(base_name), sparse_id_(sparse_id),
       dims_(m_.GetArrayDims(wpmb, false)), coarse_dims_(m_.GetArrayDims(wpmb, true)) {
   PARTHENON_REQUIRE_THROWS(m_.IsSet(Metadata::Real),
diff --git a/src/interface/variable.hpp b/src/interface/variable.hpp
index 23f641aa1637d..96b08ae13806a 100644
--- a/src/interface/variable.hpp
+++ b/src/interface/variable.hpp
@@ -13,6 +13,8 @@
 #ifndef INTERFACE_VARIABLE_HPP_
 #define INTERFACE_VARIABLE_HPP_
 
+// This file was made in part with generative AI.
+
 ///
 /// A Variable type for Placebo-K.
 /// Builds on ParArrayNDs
@@ -60,7 +62,7 @@ class Variable {
   friend class MeshBlockData<T>;
 
  public:
-  Variable(const std::string &base_name, const Metadata &metadata, int sparse_id,
+  Variable(const std::string &base_name, const Metadata &metadata, SparseID sparse_id,
            std::weak_ptr<MeshBlock> wpmb);
 
   Variable() = default;
@@ -113,7 +115,9 @@ class Variable {
   }
 
   /// Get Sparse ID (InvalidSparseID if not sparse)
-  inline int GetSparseID() const { return IsSparse() ? sparse_id_ : InvalidSparseID; }
+  inline SparseID GetSparseID() const {
+    return IsSparse() ? sparse_id_ : InvalidSparseID;
+  }
 
   inline bool IsSparse() const { return m_.IsSet(Metadata::Sparse); }
 
@@ -182,7 +186,7 @@ class Variable {
 
   Metadata m_;
   const std::string base_name_;
-  const int sparse_id_;
+  const SparseID sparse_id_;
   const std::array<int, MAX_VARIABLE_DIMENSION> dims_, coarse_dims_;
 
   // Machinery for giving each variable a unique ID that is faster to
diff --git a/src/interface/variable_pack.hpp b/src/interface/variable_pack.hpp
index 977003557f3fa..a18f9ca812419 100644
--- a/src/interface/variable_pack.hpp
+++ b/src/interface/variable_pack.hpp
@@ -13,6 +13,8 @@
 #ifndef INTERFACE_VARIABLE_PACK_HPP_
 #define INTERFACE_VARIABLE_PACK_HPP_
 
+// This file was made in part with generative AI.
+
 #include <algorithm>
 #include <array>
 #include <cstddef>
@@ -150,7 +152,7 @@ class VarListWithKeys {
   void Add(const std::shared_ptr<Variable<T>> &var,
            const std::unordered_set<int> &sparse_ids = {}) {
     if (!var->IsSparse() || sparse_ids.empty() ||
-        (sparse_ids.count(var->GetSparseID()) > 0)) {
+        (sparse_ids.count(var->GetSparseID()()) > 0)) {
       vars_.push_back(var);
       uids_.push_back(var->GetUniqueID());
       alloc_status_.push_back(var->GetAllocationStatus());
@@ -173,7 +175,8 @@ class PackIndexMap {
 
   const auto &Map() const { return map_; }
 
-  const auto &get(const std::string &base_name, int sparse_id = InvalidSparseID) const {
+  const auto &get(const std::string &base_name,
+                  SparseID sparse_id = InvalidSparseID) const {
     const auto &key = MakeVarLabel(base_name, sparse_id);
     auto itr = map_.find(key);
     if (itr == map_.end()) {
@@ -227,7 +230,7 @@ class PackIndexMap {
     return itr_shape->second;
   }
 
-  bool Has(std::string const &base_name, int sparse_id = InvalidSparseID) const {
+  bool Has(std::string const &base_name, SparseID sparse_id = InvalidSparseID) const {
     return map_.count(MakeVarLabel(base_name, sparse_id)) > 0;
   }
 
@@ -370,7 +373,7 @@ class VariablePack {
   int GetSparseIndex(const int n) const { return GetSparseID(n); }
 
   KOKKOS_FORCEINLINE_FUNCTION
-  bool IsSparse(const int n) const { return GetSparseID() != InvalidSparseID; }
+  bool IsSparse(const int n) const { return GetSparseID() != InvalidSparseID(); }
 
   KOKKOS_FORCEINLINE_FUNCTION
   int VectorComponent(const int n) const {
@@ -531,8 +534,8 @@ void AppendSparseBaseMap(const VariableVector<T> &vars, PackIndexMap *pvmap) {
     int start, stop;
     while (vi != vars.end()) {
       auto &v = *vi;
-      int sparse_id = v->GetSparseID();
-      if (sparse_id != InvalidSparseID) {
+      int sparse_id = v->GetSparseID()();
+      if (sparse_id != InvalidSparseID()) {
         std::vector<int> shape;
         auto mshape = v->metadata().Shape();
         if (mshape.size() > 0) shape.push_back(v->GetDim(4));
@@ -581,7 +584,7 @@ void FillVarView(const VariableVector<T> &vars, int vsize, bool coarse,
     for (int k = 0; k < v->GetDim(6); k++) {
       for (int j = 0; j < v->GetDim(5); j++) {
         for (int i = 0; i < v->GetDim(4); i++) {
-          host_sp(vindex) = v->GetSparseID();
+          host_sp(vindex) = v->GetSparseID()();
 
           // returns 1 for X1DIR, 2 for X2DIR, 3 for X3DIR
           // for tensors, returns flattened index.
diff --git a/src/interface/variable_state.cpp b/src/interface/variable_state.cpp
index b3d363f8613e9..250145e6ae72a 100644
--- a/src/interface/variable_state.cpp
+++ b/src/interface/variable_state.cpp
@@ -11,12 +11,14 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include "interface/variable_state.hpp"
 #include "interface/metadata.hpp"
 
 namespace parthenon {
 
-VariableState::VariableState(const Metadata &md, int sparse_id,
+VariableState::VariableState(const Metadata &md, SparseID sparse_id,
                              const std::array<int, MAX_VARIABLE_DIMENSION> &dims) {
   allocation_threshold = md.GetAllocationThreshold();
   deallocation_threshold = md.GetDeallocationThreshold();
diff --git a/src/interface/variable_state.hpp b/src/interface/variable_state.hpp
index 5e8afcd79d7bc..92d7d79071f15 100644
--- a/src/interface/variable_state.hpp
+++ b/src/interface/variable_state.hpp
@@ -13,11 +13,14 @@
 #ifndef INTERFACE_VARIABLE_STATE_HPP_
 #define INTERFACE_VARIABLE_STATE_HPP_
 
+// This file was made in part with generative AI.
+
 #include <limits>
 
 #include "basic_types.hpp"
 #include "defs.hpp"
 #include "globals.hpp"
+#include "interface/var_id.hpp"
 #include "parthenon_arrays.hpp"
 
 namespace parthenon {
@@ -25,11 +28,9 @@ namespace parthenon {
 // forward declaration
 class Metadata;
 
-static constexpr int InvalidSparseID = std::numeric_limits<int>::min();
-
 struct VariableState : public empty_state_t {
   explicit VariableState(
-      const Metadata &md, int sparse_id = InvalidSparseID,
+      const Metadata &md, SparseID sparse_id = InvalidSparseID,
       const std::array<int, MAX_VARIABLE_DIMENSION> &dims = [] {
         std::array<int, MAX_VARIABLE_DIMENSION> d;
         for (int i = 0; i < MAX_VARIABLE_DIMENSION; ++i)
@@ -39,12 +40,12 @@ struct VariableState : public empty_state_t {
 
   KOKKOS_INLINE_FUNCTION
   VariableState(Real alloc, Real dealloc, Real sparse_default_val = 0.0,
-                int sparse_id = InvalidSparseID)
+                SparseID sparse_id = InvalidSparseID)
       : allocation_threshold(alloc), deallocation_threshold(dealloc),
         sparse_default_val(sparse_default_val), sparse_id(sparse_id) {}
 
   KOKKOS_INLINE_FUNCTION
-  VariableState(Real alloc, Real dealloc, int sparse_id)
+  VariableState(Real alloc, Real dealloc, SparseID sparse_id)
       : allocation_threshold(alloc), deallocation_threshold(dealloc),
         sparse_default_val(0.0), sparse_id(sparse_id) {}
 
@@ -58,7 +59,7 @@ struct VariableState : public empty_state_t {
   Real allocation_threshold;
   Real deallocation_threshold;
   Real sparse_default_val;
-  int sparse_id;
+  SparseID sparse_id;
   int vector_component = NODIR;
   bool initialized = true;
 
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index 493f42cc49449..1b02207f74ab6 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -322,11 +322,18 @@ void MeshBlock::AllocateSparse(std::string const &label, bool only_control,
   }
 
   if (cont_set && meshblock_data.Get()->GetVarPtr(label)->IsSparse()) {
-    auto clabel = label;
-    if (!only_control) clabel = pmy_mesh->resolved_packages->GetFieldController(label);
-    const auto &var_labels = pmy_mesh->resolved_packages->GetControlledVariables(clabel);
-    for (const auto &l : var_labels)
-      AllocateVar(l);
+    if (only_control) {
+      const auto &var_labels = pmy_mesh->resolved_packages->GetControlledVariables(label);
+      for (const auto &l : var_labels)
+        AllocateVar(l);
+    } else {
+      const auto &control_group =
+          pmy_mesh->resolved_packages->GetFieldControlGroup(label);
+      const auto &var_labels = pmy_mesh->resolved_packages->GetControlledVariables(
+          control_group.begin()->label());
+      for (const auto &l : var_labels)
+        AllocateVar(l);
+    }
   } else {
     AllocateVar(label);
   }
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index 6ec3d04ab967a..a20770260d99e 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -241,7 +241,7 @@ class MeshBlock : public std::enable_shared_from_this<MeshBlock> {
   void AllocateSparse(std::string const &label, bool only_control = false,
                       bool flag_uninitialized = false);
 
-  void AllocSparseID(std::string const &base_name, const int sparse_id) {
+  void AllocSparseID(std::string const &base_name, const SparseID sparse_id) {
     AllocateSparse(MakeVarLabel(base_name, sparse_id));
   }
 
@@ -252,7 +252,8 @@ class MeshBlock : public std::enable_shared_from_this<MeshBlock> {
     return meshblock_data.Get()->IsAllocated(label);
   }
 
-  inline bool IsAllocated(std::string const &base_name, int sparse_id) const noexcept {
+  inline bool IsAllocated(std::string const &base_name,
+                          SparseID sparse_id) const noexcept {
     return IsAllocated(MakeVarLabel(base_name, sparse_id));
   }
 #else
@@ -261,7 +262,7 @@ class MeshBlock : public std::enable_shared_from_this<MeshBlock> {
   }
 
   inline constexpr bool IsAllocated(std::string const & /*base_name*/,
-                                    int /*sparse_id*/) const noexcept {
+                                    SparseID /*sparse_id*/) const noexcept {
     return true;
   }
 #endif
diff --git a/src/sparse/sparse_management.cpp b/src/sparse/sparse_management.cpp
index 5a67273a2958d..f0b0d6a478938 100644
--- a/src/sparse/sparse_management.cpp
+++ b/src/sparse/sparse_management.cpp
@@ -11,6 +11,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <string>
 #include <unordered_set>
 
@@ -139,27 +141,49 @@ void SparseDeallocOnCount(T *rc, std::size_t count,
   for (int b = 0; b < pack.GetNBlocks(); ++b) {
     auto pmbdata = rc->GetBlockDataRawPointer(b);
     auto pmb = pmbdata->GetBlockPointer();
-    for (auto &control_var : control_vars) {
-      if (exclude.count(control_var) > 0) continue;
-
-      int lo = pack.GetLowerBoundHost(b, PackIdx(packIdx[control_var]));
-      int hi = pack.GetUpperBoundHost(b, PackIdx(packIdx[control_var]));
-      if (lo <= hi) { // Check that this control variable is actually in the pack
-        auto &counter = pmbdata->Get(control_var).dealloc_count;
-        bool all_zero = true;
-        for (int iv = lo; iv <= hi; ++iv)
-          all_zero = all_zero && is_zero_h(b, iv);
-        if (all_zero) {
-          counter++;
+    // Per group, update each member's deallocation counter using the old logic.
+    // If every member in the group is ready, deallocate the whole group together.
+    for (const auto &control_group :
+         rc->GetMeshPointer()->resolved_packages->GetControlGroups()) {
+      PARTHENON_REQUIRE_THROWS(!control_group.empty(),
+                               "Encountered an empty sparse control group");
+      const auto &representative = *control_group.begin();
+
+      bool group_excluded = false;
+      bool group_has_allocated_sparse_member = false;
+      bool deallocate_group = true;
+      for (const auto &control_var : control_group) {
+        if (exclude.count(control_var.label()) > 0) {
+          group_excluded = true;
         } else {
-          counter = 0;
+          const auto pack_idx = packIdx.find(control_var.label());
+          if (pack_idx != packIdx.end()) {
+            int lo = pack.GetLowerBoundHost(b, PackIdx(pack_idx->second));
+            int hi = pack.GetUpperBoundHost(b, PackIdx(pack_idx->second));
+            if (lo <= hi) {
+              group_has_allocated_sparse_member = true;
+
+              auto &counter = pmbdata->Get(control_var.label()).dealloc_count;
+              bool all_zero = true;
+              for (int iv = lo; iv <= hi; ++iv) {
+                all_zero = all_zero && is_zero_h(b, iv);
+              }
+              if (all_zero) {
+                counter++;
+              } else {
+                counter = 0;
+              }
+              deallocate_group = deallocate_group && (counter > count);
+            }
+          }
         }
-        if (counter > count) {
-          // this variable has been flagged for deallocation deallocation_count times in
-          // a row, now deallocate it
-          counter = 0;
-          pmb->DeallocateSparse(control_var);
+      }
+
+      if (!group_excluded && group_has_allocated_sparse_member && deallocate_group) {
+        for (const auto &control_var : control_group) {
+          pmbdata->Get(control_var.label()).dealloc_count = 0;
         }
+        pmb->DeallocateSparse(representative.label());
       }
     }
   }
diff --git a/tst/regression/test_suites/restart/restart.py b/tst/regression/test_suites/restart/restart.py
index 1cb4132b17fcb..bd883134892fd 100644
--- a/tst/regression/test_suites/restart/restart.py
+++ b/tst/regression/test_suites/restart/restart.py
@@ -19,7 +19,6 @@
 import sys
 import utils.test_case
 
-
 # To prevent littering up imported folders with .pyc files or __pycache_ folder
 sys.dont_write_bytecode = True
 
diff --git a/tst/regression/test_suites/restart_fine/restart_fine.py b/tst/regression/test_suites/restart_fine/restart_fine.py
index e15452a62035a..2ffb9bf1a9d78 100644
--- a/tst/regression/test_suites/restart_fine/restart_fine.py
+++ b/tst/regression/test_suites/restart_fine/restart_fine.py
@@ -19,7 +19,6 @@
 import sys
 import utils.test_case
 
-
 # To prevent littering up imported folders with .pyc files or __pycache_ folder
 sys.dont_write_bytecode = True
 
diff --git a/tst/style/cpplint.py b/tst/style/cpplint.py
index 0ceed165b7c3f..6179b1c6c474b 100755
--- a/tst/style/cpplint.py
+++ b/tst/style/cpplint.py
@@ -2337,7 +2337,7 @@ def CloseExpression(clean_lines, linenum, pos):
         return (line, clean_lines.NumLines(), -1)
 
     # Check first line
-    (end_pos, stack) = FindEndOfExpressionInLine(line, pos, [])
+    end_pos, stack = FindEndOfExpressionInLine(line, pos, [])
     if end_pos > -1:
         return (line, linenum, end_pos)
 
@@ -2345,7 +2345,7 @@ def CloseExpression(clean_lines, linenum, pos):
     while stack and linenum < clean_lines.NumLines() - 1:
         linenum += 1
         line = clean_lines.elided[linenum]
-        (end_pos, stack) = FindEndOfExpressionInLine(line, 0, stack)
+        end_pos, stack = FindEndOfExpressionInLine(line, 0, stack)
         if end_pos > -1:
             return (line, linenum, end_pos)
 
@@ -2455,7 +2455,7 @@ def ReverseCloseExpression(clean_lines, linenum, pos):
         return (line, 0, -1)
 
     # Check last line
-    (start_pos, stack) = FindStartOfExpressionInLine(line, pos, [])
+    start_pos, stack = FindStartOfExpressionInLine(line, pos, [])
     if start_pos > -1:
         return (line, linenum, start_pos)
 
@@ -2463,7 +2463,7 @@ def ReverseCloseExpression(clean_lines, linenum, pos):
     while stack and linenum > 0:
         linenum -= 1
         line = clean_lines.elided[linenum]
-        (start_pos, stack) = FindStartOfExpressionInLine(line, len(line) - 1, stack)
+        start_pos, stack = FindStartOfExpressionInLine(line, len(line) - 1, stack)
         if start_pos > -1:
             return (line, linenum, start_pos)
 
@@ -2517,7 +2517,7 @@ def PathSplitToList(path):
     """
     lst = []
     while True:
-        (head, tail) = os.path.split(path)
+        head, tail = os.path.split(path)
         if head == path:  # absolute paths end
             lst.append(head)
             break
@@ -3364,7 +3364,7 @@ def InTemplateArgumentList(self, clean_lines, linenum, pos):
 
             # We can't be sure if we just find a single '<', and need to
             # find the matching '>'.
-            (_, end_line, end_pos) = CloseExpression(clean_lines, linenum, pos - 1)
+            _, end_line, end_pos = CloseExpression(clean_lines, linenum, pos - 1)
             if end_pos < 0:
                 # Not sure if template argument list or syntax error in file
                 return False
@@ -4362,7 +4362,7 @@ def CheckOperatorSpacing(filename, clean_lines, linenum, error):
         # space.  This is done to avoid some false positives with shifts.
         match = re.match(r"^(.*[^\s<])<[^\s=<,]", line)
         if match:
-            (_, _, end_pos) = CloseExpression(clean_lines, linenum, len(match.group(1)))
+            _, _, end_pos = CloseExpression(clean_lines, linenum, len(match.group(1)))
             if end_pos <= -1:
                 error(
                     filename,
@@ -4377,7 +4377,7 @@ def CheckOperatorSpacing(filename, clean_lines, linenum, error):
         # false positives with shifts.
         match = re.match(r"^(.*[^-\s>])>[^\s=>,]", line)
         if match:
-            (_, _, start_pos) = ReverseCloseExpression(
+            _, _, start_pos = ReverseCloseExpression(
                 clean_lines, linenum, len(match.group(1))
             )
             if start_pos <= -1:
@@ -4637,7 +4637,7 @@ def CheckBracesSpacing(filename, clean_lines, linenum, nesting_state, error):
         # spurious semicolons, e.g. "if (cond){};", but we will catch the
         # spurious semicolon with a separate check.
         leading_text = match.group(1)
-        (endline, endlinenum, endpos) = CloseExpression(
+        endline, endlinenum, endpos = CloseExpression(
             clean_lines, linenum, len(match.group(1))
         )
         trailing_text = ""
@@ -4700,7 +4700,7 @@ def IsDecltype(clean_lines, linenum, column):
     Returns:
       True if this token is decltype() expression, False otherwise.
     """
-    (text, _, start_col) = ReverseCloseExpression(clean_lines, linenum, column)
+    text, _, start_col = ReverseCloseExpression(clean_lines, linenum, column)
     if start_col < 0:
         return False
     if re.search(r"\bdecltype\s*$", text[0:start_col]):
@@ -4853,7 +4853,7 @@ def CheckBraces(filename, clean_lines, linenum, error):
         pos = line.find("else if")
         pos = line.find("(", pos)
         if pos > 0:
-            (endline, _, endpos) = CloseExpression(clean_lines, linenum, pos)
+            endline, _, endpos = CloseExpression(clean_lines, linenum, pos)
             brace_on_right = endline[endpos:].find("{") != -1
             if brace_on_left != brace_on_right:  # must be brace after if
                 error(
@@ -4923,7 +4923,7 @@ def CheckBraces(filename, clean_lines, linenum, error):
         if if_match:
             # This could be a multiline if condition, so find the end first.
             pos = if_match.end() - 1
-            (endline, endlinenum, endpos) = CloseExpression(clean_lines, linenum, pos)
+            endline, endlinenum, endpos = CloseExpression(clean_lines, linenum, pos)
         # Check for an opening brace, either directly after the if or on the next
         # line. If found, this isn't a single-statement conditional.
         if not re.match(
@@ -5128,7 +5128,7 @@ def CheckTrailingSemicolon(filename, clean_lines, linenum, error):
 
     # Check matching closing brace
     if match:
-        (endline, endlinenum, endpos) = CloseExpression(
+        endline, endlinenum, endpos = CloseExpression(
             clean_lines, linenum, len(match.group(1))
         )
         if endpos > -1 and re.match(r"^\s*;", endline[endpos:]):
@@ -5176,7 +5176,7 @@ def CheckEmptyBlockBody(filename, clean_lines, linenum, error):
     matched = re.match(r"\s*(for|while|if)\s*\(", line)
     if matched:
         # Find the end of the conditional expression.
-        (end_line, end_linenum, end_pos) = CloseExpression(
+        end_line, end_linenum, end_pos = CloseExpression(
             clean_lines, linenum, line.find("(")
         )
 
@@ -5226,7 +5226,7 @@ def CheckEmptyBlockBody(filename, clean_lines, linenum, error):
             if opening_linenum == end_linenum:
                 # We need to make opening_pos relative to the start of the entire line.
                 opening_pos += end_pos
-            (closing_line, closing_linenum, closing_pos) = CloseExpression(
+            closing_line, closing_linenum, closing_pos = CloseExpression(
                 clean_lines, opening_linenum, opening_pos
             )
             if closing_pos < 0:
@@ -5315,12 +5315,12 @@ def CheckCheck(filename, clean_lines, linenum, error):
 
     # Decide the set of replacement macros that should be suggested
     lines = clean_lines.elided
-    (check_macro, start_pos) = FindCheckMacro(lines[linenum])
+    check_macro, start_pos = FindCheckMacro(lines[linenum])
     if not check_macro:
         return
 
     # Find end of the boolean expression by matching parentheses
-    (last_line, end_line, end_pos) = CloseExpression(clean_lines, linenum, start_pos)
+    last_line, end_line, end_pos = CloseExpression(clean_lines, linenum, start_pos)
     if end_pos < 0:
         return
 
@@ -5354,7 +5354,7 @@ def CheckCheck(filename, clean_lines, linenum, error):
             if token == "(":
                 # Parenthesized operand
                 expression = matched.group(2)
-                (end, _) = FindEndOfExpressionInLine(expression, 0, ["("])
+                end, _ = FindEndOfExpressionInLine(expression, 0, ["("])
                 if end < 0:
                     return  # Unmatched parenthesis
                 lhs += "(" + expression[0:end]
@@ -6499,7 +6499,7 @@ def CheckForNonConstReference(filename, clean_lines, linenum, nesting_state, err
             # Check for templated parameter that is split across multiple lines
             endpos = line.rfind(">")
             if endpos > -1:
-                (_, startline, startpos) = ReverseCloseExpression(
+                _, startline, startpos = ReverseCloseExpression(
                     clean_lines, linenum, endpos
                 )
                 if startpos > -1 and startline < linenum:
@@ -7261,7 +7261,7 @@ def CheckRedundantVirtual(filename, clean_lines, linenum, error):
         parameter_list = re.match(r"^([^(]*)\(", line)
         if parameter_list:
             # Match parentheses to find the end of the parameter list
-            (_, end_line, end_col) = CloseExpression(
+            _, end_line, end_col = CloseExpression(
                 clean_lines, start_line, start_col + len(parameter_list.group(1))
             )
             break
@@ -7815,7 +7815,7 @@ def ParseArguments(args):
       The list of filenames to lint.
     """
     try:
-        (opts, filenames) = getopt.getopt(
+        opts, filenames = getopt.getopt(
             args,
             "",
             [
diff --git a/tst/unit/CMakeLists.txt b/tst/unit/CMakeLists.txt
index af6ab9c70e8aa..e27425b2f538b 100644
--- a/tst/unit/CMakeLists.txt
+++ b/tst/unit/CMakeLists.txt
@@ -39,6 +39,7 @@ list(APPEND unit_tests_SOURCES
     test_mesh_data.cpp
     test_output_utils.cpp
     test_pararrays.cpp
+    test_sparse_id.cpp
     test_sparse_pack.cpp
     test_parameter_input.cpp
     test_error_checking.cpp
diff --git a/tst/unit/test_meshblock_data_iterator.cpp b/tst/unit/test_meshblock_data_iterator.cpp
index d7b1ee709eeda..bfb0f320f6c30 100644
--- a/tst/unit/test_meshblock_data_iterator.cpp
+++ b/tst/unit/test_meshblock_data_iterator.cpp
@@ -14,6 +14,9 @@
 // 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.
 //========================================================================================
+
+// This file was made in part with generative AI.
+
 #include <array>
 #include <cmath>
 #include <iostream>
@@ -32,6 +35,7 @@
 #include "interface/meshblock_data.hpp"
 #include "interface/metadata.hpp"
 #include "interface/state_descriptor.hpp"
+#include "interface/var_id.hpp"
 #include "interface/variable.hpp"
 #include "interface/variable_pack.hpp"
 #include "kokkos_abstraction.hpp"
@@ -51,6 +55,7 @@ using parthenon::par_for;
 using parthenon::ParArray4D;
 using parthenon::ParArrayND;
 using parthenon::Real;
+using parthenon::SparseID;
 using parthenon::StateDescriptor;
 using parthenon::Variable;
 using parthenon::VariableVector;
@@ -357,8 +362,8 @@ TEST_CASE("Can pull variables from containers based on Metadata",
       mbd.Initialize(pkg, dummy_mb);
 
       // TODO(JL) test packs with unallocated sparse fields
-      dummy_mb->AllocSparseID("vsparse", 1);
-      dummy_mb->AllocSparseID("vsparse", 13);
+      dummy_mb->AllocSparseID("vsparse", SparseID::Scalar(1));
+      dummy_mb->AllocSparseID("vsparse", SparseID::Scalar(13));
       dummy_mb->AllocateSparse("vsparse_42");
 
       THEN("the low and high index bounds are correct as returned by PackVariables") {
@@ -399,8 +404,8 @@ TEST_CASE("Can pull variables from containers based on Metadata",
         Kokkos::parallel_reduce(
             "add correct checks", 1,
             KOKKOS_LAMBDA(const int i, int &sum) {
-              sum = (v.GetSparseID(v3first) == parthenon::InvalidSparseID);
-              sum += (v.GetSparseID(v6first) == parthenon::InvalidSparseID);
+              sum = (v.GetSparseID(v3first) == parthenon::InvalidSparseID());
+              sum += (v.GetSparseID(v6first) == parthenon::InvalidSparseID());
               sum += (v.GetSparseID(vs1) == 1);
               sum += (v.GetSparseID(vs13) == 13);
               sum += (v.GetSparseID(vs42) == 42);
diff --git a/tst/unit/test_sparse_id.cpp b/tst/unit/test_sparse_id.cpp
new file mode 100644
index 0000000000000..99cb84015c421
--- /dev/null
+++ b/tst/unit/test_sparse_id.cpp
@@ -0,0 +1,90 @@
+//========================================================================================
+// Parthenon performance portable AMR framework
+// Copyright(C) 2020-2024 The Parthenon collaboration
+// Licensed under the 3-clause BSD License, see LICENSE file for details
+//========================================================================================
+// (C) (or copyright) 2026. 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.
+//========================================================================================
+
+// This file was made in part with generative AI.
+
+#include <algorithm>
+#include <map>
+#include <string>
+#include <unordered_set>
+#include <vector>
+
+#include <catch2/catch.hpp>
+
+#include "interface/var_id.hpp"
+
+using parthenon::InvalidSparseID;
+using parthenon::MakeVarLabel;
+using parthenon::SparseID;
+using parthenon::SparseIDHasher;
+
+TEST_CASE("SparseID basics", "[SparseID]") {
+  GIVEN("Scalar, pair, and invalid sparse ids") {
+    const auto scalar = SparseID::Scalar(7);
+    const auto pair = SparseID::Pair(7, 11);
+    const auto invalid = InvalidSparseID;
+
+    THEN("Scalar and pair accessors behave as expected") {
+      REQUIRE(scalar() == 7);
+      REQUIRE(scalar(0) == 7);
+      REQUIRE(scalar(1) == parthenon::InvalidSparseIDValue);
+      REQUIRE(pair(0) == 7);
+      REQUIRE(pair(1) == 11);
+    }
+
+    THEN("Validity is based on the first component") {
+      REQUIRE(parthenon::IsValidSparseID(scalar));
+      REQUIRE(parthenon::IsValidSparseID(pair));
+      REQUIRE_FALSE(parthenon::IsValidSparseID(invalid));
+    }
+
+    THEN("Labels preserve the current scalar format") {
+      REQUIRE(MakeVarLabel("foo", invalid) == "foo");
+      REQUIRE(MakeVarLabel("foo", scalar) == "foo_7");
+      REQUIRE(MakeVarLabel("foo", pair) == "foo_7_11");
+    }
+
+    THEN("Ordering is lexicographic in the two components") {
+      std::vector<SparseID> ids = {pair, SparseID::Pair(2, 0), scalar,
+                                   SparseID::Scalar(3), SparseID::Pair(7, 1)};
+      std::sort(ids.begin(), ids.end());
+
+      std::vector<SparseID> expected = {SparseID::Pair(2, 0), SparseID::Scalar(3), scalar,
+                                        SparseID::Pair(7, 1), pair};
+      REQUIRE(ids == expected);
+    }
+
+    THEN("Hashing and equality work in associative containers") {
+      std::unordered_set<SparseID, SparseIDHasher> ids;
+      ids.insert(scalar);
+      ids.insert(pair);
+      ids.insert(SparseID::Scalar(7));
+      ids.insert(SparseID::Pair(7, 11));
+
+      REQUIRE(ids.size() == 2);
+      REQUIRE(ids.count(SparseID::Scalar(7)) == 1);
+      REQUIRE(ids.count(SparseID::Pair(7, 11)) == 1);
+
+      std::map<SparseID, std::string> labels;
+      labels[scalar] = "scalar";
+      labels[pair] = "pair";
+
+      REQUIRE(labels.at(SparseID::Scalar(7)) == "scalar");
+      REQUIRE(labels.at(SparseID::Pair(7, 11)) == "pair");
+    }
+  }
+}
diff --git a/tst/unit/test_state_descriptor.cpp b/tst/unit/test_state_descriptor.cpp
index f3336b03992d7..d913730641cb5 100644
--- a/tst/unit/test_state_descriptor.cpp
+++ b/tst/unit/test_state_descriptor.cpp
@@ -15,6 +15,8 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
+// This file was made in part with generative AI.
+
 #include <memory>
 #include <set>
 #include <string>
@@ -38,10 +40,12 @@ using parthenon::IndexRange;
 using parthenon::Metadata;
 using parthenon::MetadataFlag;
 using FC_t = parthenon::Metadata::FlagCollection;
+using parthenon::ControlGroup;
 using parthenon::Packages_t;
 using parthenon::ParArrayND;
 using parthenon::Real;
 using parthenon::ResolvePackages;
+using parthenon::SparseID;
 using parthenon::SparsePool;
 using parthenon::StateDescriptor;
 using FlagVec = std::vector<MetadataFlag>;
@@ -240,8 +244,10 @@ TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]")
           REQUIRE(!(pkg3->SparseBaseNamePresent("sparse")));
         }
         AND_THEN("The appropriate sparse metadata was added") {
-          REQUIRE(pkg3->FieldPresent("package1::sparse", sparse_ids[2]));
-          REQUIRE(pkg3->FieldPresent("package2::sparse", sparse_ids[3]));
+          REQUIRE(
+              pkg3->FieldPresent("package1::sparse", SparseID::Scalar(sparse_ids[2])));
+          REQUIRE(
+              pkg3->FieldPresent("package2::sparse", SparseID::Scalar(sparse_ids[3])));
         }
       }
     }
@@ -252,14 +258,15 @@ TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]")
         auto pkg4 = ResolvePackages(packages);
         AND_THEN("The sparse variable is present") {
           for (int i = 0; i < sparse_ids.size(); i++) {
-            REQUIRE(pkg4->FieldMetadata("sparse", sparse_ids[i]) == (m_sparse_provides));
+            REQUIRE(pkg4->FieldMetadata("sparse", SparseID::Scalar(sparse_ids[i])) ==
+                    (m_sparse_provides));
           }
         }
         AND_THEN("The sparse ids in the sparse pool are sorted") {
           auto &pool = (pkg4->GetSparsePool("sparse")).pool();
           std::vector<int> local_ids;
           for (auto &[id, m] : pool) {
-            local_ids.push_back(id);
+            local_ids.push_back(id());
           }
           REQUIRE(std::is_sorted(local_ids.begin(), local_ids.end()));
         }
@@ -333,7 +340,7 @@ TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]")
         AND_THEN("The overridable variables are retained") {
           REQUIRE(pkg4->FieldPresent("dense"));
           for (const int sid : sparse_ids) {
-            REQUIRE(pkg4->FieldPresent("sparse", sid));
+            REQUIRE(pkg4->FieldPresent("sparse", SparseID::Scalar(sid)));
           }
           REQUIRE(pkg4->SwarmPresent("myswarm"));
           REQUIRE(pkg4->SwarmValuePresent("value1", "myswarm"));
@@ -362,9 +369,9 @@ TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]")
       pkg2->AddSparsePool("sparse_c", m_sparse_provides, "sparse_b", sparse_ids);
 
       for (const int sid : sparse_ids) {
-        REQUIRE(pkg1->FieldPresent("sparse", sid));
-        REQUIRE(pkg2->FieldPresent("sparse", sid));
-        REQUIRE(pkg3->FieldPresent("sparse", sid));
+        REQUIRE(pkg1->FieldPresent("sparse", SparseID::Scalar(sid)));
+        REQUIRE(pkg2->FieldPresent("sparse", SparseID::Scalar(sid)));
+        REQUIRE(pkg3->FieldPresent("sparse", SparseID::Scalar(sid)));
       }
 
       THEN("We can safely resolve conflicts") {
@@ -388,10 +395,11 @@ TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]")
           REQUIRE(!(pkg4->SwarmValuePresent("overridable", "myswarm")));
           REQUIRE(pkg4->SparseBaseNamePresent("sparse"));
           for (const int sid : sparse_ids) {
-            REQUIRE(pkg4->FieldPresent("sparse", sid));
+            REQUIRE(pkg4->FieldPresent("sparse", SparseID::Scalar(sid)));
           }
           for (const int sid : sparse_ids) {
-            REQUIRE(pkg4->FieldMetadata("sparse", sid) == m_sparse_provides);
+            REQUIRE(pkg4->FieldMetadata("sparse", SparseID::Scalar(sid)) ==
+                    m_sparse_provides);
           }
         }
         AND_THEN("The correct sparse allocation control is resolved.") {
@@ -440,9 +448,8 @@ TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]")
                     (pkg3->RefinementFuncID(cell_funcs)));
             REQUIRE(pkg3->FieldMetadata("dense").GetRefinementFunctions() == cell_funcs);
             for (int i = 0; i < sparse_ids.size(); i++) {
-              REQUIRE(
-                  pkg3->FieldMetadata("sparse", sparse_ids[i]).GetRefinementFunctions() ==
-                  my_funcs);
+              REQUIRE(pkg3->FieldMetadata("sparse", SparseID::Scalar(sparse_ids[i]))
+                          .GetRefinementFunctions() == my_funcs);
             }
           }
         }
@@ -469,14 +476,14 @@ TEST_CASE("Test SparsePool interface", "[StateDescriptor]") {
     THEN("We can create a SparsePool with sparse metadata") {
       SparsePool pool("sparse", sparse_vec);
       AND_THEN("We can add sparse indices to the pool") {
-        const auto m2 = pool.Add(2);
+        const auto m2 = pool.Add(SparseID::Scalar(2));
         REQUIRE(m2 == sparse_vec);
         REQUIRE(m2.IsSet(Metadata::Vector));
         REQUIRE(!m2.IsSet(Metadata::Tensor));
 
         const int sparse_id = 5;
         const std::vector<int> shape = {2, 2, 4};
-        const auto m5 = pool.Add(sparse_id, shape, Metadata::Tensor);
+        const auto m5 = pool.Add(SparseID::Scalar(sparse_id), shape, Metadata::Tensor);
 
         const std::set<MetadataFlag> expected_flags{
             Metadata::Independent, Metadata::Sparse,   Metadata::Tensor,
@@ -490,12 +497,14 @@ TEST_CASE("Test SparsePool interface", "[StateDescriptor]") {
         REQUIRE(!m5.IsSet(Metadata::Vector));
         REQUIRE(m5.IsSet(Metadata::Tensor));
 
-        const auto mm17 = pool.Add(-17, {1}, Metadata::None, {"foo"});
+        const auto mm17 = pool.Add(SparseID::Scalar(-17), {1}, Metadata::None, {"foo"});
         REQUIRE(!mm17.IsSet(Metadata::Vector));
         REQUIRE(!mm17.IsSet(Metadata::Tensor));
         REQUIRE(mm17.getComponentLabels() == std::vector<std::string>{"foo"});
 
-        AND_THEN("We can't add the same sparse ID twice") { REQUIRE_THROWS(pool.Add(2)); }
+        AND_THEN("We can't add the same sparse ID twice") {
+          REQUIRE_THROWS(pool.Add(SparseID::Scalar(2)));
+        }
       }
 
       AND_THEN("We can't add InvalidSparseID") {
@@ -514,8 +523,8 @@ TEST_CASE("Test SparsePool interface", "[StateDescriptor]") {
 
     THEN("We can add sparse pools in different ways") {
       SparsePool pool1("pool1", meta_sparse);
-      pool1.Add(0);
-      pool1.Add(55);
+      pool1.Add(SparseID::Scalar(0));
+      pool1.Add(SparseID::Scalar(55));
       REQUIRE(pkg->AddSparsePool(pool1));
 
       const std::vector<int> sparse_ids_2{1, 55, 100};
@@ -529,13 +538,13 @@ TEST_CASE("Test SparsePool interface", "[StateDescriptor]") {
           "pool3", meta_sparse, sparse_ids_3, shapes,
           std::vector<MetadataFlag>{Metadata::Vector, Metadata::Tensor}));
 
-      REQUIRE(pkg->FieldPresent("pool1", 0));
-      REQUIRE(pkg->FieldPresent("pool1", 55));
-      REQUIRE(pkg->FieldPresent("pool2", 1));
-      REQUIRE(pkg->FieldPresent("pool2", 55));
-      REQUIRE(pkg->FieldPresent("pool2", 100));
-      REQUIRE(pkg->FieldPresent("pool3", 0));
-      REQUIRE(pkg->FieldPresent("pool3", 100));
+      REQUIRE(pkg->FieldPresent("pool1", SparseID::Scalar(0)));
+      REQUIRE(pkg->FieldPresent("pool1", SparseID::Scalar(55)));
+      REQUIRE(pkg->FieldPresent("pool2", SparseID::Scalar(1)));
+      REQUIRE(pkg->FieldPresent("pool2", SparseID::Scalar(55)));
+      REQUIRE(pkg->FieldPresent("pool2", SparseID::Scalar(100)));
+      REQUIRE(pkg->FieldPresent("pool3", SparseID::Scalar(0)));
+      REQUIRE(pkg->FieldPresent("pool3", SparseID::Scalar(100)));
 
       AND_THEN("We can't add a SparsePool with wrong number of Vector/Tensor flags") {
         REQUIRE_THROWS(pkg->AddSparsePool(
@@ -570,6 +579,60 @@ TEST_CASE("Test SparsePool interface", "[StateDescriptor]") {
           pkg->AddSparsePool("fake2_sparse", meta_sparse, std::vector<int>{13, 27, 9}));
     }
   }
+
+  GIVEN("Sparse pools that group control by the first sparse ID component") {
+    Packages_t packages;
+    auto pkg = std::make_shared<StateDescriptor>("pkg");
+    packages.Add(pkg);
+    Metadata meta_sparse({Metadata::Sparse});
+
+    SparsePool rho("rho", meta_sparse);
+    rho.SetControlSparseIDMode(SparsePool::ControlSparseIDMode::FirstComponent);
+    rho.Add(SparseID::Pair(0, 0));
+    rho.Add(SparseID::Pair(0, 1));
+    rho.Add(SparseID::Pair(0, 2));
+    rho.Add(SparseID::Pair(1, 0));
+    REQUIRE(pkg->AddSparsePool(rho));
+
+    SparsePool other("other", meta_sparse, "rho");
+    other.SetControlSparseIDMode(SparsePool::ControlSparseIDMode::FirstComponent);
+    other.Add(SparseID::Pair(0, 0));
+    other.Add(SparseID::Pair(0, 1));
+    other.Add(SparseID::Pair(0, 2));
+    other.Add(SparseID::Pair(1, 0));
+    REQUIRE(pkg->AddSparsePool(other));
+
+    auto resolved = ResolvePackages(packages);
+
+    THEN("The first sparse ID component is used to build controller groups") {
+      REQUIRE(resolved->GetFieldControlGroup("rho_0_1") ==
+              ControlGroup{parthenon::VarID("rho", SparseID::Pair(0, 0)),
+                           parthenon::VarID("rho", SparseID::Pair(0, 1)),
+                           parthenon::VarID("rho", SparseID::Pair(0, 2))});
+      REQUIRE(resolved->GetFieldControlGroup("other_0_2") ==
+              ControlGroup{parthenon::VarID("rho", SparseID::Pair(0, 0)),
+                           parthenon::VarID("rho", SparseID::Pair(0, 1)),
+                           parthenon::VarID("rho", SparseID::Pair(0, 2))});
+      REQUIRE(resolved->GetFieldControlGroup("rho_1_0") ==
+              ControlGroup{parthenon::VarID("rho", SparseID::Pair(1, 0))});
+    }
+
+    AND_THEN("Legacy controller access remains valid for singleton groups") {
+      REQUIRE(resolved->GetFieldController("rho_1_0") == "rho_1_0");
+      REQUIRE_THROWS(resolved->GetFieldController("rho_0_1"));
+    }
+
+    AND_THEN("All sparse fields with the same first component share a controlled set") {
+      auto controlled_0 = resolved->GetControlledVariables("rho_0_0");
+      REQUIRE(std::count(controlled_0.begin(), controlled_0.end(), "rho_0_0") == 1);
+      REQUIRE(std::count(controlled_0.begin(), controlled_0.end(), "rho_0_1") == 1);
+      REQUIRE(std::count(controlled_0.begin(), controlled_0.end(), "rho_0_2") == 1);
+      REQUIRE(std::count(controlled_0.begin(), controlled_0.end(), "other_0_0") == 1);
+      REQUIRE(std::count(controlled_0.begin(), controlled_0.end(), "other_0_1") == 1);
+      REQUIRE(std::count(controlled_0.begin(), controlled_0.end(), "other_0_2") == 1);
+      REQUIRE(controlled_0.size() == 6);
+    }
+  }
 }
 
 TEST_CASE("Test getting a vector of variable names given criteria", "[StateDescriptor]") {