From b4ea497474766c8897726aef10532836f5272cba Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Mon, 29 Jul 2024 16:30:51 -0600
Subject: [PATCH 01/40] make tensor indices accessible in packs

---
 src/interface/sparse_pack_base.cpp | 7 +++++++
 src/interface/variable_state.hpp   | 5 ++++-
 2 files changed, 11 insertions(+), 1 deletion(-)

diff --git a/src/interface/sparse_pack_base.cpp b/src/interface/sparse_pack_base.cpp
index 5c387973804e..dc0b36fc860c 100644
--- a/src/interface/sparse_pack_base.cpp
+++ b/src/interface/sparse_pack_base.cpp
@@ -249,6 +249,13 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
                   PARTHENON_REQUIRE(
                       pack.pack_h_(0, b, idx).size() > 0,
                       "Seems like this variable might not actually be allocated.");
+                  
+                  // Save the tensor indices for easy access 
+                  for (int l = 0; l < leading_dim; ++l) { 
+                    pack.pack_h_(l, b, idx).t = t;
+                    pack.pack_h_(l, b, idx).u = u;
+                    pack.pack_h_(l, b, idx).v = v;
+                  }
 
                   if (desc.flat) {
                     coords_h(idx) = pmbd->GetBlockPointer()->coords_device;
diff --git a/src/interface/variable_state.hpp b/src/interface/variable_state.hpp
index d8a2ea1ba746..f25b39d423d5 100644
--- a/src/interface/variable_state.hpp
+++ b/src/interface/variable_state.hpp
@@ -61,7 +61,10 @@ struct VariableState : public empty_state_t {
   int sparse_id;
   int vector_component = NODIR;
   bool initialized = true;
-
+  
+  int t{0};
+  int u{0};
+  int v{0};
   TopologicalType topological_type = TopologicalType::Cell;
   TopologicalElement topological_element = TopologicalElement::CC;
   std::size_t tensor_components;

From eac315a631001f661b2c2da16cb6cea2598d841e Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Mon, 29 Jul 2024 16:31:01 -0600
Subject: [PATCH 02/40] add radiation constant

---
 src/utils/constants.hpp | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/utils/constants.hpp b/src/utils/constants.hpp
index 72b706b5576f..17386cf54e61 100644
--- a/src/utils/constants.hpp
+++ b/src/utils/constants.hpp
@@ -117,6 +117,10 @@ class PhysicalConstants {
   static constexpr double stefan_boltzmann = 2.0 * M_PI * M_PI * M_PI * M_PI * M_PI * kb *
                                              kb * kb * kb / (15.0 * h * h * h * c * c);
   static constexpr double sb = stefan_boltzmann;
+  
+  // Radiation constant 
+  static constexpr double radiation_constant = 4.0 * stefan_boltzmann / speed_of_light;
+  static constexpr double ar = radiation_constant;
 
   // Faraday constant
   static constexpr double faraday_constant = 96485.33645957 * capacitance;

From d21f39096ad9de47bf0587387a14b7b8044da0b2 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 30 Jul 2024 14:53:16 -0600
Subject: [PATCH 03/40] Add some nice stuff for type lists

---
 src/CMakeLists.txt                     |  1 +
 src/interface/make_pack_descriptor.hpp | 11 ++++
 src/interface/sparse_pack.hpp          |  1 +
 src/solvers/bicgstab_solver.hpp        | 27 +++++----
 src/solvers/mg_solver.hpp              | 19 +++++--
 src/utils/type_list.hpp                | 76 ++++++++++++++++++++++++++
 6 files changed, 118 insertions(+), 17 deletions(-)
 create mode 100644 src/utils/type_list.hpp

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 708d59b5ea65..7f66585a8e8c 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -270,6 +270,7 @@ add_library(parthenon
   utils/sort.hpp
   utils/string_utils.cpp
   utils/string_utils.hpp
+  utils/type_list.hpp
   utils/unique_id.cpp
   utils/unique_id.hpp
   utils/utils.hpp
diff --git a/src/interface/make_pack_descriptor.hpp b/src/interface/make_pack_descriptor.hpp
index f3cfe17ec4f3..7620a99a0c02 100644
--- a/src/interface/make_pack_descriptor.hpp
+++ b/src/interface/make_pack_descriptor.hpp
@@ -33,6 +33,7 @@
 #include "interface/state_descriptor.hpp"
 #include "interface/variable.hpp"
 #include "mesh/mesh.hpp"
+#include "utils/type_list.hpp"
 
 namespace parthenon {
 
@@ -136,6 +137,16 @@ inline auto MakePackDescriptor(StateDescriptor *psd, const std::vector<Uid_t> &v
   return typename SparsePack<>::Descriptor(base_desc);
 }
 
+template <class... Types, class... Args> 
+inline auto MakePackDescriptorFromTypeList(TypeList<Types...>, Args&&... args) { 
+  return MakePackDescriptor<Types...>(std::forward<Args>(args)...);
+}
+
+template <class TL, class... Args> 
+inline auto MakePackDescriptorFromTypeList(Args&&... args) { 
+  return MakePackDescriptorFromTypeList(TL(), std::forward<Args>(args)...);
+}
+
 } // namespace parthenon
 
 #endif // INTERFACE_MAKE_PACK_DESCRIPTOR_HPP_
diff --git a/src/interface/sparse_pack.hpp b/src/interface/sparse_pack.hpp
index b804839fa5a6..7ef7a6fe1f41 100644
--- a/src/interface/sparse_pack.hpp
+++ b/src/interface/sparse_pack.hpp
@@ -34,6 +34,7 @@
 #include "interface/variable.hpp"
 #include "utils/concepts_lite.hpp"
 #include "utils/utils.hpp"
+#include "utils/type_list.hpp"
 
 namespace parthenon {
 
diff --git a/src/solvers/bicgstab_solver.hpp b/src/solvers/bicgstab_solver.hpp
index 9f30de206478..9fa4d81d239a 100644
--- a/src/solvers/bicgstab_solver.hpp
+++ b/src/solvers/bicgstab_solver.hpp
@@ -24,8 +24,8 @@
 #include "kokkos_abstraction.hpp"
 #include "solvers/mg_solver.hpp"
 #include "solvers/solver_utils.hpp"
-
 #include "tasks/tasks.hpp"
+#include "utils/type_list.hpp"
 
 namespace parthenon {
 
@@ -69,21 +69,25 @@ class BiCGSTABSolver {
   PARTHENON_INTERNALSOLVERVARIABLE(u, r);
   PARTHENON_INTERNALSOLVERVARIABLE(u, p);
   PARTHENON_INTERNALSOLVERVARIABLE(u, x);
+  
+  using internal_types_tl = TypeList<rhat0, v, h, s, t, r, p, x>;
+  using preconditioner_t = MGSolver<u, rhs, equations>; 
+  using all_internal_types_tl = concatenate_type_lists_t<internal_types_tl, typename preconditioner_t::internal_types_tl>;
 
   std::vector<std::string> GetInternalVariableNames() const {
-    std::vector<std::string> names{rhat0::name(), v::name(), h::name(), s::name(),
-                                   t::name(),     r::name(), p::name(), x::name()};
+    std::vector<std::string> names;
     if (params_.precondition) {
-      auto pre_names = preconditioner.GetInternalVariableNames();
-      names.insert(names.end(), pre_names.begin(), pre_names.end());
+      all_internal_types_tl::IterateTypes([&names](auto t){names.push_back(decltype(t)::name());});
+    } else {
+      internal_types_tl::IterateTypes([&names](auto t){names.push_back(decltype(t)::name());});
     }
     return names;
   }
 
   BiCGSTABSolver(StateDescriptor *pkg, BiCGSTABParams params_in,
-                 equations eq_in = equations(), std::vector<int> shape = {})
-      : preconditioner(pkg, params_in.mg_params, eq_in, shape), params_(params_in),
-        iter_counter(0), eqs_(eq_in) {
+                 equations eq_in = equations(), std::vector<int> shape = {}, const std::string &container = "base")
+      : preconditioner(pkg, params_in.mg_params, eq_in, shape, container), params_(params_in),
+        iter_counter(0), eqs_(eq_in), container_(container) {
     using namespace refinement_ops;
     auto m_no_ghost =
         Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, shape);
@@ -105,8 +109,8 @@ class BiCGSTABSolver {
   TaskID AddTasks(TaskList &tl, TaskID dependence, Mesh *pmesh, const int partition) {
     using namespace utils;
     TaskID none;
-    auto &md = pmesh->mesh_data.GetOrAdd("base", partition);
-    std::string label = "bicg_comm_" + std::to_string(partition);
+    auto &md = pmesh->mesh_data.GetOrAdd(container_, partition);
+    std::string label = container_ + "bicg_comm_" + std::to_string(partition);
     auto &md_comm =
         pmesh->mesh_data.AddShallow(label, md, std::vector<std::string>{u::name()});
     iter_counter = 0;
@@ -313,7 +317,7 @@ class BiCGSTABSolver {
   BiCGSTABParams &GetParams() { return params_; }
 
  protected:
-  MGSolver<u, rhs, equations> preconditioner;
+  preconditioner_t preconditioner;
   BiCGSTABParams params_;
   int iter_counter;
   AllReduce<Real> rtr, pAp, rhat0v, rhat0r, ts, tt, residual, rhs2;
@@ -321,6 +325,7 @@ class BiCGSTABSolver {
   equations eqs_;
   Real final_residual;
   int final_iteration;
+  std::string container_;
 };
 
 } // namespace solvers
diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
index 2a8e464503dd..4c5b7c5b9343 100644
--- a/src/solvers/mg_solver.hpp
+++ b/src/solvers/mg_solver.hpp
@@ -27,6 +27,7 @@
 #include "solvers/solver_utils.hpp"
 #include "tasks/tasks.hpp"
 #include "utils/robust.hpp"
+#include "utils/type_list.hpp"
 
 namespace parthenon {
 
@@ -76,13 +77,17 @@ class MGSolver {
   PARTHENON_INTERNALSOLVERVARIABLE(u, temp); // Temporary storage
   PARTHENON_INTERNALSOLVERVARIABLE(u, u0);   // Storage for initial solution during FAS
   PARTHENON_INTERNALSOLVERVARIABLE(u, D);    // Storage for (approximate) diagonal
+  
+  using internal_types_tl = TypeList<res_err, temp, u0, D>; 
   std::vector<std::string> GetInternalVariableNames() const {
-    return {res_err::name(), temp::name(), u0::name(), D::name()};
+    std::vector<std::string> names; 
+    internal_types_tl::IterateTypes([&names](auto t){names.push_back(decltype(t)::name());});
+    return names;
   }
 
   MGSolver(StateDescriptor *pkg, MGParams params_in, equations eq_in = equations(),
-           std::vector<int> shape = {})
-      : params_(params_in), iter_counter(0), eqs_(eq_in) {
+           std::vector<int> shape = {}, const std::string& container = "base")
+      : params_(params_in), iter_counter(0), eqs_(eq_in), container_(container) {
     using namespace parthenon::refinement_ops;
     // The ghost cells of res_err need to be filled, but this is accomplished by
     // copying res_err into u, communicating, then copying u back into res_err
@@ -128,7 +133,7 @@ class MGSolver {
     auto partitions = pmesh->GetDefaultBlockPartitions(GridIdentifier::leaf());
     if (partition >= partitions.size())
       PARTHENON_FAIL("Does not work with non-default partitioning.");
-    auto &md = pmesh->mesh_data.Add("base", partitions[partition]);
+    auto &md = pmesh->mesh_data.Add(container_, partitions[partition]);
     auto comm = AddBoundaryExchangeTasks<BoundaryType::any>(mg_finest, itl, md,
                                                             pmesh->multilevel);
     auto calc_pointwise_res = eqs_.template Ax<u, res_err>(itl, comm, md);
@@ -207,6 +212,8 @@ class MGSolver {
   equations eqs_;
   Real final_residual;
   int final_iteration;
+  std::string container_;
+
   // These functions apparently have to be public to compile with cuda since
   // they contain device side lambdas
  public:
@@ -344,7 +351,7 @@ class MGSolver {
     auto partitions =
         pmesh->GetDefaultBlockPartitions(GridIdentifier::two_level_composite(level));
     if (partition >= partitions.size()) return dependence;
-    auto &md = pmesh->mesh_data.Add("base", partitions[partition]);
+    auto &md = pmesh->mesh_data.Add(container_, partitions[partition]);
 
     auto task_out = dependence;
     if (level < max_level) {
@@ -399,7 +406,7 @@ class MGSolver {
     auto partitions =
         pmesh->GetDefaultBlockPartitions(GridIdentifier::two_level_composite(level));
     if (partition >= partitions.size()) return dependence;
-    auto &md = pmesh->mesh_data.Add("base", partitions[partition]);
+    auto &md = pmesh->mesh_data.Add(container_, partitions[partition]);
     auto &md_comm = pmesh->mesh_data.AddShallow(
         "mg_comm", md, std::vector<std::string>{u::name(), res_err::name()});
 
diff --git a/src/utils/type_list.hpp b/src/utils/type_list.hpp
new file mode 100644
index 000000000000..ede6ebdc4dc7
--- /dev/null
+++ b/src/utils/type_list.hpp
@@ -0,0 +1,76 @@
+//========================================================================================
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
+//
+// This program was produced under U.S. Government contract 89233218CNA000001 for Los
+// Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
+// for the U.S. Department of Energy/National Nuclear Security Administration. All rights
+// in the program are reserved by Triad National Security, LLC, and the U.S. Department
+// of Energy/National Nuclear Security Administration. The Government is granted for
+// itself and others acting on its behalf a nonexclusive, paid-up, irrevocable worldwide
+// license in this material to reproduce, prepare derivative works, distribute copies to
+// the public, perform publicly and display publicly, and to permit others to do so.
+//========================================================================================
+
+#ifndef UTILS_TYPE_LIST_HPP_
+#define UTILS_TYPE_LIST_HPP_
+
+#include <tuple>
+#include <type_traits>
+#include <utility>
+
+namespace parthenon {
+template <class... Args> 
+struct TypeList {
+  using types = std::tuple<Args...>;
+  
+  static constexpr std::size_t n_types{sizeof...(Args)};
+
+  template <std::size_t I>
+  using type = typename std::tuple_element<I, types>::type;
+  
+  template <std::size_t... Idxs>
+  using sublist = TypeList<type<Idxs>...>;
+
+  template <class T, std::size_t I = 0>
+  static constexpr std::size_t GetIdx(){
+    static_assert(I < n_types, "Type is not present in TypeList.");
+    if constexpr (std::is_same_v<T, type<I>>) {
+      return I;
+    } else {
+      return GetIdx<T, I + 1>();   
+    }
+  } 
+  
+  template<class F>
+  static void IterateTypes(F func) {
+    (func(Args()), ...);
+  }
+
+ private:  
+  template<std::size_t Start, std::size_t End>
+  static auto ContinuousSublist() {
+    return ContinuousSublistImpl<Start>(std::make_index_sequence<End - Start + 1>());
+  }
+  template<std::size_t Start, std::size_t... Is> 
+  static auto ContinuousSublistImpl(std::index_sequence<Is...>) {
+    return sublist<(Start + Is)...>();
+  }
+
+ public:
+  template<std::size_t Start, std::size_t End>
+  using continuous_sublist = decltype(ContinuousSublist<Start, End>());
+};
+
+namespace impl {
+  template<class... Args>
+  auto ConcatenateTypeLists(TypeList<Args...>){return TypeList<Args...>();}
+  
+  template<class... Args1, class... Args2, class... Args> 
+  auto ConcatenateTypeLists(TypeList<Args1...>, TypeList<Args2...>, Args...){ return ConcatenateTypeLists(TypeList<Args1..., Args2...>(), Args()...);} 
+}
+template <class... TLs>
+using concatenate_type_lists_t = decltype(impl::ConcatenateTypeLists(std::declval<TLs>()...));
+
+} // namespace parthenon
+
+#endif // UTILS_TYPE_LIST_HPP_
\ No newline at end of file

From 662c79224f4969608d8111beaf101893c42c2fea Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 30 Jul 2024 15:03:46 -0600
Subject: [PATCH 04/40] reorder

---
 src/solvers/bicgstab_solver.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/solvers/bicgstab_solver.hpp b/src/solvers/bicgstab_solver.hpp
index 9fa4d81d239a..a6579df0fd0b 100644
--- a/src/solvers/bicgstab_solver.hpp
+++ b/src/solvers/bicgstab_solver.hpp
@@ -91,7 +91,6 @@ class BiCGSTABSolver {
     using namespace refinement_ops;
     auto m_no_ghost =
         Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, shape);
-    pkg->AddField(x::name(), m_no_ghost);
     pkg->AddField(rhat0::name(), m_no_ghost);
     pkg->AddField(v::name(), m_no_ghost);
     pkg->AddField(h::name(), m_no_ghost);
@@ -99,6 +98,7 @@ class BiCGSTABSolver {
     pkg->AddField(t::name(), m_no_ghost);
     pkg->AddField(r::name(), m_no_ghost);
     pkg->AddField(p::name(), m_no_ghost);
+    pkg->AddField(x::name(), m_no_ghost);
   }
 
   template <class TL_t>

From d8306113dabf291b3ff7d3a4f5b9c22159959c3e Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 30 Jul 2024 20:56:34 -0600
Subject: [PATCH 05/40] add method for getting names from a typelist

---
 src/utils/type_list.hpp | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/utils/type_list.hpp b/src/utils/type_list.hpp
index ede6ebdc4dc7..ad64094a222c 100644
--- a/src/utils/type_list.hpp
+++ b/src/utils/type_list.hpp
@@ -68,9 +68,17 @@ namespace impl {
   template<class... Args1, class... Args2, class... Args> 
   auto ConcatenateTypeLists(TypeList<Args1...>, TypeList<Args2...>, Args...){ return ConcatenateTypeLists(TypeList<Args1..., Args2...>(), Args()...);} 
 }
+
 template <class... TLs>
 using concatenate_type_lists_t = decltype(impl::ConcatenateTypeLists(std::declval<TLs>()...));
 
+// Relevant only for lists of variable types 
+template <class TL> 
+auto GetNames(){
+  std::vector<std::string> names; 
+  TL::IterateTypes([&names](auto t){ names.push_back(decltype(t)::name()); });
+  return names;
+}
 } // namespace parthenon
 
 #endif // UTILS_TYPE_LIST_HPP_
\ No newline at end of file

From 090a8d968f67aa9785bb6ecd9d64156d28206b5d Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 30 Jul 2024 21:47:12 -0600
Subject: [PATCH 06/40] format and lint

---
 src/interface/make_pack_descriptor.hpp |  8 ++--
 src/interface/sparse_pack.hpp          |  2 +-
 src/interface/sparse_pack_base.cpp     |  6 +--
 src/interface/variable_state.hpp       |  8 ++--
 src/solvers/bicgstab_solver.hpp        | 21 ++++++----
 src/solvers/mg_solver.hpp              | 11 +++---
 src/utils/constants.hpp                |  4 +-
 src/utils/type_list.hpp                | 55 +++++++++++++++-----------
 8 files changed, 64 insertions(+), 51 deletions(-)

diff --git a/src/interface/make_pack_descriptor.hpp b/src/interface/make_pack_descriptor.hpp
index 7620a99a0c02..afae3c314787 100644
--- a/src/interface/make_pack_descriptor.hpp
+++ b/src/interface/make_pack_descriptor.hpp
@@ -137,13 +137,13 @@ inline auto MakePackDescriptor(StateDescriptor *psd, const std::vector<Uid_t> &v
   return typename SparsePack<>::Descriptor(base_desc);
 }
 
-template <class... Types, class... Args> 
-inline auto MakePackDescriptorFromTypeList(TypeList<Types...>, Args&&... args) { 
+template <class... Types, class... Args>
+inline auto MakePackDescriptorFromTypeList(TypeList<Types...>, Args &&...args) {
   return MakePackDescriptor<Types...>(std::forward<Args>(args)...);
 }
 
-template <class TL, class... Args> 
-inline auto MakePackDescriptorFromTypeList(Args&&... args) { 
+template <class TL, class... Args>
+inline auto MakePackDescriptorFromTypeList(Args &&...args) {
   return MakePackDescriptorFromTypeList(TL(), std::forward<Args>(args)...);
 }
 
diff --git a/src/interface/sparse_pack.hpp b/src/interface/sparse_pack.hpp
index 7ef7a6fe1f41..e939e3dc7004 100644
--- a/src/interface/sparse_pack.hpp
+++ b/src/interface/sparse_pack.hpp
@@ -33,8 +33,8 @@
 #include "interface/sparse_pack_base.hpp"
 #include "interface/variable.hpp"
 #include "utils/concepts_lite.hpp"
-#include "utils/utils.hpp"
 #include "utils/type_list.hpp"
+#include "utils/utils.hpp"
 
 namespace parthenon {
 
diff --git a/src/interface/sparse_pack_base.cpp b/src/interface/sparse_pack_base.cpp
index dc0b36fc860c..d5a66db01498 100644
--- a/src/interface/sparse_pack_base.cpp
+++ b/src/interface/sparse_pack_base.cpp
@@ -249,9 +249,9 @@ SparsePackBase SparsePackBase::Build(T *pmd, const PackDescriptor &desc,
                   PARTHENON_REQUIRE(
                       pack.pack_h_(0, b, idx).size() > 0,
                       "Seems like this variable might not actually be allocated.");
-                  
-                  // Save the tensor indices for easy access 
-                  for (int l = 0; l < leading_dim; ++l) { 
+
+                  // Save the tensor indices for easy access
+                  for (int l = 0; l < leading_dim; ++l) {
                     pack.pack_h_(l, b, idx).t = t;
                     pack.pack_h_(l, b, idx).u = u;
                     pack.pack_h_(l, b, idx).v = v;
diff --git a/src/interface/variable_state.hpp b/src/interface/variable_state.hpp
index f25b39d423d5..5e8afcd79d7b 100644
--- a/src/interface/variable_state.hpp
+++ b/src/interface/variable_state.hpp
@@ -61,10 +61,10 @@ struct VariableState : public empty_state_t {
   int sparse_id;
   int vector_component = NODIR;
   bool initialized = true;
-  
-  int t{0};
-  int u{0};
-  int v{0};
+
+  int t{-1};
+  int u{-1};
+  int v{-1};
   TopologicalType topological_type = TopologicalType::Cell;
   TopologicalElement topological_element = TopologicalElement::CC;
   std::size_t tensor_components;
diff --git a/src/solvers/bicgstab_solver.hpp b/src/solvers/bicgstab_solver.hpp
index a6579df0fd0b..caa594337b61 100644
--- a/src/solvers/bicgstab_solver.hpp
+++ b/src/solvers/bicgstab_solver.hpp
@@ -69,25 +69,30 @@ class BiCGSTABSolver {
   PARTHENON_INTERNALSOLVERVARIABLE(u, r);
   PARTHENON_INTERNALSOLVERVARIABLE(u, p);
   PARTHENON_INTERNALSOLVERVARIABLE(u, x);
-  
+
   using internal_types_tl = TypeList<rhat0, v, h, s, t, r, p, x>;
-  using preconditioner_t = MGSolver<u, rhs, equations>; 
-  using all_internal_types_tl = concatenate_type_lists_t<internal_types_tl, typename preconditioner_t::internal_types_tl>;
+  using preconditioner_t = MGSolver<u, rhs, equations>;
+  using all_internal_types_tl =
+      concatenate_type_lists_t<internal_types_tl,
+                               typename preconditioner_t::internal_types_tl>;
 
   std::vector<std::string> GetInternalVariableNames() const {
     std::vector<std::string> names;
     if (params_.precondition) {
-      all_internal_types_tl::IterateTypes([&names](auto t){names.push_back(decltype(t)::name());});
+      all_internal_types_tl::IterateTypes(
+          [&names](auto t) { names.push_back(decltype(t)::name()); });
     } else {
-      internal_types_tl::IterateTypes([&names](auto t){names.push_back(decltype(t)::name());});
+      internal_types_tl::IterateTypes(
+          [&names](auto t) { names.push_back(decltype(t)::name()); });
     }
     return names;
   }
 
   BiCGSTABSolver(StateDescriptor *pkg, BiCGSTABParams params_in,
-                 equations eq_in = equations(), std::vector<int> shape = {}, const std::string &container = "base")
-      : preconditioner(pkg, params_in.mg_params, eq_in, shape, container), params_(params_in),
-        iter_counter(0), eqs_(eq_in), container_(container) {
+                 equations eq_in = equations(), std::vector<int> shape = {},
+                 const std::string &container = "base")
+      : preconditioner(pkg, params_in.mg_params, eq_in, shape, container),
+        params_(params_in), iter_counter(0), eqs_(eq_in), container_(container) {
     using namespace refinement_ops;
     auto m_no_ghost =
         Metadata({Metadata::Cell, Metadata::Derived, Metadata::OneCopy}, shape);
diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
index 4c5b7c5b9343..c44d3c020d25 100644
--- a/src/solvers/mg_solver.hpp
+++ b/src/solvers/mg_solver.hpp
@@ -77,16 +77,17 @@ class MGSolver {
   PARTHENON_INTERNALSOLVERVARIABLE(u, temp); // Temporary storage
   PARTHENON_INTERNALSOLVERVARIABLE(u, u0);   // Storage for initial solution during FAS
   PARTHENON_INTERNALSOLVERVARIABLE(u, D);    // Storage for (approximate) diagonal
-  
-  using internal_types_tl = TypeList<res_err, temp, u0, D>; 
+
+  using internal_types_tl = TypeList<res_err, temp, u0, D>;
   std::vector<std::string> GetInternalVariableNames() const {
-    std::vector<std::string> names; 
-    internal_types_tl::IterateTypes([&names](auto t){names.push_back(decltype(t)::name());});
+    std::vector<std::string> names;
+    internal_types_tl::IterateTypes(
+        [&names](auto t) { names.push_back(decltype(t)::name()); });
     return names;
   }
 
   MGSolver(StateDescriptor *pkg, MGParams params_in, equations eq_in = equations(),
-           std::vector<int> shape = {}, const std::string& container = "base")
+           std::vector<int> shape = {}, const std::string &container = "base")
       : params_(params_in), iter_counter(0), eqs_(eq_in), container_(container) {
     using namespace parthenon::refinement_ops;
     // The ghost cells of res_err need to be filled, but this is accomplished by
diff --git a/src/utils/constants.hpp b/src/utils/constants.hpp
index 17386cf54e61..69f03d6ca716 100644
--- a/src/utils/constants.hpp
+++ b/src/utils/constants.hpp
@@ -117,8 +117,8 @@ class PhysicalConstants {
   static constexpr double stefan_boltzmann = 2.0 * M_PI * M_PI * M_PI * M_PI * M_PI * kb *
                                              kb * kb * kb / (15.0 * h * h * h * c * c);
   static constexpr double sb = stefan_boltzmann;
-  
-  // Radiation constant 
+
+  // Radiation constant
   static constexpr double radiation_constant = 4.0 * stefan_boltzmann / speed_of_light;
   static constexpr double ar = radiation_constant;
 
diff --git a/src/utils/type_list.hpp b/src/utils/type_list.hpp
index ad64094a222c..965c4e4be9bb 100644
--- a/src/utils/type_list.hpp
+++ b/src/utils/type_list.hpp
@@ -14,71 +14,78 @@
 #ifndef UTILS_TYPE_LIST_HPP_
 #define UTILS_TYPE_LIST_HPP_
 
+#include <string>
 #include <tuple>
 #include <type_traits>
 #include <utility>
+#include <vector>
 
 namespace parthenon {
-template <class... Args> 
+template <class... Args>
 struct TypeList {
   using types = std::tuple<Args...>;
-  
+
   static constexpr std::size_t n_types{sizeof...(Args)};
 
   template <std::size_t I>
   using type = typename std::tuple_element<I, types>::type;
-  
+
   template <std::size_t... Idxs>
   using sublist = TypeList<type<Idxs>...>;
 
   template <class T, std::size_t I = 0>
-  static constexpr std::size_t GetIdx(){
+  static constexpr std::size_t GetIdx() {
     static_assert(I < n_types, "Type is not present in TypeList.");
     if constexpr (std::is_same_v<T, type<I>>) {
       return I;
     } else {
-      return GetIdx<T, I + 1>();   
+      return GetIdx<T, I + 1>();
     }
-  } 
-  
-  template<class F>
+  }
+
+  template <class F>
   static void IterateTypes(F func) {
     (func(Args()), ...);
   }
 
- private:  
-  template<std::size_t Start, std::size_t End>
+ private:
+  template <std::size_t Start, std::size_t End>
   static auto ContinuousSublist() {
     return ContinuousSublistImpl<Start>(std::make_index_sequence<End - Start + 1>());
   }
-  template<std::size_t Start, std::size_t... Is> 
+  template <std::size_t Start, std::size_t... Is>
   static auto ContinuousSublistImpl(std::index_sequence<Is...>) {
     return sublist<(Start + Is)...>();
   }
 
  public:
-  template<std::size_t Start, std::size_t End>
+  template <std::size_t Start, std::size_t End>
   using continuous_sublist = decltype(ContinuousSublist<Start, End>());
 };
 
 namespace impl {
-  template<class... Args>
-  auto ConcatenateTypeLists(TypeList<Args...>){return TypeList<Args...>();}
-  
-  template<class... Args1, class... Args2, class... Args> 
-  auto ConcatenateTypeLists(TypeList<Args1...>, TypeList<Args2...>, Args...){ return ConcatenateTypeLists(TypeList<Args1..., Args2...>(), Args()...);} 
+template <class... Args>
+auto ConcatenateTypeLists(TypeList<Args...>) {
+  return TypeList<Args...>();
+}
+
+template <class... Args1, class... Args2, class... Args>
+auto ConcatenateTypeLists(TypeList<Args1...>, TypeList<Args2...>, Args...) {
+  return ConcatenateTypeLists(TypeList<Args1..., Args2...>(), Args()...);
 }
+} // namespace impl
 
 template <class... TLs>
-using concatenate_type_lists_t = decltype(impl::ConcatenateTypeLists(std::declval<TLs>()...));
+using concatenate_type_lists_t =
+    decltype(impl::ConcatenateTypeLists(std::declval<TLs>()...));
 
-// Relevant only for lists of variable types 
-template <class TL> 
-auto GetNames(){
-  std::vector<std::string> names; 
-  TL::IterateTypes([&names](auto t){ names.push_back(decltype(t)::name()); });
+// Relevant only for lists of variable types
+template <class TL>
+auto GetNames() {
+  std::vector<std::string> names;
+  TL::IterateTypes([&names](auto t) { names.push_back(decltype(t)::name()); });
   return names;
 }
 } // namespace parthenon
 
-#endif // UTILS_TYPE_LIST_HPP_
\ No newline at end of file
+#endif // UTILS_TYPE_LIST_HPP_

From 052b4d331eafea0c4aeec577021fc3b67ca18323 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 30 Jul 2024 22:03:09 -0600
Subject: [PATCH 07/40] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 9cb66c3d573f..40075a33864a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1143]](https://github.com/parthenon-hpc-lab/parthenon/pull/1143) Add tensor indices to VariableState, add radiation constant to constants, add TypeLists, allow for arbitrary containers for solvers
 - [[PR 1140]](https://github.com/parthenon-hpc-lab/parthenon/pull/1140) Allow for relative convergence tolerance in BiCGSTAB solver.
 - [[PR 1047]](https://github.com/parthenon-hpc-lab/parthenon/pull/1047) General three- and four-valent 2D forests w/ arbitrary orientations.
 - [[PR 1130]](https://github.com/parthenon-hpc-lab/parthenon/pull/1130) Enable `parthenon::par_reduce` for MD loops with Kokkos 1D Range

From 311de157c4285cd6db6ff63c13a9e85b6c4c0533 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Wed, 31 Jul 2024 11:00:46 -0600
Subject: [PATCH 08/40] Update src/interface/make_pack_descriptor.hpp

Co-authored-by: Adam C Reyes <adam.c.reyes@gmail.com>
---
 src/interface/make_pack_descriptor.hpp | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/interface/make_pack_descriptor.hpp b/src/interface/make_pack_descriptor.hpp
index afae3c314787..2a604338cdb7 100644
--- a/src/interface/make_pack_descriptor.hpp
+++ b/src/interface/make_pack_descriptor.hpp
@@ -137,8 +137,8 @@ inline auto MakePackDescriptor(StateDescriptor *psd, const std::vector<Uid_t> &v
   return typename SparsePack<>::Descriptor(base_desc);
 }
 
-template <class... Types, class... Args>
-inline auto MakePackDescriptorFromTypeList(TypeList<Types...>, Args &&...args) {
+template <template <class...> class TL, class... Types, class... Args>
+inline auto MakePackDescriptorFromTypeList(TL<Types...>, Args &&...args) {
   return MakePackDescriptor<Types...>(std::forward<Args>(args)...);
 }
 
@@ -146,7 +146,6 @@ template <class TL, class... Args>
 inline auto MakePackDescriptorFromTypeList(Args &&...args) {
   return MakePackDescriptorFromTypeList(TL(), std::forward<Args>(args)...);
 }
-
 } // namespace parthenon
 
 #endif // INTERFACE_MAKE_PACK_DESCRIPTOR_HPP_

From 0b331fe52b3d298e5b5b7a30a29f210f6f9bbe85 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Wed, 31 Jul 2024 12:09:22 -0600
Subject: [PATCH 09/40] remove some compiler warnings

---
 src/interface/update.cpp     |  5 -----
 src/mesh/mesh_refinement.cpp | 16 +---------------
 src/tasks/tasks.hpp          |  3 +++
 3 files changed, 4 insertions(+), 20 deletions(-)

diff --git a/src/interface/update.cpp b/src/interface/update.cpp
index cf79d9980cba..6282490a1073 100644
--- a/src/interface/update.cpp
+++ b/src/interface/update.cpp
@@ -156,11 +156,6 @@ TaskStatus SparseDealloc(MeshData<Real> *md) {
   auto packIdx = desc.GetMap();
 
   ParArray2D<bool> is_zero("IsZero", pack.GetNBlocks(), pack.GetMaxNumberOfVars());
-  const int Ni = ib.e + 1 - ib.s;
-  const int Nj = jb.e + 1 - jb.s;
-  const int Nk = kb.e + 1 - kb.s;
-  const int NjNi = Nj * Ni;
-  const int NkNjNi = Nk * NjNi;
   Kokkos::parallel_for(
       PARTHENON_AUTO_LABEL,
       Kokkos::TeamPolicy<>(parthenon::DevExecSpace(), pack.GetNBlocks(), Kokkos::AUTO),
diff --git a/src/mesh/mesh_refinement.cpp b/src/mesh/mesh_refinement.cpp
index 88d1da0a281a..75c9964e6a6b 100644
--- a/src/mesh/mesh_refinement.cpp
+++ b/src/mesh/mesh_refinement.cpp
@@ -98,21 +98,7 @@ void MeshRefinement::SetRefinement(AmrTag flag) {
       deref_count_ = 0;
     } else {
       deref_count_++;
-      int ec = 0, js, je, ks, ke;
-      if (!pmb->block_size.symmetry(X2DIR)) {
-        js = -1;
-        je = 1;
-      } else {
-        js = 0;
-        je = 0;
-      }
-      if (!pmb->block_size.symmetry(X3DIR)) {
-        ks = -1;
-        ke = 1;
-      } else {
-        ks = 0;
-        ke = 0;
-      }
+      int ec = 0;
       for (const auto &nb : pmb->neighbors) {
         if (nb.loc.level() > pmb->loc.level()) ec++;
       }
diff --git a/src/tasks/tasks.hpp b/src/tasks/tasks.hpp
index b495224d0ddf..21c85c317c79 100644
--- a/src/tasks/tasks.hpp
+++ b/src/tasks/tasks.hpp
@@ -296,7 +296,10 @@ class TaskList {
                          std::forward<Args>(args)...);
     } else {
       static_assert(always_false<Arg1>, "Bad signature for AddTask.");
+      return TaskID();
     }
+    // Stops some compilers from complaining
+    return TaskID();
   }
 
   template <class... Args>

From 4a3f79a2a13c5bb962bf2b1ab02b65af1ccfc15d Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Wed, 31 Jul 2024 12:13:14 -0600
Subject: [PATCH 10/40] remove more warnings

---
 src/solvers/mg_solver.hpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/solvers/mg_solver.hpp b/src/solvers/mg_solver.hpp
index c44d3c020d25..ee8cfff177ab 100644
--- a/src/solvers/mg_solver.hpp
+++ b/src/solvers/mg_solver.hpp
@@ -394,13 +394,13 @@ class MGSolver {
       PARTHENON_FAIL("Unknown solver type.");
     }
 
-    auto decorate_task_name = [partition, level](const std::string &in, auto b) {
-      return std::make_tuple(in + "(p:" + std::to_string(partition) +
-                                 ", l:" + std::to_string(level) + ")",
-                             1, b);
-    };
+//    auto decorate_task_name = [partition, level](const std::string &in, auto b) {
+//      return std::make_tuple(in + "(p:" + std::to_string(partition) +
+//                                 ", l:" + std::to_string(level) + ")",
+//                             1, b);
+//    };
 
-//#define BTF(...) decorate_task_name(TF(__VA_ARGS__))
+// #define BTF(...) decorate_task_name(TF(__VA_ARGS__))
 #define BTF(...) TF(__VA_ARGS__)
     bool multilevel = (level != min_level);
 

From f63eadbfb8fb21e4588bbbda9fb2a8efe60b5884 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 6 Aug 2024 12:12:11 -0600
Subject: [PATCH 11/40] Add a little more documentation

---
 doc/sphinx/src/concepts_lite.rst |  6 ----
 doc/sphinx/src/utilities.rst     | 48 ++++++++++++++++++++++++++++++++
 src/utils/type_list.hpp          |  6 ++++
 3 files changed, 54 insertions(+), 6 deletions(-)
 delete mode 100644 doc/sphinx/src/concepts_lite.rst
 create mode 100644 doc/sphinx/src/utilities.rst

diff --git a/doc/sphinx/src/concepts_lite.rst b/doc/sphinx/src/concepts_lite.rst
deleted file mode 100644
index fa7f01066ed3..000000000000
--- a/doc/sphinx/src/concepts_lite.rst
+++ /dev/null
@@ -1,6 +0,0 @@
-C++11 Style Concepts Implementation
-===================================
-
-*This documentation needs to be written (see issue #695), but there are
-extensive comments in src/utlils/concepts_lite.hpp and examples of
-useage in tst/unit/test_concepts_lite.hpp*
diff --git a/doc/sphinx/src/utilities.rst b/doc/sphinx/src/utilities.rst
new file mode 100644
index 000000000000..43ae4ec1ee77
--- /dev/null
+++ b/doc/sphinx/src/utilities.rst
@@ -0,0 +1,48 @@
+``TypeList``s
+=============
+
+Provides a wrapper class around a variadic pack of types to simplify
+performing compile time operations on the pack. There are templated 
+types defined giving the type at a particular index in the pack, types 
+for extracting sub-``TypeList``s of the original type list, and ``constexpr``
+functions for getting the index of the first instance of a type in the 
+pack. Additionally it provides a capability for iterating an ``auto`` lambda
+over the type list, which can be useful for calling a ``static`` function 
+defined for each of the types on each of the types. *In the future, it
+would be nice to have the capability to make a unique type list from
+another type list (i.e. the unique one only a single instance of each type
+in the original type list)*
+
+``TypeList``s have many applications and are commonly found in many 
+codebases, but in Parthenon one of the main use cases is for storing 
+lists of types associated with field variables that are used in type 
+based ``SparsePack``s.  
+
+Robust
+======
+
+Provides a number of functions for doing operations on floating point
+numbers that are bounded, avoid division by zero, etc. 
+
+C++11 Style Concepts Implementation
+===================================
+
+*This documentation needs to be written (see issue #695), but there are
+extensive comments in src/utlils/concepts_lite.hpp and examples of
+useage in tst/unit/test_concepts_lite.hpp*
+
+``Indexer``
+===========
+
+Provides functionality for iterating over an arbitrary dimensional
+hyper-rectangular index space using a flattened loop. Specific 
+instantiations, e.g. ``Indexer5D``, are provided up to eight 
+dimensions. Useage:
+
+.. code:: cpp
+
+  Indexer4D idxer({0, 3}, {1, 2}, {0, 5}, {10, 16});
+  for (int flat_idx = 0; flat_idx < idxer.size(); ++flat_idx) {
+    auto [i, j, k, l] = idxer(flat_idx);
+    // Do stuff in the 4D index space...
+  }
diff --git a/src/utils/type_list.hpp b/src/utils/type_list.hpp
index 965c4e4be9bb..a17b7b8f7063 100644
--- a/src/utils/type_list.hpp
+++ b/src/utils/type_list.hpp
@@ -21,6 +21,12 @@
 #include <vector>
 
 namespace parthenon {
+
+// Convenience struct for holding a variadic pack of types
+// and providing compile time indexing into that pack as
+// well as the ability to get the index of a given type within
+// the pack. Functions are available below for compile time
+// concatenation of TypeLists
 template <class... Args>
 struct TypeList {
   using types = std::tuple<Args...>;

From 57b7304ebe1d3687f6ac66b6172e4c69bb7d7aec Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Thu, 8 Aug 2024 21:28:37 +0200
Subject: [PATCH 12/40] Reduce memory consumption for buffer pool

---
 src/bvals/comms/build_boundary_buffers.cpp | 39 +++++++++++++++++-----
 src/utils/object_pool.hpp                  |  4 +++
 2 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/src/bvals/comms/build_boundary_buffers.cpp b/src/bvals/comms/build_boundary_buffers.cpp
index 918f4d5017c2..a4c268d39d35 100644
--- a/src/bvals/comms/build_boundary_buffers.cpp
+++ b/src/bvals/comms/build_boundary_buffers.cpp
@@ -16,10 +16,12 @@
 //========================================================================================
 
 #include <algorithm>
+#include <cstddef>
 #include <iostream> // debug
 #include <memory>
 #include <random>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 #include "bvals_in_one.hpp"
@@ -44,25 +46,46 @@ template <BoundaryType BTYPE>
 void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
                                Mesh::comm_buf_map_t &buf_map) {
   Mesh *pmesh = md->GetMeshPointer();
+  std::unordered_map<int, int> nbufs;
+
+  ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
+    // Calculate the required size of the buffer for this boundary
+    int buf_size = GetBufferSize(pmb, nb, v);
+    if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
+
+    nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
+  });
+
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
     int buf_size = GetBufferSize(pmb, nb, v);
     if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
 
     // Add a buffer pool if one does not exist for this size
+    using buf_t = buf_pool_t<Real>::base_t;
     if (pmesh->pool_map.count(buf_size) == 0) {
-      pmesh->pool_map.emplace(std::make_pair(
-          buf_size, buf_pool_t<Real>([buf_size](buf_pool_t<Real> *pool) {
-            using buf_t = buf_pool_t<Real>::base_t;
-            // TODO(LFR): Make nbuf a user settable parameter
-            const int nbuf = 200;
-            buf_t chunk("pool buffer", buf_size * nbuf);
-            for (int i = 1; i < nbuf; ++i) {
+      pmesh->pool_map.emplace(
+          buf_size, buf_pool_t<Real>([buf_size, &nbufs](buf_pool_t<Real> *pool) {
+            const auto pool_size = static_cast<int64_t>(nbufs[buf_size]) * buf_size;
+            buf_t chunk("pool buffer", pool_size);
+            for (int i = 1; i < nbufs[buf_size]; ++i) {
               pool->AddFreeObjectToPool(
                   buf_t(chunk, std::make_pair(i * buf_size, (i + 1) * buf_size)));
             }
             return buf_t(chunk, std::make_pair(0, buf_size));
-          })));
+          }));
+      // or add to existing pool (if required)
+    } else {
+      auto pool = pmesh->pool_map[buf_size];
+      const auto new_buffers_req = nbufs[buf_size] - pool.NumAvailable();
+      if (new_buffers_req > 1) {
+        const auto pool_size = static_cast<int64_t>(new_buffers_req) * buf_size;
+        buf_t chunk("pool buffer", pool_size);
+        for (int i = 1; i < new_buffers_req; ++i) {
+          pool.AddFreeObjectToPool(
+              buf_t(chunk, std::make_pair(i * buf_size, (i + 1) * buf_size)));
+        }
+      }
     }
 
     const int receiver_rank = nb.rank;
diff --git a/src/utils/object_pool.hpp b/src/utils/object_pool.hpp
index 89167334277f..f69efad3810f 100644
--- a/src/utils/object_pool.hpp
+++ b/src/utils/object_pool.hpp
@@ -62,6 +62,10 @@ class ObjectPool {
     std::cout << inuse_.size() << " used objects." << std::endl;
   }
 
+  auto NumAvailable() const {
+    return available_.size();
+  }
+
   std::uint64_t SizeInBytes() const {
     constexpr std::uint64_t datum_size = sizeof(typename base_t::value_type);
     std::uint64_t object_size = 0;

From 5522e48181f751ba2741471aa972bd41e8b632b0 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Fri, 9 Aug 2024 09:59:42 +0200
Subject: [PATCH 13/40] Only count vars for allocated variables

---
 src/bvals/comms/bnd_info.cpp               |  2 +-
 src/bvals/comms/build_boundary_buffers.cpp | 10 ++++++----
 src/utils/object_pool.hpp                  |  4 +---
 3 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/bvals/comms/bnd_info.cpp b/src/bvals/comms/bnd_info.cpp
index d30522fd6998..1505b56f956a 100644
--- a/src/bvals/comms/bnd_info.cpp
+++ b/src/bvals/comms/bnd_info.cpp
@@ -332,7 +332,7 @@ BndInfo BndInfo::GetSetBndInfo(MeshBlock *pmb, const NeighborBlock &nb,
     out.buf_allocated = false;
   } else {
     printf("%i [rank: %i] -> %i [rank: %i] (Set %s) is in state %i.\n", nb.gid, nb.rank,
-           pmb->gid, Globals::my_rank, v->label().c_str(), buf_state);
+           pmb->gid, Globals::my_rank, v->label().c_str(), static_cast<int>(buf_state));
     PARTHENON_FAIL("Buffer should be in a received state.");
   }
   return out;
diff --git a/src/bvals/comms/build_boundary_buffers.cpp b/src/bvals/comms/build_boundary_buffers.cpp
index a4c268d39d35..800443381a51 100644
--- a/src/bvals/comms/build_boundary_buffers.cpp
+++ b/src/bvals/comms/build_boundary_buffers.cpp
@@ -53,7 +53,9 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     int buf_size = GetBufferSize(pmb, nb, v);
     if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
 
-    nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
+    if (v->IsAllocated()) {
+      nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
+    }
   });
 
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
@@ -64,7 +66,7 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     // Add a buffer pool if one does not exist for this size
     using buf_t = buf_pool_t<Real>::base_t;
     if (pmesh->pool_map.count(buf_size) == 0) {
-      pmesh->pool_map.emplace(
+      pmesh->pool_map.emplace(std::make_pair(
           buf_size, buf_pool_t<Real>([buf_size, &nbufs](buf_pool_t<Real> *pool) {
             const auto pool_size = static_cast<int64_t>(nbufs[buf_size]) * buf_size;
             buf_t chunk("pool buffer", pool_size);
@@ -73,10 +75,10 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
                   buf_t(chunk, std::make_pair(i * buf_size, (i + 1) * buf_size)));
             }
             return buf_t(chunk, std::make_pair(0, buf_size));
-          }));
+          })));
       // or add to existing pool (if required)
     } else {
-      auto pool = pmesh->pool_map[buf_size];
+      auto &pool = pmesh->pool_map.at(buf_size);
       const auto new_buffers_req = nbufs[buf_size] - pool.NumAvailable();
       if (new_buffers_req > 1) {
         const auto pool_size = static_cast<int64_t>(new_buffers_req) * buf_size;
diff --git a/src/utils/object_pool.hpp b/src/utils/object_pool.hpp
index f69efad3810f..e758a7867fe2 100644
--- a/src/utils/object_pool.hpp
+++ b/src/utils/object_pool.hpp
@@ -62,9 +62,7 @@ class ObjectPool {
     std::cout << inuse_.size() << " used objects." << std::endl;
   }
 
-  auto NumAvailable() const {
-    return available_.size();
-  }
+  auto NumAvailable() const { return available_.size(); }
 
   std::uint64_t SizeInBytes() const {
     constexpr std::uint64_t datum_size = sizeof(typename base_t::value_type);

From 4bef0509faa15b201f4e248aedc517f0765adbd7 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Fri, 9 Aug 2024 11:15:43 +0200
Subject: [PATCH 14/40] printf debugging...

---
 src/bvals/comms/build_boundary_buffers.cpp | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/bvals/comms/build_boundary_buffers.cpp b/src/bvals/comms/build_boundary_buffers.cpp
index 800443381a51..44d2b1af9012 100644
--- a/src/bvals/comms/build_boundary_buffers.cpp
+++ b/src/bvals/comms/build_boundary_buffers.cpp
@@ -53,11 +53,12 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     int buf_size = GetBufferSize(pmb, nb, v);
     if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
 
-    if (v->IsAllocated()) {
-      nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
-    }
+    nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
+    std::cerr << "\nnbufs[" << buf_size << "] += 1, which is now " << nbufs[buf_size];
   });
 
+  std::cerr << "\nbefore nbufs.size = " << nbufs.size();
+
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
     int buf_size = GetBufferSize(pmb, nb, v);
@@ -66,6 +67,9 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     // Add a buffer pool if one does not exist for this size
     using buf_t = buf_pool_t<Real>::base_t;
     if (pmesh->pool_map.count(buf_size) == 0) {
+      std::cerr << "\ninside nbufs.size = " << nbufs.size();
+      std::cerr << "\nGoing to request " << nbufs[buf_size] << " buffers of size "
+                << buf_size;
       pmesh->pool_map.emplace(std::make_pair(
           buf_size, buf_pool_t<Real>([buf_size, &nbufs](buf_pool_t<Real> *pool) {
             const auto pool_size = static_cast<int64_t>(nbufs[buf_size]) * buf_size;

From 5e42a4fe53e2bdb3dc7a8ce57ec05d79f439071f Mon Sep 17 00:00:00 2001
From: "Adam M. Dempsey" <adamdemps@gmail.com>
Date: Fri, 9 Aug 2024 09:42:42 -0600
Subject: [PATCH 15/40] Return DriverStatus::timeout when appropriate (#1135)

* If the driver should KeepGoing() but finished set status to timeout instead of complete

* changelog for timeout

* Return 1 when driver fails, not 2

* Add a default error handler to the test cases. Only raise an error when the return code is 1

* copyright

* Move changelog entry to breaking section

---------

Co-authored-by: Adam Dempsey <adempsey@lanl.gov>
---
 CHANGELOG.md                      |  1 +
 example/sparse_advection/main.cpp |  6 ++---
 src/driver/driver.cpp             |  2 +-
 tst/regression/utils/test_case.py | 37 ++++++++++++++++++++++---------
 4 files changed, 31 insertions(+), 15 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 188e60294a10..89afa174910f 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -83,6 +83,7 @@
 - [[PR 1108]](https://github.com/parthenon-hpc-lab/parthenon/pull/1108) Remove NaN payload tags infrastructure
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 1135]](https://github.com/parthenon-hpc-lab/parthenon/pull/1135) Drivers now correctly return DriverStatus::timeout on hittig walltime limit
 - [[PR 1128]](https://github.com/parthenon-hpc-lab/parthenon/pull/1128) Add cycle and nbtotal to hst
 - [[PR 1108]](https://github.com/parthenon-hpc-lab/parthenon/pull/1108) Remove NaN payload tags infrastructure
 - [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
diff --git a/example/sparse_advection/main.cpp b/example/sparse_advection/main.cpp
index 3f513670596b..585d1a35ab1e 100644
--- a/example/sparse_advection/main.cpp
+++ b/example/sparse_advection/main.cpp
@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2021-2023. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2021-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -59,9 +59,9 @@ int main(int argc, char *argv[]) {
 
   if (driver_status == DriverStatus::complete) {
     return 0;
-  } else if (driver_status == DriverStatus::timeout) {
-    return 1;
   } else if (driver_status == DriverStatus::failed) {
+    return 1;
+  } else if (driver_status == DriverStatus::timeout) {
     return 2;
   } else {
     return 3;
diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index d33ee9122cf9..332468e2b83b 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -152,7 +152,7 @@ DriverStatus EvolutionDriver::Execute() {
     pmesh->UserWorkAfterLoop(pmesh, pinput, tm);
   }
 
-  DriverStatus status = DriverStatus::complete;
+  DriverStatus status = tm.KeepGoing() ? DriverStatus::timeout : DriverStatus::complete;
   // Do *not* write the "final" output, if this is analysis run.
   // The analysis output itself has already been written above before the main loop.
   if (signal != OutputSignal::analysis) {
diff --git a/tst/regression/utils/test_case.py b/tst/regression/utils/test_case.py
index 7b09d00358d5..b93bdade88b2 100644
--- a/tst/regression/utils/test_case.py
+++ b/tst/regression/utils/test_case.py
@@ -49,6 +49,10 @@ def Prepare(parameters):
     def Analyse(parameters):
         raise NotImplementedError("Every TestCase must initialse an Analyse method")
 
+    def ErrorOnNonZeroReturnCode(self, parameters, returncode):
+        # Do not raise an error when returncode /= 1
+        return returncode == 1
+
 
 class TestManager:
     def __init__(self, run_test_path, **kwargs):
@@ -244,18 +248,29 @@ def Run(self):
             print(proc.stdout.decode())
             self.parameters.stdouts.append(proc.stdout)
         except subprocess.CalledProcessError as err:
-            print("\n*****************************************************************")
-            print("Subprocess error message")
-            print("*****************************************************************\n")
-            print(str(repr(err.output)).replace("\\n", os.linesep))
-            print("\n*****************************************************************")
-            print("Error detected while running subprocess command")
-            print("*****************************************************************\n")
-            raise TestManagerError(
-                "\nReturn code {0} from command '{1}'".format(
-                    err.returncode, " ".join(err.cmd)
+            if self.test_case.ErrorOnNonZeroReturnCode(self.parameters, err.returncode):
+                print(
+                    "\n*****************************************************************"
                 )
-            )
+                print("Subprocess error message")
+                print(
+                    "*****************************************************************\n"
+                )
+                print(str(repr(err.output)).replace("\\n", os.linesep))
+                print(
+                    "\n*****************************************************************"
+                )
+                print("Error detected while running subprocess command")
+                print(
+                    "*****************************************************************\n"
+                )
+                raise TestManagerError(
+                    "\nReturn code {0} from command '{1}'".format(
+                        err.returncode, " ".join(err.cmd)
+                    )
+                )
+            else:
+                self.parameters.stdouts.append(err.stdout)
         # Reset parameters
         self.parameters.coverage_status = "only-regression"
 

From edd05ee2c796b61ce1d46ca22a3eac6c61cf918f Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Sat, 10 Aug 2024 16:45:39 -0600
Subject: [PATCH 16/40] Add explicit free after calling __cxa_demangle

---
 src/tasks/tasks.hpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/tasks/tasks.hpp b/src/tasks/tasks.hpp
index 21c85c317c79..91b9e103d6aa 100644
--- a/src/tasks/tasks.hpp
+++ b/src/tasks/tasks.hpp
@@ -503,8 +503,10 @@ class TaskList {
     if (!label.has_value()) label = "anon";
 #ifdef HAS_CXX_ABI
     int status;
-    auto signature =
-        std::string(abi::__cxa_demangle(typeid(F).name(), NULL, NULL, &status));
+    // _cxa_demangle calls malloc so we must call free
+    char *signature_name = abi::__cxa_demangle(typeid(F).name(), NULL, NULL, &status);
+    auto signature = std::string(signature_name);
+    std::free(signature_name);
 #else
     std::string signature = " (...)";
 #endif

From 3bd60362d4d4f27aecdb84977981e22536a02675 Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Sat, 10 Aug 2024 17:02:43 -0600
Subject: [PATCH 17/40] Changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 89afa174910f..c404edb308a7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -35,6 +35,7 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1152]](https://github.com/parthenon-hpc-lab/parthenon/pull/1152) Fix memory leak in task graph outputs related to `abi::__cxa_demangle`
 - [[PR 1146]](https://github.com/parthenon-hpc-lab/parthenon/pull/1146) Fix an issue outputting >4GB single variables per rank
 - [[PR 1144]](https://github.com/parthenon-hpc-lab/parthenon/pull/1144) Fix some restarts w/non-CC fields
 - [[PR 1132]](https://github.com/parthenon-hpc-lab/parthenon/pull/1132) Fix regional dependencies for iterative task lists and make solvers work for arbirtrary MeshData partitioning

From 3027cd0a34dbc8293d535b9172ff131e268e77a5 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Mon, 12 Aug 2024 15:07:17 +0200
Subject: [PATCH 18/40] Fix allocation patterns

---
 src/bvals/comms/build_boundary_buffers.cpp | 55 +++++++++++++---------
 src/utils/object_pool.hpp                  |  2 +-
 2 files changed, 33 insertions(+), 24 deletions(-)

diff --git a/src/bvals/comms/build_boundary_buffers.cpp b/src/bvals/comms/build_boundary_buffers.cpp
index 44d2b1af9012..ca49d3210d0a 100644
--- a/src/bvals/comms/build_boundary_buffers.cpp
+++ b/src/bvals/comms/build_boundary_buffers.cpp
@@ -46,7 +46,8 @@ template <BoundaryType BTYPE>
 void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
                                Mesh::comm_buf_map_t &buf_map) {
   Mesh *pmesh = md->GetMeshPointer();
-  std::unordered_map<int, int> nbufs;
+  std::unordered_map<int, int>
+      nbufs; // total (existing and new) number of buffers for given size
 
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
@@ -54,11 +55,8 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
 
     nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
-    std::cerr << "\nnbufs[" << buf_size << "] += 1, which is now " << nbufs[buf_size];
   });
 
-  std::cerr << "\nbefore nbufs.size = " << nbufs.size();
-
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
     int buf_size = GetBufferSize(pmb, nb, v);
@@ -67,30 +65,41 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     // Add a buffer pool if one does not exist for this size
     using buf_t = buf_pool_t<Real>::base_t;
     if (pmesh->pool_map.count(buf_size) == 0) {
-      std::cerr << "\ninside nbufs.size = " << nbufs.size();
-      std::cerr << "\nGoing to request " << nbufs[buf_size] << " buffers of size "
-                << buf_size;
-      pmesh->pool_map.emplace(std::make_pair(
-          buf_size, buf_pool_t<Real>([buf_size, &nbufs](buf_pool_t<Real> *pool) {
-            const auto pool_size = static_cast<int64_t>(nbufs[buf_size]) * buf_size;
+      std::cerr << "Setting up a new pool for buffers of size " << buf_size << "\n";
+      // Might be worth discussing what a good default is.
+      // Using the number of packs, assumes that all blocks in a pack have fairly similar
+      // buffer configurations, which may or may not be a good approximation.
+      // An alternative would be "1", which would reduce the memory footprint, but
+      // increase the number of individual memory allocations.
+      const int64_t nbuf = pmesh->DefaultNumPartitions();
+      pmesh->pool_map.emplace(
+          buf_size, buf_pool_t<Real>([buf_size, nbuf](buf_pool_t<Real> *pool) {
+            std::cerr << "Dynamically adding " << nbuf << " buffers of size " << buf_size
+                      << " to a pool with current size " << pool->NumBuffersInPool()
+                      << " and future size " << pool->NumBuffersInPool() + nbuf << "\n";
+
+            const auto pool_size = nbuf * buf_size;
             buf_t chunk("pool buffer", pool_size);
-            for (int i = 1; i < nbufs[buf_size]; ++i) {
+            for (int i = 1; i < nbuf; ++i) {
               pool->AddFreeObjectToPool(
                   buf_t(chunk, std::make_pair(i * buf_size, (i + 1) * buf_size)));
             }
             return buf_t(chunk, std::make_pair(0, buf_size));
-          })));
-      // or add to existing pool (if required)
-    } else {
-      auto &pool = pmesh->pool_map.at(buf_size);
-      const auto new_buffers_req = nbufs[buf_size] - pool.NumAvailable();
-      if (new_buffers_req > 1) {
-        const auto pool_size = static_cast<int64_t>(new_buffers_req) * buf_size;
-        buf_t chunk("pool buffer", pool_size);
-        for (int i = 1; i < new_buffers_req; ++i) {
-          pool.AddFreeObjectToPool(
-              buf_t(chunk, std::make_pair(i * buf_size, (i + 1) * buf_size)));
-        }
+          }));
+    }
+    // Now that the pool is guaranteed to exist we can add free objects of the required
+    // amount.
+    auto &pool = pmesh->pool_map.at(buf_size);
+    const std::int64_t new_buffers_req = nbufs.at(buf_size) - pool.NumBuffersInPool();
+    if (new_buffers_req > 0) {
+      std::cerr << "Reserving " << new_buffers_req << " new buffers of size " << buf_size
+                << " to pool with " << pool.NumBuffersInPool() << " buffers because "
+                << nbufs.at(buf_size) << " are required in total.\n";
+      const auto pool_size = new_buffers_req * buf_size;
+      buf_t chunk("pool buffer", pool_size);
+      for (int i = 0; i < new_buffers_req; ++i) {
+        pool.AddFreeObjectToPool(
+            buf_t(chunk, std::make_pair(i * buf_size, (i + 1) * buf_size)));
       }
     }
 
diff --git a/src/utils/object_pool.hpp b/src/utils/object_pool.hpp
index e758a7867fe2..c7452499f126 100644
--- a/src/utils/object_pool.hpp
+++ b/src/utils/object_pool.hpp
@@ -62,7 +62,7 @@ class ObjectPool {
     std::cout << inuse_.size() << " used objects." << std::endl;
   }
 
-  auto NumAvailable() const { return available_.size(); }
+  auto NumBuffersInPool() const { return inuse_.size() + available_.size(); }
 
   std::uint64_t SizeInBytes() const {
     constexpr std::uint64_t datum_size = sizeof(typename base_t::value_type);

From 9aeec293698c4de4613e48b137c048bf452abf09 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Mon, 12 Aug 2024 21:08:33 -0600
Subject: [PATCH 19/40] segfaulting in the first step

---
 src/bvals/comms/boundary_communication.cpp | 14 +++--
 src/coordinates/uniform_cartesian.hpp      |  2 +-
 src/interface/mesh_data.hpp                |  1 +
 src/mesh/mesh-amr_loadbalance.cpp          | 64 ++++++++++++----------
 src/mesh/mesh.cpp                          | 29 +++++-----
 src/mesh/mesh.hpp                          |  2 +-
 src/mesh/meshblock.cpp                     | 42 ++++++++------
 src/mesh/meshblock.hpp                     |  2 +
 src/outputs/parthenon_hdf5.cpp             | 13 ++---
 src/utils/communication_buffer.hpp         |  1 +
 src/utils/partition_stl_containers.hpp     |  4 +-
 11 files changed, 98 insertions(+), 76 deletions(-)

diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp
index bd2dedd377aa..01e1b05233b8 100644
--- a/src/bvals/comms/boundary_communication.cpp
+++ b/src/bvals/comms/boundary_communication.cpp
@@ -79,10 +79,12 @@ TaskStatus SendBoundBufs(std::shared_ptr<MeshData<Real>> &md) {
     }
   }
   // Restrict
-  auto pmb = md->GetBlockData(0)->GetBlockPointer();
-  StateDescriptor *resolved_packages = pmb->resolved_packages.get();
-  refinement::Restrict(resolved_packages, cache.prores_cache, pmb->cellbounds,
-                       pmb->c_cellbounds);
+  if (md->NumBlocks() > 0) {
+    auto pmb = md->GetBlockData(0)->GetBlockPointer();
+    StateDescriptor *resolved_packages = pmb->resolved_packages.get();
+    refinement::Restrict(resolved_packages, cache.prores_cache, pmb->cellbounds,
+                         pmb->c_cellbounds);
+  }
 
   // Load buffer data
   auto &bnd_info = cache.bnd_info;
@@ -335,7 +337,7 @@ TaskStatus SetBounds(std::shared_ptr<MeshData<Real>> &md) {
 #endif
   std::for_each(std::begin(cache.buf_vec), std::end(cache.buf_vec),
                 [](auto pbuf) { pbuf->Stale(); });
-  if (nbound > 0 && pmesh->multilevel) {
+  if (nbound > 0 && pmesh->multilevel && md->NumBlocks() > 0) {
     // Restrict
     auto pmb = md->GetBlockData(0)->GetBlockPointer();
     StateDescriptor *resolved_packages = pmb->resolved_packages.get();
@@ -377,7 +379,7 @@ TaskStatus ProlongateBounds(std::shared_ptr<MeshData<Real>> &md) {
     }
   }
 
-  if (nbound > 0 && pmesh->multilevel) {
+  if (nbound > 0 && pmesh->multilevel  && md->NumBlocks() > 0) {
     auto pmb = md->GetBlockData(0)->GetBlockPointer();
     StateDescriptor *resolved_packages = pmb->resolved_packages.get();
 
diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp
index 04cec9d41974..2b37fcd7cc52 100644
--- a/src/coordinates/uniform_cartesian.hpp
+++ b/src/coordinates/uniform_cartesian.hpp
@@ -318,7 +318,7 @@ class UniformCartesian {
   const std::array<Real, 3> &GetXmin() const { return xmin_; }
   const std::array<int, 3> &GetStartIndex() const { return istart_; }
   const char *Name() const { return name_; }
-
+  static const char *StaticName() { return name_; }
  private:
   std::array<int, 3> istart_;
   std::array<Real, 3> xmin_, dx_, area_;
diff --git a/src/interface/mesh_data.hpp b/src/interface/mesh_data.hpp
index e00124132c81..621ffce9763f 100644
--- a/src/interface/mesh_data.hpp
+++ b/src/interface/mesh_data.hpp
@@ -519,6 +519,7 @@ class MeshData {
 
 template <typename T, typename... Args>
 std::vector<Uid_t> UidIntersection(MeshData<T> *md1, MeshData<T> *md2, Args &&...args) {
+  if (md1->NumBlocks() == 0 || md2->NumBlocks() == 0) return std::vector<Uid_t>();
   return UidIntersection(md1->GetBlockData(0).get(), md2->GetBlockData(0).get(),
                          std::forward<Args>(args)...);
 }
diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp
index 52749c2c18cf..22ca36c7772c 100644
--- a/src/mesh/mesh-amr_loadbalance.cpp
+++ b/src/mesh/mesh-amr_loadbalance.cpp
@@ -388,18 +388,16 @@ void AssignBlocks(std::vector<double> const &costlist, std::vector<int> &ranklis
 
 void UpdateBlockList(std::vector<int> const &ranklist, std::vector<int> &nslist,
                      std::vector<int> &nblist) {
-  nslist.resize(Globals::nranks);
-  nblist.resize(Globals::nranks);
-
-  nslist[0] = 0;
-  int rank = 0;
-  for (int block_id = 1; block_id < ranklist.size(); block_id++) {
-    if (ranklist[block_id] != ranklist[block_id - 1]) {
-      nblist[rank] = block_id - nslist[rank];
-      nslist[++rank] = block_id;
-    }
-  }
-  nblist[rank] = ranklist.size() - nslist[rank];
+  // First count the number of blocks on each rank
+  nblist = std::vector<int>(Globals::nranks, 0);
+  for (int b = ranklist.size() - 1; b >= 0; --b)
+    nblist[ranklist[b]]++; 
+  
+  // Then find the starting gid of the blocks, assuming they 
+  // are apportioned in increasing order
+  nslist = std::vector<int>(Globals::nranks, 0);
+  for (int b = 1; b < nslist.size(); ++b)
+    nslist[b] = nslist[b - 1] + nblist[b - 1]; 
 }
 } // namespace
 
@@ -424,6 +422,10 @@ void Mesh::CalculateLoadBalance(std::vector<double> const &costlist,
   // on each rank.
   UpdateBlockList(ranklist, nslist, nblist);
 
+  for (int i=0; i<ranklist.size(); ++i) { 
+    printf("%i %i\n", i, ranklist[i]);
+  }
+
 #ifdef MPI_PARALLEL
   if (total_blocks % (Globals::nranks) != 0 && !adaptive && !lb_flag_ &&
       maxcost == mincost && Globals::my_rank == 0) {
@@ -912,23 +914,25 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
       auto pmb = FindMeshBlock(nn);
       if (newloc[nn].level() > loclist[on].level()) nprolong += pmb->vars_cc_.size();
     }
-    prolongation_cache.Initialize(nprolong, resolved_packages.get());
-    int iprolong = 0;
-    for (int nn = nbs; nn <= nbe; nn++) {
-      int on = newtoold[nn];
-      if (newloc[nn].level() > loclist[on].level()) {
-        auto pmb = FindMeshBlock(nn);
-        for (auto &var : pmb->vars_cc_) {
-          prolongation_cache.RegisterRegionHost(
-              iprolong++,
-              ProResInfo::GetInteriorProlongate(pmb.get(), NeighborBlock(), var),
-              var.get(), resolved_packages.get());
+    if (nprolong > 0) {
+      prolongation_cache.Initialize(nprolong, resolved_packages.get());
+      int iprolong = 0;
+      for (int nn = nbs; nn <= nbe; nn++) {
+        int on = newtoold[nn];
+        if (newloc[nn].level() > loclist[on].level()) {
+          auto pmb = FindMeshBlock(nn);
+          for (auto &var : pmb->vars_cc_) {
+            prolongation_cache.RegisterRegionHost(
+                iprolong++,
+                ProResInfo::GetInteriorProlongate(pmb.get(), NeighborBlock(), var),
+                var.get(), resolved_packages.get());
+          }
         }
+        prolongation_cache.CopyToDevice();
       }
-      prolongation_cache.CopyToDevice();
+      refinement::ProlongateShared(resolved_packages.get(), prolongation_cache,
+                                   block_list[0]->cellbounds, block_list[0]->c_cellbounds);
     }
-    refinement::ProlongateShared(resolved_packages.get(), prolongation_cache,
-                                 block_list[0]->cellbounds, block_list[0]->c_cellbounds);
 
     // update the lists
     loclist = std::move(newloc);
@@ -978,9 +982,11 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
                                     // ghosts of non-cell centered vars may get some junk
       // Now there is the correct data for prolongating on un-shared topological elements
       // on the new fine blocks
-      refinement::ProlongateInternal(resolved_packages.get(), prolongation_cache,
-                                     block_list[0]->cellbounds,
-                                     block_list[0]->c_cellbounds);
+      if (block_list.size() > 0) {
+        refinement::ProlongateInternal(resolved_packages.get(), prolongation_cache,
+                                       block_list[0]->cellbounds,
+                                       block_list[0]->c_cellbounds);
+      }
     }
 
     // Rebuild just the ownership model, this time weighting the "new" fine blocks just
diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index cf2ee90bf2bb..cfe41af9351a 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -330,7 +330,7 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in,
       msg << "### FATAL ERROR in Mesh constructor" << std::endl
           << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks ("
           << Globals::nranks << ")" << std::endl;
-      PARTHENON_FAIL(msg);
+      PARTHENON_WARN(msg);
     } else { // test
       std::cout << "### Warning in Mesh constructor" << std::endl
                 << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks ("
@@ -356,6 +356,9 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in,
   // create MeshBlock list for this process
   int nbs = nslist[Globals::my_rank];
   int nbe = nbs + nblist[Globals::my_rank] - 1;
+
+  printf("rank: %i nbs: %i nbe: %i\n", Globals::my_rank, nbs, nbe);
+
   // create MeshBlock list for this process
   block_list.clear();
   block_list.resize(nbe - nbs + 1);
@@ -399,8 +402,8 @@ void Mesh::BuildBlockPartitions(GridIdentifier grid) {
       grid.type == GridType::leaf ? block_list : gmg_block_lists[grid.logical_level],
       DefaultPackSize());
   // Account for possibly empty block_list
-  if (partition_blocklists.size() == 0)
-    partition_blocklists = std::vector<BlockList_t>(1);
+  //if (partition_blocklists.size() == 0)
+  //  partition_blocklists = std::vector<BlockList_t>(1);
   std::vector<std::shared_ptr<BlockListPartition>> out;
   int id = 0;
   for (auto &part_bl : partition_blocklists)
@@ -725,6 +728,9 @@ void Mesh::PreCommFillDerived() {
   }
   for (auto &partition : GetDefaultBlockPartitions()) {
     auto &md = mesh_data.Add("base", partition);
+    printf("partition size = %i\n", partition->block_list.size());
+    PARTHENON_REQUIRE(partition->pmesh == this, "Bad partition mesh pointer");
+    PARTHENON_REQUIRE(md->GetParentPointer() == this, "Bad mesh pointer");
     Update::PreCommFillDerived(md.get());
   }
 }
@@ -853,6 +859,7 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput *
 
 /// Finds location of a block with ID `tgid`.
 std::shared_ptr<MeshBlock> Mesh::FindMeshBlock(int tgid) const {
+  PARTHENON_REQUIRE(block_list.size() > 0, "Trying to call FindMeshBlock with empty block list");
   // Attempt to simply index into the block list.
   const int nbs = block_list[0]->gid;
   const int i = tgid - nbs;
@@ -878,24 +885,20 @@ bool Mesh::SetBlockSizeAndBoundaries(LogicalLocation loc, RegionSize &block_size
 }
 
 std::int64_t Mesh::GetTotalCells() {
-  auto &pmb = block_list.front();
-  return static_cast<std::int64_t>(nbtotal) * pmb->block_size.nx(X1DIR) *
-         pmb->block_size.nx(X2DIR) * pmb->block_size.nx(X3DIR);
+  return static_cast<std::int64_t>(nbtotal) * GetNumberOfMeshBlockCells();
 }
 
-// TODO(JMM): Move block_size into mesh.
 int Mesh::GetNumberOfMeshBlockCells() const {
-  return block_list.front()->GetNumberOfMeshBlockCells();
+  return base_block_size.nx(X1DIR) * base_block_size.nx(X2DIR) * base_block_size.nx(X3DIR);
 }
 
-const IndexShape &Mesh::GetLeafBlockCellBounds(CellLevel level) const {
-  MeshBlock *pmb = block_list[0].get();
+const IndexShape Mesh::GetLeafBlockCellBounds(CellLevel level) const {
   if (level == CellLevel::same) {
-    return pmb->cellbounds;
+    return std::get<0>(GetIndexShapes({base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)}, multilevel, this));
   } else if (level == CellLevel::fine) {
-    return pmb->f_cellbounds;
+    return std::get<1>(GetIndexShapes({base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)}, multilevel, this));
   } else { // if (level == CellLevel::coarse) {
-    return pmb->c_cellbounds;
+    return std::get<2>(GetIndexShapes({base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)}, multilevel, this));
   }
 }
 
diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp
index 2872a90cf864..ee7c06140b26 100644
--- a/src/mesh/mesh.hpp
+++ b/src/mesh/mesh.hpp
@@ -106,7 +106,7 @@ class Mesh {
   RegionSize GetBlockSize(const LogicalLocation &loc) const {
     return forest.GetBlockDomain(loc);
   }
-  const IndexShape &GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const;
+  const IndexShape GetLeafBlockCellBounds(CellLevel level = CellLevel::same) const;
 
   const forest::Forest &Forest() const { return forest; }
 
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index b36fb6489e91..ad1831121b56 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -201,27 +201,33 @@ void MeshBlock::Initialize(int igid, int ilid, LogicalLocation iloc,
 
 MeshBlock::~MeshBlock() = default;
 
+std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel, const Mesh *pmesh) {
+  IndexShape cellbounds(nx[2], nx[1], nx[0], Globals::nghost);
+  IndexShape f_cellbounds(2 * nx[2], 2 * nx[1], 2 * nx[0], Globals::nghost);
+  IndexShape c_cellbounds(nx[2] / 2, nx[1] / 2, nx[0] / 2, 0);
+  if (multilevel) {
+    // Prevent the coarse bounds from going to zero
+      int cnx1 = nx[0] / 2;
+      int cnx2 = nx[1] / 2;
+      int cnx3 = nx[2] / 2;
+      if (pmesh != nullptr) {
+        cnx1 = pmesh->mesh_size.symmetry(X1DIR) ? 0 : std::max(1, nx[0] / 2);
+        cnx2 = pmesh->mesh_size.symmetry(X2DIR) ? 0 : std::max(1, nx[1] / 2);
+        cnx3 = pmesh->mesh_size.symmetry(X3DIR) ? 0 : std::max(1, nx[2] / 2);
+      }
+      c_cellbounds = IndexShape(cnx3, cnx2, cnx1, Globals::nghost); 
+  }
+  return {cellbounds, f_cellbounds, c_cellbounds};
+}
+
 void MeshBlock::InitializeIndexShapesImpl(const int nx1, const int nx2, const int nx3,
                                           bool init_coarse, bool multilevel) {
-  cellbounds = IndexShape(nx3, nx2, nx1, Globals::nghost);
-  f_cellbounds = IndexShape(2 * nx3, 2 * nx2, 2 * nx1, Globals::nghost);
-
+  auto [cb, fcb, ccb] = GetIndexShapes({nx1, nx2, nx3}, multilevel, pmy_mesh); 
+  cellbounds = cb; 
+  f_cellbounds = fcb; 
   if (init_coarse) {
-    if (multilevel) {
-      // Prevent the coarse bounds from going to zero
-      int cnx1 = nx1 / 2;
-      int cnx2 = nx2 / 2;
-      int cnx3 = nx3 / 2;
-      if (pmy_mesh != nullptr) {
-        cnx1 = pmy_mesh->mesh_size.symmetry(X1DIR) ? 0 : std::max(1, nx1 / 2);
-        cnx2 = pmy_mesh->mesh_size.symmetry(X2DIR) ? 0 : std::max(1, nx2 / 2);
-        cnx3 = pmy_mesh->mesh_size.symmetry(X3DIR) ? 0 : std::max(1, nx3 / 2);
-      }
-      cnghost = (Globals::nghost + 1) / 2 + 1;
-      c_cellbounds = IndexShape(cnx3, cnx2, cnx1, Globals::nghost);
-    } else {
-      c_cellbounds = IndexShape(nx3 / 2, nx2 / 2, nx1 / 2, 0);
-    }
+    cnghost = (Globals::nghost + 1) / 2 + 1;
+    c_cellbounds = ccb;
   }
 }
 
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index 41ae52ea5930..1ed68455cb7e 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -63,6 +63,8 @@ KOKKOS_FORCEINLINE_FUNCTION void par_for_inner(const team_mbr_t &team_member,
   parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, team_member, il, iu, function);
 }
 
+std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel, const Mesh *pmesh);
+
 //----------------------------------------------------------------------------------------
 //! \class MeshBlock
 //! \brief data/functions associated with a single block
diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp
index 7872e2ad4d03..ed3818b11af1 100644
--- a/src/outputs/parthenon_hdf5.cpp
+++ b/src/outputs/parthenon_hdf5.cpp
@@ -82,10 +82,8 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
   const IndexDomain theDomain =
       (output_params.include_ghost_zones ? IndexDomain::entire : IndexDomain::interior);
 
-  auto const &first_block = *(pm->block_list.front());
-
-  const auto &cellbounds = first_block.cellbounds;
-  const auto &f_cellbounds = first_block.f_cellbounds;
+  const auto &cellbounds = pm->GetLeafBlockCellBounds(CellLevel::same);
+  const auto &f_cellbounds = pm->GetLeafBlockCellBounds(CellLevel::coarse);
 
   auto const nx1 = cellbounds.ncellsi(theDomain);
   auto const nx2 = cellbounds.ncellsj(theDomain);
@@ -153,7 +151,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
                        info_group);
     // write number of ghost cells in simulation
     HDF5WriteAttribute("NGhost", Globals::nghost, info_group);
-    HDF5WriteAttribute("Coordinates", std::string(first_block.coords.Name()).c_str(),
+    HDF5WriteAttribute("Coordinates", std::string(Coordinates_t::StaticName()).c_str(),
                        info_group);
 
     // restart info, write always
@@ -267,8 +265,9 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
 
   // get list of all vars, just use the first block since the list is
   // the same for all blocks
-  auto all_vars_info =
-      VarInfo::GetAll(get_vars(pm->block_list.front()), cellbounds, f_cellbounds);
+  auto all_vars_info = 
+      VarInfo::GetAll(pm->block_list.size() > 0 ? get_vars(pm->block_list.front()) : VariableVector<Real>(), cellbounds, f_cellbounds);
+
 
   // We need to add information about the sparse variables to the HDF5 file, namely:
   // 1) Which variables are sparse
diff --git a/src/utils/communication_buffer.hpp b/src/utils/communication_buffer.hpp
index edfc04f11c77..9f5e32bf050f 100644
--- a/src/utils/communication_buffer.hpp
+++ b/src/utils/communication_buffer.hpp
@@ -275,6 +275,7 @@ bool CommBuffer<T>::IsAvailableForWrite() {
   } else if (*comm_type_ == BuffCommType::both) {
     return (*state_ == BufferState::stale);
   } else {
+    printf("send_rank: %i recv_rank: %i\n", send_rank_, recv_rank_);
     PARTHENON_FAIL("Receiving buffer is never available for write.");
   }
 }
diff --git a/src/utils/partition_stl_containers.hpp b/src/utils/partition_stl_containers.hpp
index 8c943c2368e2..4057c5170001 100644
--- a/src/utils/partition_stl_containers.hpp
+++ b/src/utils/partition_stl_containers.hpp
@@ -47,8 +47,10 @@ auto ToSizeN(Container_t<T, extra...> &container, const int N) {
   using std::to_string;
   using namespace partition_impl;
 
+  if (N == 0) return Partition_t<T>();
+  
   PARTHENON_REQUIRE_THROWS(N > 0, "Your partition must be at least size 1");
-
+  
   int nelements = container.size();
   int npartitions = IntCeil(nelements, N);
 

From 65bdeaae4230c7ecb08f16982872433af046229d Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 11:10:48 -0600
Subject: [PATCH 20/40] allow for initialization with fewer ranks

---
 src/mesh/mesh-amr_loadbalance.cpp | 45 +++++++++++++++----------------
 src/mesh/mesh.cpp                 | 24 ++++++++---------
 2 files changed, 33 insertions(+), 36 deletions(-)

diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp
index 22ca36c7772c..983ff940bcb0 100644
--- a/src/mesh/mesh-amr_loadbalance.cpp
+++ b/src/mesh/mesh-amr_loadbalance.cpp
@@ -422,10 +422,6 @@ void Mesh::CalculateLoadBalance(std::vector<double> const &costlist,
   // on each rank.
   UpdateBlockList(ranklist, nslist, nblist);
 
-  for (int i=0; i<ranklist.size(); ++i) { 
-    printf("%i %i\n", i, ranklist[i]);
-  }
-
 #ifdef MPI_PARALLEL
   if (total_blocks % (Globals::nranks) != 0 && !adaptive && !lb_flag_ &&
       maxcost == mincost && Globals::my_rank == 0) {
@@ -737,30 +733,31 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
   int nbe = nbs + nblist[Globals::my_rank] - 1;
 
   // Restrict fine to coarse buffers
-  ProResCache_t restriction_cache;
-  int nrestrict = 0;
-  for (int on = onbs; on <= onbe; on++) {
-    int nn = oldtonew[on];
-    auto pmb = FindMeshBlock(on);
-    if (newloc[nn].level() < loclist[on].level()) nrestrict += pmb->vars_cc_.size();
-  }
-  restriction_cache.Initialize(nrestrict, resolved_packages.get());
-  int irestrict = 0;
-  for (int on = onbs; on <= onbe; on++) {
-    int nn = oldtonew[on];
-    if (newloc[nn].level() < loclist[on].level()) {
+  if (block_list.size() > 0) {
+    ProResCache_t restriction_cache;
+    int nrestrict = 0;
+    for (int on = onbs; on <= onbe; on++) {
+      int nn = oldtonew[on];
       auto pmb = FindMeshBlock(on);
-      for (auto &var : pmb->vars_cc_) {
-        restriction_cache.RegisterRegionHost(
-            irestrict++, ProResInfo::GetInteriorRestrict(pmb.get(), NeighborBlock(), var),
-            var.get(), resolved_packages.get());
+      if (newloc[nn].level() < loclist[on].level()) nrestrict += pmb->vars_cc_.size();
+    }
+    restriction_cache.Initialize(nrestrict, resolved_packages.get());
+    int irestrict = 0;
+    for (int on = onbs; on <= onbe; on++) {
+      int nn = oldtonew[on];
+      if (newloc[nn].level() < loclist[on].level()) {
+        auto pmb = FindMeshBlock(on);
+        for (auto &var : pmb->vars_cc_) {
+          restriction_cache.RegisterRegionHost(
+              irestrict++, ProResInfo::GetInteriorRestrict(pmb.get(), NeighborBlock(), var),
+              var.get(), resolved_packages.get());
+        }
       }
     }
+    restriction_cache.CopyToDevice();
+    refinement::Restrict(resolved_packages.get(), restriction_cache,
+                         block_list[0]->cellbounds, block_list[0]->c_cellbounds);
   }
-  restriction_cache.CopyToDevice();
-  refinement::Restrict(resolved_packages.get(), restriction_cache,
-                       block_list[0]->cellbounds, block_list[0]->c_cellbounds);
-
   Kokkos::fence();
 
 #ifdef MPI_PARALLEL
diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index cfe41af9351a..1ed2a5b700a5 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -326,12 +326,7 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in,
 #ifdef MPI_PARALLEL
   // check if there are sufficient blocks
   if (nbtotal < Globals::nranks) {
-    if (mesh_test == 0) {
-      msg << "### FATAL ERROR in Mesh constructor" << std::endl
-          << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks ("
-          << Globals::nranks << ")" << std::endl;
-      PARTHENON_WARN(msg);
-    } else { // test
+    if (mesh_test != 0) {
       std::cout << "### Warning in Mesh constructor" << std::endl
                 << "Too few mesh blocks: nbtotal (" << nbtotal << ") < nranks ("
                 << Globals::nranks << ")" << std::endl;
@@ -357,8 +352,6 @@ void Mesh::BuildBlockList(ParameterInput *pin, ApplicationInput *app_in,
   int nbs = nslist[Globals::my_rank];
   int nbe = nbs + nblist[Globals::my_rank] - 1;
 
-  printf("rank: %i nbs: %i nbe: %i\n", Globals::my_rank, nbs, nbe);
-
   // create MeshBlock list for this process
   block_list.clear();
   block_list.resize(nbe - nbs + 1);
@@ -401,9 +394,6 @@ void Mesh::BuildBlockPartitions(GridIdentifier grid) {
   auto partition_blocklists = partition::ToSizeN(
       grid.type == GridType::leaf ? block_list : gmg_block_lists[grid.logical_level],
       DefaultPackSize());
-  // Account for possibly empty block_list
-  //if (partition_blocklists.size() == 0)
-  //  partition_blocklists = std::vector<BlockList_t>(1);
   std::vector<std::shared_ptr<BlockListPartition>> out;
   int id = 0;
   for (auto &part_bl : partition_blocklists)
@@ -728,7 +718,6 @@ void Mesh::PreCommFillDerived() {
   }
   for (auto &partition : GetDefaultBlockPartitions()) {
     auto &md = mesh_data.Add("base", partition);
-    printf("partition size = %i\n", partition->block_list.size());
     PARTHENON_REQUIRE(partition->pmesh == this, "Bad partition mesh pointer");
     PARTHENON_REQUIRE(md->GetParentPointer() == this, "Bad mesh pointer");
     Update::PreCommFillDerived(md.get());
@@ -853,6 +842,17 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput *
     }
   } while (!init_done);
 
+#ifdef MPI_PARALLEL
+  // check if there are sufficient blocks
+  if (nbtotal < Globals::nranks) {
+    std::stringstream msg;
+    msg << "### FATAL ERROR in Mesh Initialize" << std::endl
+        << "Too few mesh blocks after initialization: nbtotal (" << nbtotal << ") < nranks ("
+        << Globals::nranks << ")" << std::endl;
+    PARTHENON_FAIL(msg);
+  }
+#endif
+
   // Initialize the "base" MeshData object
   mesh_data.Get()->Initialize(block_list, this);
 }

From 5cf2ab8b9c94a9b6239afb4f7f08f5c0b7285075 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 11:17:19 -0600
Subject: [PATCH 21/40] format and lint

---
 src/bvals/comms/boundary_communication.cpp |  2 +-
 src/coordinates/uniform_cartesian.hpp      |  1 +
 src/mesh/mesh-amr_loadbalance.cpp          | 16 +++++++------
 src/mesh/mesh.cpp                          | 22 ++++++++++++------
 src/mesh/meshblock.cpp                     | 27 +++++++++++-----------
 src/mesh/meshblock.hpp                     |  3 ++-
 src/outputs/parthenon_hdf5.cpp             |  7 +++---
 src/utils/partition_stl_containers.hpp     |  4 ++--
 8 files changed, 48 insertions(+), 34 deletions(-)

diff --git a/src/bvals/comms/boundary_communication.cpp b/src/bvals/comms/boundary_communication.cpp
index 01e1b05233b8..1f238b1fdbf7 100644
--- a/src/bvals/comms/boundary_communication.cpp
+++ b/src/bvals/comms/boundary_communication.cpp
@@ -379,7 +379,7 @@ TaskStatus ProlongateBounds(std::shared_ptr<MeshData<Real>> &md) {
     }
   }
 
-  if (nbound > 0 && pmesh->multilevel  && md->NumBlocks() > 0) {
+  if (nbound > 0 && pmesh->multilevel && md->NumBlocks() > 0) {
     auto pmb = md->GetBlockData(0)->GetBlockPointer();
     StateDescriptor *resolved_packages = pmb->resolved_packages.get();
 
diff --git a/src/coordinates/uniform_cartesian.hpp b/src/coordinates/uniform_cartesian.hpp
index 2b37fcd7cc52..90e98447f0df 100644
--- a/src/coordinates/uniform_cartesian.hpp
+++ b/src/coordinates/uniform_cartesian.hpp
@@ -319,6 +319,7 @@ class UniformCartesian {
   const std::array<int, 3> &GetStartIndex() const { return istart_; }
   const char *Name() const { return name_; }
   static const char *StaticName() { return name_; }
+
  private:
   std::array<int, 3> istart_;
   std::array<Real, 3> xmin_, dx_, area_;
diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp
index 983ff940bcb0..d3a4e2d9ed4a 100644
--- a/src/mesh/mesh-amr_loadbalance.cpp
+++ b/src/mesh/mesh-amr_loadbalance.cpp
@@ -391,13 +391,13 @@ void UpdateBlockList(std::vector<int> const &ranklist, std::vector<int> &nslist,
   // First count the number of blocks on each rank
   nblist = std::vector<int>(Globals::nranks, 0);
   for (int b = ranklist.size() - 1; b >= 0; --b)
-    nblist[ranklist[b]]++; 
-  
-  // Then find the starting gid of the blocks, assuming they 
+    nblist[ranklist[b]]++;
+
+  // Then find the starting gid of the blocks, assuming they
   // are apportioned in increasing order
   nslist = std::vector<int>(Globals::nranks, 0);
   for (int b = 1; b < nslist.size(); ++b)
-    nslist[b] = nslist[b - 1] + nblist[b - 1]; 
+    nslist[b] = nslist[b - 1] + nblist[b - 1];
 }
 } // namespace
 
@@ -749,8 +749,9 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
         auto pmb = FindMeshBlock(on);
         for (auto &var : pmb->vars_cc_) {
           restriction_cache.RegisterRegionHost(
-              irestrict++, ProResInfo::GetInteriorRestrict(pmb.get(), NeighborBlock(), var),
-              var.get(), resolved_packages.get());
+              irestrict++,
+              ProResInfo::GetInteriorRestrict(pmb.get(), NeighborBlock(), var), var.get(),
+              resolved_packages.get());
         }
       }
     }
@@ -928,7 +929,8 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
         prolongation_cache.CopyToDevice();
       }
       refinement::ProlongateShared(resolved_packages.get(), prolongation_cache,
-                                   block_list[0]->cellbounds, block_list[0]->c_cellbounds);
+                                   block_list[0]->cellbounds,
+                                   block_list[0]->c_cellbounds);
     }
 
     // update the lists
diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index 1ed2a5b700a5..0675d648bac2 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -847,8 +847,8 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput *
   if (nbtotal < Globals::nranks) {
     std::stringstream msg;
     msg << "### FATAL ERROR in Mesh Initialize" << std::endl
-        << "Too few mesh blocks after initialization: nbtotal (" << nbtotal << ") < nranks ("
-        << Globals::nranks << ")" << std::endl;
+        << "Too few mesh blocks after initialization: nbtotal (" << nbtotal
+        << ") < nranks (" << Globals::nranks << ")" << std::endl;
     PARTHENON_FAIL(msg);
   }
 #endif
@@ -859,7 +859,8 @@ void Mesh::Initialize(bool init_problem, ParameterInput *pin, ApplicationInput *
 
 /// Finds location of a block with ID `tgid`.
 std::shared_ptr<MeshBlock> Mesh::FindMeshBlock(int tgid) const {
-  PARTHENON_REQUIRE(block_list.size() > 0, "Trying to call FindMeshBlock with empty block list");
+  PARTHENON_REQUIRE(block_list.size() > 0,
+                    "Trying to call FindMeshBlock with empty block list");
   // Attempt to simply index into the block list.
   const int nbs = block_list[0]->gid;
   const int i = tgid - nbs;
@@ -889,16 +890,23 @@ std::int64_t Mesh::GetTotalCells() {
 }
 
 int Mesh::GetNumberOfMeshBlockCells() const {
-  return base_block_size.nx(X1DIR) * base_block_size.nx(X2DIR) * base_block_size.nx(X3DIR);
+  return base_block_size.nx(X1DIR) * base_block_size.nx(X2DIR) *
+         base_block_size.nx(X3DIR);
 }
 
 const IndexShape Mesh::GetLeafBlockCellBounds(CellLevel level) const {
   if (level == CellLevel::same) {
-    return std::get<0>(GetIndexShapes({base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)}, multilevel, this));
+    return std::get<0>(GetIndexShapes(
+        {base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)},
+        multilevel, this));
   } else if (level == CellLevel::fine) {
-    return std::get<1>(GetIndexShapes({base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)}, multilevel, this));
+    return std::get<1>(GetIndexShapes(
+        {base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)},
+        multilevel, this));
   } else { // if (level == CellLevel::coarse) {
-    return std::get<2>(GetIndexShapes({base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)}, multilevel, this));
+    return std::get<2>(GetIndexShapes(
+        {base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)},
+        multilevel, this));
   }
 }
 
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index ad1831121b56..935405ae614b 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -201,30 +201,31 @@ void MeshBlock::Initialize(int igid, int ilid, LogicalLocation iloc,
 
 MeshBlock::~MeshBlock() = default;
 
-std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel, const Mesh *pmesh) {
+std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel,
+                                         const Mesh *pmesh) {
   IndexShape cellbounds(nx[2], nx[1], nx[0], Globals::nghost);
   IndexShape f_cellbounds(2 * nx[2], 2 * nx[1], 2 * nx[0], Globals::nghost);
   IndexShape c_cellbounds(nx[2] / 2, nx[1] / 2, nx[0] / 2, 0);
   if (multilevel) {
     // Prevent the coarse bounds from going to zero
-      int cnx1 = nx[0] / 2;
-      int cnx2 = nx[1] / 2;
-      int cnx3 = nx[2] / 2;
-      if (pmesh != nullptr) {
-        cnx1 = pmesh->mesh_size.symmetry(X1DIR) ? 0 : std::max(1, nx[0] / 2);
-        cnx2 = pmesh->mesh_size.symmetry(X2DIR) ? 0 : std::max(1, nx[1] / 2);
-        cnx3 = pmesh->mesh_size.symmetry(X3DIR) ? 0 : std::max(1, nx[2] / 2);
-      }
-      c_cellbounds = IndexShape(cnx3, cnx2, cnx1, Globals::nghost); 
+    int cnx1 = nx[0] / 2;
+    int cnx2 = nx[1] / 2;
+    int cnx3 = nx[2] / 2;
+    if (pmesh != nullptr) {
+      cnx1 = pmesh->mesh_size.symmetry(X1DIR) ? 0 : std::max(1, nx[0] / 2);
+      cnx2 = pmesh->mesh_size.symmetry(X2DIR) ? 0 : std::max(1, nx[1] / 2);
+      cnx3 = pmesh->mesh_size.symmetry(X3DIR) ? 0 : std::max(1, nx[2] / 2);
+    }
+    c_cellbounds = IndexShape(cnx3, cnx2, cnx1, Globals::nghost);
   }
   return {cellbounds, f_cellbounds, c_cellbounds};
 }
 
 void MeshBlock::InitializeIndexShapesImpl(const int nx1, const int nx2, const int nx3,
                                           bool init_coarse, bool multilevel) {
-  auto [cb, fcb, ccb] = GetIndexShapes({nx1, nx2, nx3}, multilevel, pmy_mesh); 
-  cellbounds = cb; 
-  f_cellbounds = fcb; 
+  auto [cb, fcb, ccb] = GetIndexShapes({nx1, nx2, nx3}, multilevel, pmy_mesh);
+  cellbounds = cb;
+  f_cellbounds = fcb;
   if (init_coarse) {
     cnghost = (Globals::nghost + 1) / 2 + 1;
     c_cellbounds = ccb;
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index 1ed68455cb7e..4c36f8396d58 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -63,7 +63,8 @@ KOKKOS_FORCEINLINE_FUNCTION void par_for_inner(const team_mbr_t &team_member,
   parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, team_member, il, iu, function);
 }
 
-std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel, const Mesh *pmesh);
+std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel,
+                                         const Mesh *pmesh);
 
 //----------------------------------------------------------------------------------------
 //! \class MeshBlock
diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp
index ed3818b11af1..a1afb133d18a 100644
--- a/src/outputs/parthenon_hdf5.cpp
+++ b/src/outputs/parthenon_hdf5.cpp
@@ -265,9 +265,10 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
 
   // get list of all vars, just use the first block since the list is
   // the same for all blocks
-  auto all_vars_info = 
-      VarInfo::GetAll(pm->block_list.size() > 0 ? get_vars(pm->block_list.front()) : VariableVector<Real>(), cellbounds, f_cellbounds);
-
+  auto all_vars_info =
+      VarInfo::GetAll(pm->block_list.size() > 0 ? get_vars(pm->block_list.front())
+                                                : VariableVector<Real>(),
+                      cellbounds, f_cellbounds);
 
   // We need to add information about the sparse variables to the HDF5 file, namely:
   // 1) Which variables are sparse
diff --git a/src/utils/partition_stl_containers.hpp b/src/utils/partition_stl_containers.hpp
index 4057c5170001..8f1e1f8e90cc 100644
--- a/src/utils/partition_stl_containers.hpp
+++ b/src/utils/partition_stl_containers.hpp
@@ -48,9 +48,9 @@ auto ToSizeN(Container_t<T, extra...> &container, const int N) {
   using namespace partition_impl;
 
   if (N == 0) return Partition_t<T>();
-  
+
   PARTHENON_REQUIRE_THROWS(N > 0, "Your partition must be at least size 1");
-  
+
   int nelements = container.size();
   int npartitions = IntCeil(nelements, N);
 

From 480df064b4565fd4f93f7f6b57da0e2f2d5c96fd Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 11:22:13 -0600
Subject: [PATCH 22/40] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c404edb308a7..036f336c58a7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -29,6 +29,7 @@
 - [[PR 1019]](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Enable output for non-cell-centered variables
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 1153]](https://github.com/parthenon-hpc-lab/parthenon/pull/1153) Allow base grid with fewer blocks than ranks before initial AMR
 - [[PR 1105]](https://github.com/parthenon-hpc-lab/parthenon/pull/1105) Refactor parameter input for linear solvers
 - [[PR 1078]](https://github.com/parthenon-hpc-lab/parthenon/pull/1078) Add reduction fallback in 1D. Add IndexRange overload for 1D par loops
 - [[PR 1024]](https://github.com/parthenon-hpc-lab/parthenon/pull/1024) Add .outN. to history output filenames

From 60c6ac27d11cd386a871b24ec3bf2c198e0e56c7 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 11:24:37 -0600
Subject: [PATCH 23/40] small

---
 src/utils/communication_buffer.hpp | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/utils/communication_buffer.hpp b/src/utils/communication_buffer.hpp
index 9f5e32bf050f..edfc04f11c77 100644
--- a/src/utils/communication_buffer.hpp
+++ b/src/utils/communication_buffer.hpp
@@ -275,7 +275,6 @@ bool CommBuffer<T>::IsAvailableForWrite() {
   } else if (*comm_type_ == BuffCommType::both) {
     return (*state_ == BufferState::stale);
   } else {
-    printf("send_rank: %i recv_rank: %i\n", send_rank_, recv_rank_);
     PARTHENON_FAIL("Receiving buffer is never available for write.");
   }
 }

From c90a2e2f8da8317b82969865685c87da0cb7a9ed Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 12:55:04 -0600
Subject: [PATCH 24/40] revert hdf5 output

---
 src/outputs/parthenon_hdf5.cpp | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp
index a1afb133d18a..7872e2ad4d03 100644
--- a/src/outputs/parthenon_hdf5.cpp
+++ b/src/outputs/parthenon_hdf5.cpp
@@ -82,8 +82,10 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
   const IndexDomain theDomain =
       (output_params.include_ghost_zones ? IndexDomain::entire : IndexDomain::interior);
 
-  const auto &cellbounds = pm->GetLeafBlockCellBounds(CellLevel::same);
-  const auto &f_cellbounds = pm->GetLeafBlockCellBounds(CellLevel::coarse);
+  auto const &first_block = *(pm->block_list.front());
+
+  const auto &cellbounds = first_block.cellbounds;
+  const auto &f_cellbounds = first_block.f_cellbounds;
 
   auto const nx1 = cellbounds.ncellsi(theDomain);
   auto const nx2 = cellbounds.ncellsj(theDomain);
@@ -151,7 +153,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
                        info_group);
     // write number of ghost cells in simulation
     HDF5WriteAttribute("NGhost", Globals::nghost, info_group);
-    HDF5WriteAttribute("Coordinates", std::string(Coordinates_t::StaticName()).c_str(),
+    HDF5WriteAttribute("Coordinates", std::string(first_block.coords.Name()).c_str(),
                        info_group);
 
     // restart info, write always
@@ -266,9 +268,7 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
   // get list of all vars, just use the first block since the list is
   // the same for all blocks
   auto all_vars_info =
-      VarInfo::GetAll(pm->block_list.size() > 0 ? get_vars(pm->block_list.front())
-                                                : VariableVector<Real>(),
-                      cellbounds, f_cellbounds);
+      VarInfo::GetAll(get_vars(pm->block_list.front()), cellbounds, f_cellbounds);
 
   // We need to add information about the sparse variables to the HDF5 file, namely:
   // 1) Which variables are sparse

From 3338dc21b138afae66045d99439bbfdf667db6ee Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 14:20:59 -0600
Subject: [PATCH 25/40] fix bug

---
 src/mesh/mesh-amr_loadbalance.cpp |  2 +-
 src/mesh/mesh.cpp                 | 15 ++++++---------
 src/mesh/meshblock.cpp            | 24 ++++++++++++------------
 src/mesh/meshblock.hpp            |  4 ++--
 4 files changed, 21 insertions(+), 24 deletions(-)

diff --git a/src/mesh/mesh-amr_loadbalance.cpp b/src/mesh/mesh-amr_loadbalance.cpp
index d3a4e2d9ed4a..fdb454fb2a34 100644
--- a/src/mesh/mesh-amr_loadbalance.cpp
+++ b/src/mesh/mesh-amr_loadbalance.cpp
@@ -981,7 +981,7 @@ void Mesh::RedistributeAndRefineMeshBlocks(ParameterInput *pin, ApplicationInput
                                     // ghosts of non-cell centered vars may get some junk
       // Now there is the correct data for prolongating on un-shared topological elements
       // on the new fine blocks
-      if (block_list.size() > 0) {
+      if (nprolong > 0) {
         refinement::ProlongateInternal(resolved_packages.get(), prolongation_cache,
                                        block_list[0]->cellbounds,
                                        block_list[0]->c_cellbounds);
diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index 0675d648bac2..2672155e83dc 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -895,18 +895,15 @@ int Mesh::GetNumberOfMeshBlockCells() const {
 }
 
 const IndexShape Mesh::GetLeafBlockCellBounds(CellLevel level) const {
+  auto shapes = GetIndexShapes(
+      ndim > 0 ? base_block_size.nx(X1DIR) : 0, ndim > 1 ? base_block_size.nx(X2DIR) : 0,
+      ndim > 2 ? base_block_size.nx(X3DIR) : 0, multilevel, this);
   if (level == CellLevel::same) {
-    return std::get<0>(GetIndexShapes(
-        {base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)},
-        multilevel, this));
+    return shapes[0];
   } else if (level == CellLevel::fine) {
-    return std::get<1>(GetIndexShapes(
-        {base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)},
-        multilevel, this));
+    return shapes[1];
   } else { // if (level == CellLevel::coarse) {
-    return std::get<2>(GetIndexShapes(
-        {base_block_size.nx(X1DIR), base_block_size.nx(X2DIR), base_block_size.nx(X3DIR)},
-        multilevel, this));
+    return shapes[2];
   }
 }
 
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index 935405ae614b..3f84ce48293b 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -201,20 +201,20 @@ void MeshBlock::Initialize(int igid, int ilid, LogicalLocation iloc,
 
 MeshBlock::~MeshBlock() = default;
 
-std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel,
-                                         const Mesh *pmesh) {
-  IndexShape cellbounds(nx[2], nx[1], nx[0], Globals::nghost);
-  IndexShape f_cellbounds(2 * nx[2], 2 * nx[1], 2 * nx[0], Globals::nghost);
-  IndexShape c_cellbounds(nx[2] / 2, nx[1] / 2, nx[0] / 2, 0);
+std::array<IndexShape, 3> GetIndexShapes(const int nx1, const int nx2, const int nx3,
+                                         bool multilevel, const Mesh *pmesh) {
+  IndexShape cellbounds(nx3, nx2, nx1, Globals::nghost);
+  IndexShape f_cellbounds(2 * nx3, 2 * nx2, 2 * nx1, Globals::nghost);
+  IndexShape c_cellbounds(nx3 / 2, nx2 / 2, nx1 / 2, 0);
   if (multilevel) {
     // Prevent the coarse bounds from going to zero
-    int cnx1 = nx[0] / 2;
-    int cnx2 = nx[1] / 2;
-    int cnx3 = nx[2] / 2;
+    int cnx1 = nx1 / 2;
+    int cnx2 = nx2 / 2;
+    int cnx3 = nx3 / 2;
     if (pmesh != nullptr) {
-      cnx1 = pmesh->mesh_size.symmetry(X1DIR) ? 0 : std::max(1, nx[0] / 2);
-      cnx2 = pmesh->mesh_size.symmetry(X2DIR) ? 0 : std::max(1, nx[1] / 2);
-      cnx3 = pmesh->mesh_size.symmetry(X3DIR) ? 0 : std::max(1, nx[2] / 2);
+      cnx1 = pmesh->mesh_size.symmetry(X1DIR) ? 0 : std::max(1, nx1 / 2);
+      cnx2 = pmesh->mesh_size.symmetry(X2DIR) ? 0 : std::max(1, nx2 / 2);
+      cnx3 = pmesh->mesh_size.symmetry(X3DIR) ? 0 : std::max(1, nx3 / 2);
     }
     c_cellbounds = IndexShape(cnx3, cnx2, cnx1, Globals::nghost);
   }
@@ -223,7 +223,7 @@ std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel,
 
 void MeshBlock::InitializeIndexShapesImpl(const int nx1, const int nx2, const int nx3,
                                           bool init_coarse, bool multilevel) {
-  auto [cb, fcb, ccb] = GetIndexShapes({nx1, nx2, nx3}, multilevel, pmy_mesh);
+  auto [cb, fcb, ccb] = GetIndexShapes(nx1, nx2, nx3, multilevel, pmy_mesh);
   cellbounds = cb;
   f_cellbounds = fcb;
   if (init_coarse) {
diff --git a/src/mesh/meshblock.hpp b/src/mesh/meshblock.hpp
index 4c36f8396d58..60ce4d86282d 100644
--- a/src/mesh/meshblock.hpp
+++ b/src/mesh/meshblock.hpp
@@ -63,8 +63,8 @@ KOKKOS_FORCEINLINE_FUNCTION void par_for_inner(const team_mbr_t &team_member,
   parthenon::par_for_inner(DEFAULT_INNER_LOOP_PATTERN, team_member, il, iu, function);
 }
 
-std::array<IndexShape, 3> GetIndexShapes(std::array<int, 3> nx, bool multilevel,
-                                         const Mesh *pmesh);
+std::array<IndexShape, 3> GetIndexShapes(const int nx1, const int nx2, const int nx3,
+                                         bool multilevel, const Mesh *pmesh);
 
 //----------------------------------------------------------------------------------------
 //! \class MeshBlock

From e0354f6675d53c6420f2ad4b11f76b6420d296bc Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 13 Aug 2024 14:31:46 -0600
Subject: [PATCH 26/40] fix partitioning tests

---
 src/mesh/mesh.hpp                      | 3 ++-
 src/utils/partition_stl_containers.hpp | 2 --
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp
index ee7c06140b26..ca221d3f6e56 100644
--- a/src/mesh/mesh.hpp
+++ b/src/mesh/mesh.hpp
@@ -143,7 +143,8 @@ class Mesh {
   void LoadBalancingAndAdaptiveMeshRefinement(ParameterInput *pin,
                                               ApplicationInput *app_in);
   int DefaultPackSize() {
-    return default_pack_size_ < 1 ? block_list.size() : default_pack_size_;
+    return default_pack_size_ < 1 ? std::max(static_cast<int>(block_list.size()), 1)
+                                  : default_pack_size_;
   }
   int DefaultNumPartitions() {
     return partition::partition_impl::IntCeil(block_list.size(), DefaultPackSize());
diff --git a/src/utils/partition_stl_containers.hpp b/src/utils/partition_stl_containers.hpp
index 8f1e1f8e90cc..8c943c2368e2 100644
--- a/src/utils/partition_stl_containers.hpp
+++ b/src/utils/partition_stl_containers.hpp
@@ -47,8 +47,6 @@ auto ToSizeN(Container_t<T, extra...> &container, const int N) {
   using std::to_string;
   using namespace partition_impl;
 
-  if (N == 0) return Partition_t<T>();
-
   PARTHENON_REQUIRE_THROWS(N > 0, "Your partition must be at least size 1");
 
   int nelements = container.size();

From 9e2c7eafd85b8efb1ba40089a5b10da2c4f15e12 Mon Sep 17 00:00:00 2001
From: Ben Prather <bprather@lanl.gov>
Date: Wed, 14 Aug 2024 08:54:07 -0600
Subject: [PATCH 27/40] Add `GetPackDimension` to `StateDescriptor` (#1148)

* Add GetPackDimension to StateDescriptor

Sometimes when initializing a package, you want to know how many
variables another package has declared (e.g., a `Flux` package
might need the total number of `WithFluxes` variables declared
so far).
The `Mesh` doesn't exist yet, so you have to ask the other packages
directly what they've declared.  Thus, `GetPackDimension`.

* Changelog

* Fix GitPackDim for tensors and add unit test

* GetPackDim: Shorter general implementation from Jonah

---------

Co-authored-by: Jonah Miller <jonah.maxwell.miller@gmail.com>
---
 CHANGELOG.md                       |  1 +
 src/interface/state_descriptor.cpp | 32 ++++++++++++++++++++++++++++++
 src/interface/state_descriptor.hpp | 10 ++++++++++
 tst/unit/test_state_descriptor.cpp | 16 +++++++++++++++
 4 files changed, 59 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c404edb308a7..2cfc0df66b46 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1148]](https://github.com/parthenon-hpc-lab/parthenon/pull/1148) Add `GetPackDimension` to `StateDescriptor` for calculating pack sizes before `Mesh` initialization
 - [[PR 1143]](https://github.com/parthenon-hpc-lab/parthenon/pull/1143) Add tensor indices to VariableState, add radiation constant to constants, add TypeLists, allow for arbitrary containers for solvers
 - [[PR 1140]](https://github.com/parthenon-hpc-lab/parthenon/pull/1140) Allow for relative convergence tolerance in BiCGSTAB solver.
 - [[PR 1047]](https://github.com/parthenon-hpc-lab/parthenon/pull/1047) General three- and four-valent 2D forests w/ arbitrary orientations.
diff --git a/src/interface/state_descriptor.cpp b/src/interface/state_descriptor.cpp
index 308ba14a7b67..aa3fa5e6af27 100644
--- a/src/interface/state_descriptor.cpp
+++ b/src/interface/state_descriptor.cpp
@@ -13,6 +13,7 @@
 
 #include <iomanip>
 #include <iostream>
+#include <numeric>
 #include <sstream>
 #include <string>
 #include <unordered_map>
@@ -566,4 +567,35 @@ StateDescriptor::GetVariableNames(const Metadata::FlagCollection &flags) {
   return GetVariableNames({}, flags, {});
 }
 
+// Get the total length of this StateDescriptor's variables when packed
+int StateDescriptor::GetPackDimension(const std::vector<std::string> &req_names,
+                                      const Metadata::FlagCollection &flags,
+                                      const std::vector<int> &sparse_ids) {
+  std::vector<std::string> names = GetVariableNames(req_names, flags, sparse_ids);
+  int dimension = 0;
+  for (auto name : names) {
+    const auto &meta = metadataMap_[VarID(name)];
+    // if meta.Shape().size() < 1, then 'accumulate' will return the initialization value,
+    // which is 1. Otherwise, this multiplies all elements present in 'Shape' to obtain
+    // total length
+    dimension += std::accumulate(meta.Shape().begin(), meta.Shape().end(), 1,
+                                 [](auto a, auto b) { return a * b; });
+  }
+  return dimension;
+}
+int StateDescriptor::GetPackDimension(const std::vector<std::string> &req_names,
+                                      const std::vector<int> &sparse_ids) {
+  return GetPackDimension(req_names, Metadata::FlagCollection(), sparse_ids);
+}
+int StateDescriptor::GetPackDimension(const Metadata::FlagCollection &flags,
+                                      const std::vector<int> &sparse_ids) {
+  return GetPackDimension({}, flags, sparse_ids);
+}
+int StateDescriptor::GetPackDimension(const std::vector<std::string> &req_names) {
+  return GetPackDimension(req_names, Metadata::FlagCollection(), {});
+}
+int StateDescriptor::GetPackDimension(const Metadata::FlagCollection &flags) {
+  return GetPackDimension({}, flags, {});
+}
+
 } // namespace parthenon
diff --git a/src/interface/state_descriptor.hpp b/src/interface/state_descriptor.hpp
index 47f1a9138d09..c488d53a1bcd 100644
--- a/src/interface/state_descriptor.hpp
+++ b/src/interface/state_descriptor.hpp
@@ -276,6 +276,16 @@ class StateDescriptor {
   std::vector<std::string> GetVariableNames(const std::vector<std::string> &req_names);
   std::vector<std::string> GetVariableNames(const Metadata::FlagCollection &flags);
 
+  int GetPackDimension(const std::vector<std::string> &req_names,
+                       const Metadata::FlagCollection &flags,
+                       const std::vector<int> &sparse_ids);
+  int GetPackDimension(const std::vector<std::string> &req_names,
+                       const std::vector<int> &sparse_ids);
+  int GetPackDimension(const Metadata::FlagCollection &flags,
+                       const std::vector<int> &sparse_ids);
+  int GetPackDimension(const std::vector<std::string> &req_names);
+  int GetPackDimension(const Metadata::FlagCollection &flags);
+
   std::size_t
   RefinementFuncID(const refinement::RefinementFunctions_t &funcs) const noexcept {
     return refinementFuncMaps_.funcs_to_ids.at(funcs);
diff --git a/tst/unit/test_state_descriptor.cpp b/tst/unit/test_state_descriptor.cpp
index 706f4c73c74d..95890f9caf11 100644
--- a/tst/unit/test_state_descriptor.cpp
+++ b/tst/unit/test_state_descriptor.cpp
@@ -121,6 +121,22 @@ TEST_CASE("Test Associate in StateDescriptor", "[StateDescriptor]") {
   }
 }
 
+TEST_CASE("Test GetPackDimension in StateDescriptor", "[StateDescriptor]") {
+  GIVEN("Some flags and state descriptors") {
+    StateDescriptor state("state");
+    WHEN("We add some fields with various shapes and total size") {
+      state.AddField("foo", Metadata(std::vector<MetadataFlag>{}, std::vector<int>{4}));
+      state.AddField("bar",
+                     Metadata(std::vector<MetadataFlag>{}, std::vector<int>{4, 4}));
+      state.AddField("baz",
+                     Metadata(std::vector<MetadataFlag>{}, std::vector<int>{4, 4, 4}));
+      THEN("The total length is identified correctly") {
+        REQUIRE(state.GetPackDimension(Metadata::GetUserFlag("state")) == 84);
+      }
+    }
+  }
+}
+
 TEST_CASE("Test dependency resolution in StateDescriptor", "[StateDescriptor]") {
   GIVEN("Some empty state descriptors and metadata") {
     // metadata

From 6f3023967e9c18cc16fa44eb4b1775e483b41739 Mon Sep 17 00:00:00 2001
From: Ben Prather <bprather@lanl.gov>
Date: Wed, 21 Aug 2024 16:04:44 -0600
Subject: [PATCH 28/40] Update driver.hpp

---
 src/driver/driver.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/driver/driver.hpp b/src/driver/driver.hpp
index 19ab2821ec16..4991d954133c 100644
--- a/src/driver/driver.hpp
+++ b/src/driver/driver.hpp
@@ -79,8 +79,8 @@ class EvolutionDriver : public Driver {
     pouts = std::make_unique<Outputs>(pmesh, pinput, &tm);
   }
   DriverStatus Execute() override;
-  void SetGlobalTimeStep();
-  void OutputCycleDiagnostics();
+  virtual void SetGlobalTimeStep();
+  virtual void OutputCycleDiagnostics();
   void DumpInputParameters();
 
   virtual TaskListStatus Step() = 0;

From ae73a4a51d193a7d29f723ed77ce7359cad289cc Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Thu, 22 Aug 2024 16:27:17 +0200
Subject: [PATCH 29/40] Cleanup PR and include changelog

---
 CHANGELOG.md                               |  2 ++
 src/bvals/comms/build_boundary_buffers.cpp | 13 +++++--------
 2 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2cfc0df66b46..08ffd596dbdf 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -36,6 +36,8 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1150]](https://github.com/parthenon-hpc-lab/parthenon/pull/1150) Reduce memory consumption for buffer pool
+- [[PR 1146]](https://github.com/parthenon-hpc-lab/parthenon/pull/1146) Fix an issue outputting >4GB single variables per rank
 - [[PR 1152]](https://github.com/parthenon-hpc-lab/parthenon/pull/1152) Fix memory leak in task graph outputs related to `abi::__cxa_demangle`
 - [[PR 1146]](https://github.com/parthenon-hpc-lab/parthenon/pull/1146) Fix an issue outputting >4GB single variables per rank
 - [[PR 1144]](https://github.com/parthenon-hpc-lab/parthenon/pull/1144) Fix some restarts w/non-CC fields
diff --git a/src/bvals/comms/build_boundary_buffers.cpp b/src/bvals/comms/build_boundary_buffers.cpp
index ca49d3210d0a..aac532d037e6 100644
--- a/src/bvals/comms/build_boundary_buffers.cpp
+++ b/src/bvals/comms/build_boundary_buffers.cpp
@@ -52,6 +52,10 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
     int buf_size = GetBufferSize(pmb, nb, v);
+    //  LR: Multigrid logic requires blocks sending messages to themselves (since the same
+    //  block can show up on two multigrid levels). This doesn't require any data
+    //  transfer, so the message size can be zero. It is essentially just a flag to show
+    //  that the block is done being used on one level and can be used on the next level.
     if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
 
     nbufs[buf_size] += 1; // relying on value init of int to 0 for initial entry
@@ -60,12 +64,12 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
   ForEachBoundary<BTYPE>(md, [&](auto pmb, sp_mbd_t /*rc*/, nb_t &nb, const sp_cv_t v) {
     // Calculate the required size of the buffer for this boundary
     int buf_size = GetBufferSize(pmb, nb, v);
+    // See comment above on the same logic.
     if (pmb->gid == nb.gid && nb.offsets.IsCell()) buf_size = 0;
 
     // Add a buffer pool if one does not exist for this size
     using buf_t = buf_pool_t<Real>::base_t;
     if (pmesh->pool_map.count(buf_size) == 0) {
-      std::cerr << "Setting up a new pool for buffers of size " << buf_size << "\n";
       // Might be worth discussing what a good default is.
       // Using the number of packs, assumes that all blocks in a pack have fairly similar
       // buffer configurations, which may or may not be a good approximation.
@@ -74,10 +78,6 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
       const int64_t nbuf = pmesh->DefaultNumPartitions();
       pmesh->pool_map.emplace(
           buf_size, buf_pool_t<Real>([buf_size, nbuf](buf_pool_t<Real> *pool) {
-            std::cerr << "Dynamically adding " << nbuf << " buffers of size " << buf_size
-                      << " to a pool with current size " << pool->NumBuffersInPool()
-                      << " and future size " << pool->NumBuffersInPool() + nbuf << "\n";
-
             const auto pool_size = nbuf * buf_size;
             buf_t chunk("pool buffer", pool_size);
             for (int i = 1; i < nbuf; ++i) {
@@ -92,9 +92,6 @@ void BuildBoundaryBufferSubset(std::shared_ptr<MeshData<Real>> &md,
     auto &pool = pmesh->pool_map.at(buf_size);
     const std::int64_t new_buffers_req = nbufs.at(buf_size) - pool.NumBuffersInPool();
     if (new_buffers_req > 0) {
-      std::cerr << "Reserving " << new_buffers_req << " new buffers of size " << buf_size
-                << " to pool with " << pool.NumBuffersInPool() << " buffers because "
-                << nbufs.at(buf_size) << " are required in total.\n";
       const auto pool_size = new_buffers_req * buf_size;
       buf_t chunk("pool buffer", pool_size);
       for (int i = 0; i < new_buffers_req; ++i) {

From 6968aa625146f4703a98fcee7fd0a4c7b1a4b569 Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Fri, 23 Aug 2024 10:11:18 -0600
Subject: [PATCH 30/40] Add additional timestep controllers and safeties

---
 src/driver/driver.cpp | 57 ++++++++++++++++++++++++++++++++++---------
 src/driver/driver.hpp | 18 ++++++++++++++
 2 files changed, 64 insertions(+), 11 deletions(-)

diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index 332468e2b83b..bf840876981e 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -205,15 +205,27 @@ void EvolutionDriver::InitializeBlockTimeSteps() {
 // \brief function that loops over all MeshBlocks and find new timestep
 
 void EvolutionDriver::SetGlobalTimeStep() {
-  // don't allow dt to grow by more than 2x
-  // consider making this configurable in the input
-  if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
-    tm.dt *= 2.0;
-  }
-  Real big = std::numeric_limits<Real>::max();
-  for (auto const &pmb : pmesh->block_list) {
-    tm.dt = std::min(tm.dt, pmb->NewDt());
-    pmb->SetAllowedDt(big);
+  // Check if user wants to force the value
+  if (dt_force > 0.0) {
+    tm.dt = dt_force;
+  } else {
+    // Check for special first cycle value
+    if (tm.ncycle == 0) {
+      tm.dt = std::min(tm.dt, dt_init);
+    }
+    // don't allow dt to grow by more than 2x
+    // consider making this configurable in the input
+    if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
+      tm.dt *= 2.0;
+    }
+    // Allow the meshblocks to vote
+    Real big = std::numeric_limits<Real>::max();
+    for (auto const &pmb : pmesh->block_list) {
+      tm.dt = std::min(tm.dt, pmb->NewDt());
+      pmb->SetAllowedDt(big);
+    }
+    // Allow the user to vote
+    tm.dt = std::min(tm.dt, dt_user);
   }
 
 #ifdef MPI_PARALLEL
@@ -221,9 +233,32 @@ void EvolutionDriver::SetGlobalTimeStep() {
                                     MPI_COMM_WORLD));
 #endif
 
-  if (tm.time < tm.tlim &&
-      (tm.tlim - tm.time) < tm.dt) // timestep would take us past desired endpoint
+  // Check that we have not gone off the rails
+  if (tm.dt <= dt_min) {
+    if (++dt_min_count >= dt_min_count_max) {
+      std::stringstream msg;
+      msg << "Timesetep has fallen bellow minimum (parthenon/time/dt_min=" << dt_min
+          << ") for more than " << dt_min_count_max << " steps";
+      PARTHENON_FAIL(msg);
+    } else {
+      dt_min_count = 0;
+    }
+  }
+  if (tm.dt >= dt_max) {
+    if (++dt_max_count >= dt_max_count_max) {
+      std::stringstream msg;
+      msg << "Timesetep has risen above maximum (parthenon/time/dt_max=" << dt_max
+          << ") for more than " << dt_max_count_max << " steps";
+      PARTHENON_FAIL(msg);
+    } else {
+      dt_max_count = 0;
+    }
+  }
+
+  // Limit timestep if it would take us past desired endpoint
+  if (tm.time < tm.tlim && (tm.tlim - tm.time) < tm.dt) {
     tm.dt = tm.tlim - tm.time;
+  }
 }
 
 void EvolutionDriver::DumpInputParameters() {
diff --git a/src/driver/driver.hpp b/src/driver/driver.hpp
index 4991d954133c..b05a798d0b05 100644
--- a/src/driver/driver.hpp
+++ b/src/driver/driver.hpp
@@ -69,6 +69,21 @@ class EvolutionDriver : public Driver {
                                       std::numeric_limits<Real>::infinity());
     Real dt =
         pinput->GetOrAddReal("parthenon/time", "dt", std::numeric_limits<Real>::max());
+    dt_min = pinput->GetOrAddReal("parthenon/time", "dt_min",
+                                  std::numeric_limits<Real>::min());
+    dt_max = pinput->GetOrAddReal("parthenon/time", "dt_max",
+                                  std::numeric_limits<Real>::max());
+    dt_init = pinput->GetOrAddReal("parthenon/time", "dt_init",
+                                   std::numeric_limits<Real>::max());
+    dt_user = pinput->GetOrAddReal("parthenon/time", "dt_user",
+                                   std::numeric_limits<Real>::max());
+    dt_force = pinput->GetOrAddReal("parthenon/time", "dt_force",
+                                    std::numeric_limits<Real>::lowest());
+    dt_min_count_max =
+        pinput->GetOrAddInteger("parthenon/time", "dt_min_cycle_limit", 10);
+    dt_max_count_max = pinput->GetOrAddInteger("parthenon/time", "dt_max_cycle_limit", 1);
+    dt_min_count = 0;
+    dt_max_count = 0;
     const auto ncycle = pinput->GetOrAddInteger("parthenon/time", "ncycle", 0);
     const auto nmax = pinput->GetOrAddInteger("parthenon/time", "nlim", -1);
     const auto nout = pinput->GetOrAddInteger("parthenon/time", "ncycle_out", 1);
@@ -88,6 +103,9 @@ class EvolutionDriver : public Driver {
 
  protected:
   void PostExecute(DriverStatus status) override;
+  Real dt_user, dt_force, dt_init, dt_min, dt_max;
+  int dt_min_count, dt_max_count;
+  int dt_min_count_max, dt_max_count_max;
 
  private:
   void InitializeBlockTimeSteps();

From 0998f8bf8140882f8dfcdba322b669256630731f Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Fri, 23 Aug 2024 12:21:09 -0600
Subject: [PATCH 31/40] Fix count increments

---
 src/driver/driver.cpp | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index bf840876981e..34d07a5dee28 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -240,9 +240,9 @@ void EvolutionDriver::SetGlobalTimeStep() {
       msg << "Timesetep has fallen bellow minimum (parthenon/time/dt_min=" << dt_min
           << ") for more than " << dt_min_count_max << " steps";
       PARTHENON_FAIL(msg);
-    } else {
-      dt_min_count = 0;
     }
+  } else {
+    dt_min_count = 0;
   }
   if (tm.dt >= dt_max) {
     if (++dt_max_count >= dt_max_count_max) {
@@ -250,9 +250,9 @@ void EvolutionDriver::SetGlobalTimeStep() {
       msg << "Timesetep has risen above maximum (parthenon/time/dt_max=" << dt_max
           << ") for more than " << dt_max_count_max << " steps";
       PARTHENON_FAIL(msg);
-    } else {
-      dt_max_count = 0;
     }
+  } else {
+    dt_max_count = 0;
   }
 
   // Limit timestep if it would take us past desired endpoint

From abc0c6c63232febc584f43a25c36f235a918cb09 Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Fri, 23 Aug 2024 13:13:15 -0600
Subject: [PATCH 32/40] Add more dt controllers. Cleanup a bit

---
 src/driver/driver.cpp | 29 +++++++++++++++++------------
 src/driver/driver.hpp | 13 ++++++++++++-
 2 files changed, 29 insertions(+), 13 deletions(-)

diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index 34d07a5dee28..54e57485c893 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -205,33 +205,36 @@ void EvolutionDriver::InitializeBlockTimeSteps() {
 // \brief function that loops over all MeshBlocks and find new timestep
 
 void EvolutionDriver::SetGlobalTimeStep() {
-  // Check if user wants to force the value
   if (dt_force > 0.0) {
+    // Check if user wants to force the value
     tm.dt = dt_force;
+  } else if ((tm.ncycle == 0) && (dt_init_force) && (dt_init > 0.0)) {
+    // Check if user wants to force the first cycle value
+    tm.dt = dt_init;
   } else {
+    if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
+      // don't allow dt to grow by more than 2x
+      // consider making this configurable in the input
+      tm.dt *= dt_factor;
+    }
     // Check for special first cycle value
     if (tm.ncycle == 0) {
       tm.dt = std::min(tm.dt, dt_init);
     }
-    // don't allow dt to grow by more than 2x
-    // consider making this configurable in the input
-    if (tm.dt < 0.1 * std::numeric_limits<Real>::max()) {
-      tm.dt *= 2.0;
-    }
     // Allow the meshblocks to vote
-    Real big = std::numeric_limits<Real>::max();
     for (auto const &pmb : pmesh->block_list) {
       tm.dt = std::min(tm.dt, pmb->NewDt());
-      pmb->SetAllowedDt(big);
+      pmb->SetAllowedDt(std::numeric_limits<Real>::max());
     }
     // Allow the user to vote
     tm.dt = std::min(tm.dt, dt_user);
-  }
-
+    // Force timestep to be in the allowable range
+    tm.dt = std::max(dt_floor, std::min(tm.dt, dt_ceil));
 #ifdef MPI_PARALLEL
-  PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &tm.dt, 1, MPI_PARTHENON_REAL, MPI_MIN,
-                                    MPI_COMM_WORLD));
+    PARTHENON_MPI_CHECK(MPI_Allreduce(MPI_IN_PLACE, &tm.dt, 1, MPI_PARTHENON_REAL,
+                                      MPI_MIN, MPI_COMM_WORLD));
 #endif
+  }
 
   // Check that we have not gone off the rails
   if (tm.dt <= dt_min) {
@@ -256,6 +259,8 @@ void EvolutionDriver::SetGlobalTimeStep() {
   }
 
   // Limit timestep if it would take us past desired endpoint
+  // This comes after the bounds checking so that we don't fail when epsilon
+  // away from tlim
   if (tm.time < tm.tlim && (tm.tlim - tm.time) < tm.dt) {
     tm.dt = tm.tlim - tm.time;
   }
diff --git a/src/driver/driver.hpp b/src/driver/driver.hpp
index b05a798d0b05..691dd17563f3 100644
--- a/src/driver/driver.hpp
+++ b/src/driver/driver.hpp
@@ -75,15 +75,24 @@ class EvolutionDriver : public Driver {
                                   std::numeric_limits<Real>::max());
     dt_init = pinput->GetOrAddReal("parthenon/time", "dt_init",
                                    std::numeric_limits<Real>::max());
+    dt_init_force = pinput->GetOrAddBoolean("parthenon/time", "dt_init_force", false);
+
     dt_user = pinput->GetOrAddReal("parthenon/time", "dt_user",
                                    std::numeric_limits<Real>::max());
     dt_force = pinput->GetOrAddReal("parthenon/time", "dt_force",
                                     std::numeric_limits<Real>::lowest());
+    dt_floor = pinput->GetOrAddReal("parthenon/time", "dt_floor",
+                                    std::numeric_limits<Real>::min());
+    dt_ceil = pinput->GetOrAddReal("parthenon/time", "dt_ceil",
+                                   std::numeric_limits<Real>::max());
     dt_min_count_max =
         pinput->GetOrAddInteger("parthenon/time", "dt_min_cycle_limit", 10);
     dt_max_count_max = pinput->GetOrAddInteger("parthenon/time", "dt_max_cycle_limit", 1);
     dt_min_count = 0;
     dt_max_count = 0;
+
+    dt_factor = pinput->GetOrAddReal("parthenon/time", "dt_factor", 2.0);
+
     const auto ncycle = pinput->GetOrAddInteger("parthenon/time", "ncycle", 0);
     const auto nmax = pinput->GetOrAddInteger("parthenon/time", "nlim", -1);
     const auto nout = pinput->GetOrAddInteger("parthenon/time", "ncycle_out", 1);
@@ -103,7 +112,9 @@ class EvolutionDriver : public Driver {
 
  protected:
   void PostExecute(DriverStatus status) override;
-  Real dt_user, dt_force, dt_init, dt_min, dt_max;
+  Real dt_user, dt_force, dt_init, dt_min, dt_max, dt_floor, dt_ceil;
+  Real dt_factor;
+  bool dt_init_force;
   int dt_min_count, dt_max_count;
   int dt_min_count_max, dt_max_count_max;
 

From 21bed3ab4cd0010f94f77e523c79de9deea528fa Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Fri, 23 Aug 2024 13:14:55 -0600
Subject: [PATCH 33/40] Not a vote

---
 src/driver/driver.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/driver/driver.cpp b/src/driver/driver.cpp
index 54e57485c893..85327ecbae3b 100644
--- a/src/driver/driver.cpp
+++ b/src/driver/driver.cpp
@@ -226,7 +226,7 @@ void EvolutionDriver::SetGlobalTimeStep() {
       tm.dt = std::min(tm.dt, pmb->NewDt());
       pmb->SetAllowedDt(std::numeric_limits<Real>::max());
     }
-    // Allow the user to vote
+    // Allow the user to enforce maximum timestep
     tm.dt = std::min(tm.dt, dt_user);
     // Force timestep to be in the allowable range
     tm.dt = std::max(dt_floor, std::min(tm.dt, dt_ceil));

From f0c49f2b1742a684d4052ef7c460b6821a37172b Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Fri, 23 Aug 2024 14:56:23 -0600
Subject: [PATCH 34/40] Update parthenon/time parameters

---
 doc/sphinx/src/inputs.rst | 25 ++++++++++++++++++-------
 1 file changed, 18 insertions(+), 7 deletions(-)

diff --git a/doc/sphinx/src/inputs.rst b/doc/sphinx/src/inputs.rst
index c68c89c7b95a..6d83ffceed7e 100644
--- a/doc/sphinx/src/inputs.rst
+++ b/doc/sphinx/src/inputs.rst
@@ -21,8 +21,8 @@ General parthenon options such as problem name and parameter handling.
 +---------------------+---------+---------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | Option              | Default | Type    | Description                                                                                                                                                                                            |
 +=====================+=========+=========+========================================================================================================================================================================================================+
-|| name               || none   || string || Name of this problem or initialization, prefixed to output files.                                                                                                                                     |
 || archive_parameters || false  || string || Produce a parameter file containing all parameters known to Parthenon. Set to `true` for an output file named `parthinput.archive`. Set to `timestamp` for a file with a name containing a timestamp. |
+|| name               || none   || string || Name of this problem or initialization, prefixed to output files.                                                                                                                                     |
 +---------------------+---------+---------+--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 
 
@@ -34,12 +34,23 @@ Options related to time-stepping and printing of diagnostic data.
 +------------------------------+---------+--------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 | Option                       | Default | Type   | Description                                                                                                                                                           |
 +==============================+=========+========+=======================================================================================================================================================================+
-|| tlim                        || none   || float || Stop criterion on simulation time.                                                                                                                                   |
-|| nlim                        || -1     || int   || Stop criterion on total number of steps taken. Ignored if < 0.                                                                                                       |
-|| perf_cycle_offset           || 0      || int   || Skip the first N cycles when calculating the final performance (e.g., zone-cycles/wall_second). Allows to hide the initialization overhead in Parthenon.             |
+|| dt_ceil                     || none   || Real  || The maximum allowed timestep.                                                                                                                                        |
+|| dt_factor                   || 2.0    || Real  || The maximum allowed relative increase of the timestep over the previous value.                                                                                       |
+|| dt_floor                    || none   || Real  || The minimum allowed timestep.                                                                                                                                        |
+|| dt_force                    || none   || Real  || Force the timestep to this value, ignoring all other conditions.                                                                                                     |
+|| dt_init                     || none   || Real  || The maximum allowed timestep during the first cycle.                                                                                                                 |
+|| dt_init_force               || none   || bool  || If set to true, force the first cycle's timestep to the value given by dt_init.                                                                                      |
+|| dt_min                      || none   || Real  || If the timestep falls below dt_min for dt_min_cycle_limit cycles, Parthenon fatals.                                                                                  |
+|| dt_min_cycle_limit          || 10     || int   || The maximum number of cycles the timestep can be below dt_min.                                                                                                       |
+|| dt_max                      || none   || Real  || If the timestep falls below dt_max for dt_max_cycle_limit cycles, Parthenon fatals.                                                                                  |
+|| dt_max_cycle_limit          || 1      || int   || The maximum number of cycles the timestep an be above dt_max.                                                                                                        |
+|| dt_user                     || none   || Real  || Set a global timestep limit.                                                                                                                                         |
+|| ncrecv_bdry_buf_timeout_sec || -1.0   || Real  || Timeout in seconds for the `ReceiveBoundaryBuffers` tasks. Disabed (negative) by default. Typically no need in production runs. Useful for debugging MPI calls.      |
 || ncycle_out                  || 1      || int   || Number of cycles between short diagnostic output to standard out containing, e.g., current time, dt, zone-update/wsec. Default: 1 (i.e, every cycle).                |
 || ncycle_out_mesh             || 0      || int   || Number of cycles between printing the mesh structure to standard out. Use a negative number to also print every time the mesh was modified. Default: 0 (i.e, off).   |
-|| ncrecv_bdry_buf_timeout_sec || -1.0   || Real  || Timeout in seconds for the `ReceiveBoundaryBuffers` tasks. Disabed (negative) by default. Typically no need in production runs. Useful for debugging MPI calls.      |
+|| nlim                        || -1     || int   || Stop criterion on total number of steps taken. Ignored if < 0.                                                                                                       |
+|| perf_cycle_offset           || 0      || int   || Skip the first N cycles when calculating the final performance (e.g., zone-cycles/wall_second). Allows to hide the initialization overhead in Parthenon.             |
+|| tlim                        || none   || Real  || Stop criterion on simulation time.                                                                                                                                   |
 +------------------------------+---------+--------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------+
 
 
@@ -64,9 +75,9 @@ See the :ref:`sparse impl` documentation for details.
 +--------------------+---------+--------+----------------------------------------------------------------------------------------------------------------------------------------------+
 | Option             | Default | Type   | Description                                                                                                                                  |
 +====================+=========+========+==============================================================================================================================================+
-|| enable_sparse     || `true` || bool  || If set to false, sparse variables will always be allocated, see also :ref:`sparse run-time`                                                 |
 || alloc_threshold   || 1e-12  || float || Global (for all sparse variables) threshold to trigger allocation of a variable if cells in the receiving ghost cells are above this value. |
-|| dealloc_threshold || 1e-14  || float || Global (for all sparse variables) threshold to trigger deallocation if all active cells of a variable in a block are below this value.      |
 || dealloc_count     || 5      || int   || First deallocate a sparse variable if the `dealloc_threshold` has been met in this number of consecutive cycles.                            |
+|| dealloc_threshold || 1e-14  || float || Global (for all sparse variables) threshold to trigger deallocation if all active cells of a variable in a block are below this value.      |
+|| enable_sparse     || `true` || bool  || If set to false, sparse variables will always be allocated, see also :ref:`sparse run-time`                                                 |
 +--------------------+---------+--------+----------------------------------------------------------------------------------------------------------------------------------------------+
 

From 0fc6a9852d8d0a6c59213da882abfa3057d107cd Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Fri, 23 Aug 2024 14:58:00 -0600
Subject: [PATCH 35/40] changelog for timestep controllers

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a86e3b20889c..7b3fa70d917e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1159]](https://github.com/parthenon-hpc-lab/parthenon/pull/1159) Add additional timestep controllers in parthenon/time.
 - [[PR 1148]](https://github.com/parthenon-hpc-lab/parthenon/pull/1148) Add `GetPackDimension` to `StateDescriptor` for calculating pack sizes before `Mesh` initialization
 - [[PR 1143]](https://github.com/parthenon-hpc-lab/parthenon/pull/1143) Add tensor indices to VariableState, add radiation constant to constants, add TypeLists, allow for arbitrary containers for solvers
 - [[PR 1140]](https://github.com/parthenon-hpc-lab/parthenon/pull/1140) Allow for relative convergence tolerance in BiCGSTAB solver.

From 7529be3f4790424b13d94ac55e349527396044b1 Mon Sep 17 00:00:00 2001
From: Ben Prather <bprather@lanl.gov>
Date: Mon, 26 Aug 2024 04:23:36 -0600
Subject: [PATCH 36/40] Update driver.rst (#1158)

---
 doc/sphinx/src/driver.rst | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/doc/sphinx/src/driver.rst b/doc/sphinx/src/driver.rst
index a7c400bc9086..2587a806c070 100644
--- a/doc/sphinx/src/driver.rst
+++ b/doc/sphinx/src/driver.rst
@@ -32,7 +32,9 @@ The ``EvolutionDriver`` class derives from ``Driver``, defining the
 
 loop, including periodic outputs. It has a single pure virtual member
 function called ``Step`` which a derived class must define and which
-will be called during each pass of the loop above.
+will be called during each pass of the loop above. The
+``SetGlobalTimeStep`` and ``OutputCycleDiagnostics`` functions have
+default implementations, but can be overridden for flexibility.
 
 MultiStageDriver
 ----------------

From 83aff4ce6dbbdce870693df57a105a84c0a61594 Mon Sep 17 00:00:00 2001
From: Ben Prather <bprather@lanl.gov>
Date: Mon, 26 Aug 2024 07:43:53 -0600
Subject: [PATCH 37/40] Add `par_reduce_inner` functions (#1147)

* Add par_reduce_inner functions

* Formatting, changelog

* Be explicit that existing inner reductions are TeamThreadRange
---
 CHANGELOG.md               |  1 +
 src/kokkos_abstraction.hpp | 59 ++++++++++++++++++++++++++++++++++++++
 2 files changed, 60 insertions(+)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7b3fa70d917e..a29012d8ffb6 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1147]](https://github.com/parthenon-hpc-lab/parthenon/pull/1147) Add `par_reduce_inner` functions
 - [[PR 1159]](https://github.com/parthenon-hpc-lab/parthenon/pull/1159) Add additional timestep controllers in parthenon/time.
 - [[PR 1148]](https://github.com/parthenon-hpc-lab/parthenon/pull/1148) Add `GetPackDimension` to `StateDescriptor` for calculating pack sizes before `Mesh` initialization
 - [[PR 1143]](https://github.com/parthenon-hpc-lab/parthenon/pull/1143) Add tensor indices to VariableState, add radiation constant to constants, add TypeLists, allow for arbitrary containers for solvers
diff --git a/src/kokkos_abstraction.hpp b/src/kokkos_abstraction.hpp
index ca8c59ffe12e..8fa89f82e95e 100644
--- a/src/kokkos_abstraction.hpp
+++ b/src/kokkos_abstraction.hpp
@@ -976,6 +976,65 @@ KOKKOS_FORCEINLINE_FUNCTION void par_for_inner(team_mbr_t team_member, Args &&..
   par_for_inner(DEFAULT_INNER_LOOP_PATTERN, team_member, std::forward<Args>(args)...);
 }
 
+// Inner reduction loops
+template <typename Function, typename T>
+KOKKOS_FORCEINLINE_FUNCTION void
+par_reduce_inner(InnerLoopPatternTTR, team_mbr_t team_member, const int kl, const int ku,
+                 const int jl, const int ju, const int il, const int iu,
+                 const Function &function, T reduction) {
+  const int Nk = ku - kl + 1;
+  const int Nj = ju - jl + 1;
+  const int Ni = iu - il + 1;
+  const int NkNjNi = Nk * Nj * Ni;
+  const int NjNi = Nj * Ni;
+  Kokkos::parallel_reduce(
+      Kokkos::TeamThreadRange(team_member, NkNjNi),
+      [&](const int &idx, typename T::value_type &lreduce) {
+        int k = idx / NjNi;
+        int j = (idx - k * NjNi) / Ni;
+        int i = idx - k * NjNi - j * Ni;
+        k += kl;
+        j += jl;
+        i += il;
+        function(k, j, i, lreduce);
+      },
+      reduction);
+}
+
+template <typename Function, typename T>
+KOKKOS_FORCEINLINE_FUNCTION void
+par_reduce_inner(InnerLoopPatternTTR, team_mbr_t team_member, const int jl, const int ju,
+                 const int il, const int iu, const Function &function, T reduction) {
+  const int Nj = ju - jl + 1;
+  const int Ni = iu - il + 1;
+  const int NjNi = Nj * Ni;
+  Kokkos::parallel_reduce(
+      Kokkos::TeamThreadRange(team_member, NjNi),
+      [&](const int &idx, typename T::value_type &lreduce) {
+        int j = idx / Ni;
+        int i = idx - j * Ni;
+        j += jl;
+        i += il;
+        function(j, i, lreduce);
+      },
+      reduction);
+}
+
+template <typename Function, typename T>
+KOKKOS_FORCEINLINE_FUNCTION void
+par_reduce_inner(InnerLoopPatternTTR, team_mbr_t team_member, const int il, const int iu,
+                 const Function &function, T reduction) {
+  const int Ni = iu - il + 1;
+  Kokkos::parallel_reduce(
+      Kokkos::TeamThreadRange(team_member, Ni),
+      [&](const int &idx, typename T::value_type &lreduce) {
+        int i = idx;
+        i += il;
+        function(i, lreduce);
+      },
+      reduction);
+}
+
 // reused from kokoks/core/perf_test/PerfTest_ExecSpacePartitioning.cpp
 // commit a0d011fb30022362c61b3bb000ae3de6906cb6a7
 template <class ExecSpace>

From e0bd43bf0bdd899b99e891610c61a4434f83d086 Mon Sep 17 00:00:00 2001
From: Peter Brady <ptb@lanl.gov>
Date: Mon, 26 Aug 2024 11:50:32 -0500
Subject: [PATCH 38/40] Add time offset for LowStorageIntegrator (#1151)

* add time offset for LowStorageIntegrator

* Add docs for `c` in `LowStorageIntegrator`

* Fix math expressions in docs

* add notebook for sparse time levels

* Remove binary content from notebook

* Oops. Fixed mixed files.
---
 CHANGELOG.md                                  |   1 +
 doc/sphinx/src/integrators.rst                |  20 +-
 scripts/mathematica/sparse_integrators.nb     | 232 ++++++++++++++++++
 .../low_storage_integrator.cpp                |  19 ++
 src/time_integration/staged_integrator.hpp    |   2 +
 5 files changed, 268 insertions(+), 6 deletions(-)
 create mode 100644 scripts/mathematica/sparse_integrators.nb

diff --git a/CHANGELOG.md b/CHANGELOG.md
index a29012d8ffb6..4f474633ce8e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1151]]((https://github.com/parthenon-hpc-lab/parthenon/pull/1151)) Add time offset `c` to LowStorageIntegrator
 - [[PR 1147]](https://github.com/parthenon-hpc-lab/parthenon/pull/1147) Add `par_reduce_inner` functions
 - [[PR 1159]](https://github.com/parthenon-hpc-lab/parthenon/pull/1159) Add additional timestep controllers in parthenon/time.
 - [[PR 1148]](https://github.com/parthenon-hpc-lab/parthenon/pull/1148) Add `GetPackDimension` to `StateDescriptor` for calculating pack sizes before `Mesh` initialization
diff --git a/doc/sphinx/src/integrators.rst b/doc/sphinx/src/integrators.rst
index 8e035832f473..53b67b7a19e0 100644
--- a/doc/sphinx/src/integrators.rst
+++ b/doc/sphinx/src/integrators.rst
@@ -32,15 +32,17 @@ described in `Ketchson (2010)`_. These integrators are of the classic
 .. math::
 
    u^{(0)} &= u^n \\
-   u^{(i)} &= \sum_{k=0}^{i-1} (\alpha_{i,k} u^{(k)} + \Delta t \beta_{i, k} F(u^{(k)})\\
+   u^{(i)} &= \sum_{k=0}^{i-1} (\alpha_{i,k} u^{(k)} + \Delta t \beta_{i, k} F(t^n+c_k \Delta t, u^{(k)}))\\
    u^{n+1} &= u^{(m)}
 
 where superscripts in parentheses mean subcycles in a Runge-Kutta
-integration and :math:`F` is the right-hand-side of ODE system. The
+integration and :math:`F` is the right-hand-side of ODE system. Note
+that the time dependence of :math:`F` is explicitly included in the above
+formulation. The
 difference between these low-storage methods and the classic SSPK
 methods is that the low-storage methods typically have sparse
 :math:`\alpha` and :math:`\beta` matrices, which are replaced by
-diagonal termes, named :math:`\gamma_0` and :math:`\gamma_1`
+diagonal terms, named :math:`\gamma_0` and :math:`\gamma_1`
 respectively. 
 
 These methods can be generalized to support more general methods with
@@ -52,11 +54,17 @@ The full update then takes the form:
 .. math::
 
    u^{(1)} &:= u^{(1)} + \delta_s u^{(0)} \\
-   u^{(0)} &:= \gamma_{s0} u^{(0)} + \gamma_{s1} u^{(1)} + \beta_{s,s-1} \Delta t F(u^{(0)})
+   u^{(0)} &:= \gamma_{s0} u^{(0)} + \gamma_{s1} u^{(1)} + \beta_{s,s-1} \Delta t F(t^n+c_s\Delta t, u^{(0)})
 
 where here :math:`u^{(0)}` and :math:`u^{(1)}` are the two storage
 buffers required to compute the update for a given Runge-Kutta stage
-:math:`s`.
+:math:`s`.  While the :math:`\delta`, :math:`\beta`, :math:`\gamma_0` and :math:`\gamma_1`
+associated with a particular scheme are published in the literature, :math:`c` is not.
+Instead, :math:`c` is computed following the procedure outlined in 
+`Ketchson (2010)`_ for obtaining the Butcher coefficients from their low-storage
+counterparts.
+A Mathematica notebook to calculate :math:`c` is provided
+`here <https://github.com/parthenon-hpc-lab/parthenon/blob/develop/scripts/mathematica/sparse_integrators.nb>`__.
 
 .. _Ketchson (2010): https://doi.org/10.1016/j.jcp.2009.11.006
 
@@ -65,7 +73,7 @@ buffers required to compute the update for a given Runge-Kutta stage
 .. _Athena++ paper: https://doi.org/10.3847/1538-4365/ab929b
 
 The ``LowStorageIntegrator`` contains arrays for ``delta``, ``beta``,
-``gam0``, and ``gam1``. Available integration methods are:
+``gam0``, ``gam1``, and ``c``. Available integration methods are:
 
 * ``RK1``, which is simply forward Euler.
 
diff --git a/scripts/mathematica/sparse_integrators.nb b/scripts/mathematica/sparse_integrators.nb
new file mode 100644
index 000000000000..9d7acd9f8bae
--- /dev/null
+++ b/scripts/mathematica/sparse_integrators.nb
@@ -0,0 +1,232 @@
+Notebook[{Cell[
+CellGroupData[{Cell[
+"\<Computing Butcher Coefficients\>","Chapter",
+CellChangeTimes -> {{3933340025.507483`,3933340039.6422744`}},ExpressionUUID -> "da94eb33-72d4-4c5a-bfbb-d9684d02fbee"],Cell[
+TextData[
+{"The goal is to compute the time offset, \[ScriptC], associated with the stages of a low storage RK integrator.  This is done by computing the Butcher table coefficients from the low storage coefficients following the algorithm in   ",ButtonBox[
+"https://doi.org/10.1016/j.jcp.2009.11.006",BaseStyle -> "Hyperlink",ButtonData -> {URL[
+"https://doi.org/10.1016/j.jcp.2009.11.006"],None},ButtonNote -> "https://doi.org/10.1016/j.jcp.2009.11.006"]}],
+"Text",CellChangeTimes -> {{3933340050.500537`,3933340155.1628`},{3933340346.592301`,3933340396.448125`},{3933340431.4327593`,3933340466.733057`},{3933340534.786607`,3933340550.3830557`}},
+ExpressionUUID -> "322d54b1-bb59-4650-9324-a3f721acdac0"]},Open],ExpressionUUID -> "29cf4e4d-2ee9-462d-827b-1aff875fd6c3"],Cell[
+CellGroupData[
+{Cell[
+"\<2S Method - \[Alpha], \[Beta] matrices\>","Chapter",CellChangeTimes -> {{3933340715.932539`,3933340736.7708607`}},
+ExpressionUUID -> "1a30b021-4730-4bd2-a854-e1f56039de6b"],Cell[
+"\<The sparse \[Alpha], and \[Beta] matrices are constructed according to Eqs 7 and 8 in the paper as well as the 3 recurrence relations specific to the 2S methods\>",
+"Text",CellChangeTimes -> {{3933340749.2902527`,3933340759.7757583`},{3933340811.433118`,3933340855.975274`},{3933340979.951785`,3933341019.3198247`}},
+ExpressionUUID -> "6b67523b-41f1-41dd-a447-8333be1a45e4"],Cell[
+BoxData[{RowBox[{RowBox[
+{"\[Beta]ij","[",RowBox[{RowBox[{"i_","/;",RowBox[{"i","<=","1"}]}],",","j_"}],"]"}],":=","0"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Alpha]ij","[",RowBox[{RowBox[{"i_","/;",RowBox[{"i","<=","1"}]}],",","j_"}],"]"}],":=","0"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Delta]","[","1","]"}],":=","1"}],"\[IndentingNewLine]",RowBox[{RowBox[
+{"\[Gamma]","[",RowBox[{RowBox[{"i_","/;",RowBox[{"i",">=","2"}]}],",","1"}],"]"}],":="," ",RowBox[
+{"1","-",RowBox[{RowBox[{"\[Gamma]","[",RowBox[{"i",",","2"}],"]"}],RowBox[{"Sum","[",RowBox[
+{RowBox[{"\[Delta]","[","j","]"}],",",RowBox[{"{",RowBox[{"j",",",RowBox[{"i","-","1"}]}],"}"}]}],"]"}]}]}]}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"2",",","2"}],"]"}],":=","1"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{RowBox[{"i_","/;",RowBox[{"i",">","2"}]}],","," ","j_"}],"]"}],":=",RowBox[
+{RowBox[{RowBox[{"-",RowBox[{"\[Gamma]","[",RowBox[{"i",",","2"}],"]"}]}],RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{RowBox[{"i","-","1"}],",","j"}],"]"}],"/",RowBox[
+{"\[Gamma]","[",RowBox[{RowBox[{"i","-","1"}],",","2"}],"]"}]}]}],"/;",RowBox[{"j","==",RowBox[
+{"i","-","2"}]}]}]}],"\[IndentingNewLine]",RowBox[{RowBox[{"\[Alpha]ij","[",RowBox[
+{RowBox[{"i_","/;",RowBox[{"i",">","2"}]}],",","j_"}],"]"}],":=",RowBox[{RowBox[{RowBox[
+{"-",RowBox[{"\[Gamma]","[",RowBox[{"i",",","2"}],"]"}]}],RowBox[{RowBox[{"\[Gamma]","[",RowBox[
+{RowBox[{"i","-","1"}],",","1"}],"]"}],"/",RowBox[{"\[Gamma]","[",RowBox[{RowBox[
+{"i","-","1"}],",","2"}],"]"}]}]}],"/;",RowBox[{"j","==",RowBox[{"i","-","2"}]}]}]}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Alpha]ij","[",RowBox[{RowBox[{"i_","/;",RowBox[{"i",">","2"}]}],",","j_"}],"]"}],":=",RowBox[
+{RowBox[{"1","-",RowBox[{"\[Alpha]ij","[",RowBox[{"i",",",RowBox[{"i","-","2"}]}],"]"}]}],"/;",RowBox[
+{"j","==",RowBox[{"i","-","1"}]}]}]}]}],"Input",CellChangeTimes -> {{3933341754.0773373`,3933341778.301538`},{3933342058.0719028`,3933342063.8268886`},{3933342158.92589`,3933342180.880818`},3933342222.003084`,{3933342289.2385783`,3933342462.3291593`},{3933342529.645609`,3933342536.171939`},{3933342681.3267727`,3933342805.642809`},{3933342842.879986`,3933342844.6210423`},3933342921.143648`,{3933343042.85424`,3933343059.4959207`},{3933343099.0392666`,3933343170.879775`},{3933343202.586958`,3933343276.1652203`},{3933343354.2191763`,3933343392.734193`},{3933343511.178874`,3933343550.276238`}},
+CellLabel -> "In[2]:=",ExpressionUUID -> "43cebf24-b6eb-46d1-8262-c6c28e220a8b"]},
+Open],ExpressionUUID -> "4f9d7fa0-9353-414b-a302-36642d9e0b06"],Cell[
+CellGroupData[
+{Cell[
+"\<Butcher table from \[Alpha], \[Beta]\>","Chapter",CellChangeTimes -> {{3933431476.5859537`,3933431493.3749733`}},
+ExpressionUUID -> "22b11d87-dc5a-4a2a-96b3-fefc4ad5cdd2"],Cell[
+BoxData[RowBox[{RowBox[
+{RowBox[{"butcher2S","[","m_","]"}],":=",RowBox[{"Module","[",RowBox[{RowBox[{"{","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Alpha]","=",RowBox[{"SparseArray","[",RowBox[{RowBox[{"{","\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{RowBox[{"{",RowBox[{"i_",",","j_"}],"}"}],"/;",RowBox[{"i","==",RowBox[
+{"j","+","2"}]}]}],"->",RowBox[{"\[Alpha]ij","[",RowBox[{"i",",","j"}],"]"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{"{",RowBox[{"i_",",","j_"}],"}"}],"/;"," ",RowBox[{RowBox[{"i",">","2"}],"&&",RowBox[
+{"i","==",RowBox[{"j","+","1"}]}]}]}],"->",RowBox[{"\[Alpha]ij","[",RowBox[{"i",",","j"}],"]"}]}]}],"\[IndentingNewLine]","}"}],","," ",RowBox[
+{"{",RowBox[{RowBox[{"m","+","1"}],","," ","m"}],"}"}]}],"]"}]}],",","\[IndentingNewLine]",RowBox[
+{"\[Beta]","=",RowBox[{"SparseArray","[",RowBox[{RowBox[{"{","\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{RowBox[{"{",RowBox[{"i_",",","j_"}],"}"}],"/;",RowBox[{"i","==",RowBox[
+{"j","+","1"}]}]}],"->",RowBox[{"\[Beta]ij","[",RowBox[{"i",",","j"}],"]"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{"{",RowBox[{"i_",",","j_"}],"}"}],"/;",RowBox[{"i","==",RowBox[{"j","+","2"}]}]}],"->",RowBox[
+{"\[Beta]ij","[",RowBox[{"i",",","j"}],"]"}]}]}],"\[IndentingNewLine]","}"}],",",RowBox[
+{"{",RowBox[{RowBox[{"m","+","1"}],",","m"}],"}"}]}],"]"}]}],",","\[IndentingNewLine]","\[Alpha]0",",","\[Alpha]1",",","\[Beta]0",",","\[Beta]1",",","A",",","b",",","c"}],"}"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{"\[Alpha]0","=",RowBox[{"\[Alpha]","[",RowBox[{"[",RowBox[{RowBox[
+{";;","m"}],",","All"}],"]"}],"]"}]}],";","\[IndentingNewLine]",RowBox[{"\[Beta]0","=",RowBox[
+{"\[Beta]","[",RowBox[{"[",RowBox[{RowBox[{";;","m"}],",","All"}],"]"}],"]"}]}],";","\[IndentingNewLine]",RowBox[
+{"\[Alpha]1","=",RowBox[{"\[Alpha]","[",RowBox[{"[",RowBox[{RowBox[{RowBox[{"m","+","1"}],";;"}],",","All"}],"]"}],"]"}]}],";","\[IndentingNewLine]",RowBox[
+{"\[Beta]1","=",RowBox[{"\[Beta]","[",RowBox[{"[",RowBox[{RowBox[{RowBox[{"m","+","1"}],";;"}],",","All"}],"]"}],"]"}]}],";"}],"\[IndentingNewLine]",RowBox[
+{"(*"," ",RowBox[{"Equations"," ","9","a"," ","and"," ","9","b"}]," ","*)"}],";","\[IndentingNewLine]",RowBox[
+{"A","=",RowBox[{RowBox[{"Inverse","[",RowBox[{RowBox[{"IdentityMatrix","[","m","]"}],"-","\[Alpha]0"}],"]"}],".","\[Beta]0"}]}],";","\[IndentingNewLine]",RowBox[
+{"b","=",RowBox[{RowBox[{"Transpose","[",RowBox[{"\[Beta]1","+",RowBox[{"\[Alpha]1",".","A"}]}],"]"}],"[",RowBox[
+{"[",RowBox[{"All",",","1"}],"]"}],"]"}]}],";","\[IndentingNewLine]",RowBox[{"c","=",RowBox[
+{"Plus","@@@","A"}]}],";","\[IndentingNewLine]",RowBox[{"{",RowBox[{"A",",","b",",","c"}],"}"}]}]}],"]"}]}],"\[IndentingNewLine]"}]],
+"Input",CellChangeTimes -> {{3933342544.173015`,3933342557.6832466`},{3933342926.48391`,3933342963.9052906`},{3933343071.138041`,3933343072.561965`},{3933343286.095475`,3933343336.691963`},3933343397.5969543`,{3933343433.2680583`,3933343468.636072`},{3933343561.2381763`,3933343605.468481`},{3933345349.796748`,3933345450.195347`},{3933345485.539579`,3933345698.9406357`},{3933345742.86693`,3933345778.218711`},{3933345924.340559`,3933345940.626734`},{3933346023.1474338`,3933346030.219181`},{3933346656.9510098`,3933346662.15865`},{3933347097.005034`,3933347107.406505`}},
+CellLabel -> "In[10]:=",ExpressionUUID -> "01548abc-a98a-4f8b-90d2-8637f55a6c26"]},
+Open],ExpressionUUID -> "5471d1fd-8676-4097-9cdf-e59ba727d56a"],Cell[
+CellGroupData[
+{Cell[
+"\<Low Storage coefficients for different schemes\>","Chapter",CellChangeTimes -> {{3933349253.0424547`,3933349274.8627`}},
+ExpressionUUID -> "4af157fc-17f1-46e5-a272-f102fe2d6825"],Cell[
+BoxData[RowBox[{RowBox[
+{"scheme","=",RowBox[{"<|","\[IndentingNewLine]",RowBox[{RowBox[{"rk2","->"," ",RowBox[
+{"{","\[IndentingNewLine]",RowBox[{RowBox[{"m","->","2"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"2",",","1"}],"]"}],"->","1"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"3",",","2"}],"]"}],"->",RowBox[{"1","/","2"}]}],",","\n","          ",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"3",",","2"}],"]"}],"->",RowBox[{"1","/","2"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Delta]","[","2","]"}],"->","0"}]}],"\[IndentingNewLine]","}"}]}],",","\[IndentingNewLine]",RowBox[
+{"vl2","->",RowBox[{"{","\[IndentingNewLine]",RowBox[{RowBox[{"m","->","2"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"2",",","1"}],"]"}],"->",RowBox[{"1","/","2"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"3",",","2"}],"]"}],"->","1"}],",","\n","\t",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"3",",","2"}],"]"}],"->","1"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Delta]","[","2","]"}],"->","0"}]}],"\[IndentingNewLine]","}"}]}],",","\[IndentingNewLine]",RowBox[
+{"rk3","->",RowBox[{"{","\[IndentingNewLine]",RowBox[{RowBox[{"m","->"," ","3"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Delta]","[","2","]"}],"->","0"}],",","\[IndentingNewLine]",RowBox[{RowBox[
+{"\[Delta]","[","3","]"}],"->","0"}],",","\[IndentingNewLine]",RowBox[{RowBox[{"\[Beta]ij","[",RowBox[
+{"2",",","1"}],"]"}],"->","1"}],",","\[IndentingNewLine]",RowBox[{RowBox[{"\[Beta]ij","[",RowBox[
+{"3",",","2"}],"]"}],"->",RowBox[{"1","/","4"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"4",",","3"}],"]"}],"->",RowBox[{"2","/","3"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"3",",","2"}],"]"}],"->",RowBox[{"3","/","4"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"4",",","2"}],"]"}],"->",RowBox[{"1","/","3"}]}]}],"\[IndentingNewLine]","}"}]}],",","\[IndentingNewLine]",RowBox[
+{"rk45","->",RowBox[{"{","\[IndentingNewLine]",RowBox[{RowBox[{"m","->","5"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"2",",","2"}],"]"}],"->","1"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"3",",","2"}],"]"}],"->","4.666545952121251"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"4",",","2"}],"]"}],"->","0.964197464041912"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"5",",","2"}],"]"}],"->"," ",RowBox[{"-","3.398279365655790"}]}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Gamma]","[",RowBox[{"6",",","2"}],"]"}],"->","0.229588412671583"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"2",",","1"}],"]"}],"->","0.357534921136978"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"3",",","2"}],"]"}],"->","2.364680399061355"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"4",",","3"}],"]"}],"->","0.016239790859612"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"5",",","4"}],"]"}],"->","0.498173799587251"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"6",",","5"}],"]"}],"->","0.433334235669763"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Delta]","[","2","]"}],"->","0"}],",","\[IndentingNewLine]",RowBox[{RowBox[
+{"\[Delta]","[","3","]"}],"->","0"}],",","\[IndentingNewLine]",RowBox[{RowBox[{"\[Delta]","[","4","]"}],"->","0"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"\[Delta]","[","5","]"}],"->","0"}],",","\[IndentingNewLine]",RowBox[{RowBox[
+{"\[Delta]","[","6","]"}],"->","0"}]}],"\[IndentingNewLine]","}"}]}]}],"\[IndentingNewLine]","|>"}]}],";"}]],
+"Input",CellChangeTimes -> {{3933349298.7741137`,3933349462.4821453`},{3933349492.6990557`,3933349514.3585224`},{3933349547.904687`,3933349616.3705506`},{3933349679.864212`,3933349681.8224382`},{3933351106.428812`,3933351143.1348677`},{3933351189.140523`,3933351377.140004`},{3933351436.3674917`,3933351545.6507707`},{3933351603.7029676`,3933351708.463834`},3933352267.851912`,3933352782.149241`},
+CellLabel -> "In[11]:=",ExpressionUUID -> "6e0706a9-34d0-45d8-be04-7222d8a1a077"]},
+Open],ExpressionUUID -> "3d6768ad-4a9b-4ba2-b4c9-4ed6f208a0da"],Cell[
+CellGroupData[
+{Cell[
+"\<Make a Low-Storage RK method from coefficients\>","Chapter",CellChangeTimes -> {{3933431529.519893`,3933431551.03344`}},
+ExpressionUUID -> "d3c41edf-8536-48a2-beda-277d21e26c81"],Cell[
+BoxData[{RowBox[{RowBox[
+{"makeIntegrator","[","key_","]"}],":=",RowBox[{"Module","[",RowBox[{RowBox[{"{","\[IndentingNewLine]",RowBox[
+{RowBox[{"m","=",RowBox[{RowBox[{RowBox[{RowBox[{"Key","[","key","]"}],"[","scheme","]"}],"[",RowBox[
+{"[","1","]"}],"]"}],"[",RowBox[{"[","2","]"}],"]"}]}],",","\[IndentingNewLine]",RowBox[
+{"subs","=",RowBox[{RowBox[{RowBox[{"Key","[","key","]"}],"[","scheme","]"}],"[",RowBox[
+{"[",RowBox[{"2",";;"}],"]"}],"]"}]}],",","\[IndentingNewLine]","butcher"}],"}"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"butcher","=",RowBox[{RowBox[{"butcher2S","[","m","]"}],"/.","subs"}]}],";","\[IndentingNewLine]",RowBox[
+{"<|","\[IndentingNewLine]",RowBox[{RowBox[{RowBox[{"beta","->",RowBox[{"Table","[",RowBox[
+{RowBox[{"\[Beta]ij","[",RowBox[{"i",",",RowBox[{"i","-","1"}]}],"]"}],",",RowBox[
+{"{",RowBox[{"i",",","2",",",RowBox[{"m","+","1"}]}],"}"}]}],"]"}]}],"/.","subs"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"delta","->",RowBox[{"Table","[",RowBox[{RowBox[{"\[Delta]","[",RowBox[{"i","-","1"}],"]"}],",",RowBox[
+{"{",RowBox[{"i",",","2",",",RowBox[{"m","+","1"}]}],"}"}]}],"]"}]}],"/.","subs"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"gam0","->",RowBox[{"Table","[",RowBox[{RowBox[{"\[Gamma]","[",RowBox[{"i",",","1"}],"]"}],",",RowBox[
+{"{",RowBox[{"i",",","2",",",RowBox[{"m","+","1"}]}],"}"}]}],"]"}]}],"/.","subs"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{"gam1","->",RowBox[{"Table","[",RowBox[{RowBox[{"\[Gamma]","[",RowBox[{"i",",","2"}],"]"}],",",RowBox[
+{"{",RowBox[{"i",",","2",",",RowBox[{"m","+","1"}]}],"}"}]}],"]"}]}],"/.","subs"}],",","\[IndentingNewLine]",RowBox[
+{"c","->",RowBox[{"butcher","[",RowBox[{"[","3","]"}],"]"}]}]}],"\[IndentingNewLine]","|>"}]}]}],"\[IndentingNewLine]","]"}]}],"\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{"rkInteriorStep","[",RowBox[{RowBox[{"{",RowBox[{"s1_",",","s2_",",","time_",","," ","dt_",",","rhs_"}],"}"}],","," ",RowBox[
+{"{",RowBox[{"beta_",",","gam0_",",","gam1_",",","delta_",",","c_"}],"}"}]}],"]"}],":=","\[IndentingNewLine]",RowBox[
+{"Module","[",RowBox[{RowBox[{"{",RowBox[{"s2buf","=",RowBox[{"s2","+",RowBox[{"delta"," ","s1"}]}]}],"}"}],",","\[IndentingNewLine]",RowBox[
+{"{",RowBox[{RowBox[{RowBox[{"gam0"," ","s1"}],"+",RowBox[{"gam1"," ","s2buf"}]," ","+",RowBox[
+{"beta"," ","dt"," ",RowBox[{"rhs","[",RowBox[{"s1",",",RowBox[{"time","+",RowBox[
+{"c"," ","dt"}]}]}]," ","]"}]}]}],","," ","s2buf",","," ","time",","," ","dt",","," ","rhs"}],"}"}]}],"]"}]}],"\[IndentingNewLine]"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{"rkStep","[",RowBox[{"yn_",",","time_",",","dt_",","," ","rhs_",",","beta_",",","gam0_",",","gam1_",",","delta_",",","c_"}]," ","]"}],":=","\[IndentingNewLine]",RowBox[
+{"Fold","[",RowBox[{"rkInteriorStep",","," ",RowBox[{"{",RowBox[{"yn",",","0",",","time",",","dt",",","rhs"}],"}"}],","," ",RowBox[
+{"Transpose","[",RowBox[{"List","[",RowBox[{"beta",",","gam0",",","gam1",",","delta",",","c"}],"]"}],"]"}]}],"]"}]}],"\[IndentingNewLine]"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"rkIntegrate","[",RowBox[{"y0_",",","t0_",",","tn_",",","dt_",",","rhs_",",","rkInfo_"}],"]"}],":=",RowBox[
+{"Module","[",RowBox[{RowBox[{"{","\[IndentingNewLine]",RowBox[{"rkFunc","=",RowBox[
+{"Function","[",RowBox[{RowBox[{"{",RowBox[{"y",",","t"}],"}"}],",",RowBox[{"First","[",RowBox[
+{"rkStep","[",RowBox[{"y",",","t",",","dt",",","rhs",",",RowBox[{RowBox[{"Key","[","beta","]"}],"[","rkInfo","]"}],",",RowBox[
+{RowBox[{"Key","[","gam0","]"}],"[","rkInfo","]"}],",",RowBox[{RowBox[{"Key","[","gam1","]"}],"[","rkInfo","]"}],",",RowBox[
+{RowBox[{"Key","[","delta","]"}],"[","rkInfo","]"}],",",RowBox[{RowBox[{"Key","[","c","]"}],"[","rkInfo","]"}]}],"]"}],"]"}]}],"]"}]}],"\[IndentingNewLine]","}"}],",","\[IndentingNewLine]",RowBox[
+{"Transpose","[",RowBox[{"List","[",RowBox[{RowBox[{"Range","[",RowBox[{"t0",",","tn",",","dt"}],"]"}],",",RowBox[
+{"FoldList","[",RowBox[{"rkFunc",",","y0",",",RowBox[{"Range","[",RowBox[{"t0",",",RowBox[
+{"tn","-","dt"}],",","dt"}],"]"}]}],"]"}]}],"]"}],"]"}]}],"\[IndentingNewLine]","]"}]}]}],
+"Input",CellChangeTimes -> {{3933352274.2790194`,3933352496.842067`},{3933352610.3748407`,3933352663.103447`},{3933352717.881524`,3933352751.332842`},{3933352828.5289087`,3933352834.065763`},{3933352907.569047`,3933352908.57596`},{3933352963.7358217`,3933352968.7869663`},{3933415139.4766397`,3933415319.8521843`},{3933415363.8431273`,3933415385.3331003`},{3933431555.302792`,3933431583.660833`},{3933432296.705492`,3933432305.2384553`},{3933432337.6400623`,3933432358.49339`},{3933432427.7060084`,3933432429.327286`}},
+CellLabel -> "In[22]:=",ExpressionUUID -> "6cd24dc7-b01a-4b10-a7ae-b14670710a1e"],Cell[
+
+"\<Example:\>","Text",CellChangeTimes -> {{3933432251.855978`,3933432252.864431`}},
+ExpressionUUID -> "ca3fc76b-a4e0-47eb-be32-7bba83b6a8ef"],Cell[
+BoxData[RowBox[{"makeIntegrator","[","rk45","]"}]],
+"Input",CellChangeTimes -> {{3933432257.0367227`,3933432266.51719`},{3933432363.259603`,3933432365.329129`}},
+CellLabel -> "In[26]:=",ExpressionUUID -> "246e534d-a354-4595-8207-8efbef496a07"]},
+Open],ExpressionUUID -> "fe0f079d-6baa-4a6e-95e7-b8f76edef3e3"],Cell[
+CellGroupData[
+{Cell[
+"\<Sanity Check\>","Chapter",CellChangeTimes -> {{3933431595.133444`,3933431601.2454133`}},
+ExpressionUUID -> "6860462f-ca26-43a7-b875-4e3cfcb16c41"],Cell[
+BoxData[{RowBox[{RowBox[
+{"s","=",RowBox[{"NDSolve","[",RowBox[{RowBox[{"{",RowBox[{RowBox[{RowBox[{RowBox[
+{"y","'"}],"[","x","]"}],"==",RowBox[{RowBox[{"y","[","x","]"}]," ",RowBox[{"Cos","[",RowBox[
+{RowBox[{"y","[","x","]"}],"+","x"}],"]"}]}]}],",",RowBox[{RowBox[{"y","[","0","]"}],"==","1"}]}],"}"}],",","y",",",RowBox[
+{"{",RowBox[{"x",",","0",",","30"}],"}"}]}],"]"}]}],";"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"p1"," ","="," ",RowBox[{"Plot","[",RowBox[{RowBox[{"Evaluate","[",RowBox[
+{RowBox[{"y","[","x","]"}],"/.","s"}],"]"}],",",RowBox[{"{",RowBox[{"x",",","0",",","30"}],"}"}],",",RowBox[
+{"PlotRange","->","All"}],","," ",RowBox[{"PlotStyle","->","Red"}]}],"]"}]}],";"}],"\[IndentingNewLine]",RowBox[
+{RowBox[{"p2"," ","="," ",RowBox[{"ListPlot","[",RowBox[{RowBox[{"rkIntegrate","[",RowBox[
+{"1",",","0",",","30",",","0.5",",",RowBox[{"Function","[",RowBox[{RowBox[{"{",RowBox[
+{"y",",","t"}],"}"}],",",RowBox[{"y"," ",RowBox[{"Cos","[",RowBox[{"t","+","y"}],"]"}]}]}],"]"}],",",RowBox[
+{"makeIntegrator","[","vl2","]"}]}],"]"}],",",RowBox[{"PlotStyle","->","Black"}]}],"]"}]}],";"}],"\[IndentingNewLine]",RowBox[
+{"Show","[",RowBox[{"p1",",","p2"}],"]"}]}],"Input",CellChangeTimes -> {{3932132868.5792117`,3932132868.58528`},{3932132948.5663033`,3932132949.082738`},{3932133314.530981`,3932133324.9773607`},{3932207392.077511`,3932207394.663705`},{3933431608.965467`,3933431707.9118237`},{3933432438.919209`,3933432440.836319`}},
+CellLabel -> "In[27]:=",ExpressionUUID -> "bba2b8fa-1ef3-4355-ae5c-d1a48ea30cd9"]},
+Open],ExpressionUUID -> "e9a6bb6e-848f-40b0-94cf-59fbb3a3b072"],Cell[
+CellGroupData[
+{Cell[
+"\<Manufactured Solution\>","Chapter",CellChangeTimes -> {{3933431733.18156`,3933431737.692576`}},
+ExpressionUUID -> "cd56dc27-c307-4105-b857-f67b42a299e1"],Cell[
+"\<Chose a simple time-varying MS.  Chose an ODE whose r.h.s. is not a function of time.  This way all of the time information comes through the MS source term.  This seems like it should eliminate any fortuitous error cancellation from having the time dependence in multiple places?\>",
+"Text",CellChangeTimes -> {{3933431841.6903787`,3933431950.231882`}},ExpressionUUID -> "9792a976-225f-4011-86dd-a425d96996da"],Cell[
+BoxData[
+RowBox[{RowBox[{"(*"," ",RowBox[{RowBox[{"y","'"}],"="," ",RowBox[{"rhs","[",RowBox[
+{"y",",","t"}],"]"}]}]," ","*)"}],"\[IndentingNewLine]",RowBox[{RowBox[{RowBox[{"rhs","[",RowBox[
+{"y_",",","t_"}],"]"}],":=",RowBox[{"y"," ",RowBox[{"Cos","[","y","]"}]}]}]," ","\[IndentingNewLine]",RowBox[
+{RowBox[{"yMs","[","t_","]"}],":=",RowBox[{"1","+",RowBox[{RowBox[{"Sin","[","t","]"}],"/","2"}]}]}]}]}]],
+"Input",CellChangeTimes -> {{3932207456.145615`,3932207468.0975094`},{3932208364.561344`,3932208471.9323378`},3932208984.1193123`,{3932209233.476224`,3932209258.633163`},{3933433229.070273`,3933433284.047635`}},
+CellLabel -> "In[48]:=",ExpressionUUID -> "2de5960c-8a48-4e20-8369-557b18e8a543"],Cell[
+BoxData[
+RowBox[{RowBox[{"rkError","[",RowBox[{"dt_",",","rkInfo_",",",RowBox[{"levels_",":","3"}]}],"]"}],":=",RowBox[
+{"Module","[",RowBox[{RowBox[{"{","\[IndentingNewLine]",RowBox[{RowBox[{"results","=",RowBox[
+{"Table","[",RowBox[{RowBox[{"rkIntegrate","[",RowBox[{RowBox[{"yMs","[","0","]"}],",","0",",","200",",",RowBox[
+{"dt","/",RowBox[{"2","^","i"}]}],",",RowBox[{"Function","[",RowBox[{RowBox[{"{",RowBox[
+{"y",",","t"}],"}"}],",",RowBox[{RowBox[{"(",RowBox[{RowBox[{"D","[",RowBox[{RowBox[
+{"yMs","[","x","]"}],",","x"}],"]"}],"/.",RowBox[{"x","->","t"}]}],")"}],"-",RowBox[
+{"rhs","[",RowBox[{RowBox[{"yMs","[","t","]"}],",","t"}],"]"}],"+",RowBox[{"rhs","[",RowBox[
+{"y",",","t"}],"]"}]}]}],"]"}],",","rkInfo"}]," ","]"}],",",RowBox[{"{",RowBox[{"i",",","0",",",RowBox[
+{"levels","-","1"}]}],"}"}]}],"]"}]}],",","\[IndentingNewLine]","errors",",","\[IndentingNewLine]","linf"}],"\[IndentingNewLine]","}"}],",","\[IndentingNewLine]",RowBox[
+{RowBox[{RowBox[{"errors","=",RowBox[{"Map","[",RowBox[{RowBox[{"Function","[",RowBox[
+{RowBox[{"{","args","}"}],",",RowBox[{"{",RowBox[{RowBox[{"args","[",RowBox[{"[","1","]"}],"]"}],",",RowBox[
+{"Abs","[",RowBox[{RowBox[{"args","[",RowBox[{"[","2","]"}],"]"}],"-",RowBox[{"yMs","[",RowBox[
+{"args","[",RowBox[{"[","1","]"}],"]"}],"]"}]}],"]"}]}],"}"}]}],"]"}],",","results",",",RowBox[
+{"{","2","}"}]}],"]"}]}],";","\[IndentingNewLine]",RowBox[{"Print","[",RowBox[{"ListPlot","/@","errors"}],"]"}],";","\[IndentingNewLine]",RowBox[
+{"linf","=",RowBox[{RowBox[{"Map","[",RowBox[{"Function","[",RowBox[{RowBox[{"{","err","}"}],",",RowBox[
+{"Max","[",RowBox[{"err","[",RowBox[{"[",RowBox[{"All",",","2"}],"]"}],"]"}],"]"}]}],"]"}],"]"}],"[","errors","]"}]}],";"}],"\[IndentingNewLine]",RowBox[
+{"(*"," ",RowBox[{"Order"," ","of"," ","temporal"," ","convergence"}]," ","*)"}],";","\[IndentingNewLine]",RowBox[
+{RowBox[{"Map","[",RowBox[{RowBox[{"Log","[",RowBox[{"2",",","#"}],"]"}],"&"}],"]"}],"[",RowBox[
+{RowBox[{"linf","[",RowBox[{"[",RowBox[{"1",";;",RowBox[{"levels","-","1"}]}],"]"}],"]"}],"/",RowBox[
+{"linf","[",RowBox[{"[",RowBox[{"2",";;","levels"}],"]"}],"]"}]}],"]"}]}]}],"\[IndentingNewLine]","]"}]}]],
+"Input",CellChangeTimes -> {{3932217127.038363`,3932217252.61661`},{3932217291.104064`,3932217315.840137`},{3932217436.2460036`,3932217466.72156`},{3932217512.8398857`,3932217513.256936`},{3932217615.848021`,3932217629.033698`},{3933431967.554924`,3933432160.9799633`},{3933432468.907331`,3933432533.36943`},{3933433271.819363`,3933433301.1822567`}},
+CellLabel -> "In[51]:=",ExpressionUUID -> "961eaaf4-9d04-4944-b92c-0111a205f0b5"],Cell[
+BoxData[
+RowBox[{"rkError","[",RowBox[{"0.5",",",RowBox[{"makeIntegrator","[","rk45","]"}],",","4"}],"]"}]],
+"Input",CellChangeTimes -> {{3933432155.991725`,3933432203.900763`},{3933432452.673057`,3933432454.930271`},{3933432511.084322`,3933432511.4058857`},{3933432550.377314`,3933432569.4454384`}},
+CellLabel -> "In[52]:=",ExpressionUUID -> "2e3aadee-9f34-497d-9f42-10bd8ae2d7bb"],Cell[
+BoxData[
+RowBox[{"rkError","[",RowBox[{"0.5",",",RowBox[{"makeIntegrator","[","rk3","]"}],",","4"}],"]"}]],
+"Input",CellChangeTimes -> {{3933432956.625321`,3933432957.123412`}},CellLabel -> "In[53]:=",
+ExpressionUUID -> "2499f57e-ebb6-43c5-94d0-614c0d6fd8bf"],Cell[
+BoxData[RowBox[{"rkError","[",RowBox[
+{"0.5",",",RowBox[{"makeIntegrator","[","rk2","]"}],",","4"}],"]"}]],"Input",CellChangeTimes -> {{3933432968.130974`,3933432968.3723774`}},
+CellLabel -> "In[54]:=",ExpressionUUID -> "b375d17b-9736-456e-a790-fd697b9cff4c"],Cell[
+BoxData[
+RowBox[{"rkError","[",RowBox[{"0.5",",",RowBox[{"makeIntegrator","[","vl2","]"}],",","4"}],"]"}]],
+"Input",CellChangeTimes -> {{3933432979.809408`,3933432980.5701337`}},CellLabel -> "In[55]:=",
+ExpressionUUID -> "00aee915-7b67-4859-b1e6-0c237e1859b6"]},Open],ExpressionUUID -> "4a75b700-ce54-4703-af7e-e192431e418a"]},
+StyleDefinitions -> "Default.nb",WindowSize -> {1267,1253},WindowMargins -> {{464,Automatic},{Automatic,0}},
+FrontEndVersion -> "14.1 for Wolfram Cloud 1.69 (July 16, 2024)",ExpressionUUID -> "c844a8dd-7ee1-484c-b825-8ff4cd1b0d06"]
diff --git a/src/time_integration/low_storage_integrator.cpp b/src/time_integration/low_storage_integrator.cpp
index d833ddd8ad20..4d30e70fa97c 100644
--- a/src/time_integration/low_storage_integrator.cpp
+++ b/src/time_integration/low_storage_integrator.cpp
@@ -57,11 +57,13 @@ LowStorageIntegrator::LowStorageIntegrator(const std::string &name)
     beta.resize(nstages);
     gam0.resize(nstages);
     gam1.resize(nstages);
+    c.resize(nstages);
 
     delta[0] = 1.0;
     beta[0] = 1.0;
     gam0[0] = 0.0;
     gam1[0] = 1.0;
+    c[0] = 0.0;
   } else if (name_ == "rk2") {
     nstages = 2;
     nbuffers = 2;
@@ -69,16 +71,19 @@ LowStorageIntegrator::LowStorageIntegrator(const std::string &name)
     beta.resize(nstages);
     gam0.resize(nstages);
     gam1.resize(nstages);
+    c.resize(nstages);
 
     delta[0] = 1.0;
     beta[0] = 1.0;
     gam0[0] = 0.0;
     gam1[0] = 1.0;
+    c[0] = 0.0;
 
     delta[1] = 0.0;
     beta[1] = 0.5;
     gam0[1] = 0.5;
     gam1[1] = 0.5;
+    c[1] = 1.0;
   } else if (name_ == "vl2") {
     nstages = 2;
     nbuffers = 2;
@@ -86,16 +91,19 @@ LowStorageIntegrator::LowStorageIntegrator(const std::string &name)
     beta.resize(nstages);
     gam0.resize(nstages);
     gam1.resize(nstages);
+    c.resize(nstages);
 
     delta[0] = 1.0;
     beta[0] = 0.5;
     gam0[0] = 0.0;
     gam1[0] = 1.0;
+    c[0] = 0.0;
 
     delta[1] = 0.0;
     beta[1] = 1.0;
     gam0[1] = 0.0;
     gam1[1] = 1.0;
+    c[1] = 0.5;
   } else if (name_ == "rk3") {
     nstages = 3;
     nbuffers = 2;
@@ -103,21 +111,25 @@ LowStorageIntegrator::LowStorageIntegrator(const std::string &name)
     beta.resize(nstages);
     gam0.resize(nstages);
     gam1.resize(nstages);
+    c.resize(nstages);
 
     delta[0] = 1.0;
     beta[0] = 1.0;
     gam0[0] = 0.0;
     gam1[0] = 1.0;
+    c[0] = 0.0;
 
     delta[1] = 0.0;
     beta[1] = 0.25;
     gam0[1] = 0.25;
     gam1[1] = 0.75;
+    c[1] = 1.0;
 
     delta[2] = 0.0;
     beta[2] = 2.0 / 3.0;
     gam0[2] = 2.0 / 3.0;
     gam1[2] = 1.0 / 3.0;
+    c[2] = 0.5;
   } else if (name_ == "rk4") {
     // Classic 5-stage SSPRK(5)4 in low-storage form
     // ceff = 0.377
@@ -128,31 +140,38 @@ LowStorageIntegrator::LowStorageIntegrator(const std::string &name)
     beta.resize(nstages);
     gam0.resize(nstages);
     gam1.resize(nstages);
+    c.resize(nstages);
 
     delta[0] = 1.0;
     beta[0] = 0.357534921136978;
     gam0[0] = 0.0;
     gam1[0] = 1.0;
+    c[0] = 0.0;
 
     delta[1] = 0.0;
     beta[1] = 2.364680399061355;
     gam0[1] = -3.666545952121251;
     gam1[1] = 4.666545952121251;
+    c[1] = 0.357534921136978;
 
     delta[2] = 0.0;
     beta[2] = 0.016239790859612;
     gam0[2] = 0.035802535958088;
     gam1[2] = 0.964197464041912;
+    c[2] = 1.0537621812245777;
 
     delta[3] = 0.0;
     beta[3] = 0.498173799587251;
     gam0[3] = 4.398279365655791;
     gam1[3] = -3.398279365655790;
+    c[3] = 0.05396714924417825;
 
     delta[4] = 0.0;
     beta[4] = 0.433334235669763;
     gam0[4] = 0.770411587328417;
     gam1[4] = 0.229588412671583;
+    c[4] = 0.7355363985311864;
+
   } else {
     throw std::invalid_argument("Invalid selection for the time integrator: " + name_);
   }
diff --git a/src/time_integration/staged_integrator.hpp b/src/time_integration/staged_integrator.hpp
index 329f69bc294b..e2dbe570150a 100644
--- a/src/time_integration/staged_integrator.hpp
+++ b/src/time_integration/staged_integrator.hpp
@@ -55,6 +55,8 @@ class LowStorageIntegrator : public StagedIntegrator {
   std::vector<Real> beta;
   std::vector<Real> gam0;
   std::vector<Real> gam1;
+  // Time offset of each stage. See scripts/mathematica/sparse_integrators.nb
+  std::vector<Real> c;
 };
 
 // TODO(JMM): Should this be named ButcherTableauIntegrator?

From 517bf9211217b2a96834cac3d061ea103aef8594 Mon Sep 17 00:00:00 2001
From: Ben Ryan <brryan@lanl.gov>
Date: Wed, 28 Aug 2024 10:43:17 -0600
Subject: [PATCH 39/40] Fix remaining swarm D->H->D copies (#1145)

* working

* Clean up implicit captures

* OK something isnt working with updating empty_indices

* Creating indices works but defragmentation is broken

* Defrag works

* Switch to persistent scratch memory

* Clean up

* Fix GPU issues

* Fix compile error by cleaning up code

* Remove unnecessary check against non-null user swarm BCs

* Remove unused function

* Formatting

* Fiddle with send logic

* Perform swarm boundary logic on device (#1154)

* First kernel down

* Further cleanup

* Notes

* Send seems to work

* Need to promote particles example to latest pattern

* May have fixed particles example

* cycles...

* iterative tasking is only iterating on one meshblock in particles example...

* Still not working...

* New loop seems to work

* Cleaned up some printfs/marked unused code for deletion

* New algorithm is cycling, time to clean up

* Fixed indexing bug...

* Clean up

* finished_transport shouldnt be provided by swarm

* Reverting to old manual iterative tasking

* Still working...

* Starting to make progress...

* Cleaned up

* format

* A few leftover print statements

* Oops ParArray1D isn't a host array when compiled for device

* implicit this->

* bug in nrecvd particles with 1 particle received...

* Found the bug

* Fixed bug, cleaned up

* kokkos parallel_scan -> par_scan

* Fix types, clean up code

* Add warning if using swarms but Real != double

* Cleanup from debugging

---------

Co-authored-by: Luke Roberts <lfroberts@lanl.gov>
---
 CHANGELOG.md                           |   1 +
 example/particles/parthinput.particles |   6 +-
 example/particles/particles.cpp        | 180 ++++++++---------
 example/particles/particles.hpp        |   4 +-
 src/interface/state_descriptor.cpp     |   7 +
 src/interface/swarm.cpp                | 259 +++++++++++++-----------
 src/interface/swarm.hpp                |  38 ++--
 src/interface/swarm_comms.cpp          | 268 ++++++++-----------------
 src/interface/swarm_container.cpp      |  30 ---
 src/interface/swarm_container.hpp      |   2 -
 src/mesh/forest/tree.cpp               |   2 -
 src/mesh/meshblock.cpp                 |   3 -
 src/tasks/tasks.hpp                    |   2 +
 tst/unit/test_swarm.cpp                |   4 +-
 14 files changed, 338 insertions(+), 468 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4f474633ce8e..af70ab73376b 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -40,6 +40,7 @@
 - [[PR 1004]](https://github.com/parthenon-hpc-lab/parthenon/pull/1004) Allow parameter modification from an input file for restarts
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1145]](https://github.com/parthenon-hpc-lab/parthenon/pull/1145) Fix remaining swarm D->H->D copies
 - [[PR 1150]](https://github.com/parthenon-hpc-lab/parthenon/pull/1150) Reduce memory consumption for buffer pool
 - [[PR 1146]](https://github.com/parthenon-hpc-lab/parthenon/pull/1146) Fix an issue outputting >4GB single variables per rank
 - [[PR 1152]](https://github.com/parthenon-hpc-lab/parthenon/pull/1152) Fix memory leak in task graph outputs related to `abi::__cxa_demangle`
diff --git a/example/particles/parthinput.particles b/example/particles/parthinput.particles
index bca9ec93292e..a79f774c6850 100644
--- a/example/particles/parthinput.particles
+++ b/example/particles/parthinput.particles
@@ -24,10 +24,8 @@ refinement = none
 nx1 = 16
 x1min = -0.5
 x1max = 0.5
-ix1_bc = user
-ox1_bc = user
-# ix1_bc = periodic # Optionally use periodic boundary conditions everywhere
-# ox1_bc = periodic
+ix1_bc = periodic
+ox1_bc = periodic
 
 nx2 = 16
 x2min = -0.5
diff --git a/example/particles/particles.cpp b/example/particles/particles.cpp
index 8a8d47ef111a..3f99407ded54 100644
--- a/example/particles/particles.cpp
+++ b/example/particles/particles.cpp
@@ -340,8 +340,7 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
   return TaskStatus::complete;
 }
 
-TaskStatus TransportParticles(MeshBlock *pmb, const StagedIntegrator *integrator,
-                              const double t0) {
+TaskStatus TransportParticles(MeshBlock *pmb, const double t0, const double dt) {
   PARTHENON_INSTRUMENT
 
   auto swarm = pmb->meshblock_data.Get()->GetSwarmData()->Get("my_particles");
@@ -350,8 +349,6 @@ TaskStatus TransportParticles(MeshBlock *pmb, const StagedIntegrator *integrator
 
   int max_active_index = swarm->GetMaxActiveIndex();
 
-  Real dt = integrator->dt;
-
   auto &t = swarm->Get<Real>("t").Get();
   auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
   auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
@@ -469,97 +466,31 @@ TaskStatus TransportParticles(MeshBlock *pmb, const StagedIntegrator *integrator
 // Custom step function to allow for looping over MPI-related tasks until complete
 TaskListStatus ParticleDriver::Step() {
   TaskListStatus status;
-  integrator.dt = tm.dt;
+
+  PARTHENON_REQUIRE(integrator.nstages == 1,
+                    "Only first order time integration supported!");
 
   BlockList_t &blocks = pmesh->block_list;
   auto num_task_lists_executed_independently = blocks.size();
 
   // Create all the particles that will be created during the step
   status = MakeParticlesCreationTaskCollection().Execute();
+  PARTHENON_REQUIRE(status == TaskListStatus::complete,
+                    "ParticlesCreation task list failed!");
 
-  // Loop over repeated MPI calls until every particle is finished. This logic is
-  // required because long-distance particle pushes can lead to a large, unpredictable
-  // number of MPI sends and receives.
-  bool particles_update_done = false;
-  while (!particles_update_done) {
-    status = MakeParticlesUpdateTaskCollection().Execute();
-
-    particles_update_done = true;
-    for (auto &block : blocks) {
-      // TODO(BRR) Despite this "my_particles"-specific call, this function feels like it
-      // should be generalized
-      auto swarm = block->meshblock_data.Get()->GetSwarmData()->Get("my_particles");
-      if (!swarm->finished_transport) {
-        particles_update_done = false;
-      }
-    }
-  }
+  // Transport particles iteratively until all particles reach final time
+  status = IterativeTransport();
+  // status = MakeParticlesTransportTaskCollection().Execute();
+  PARTHENON_REQUIRE(status == TaskListStatus::complete,
+                    "IterativeTransport task list failed!");
 
   // Use a more traditional task list for predictable post-MPI evaluations.
   status = MakeFinalizationTaskCollection().Execute();
+  PARTHENON_REQUIRE(status == TaskListStatus::complete, "Finalization task list failed!");
 
   return status;
 }
 
-// TODO(BRR) This should really be in parthenon/src... but it can't just live in Swarm
-// because of the loop over blocks
-TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
-  PARTHENON_INSTRUMENT
-
-  int num_sent_local = 0;
-  for (auto &block : blocks) {
-    auto sc = block->meshblock_data.Get()->GetSwarmData();
-    auto swarm = sc->Get("my_particles");
-    swarm->finished_transport = false;
-    num_sent_local += swarm->num_particles_sent_;
-  }
-
-  int num_sent_global = num_sent_local; // potentially overwritten by following Allreduce
-#ifdef MPI_PARALLEL
-  for (auto &block : blocks) {
-    auto swarm = block->meshblock_data.Get()->GetSwarmData()->Get("my_particles");
-    for (int n = 0; n < block->neighbors.size(); n++) {
-      NeighborBlock &nb = block->neighbors[n];
-      // TODO(BRR) May want logic like this if we have non-blocking TaskRegions
-      // if (nb.snb.rank != Globals::my_rank) {
-      //  if (swarm->vbswarm->bd_var_.flag[nb.bufid] != BoundaryStatus::completed) {
-      //    return TaskStatus::incomplete;
-      //  }
-      //}
-
-      // TODO(BRR) May want to move this logic into a per-cycle initialization call
-      if (swarm->vbswarm->bd_var_.flag[nb.bufid] == BoundaryStatus::completed) {
-        swarm->vbswarm->bd_var_.req_send[nb.bufid] = MPI_REQUEST_NULL;
-      }
-    }
-  }
-
-  MPI_Allreduce(&num_sent_local, &num_sent_global, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
-#endif // MPI_PARALLEL
-
-  if (num_sent_global == 0) {
-    for (auto &block : blocks) {
-      auto &pmb = block;
-      auto sc = pmb->meshblock_data.Get()->GetSwarmData();
-      auto swarm = sc->Get("my_particles");
-      swarm->finished_transport = true;
-    }
-  }
-
-  // Reset boundary statuses
-  for (auto &block : blocks) {
-    auto &pmb = block;
-    auto sc = pmb->meshblock_data.Get()->GetSwarmData();
-    auto swarm = sc->Get("my_particles");
-    for (int n = 0; n < pmb->neighbors.size(); n++) {
-      auto &nb = block->neighbors[n];
-      swarm->vbswarm->bd_var_.flag[nb.bufid] = BoundaryStatus::waiting;
-    }
-  }
-
-  return TaskStatus::complete;
-}
-
 TaskCollection ParticleDriver::MakeParticlesCreationTaskCollection() const {
   TaskCollection tc;
   TaskID none(0);
@@ -577,40 +508,93 @@ TaskCollection ParticleDriver::MakeParticlesCreationTaskCollection() const {
   return tc;
 }
 
-TaskCollection ParticleDriver::MakeParticlesUpdateTaskCollection() const {
+TaskStatus CountNumSent(const BlockList_t &blocks, const double tf_, bool *done) {
+  int num_unfinished = 0;
+  for (auto &block : blocks) {
+    auto sc = block->meshblock_data.Get()->GetSwarmData();
+    auto swarm = sc->Get("my_particles");
+    int max_active_index = swarm->GetMaxActiveIndex();
+
+    auto &t = swarm->Get<Real>("t").Get();
+
+    auto swarm_d = swarm->GetDeviceContext();
+
+    const auto &tf = tf_;
+
+    parthenon::par_reduce(
+        PARTHENON_AUTO_LABEL, 0, max_active_index,
+        KOKKOS_LAMBDA(const int n, int &num_unfinished) {
+          if (swarm_d.IsActive(n)) {
+            if (t(n) < tf) {
+              num_unfinished++;
+            }
+          }
+        },
+        Kokkos::Sum<int>(num_unfinished));
+  }
+
+#ifdef MPI_PARALLEL
+  MPI_Allreduce(MPI_IN_PLACE, &num_unfinished, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
+#endif // MPI_PARALLEL
+
+  if (num_unfinished > 0) {
+    *done = false;
+  } else {
+    *done = true;
+  }
+
+  return TaskStatus::complete;
+}
+
+TaskCollection ParticleDriver::IterativeTransportTaskCollection(bool *done) const {
   TaskCollection tc;
   TaskID none(0);
-  const double t0 = tm.time;
   const BlockList_t &blocks = pmesh->block_list;
+  const int nblocks = blocks.size();
+  const double t0 = tm.time;
+  const double dt = tm.dt;
 
-  auto num_task_lists_executed_independently = blocks.size();
-
-  TaskRegion &async_region0 = tc.AddRegion(num_task_lists_executed_independently);
-  for (int i = 0; i < blocks.size(); i++) {
+  TaskRegion &async_region = tc.AddRegion(nblocks);
+  for (int i = 0; i < nblocks; i++) {
     auto &pmb = blocks[i];
-
     auto &sc = pmb->meshblock_data.Get()->GetSwarmData();
+    auto &tl = async_region[i];
 
-    auto &tl = async_region0[i];
-
-    auto transport_particles =
-        tl.AddTask(none, TransportParticles, pmb.get(), &integrator, t0);
-
-    auto send = tl.AddTask(transport_particles, &SwarmContainer::Send, sc.get(),
-                           BoundaryCommSubset::all);
+    auto transport = tl.AddTask(none, TransportParticles, pmb.get(), t0, dt);
+    auto reset_comms =
+        tl.AddTask(transport, &SwarmContainer::ResetCommunication, sc.get());
+    auto send =
+        tl.AddTask(reset_comms, &SwarmContainer::Send, sc.get(), BoundaryCommSubset::all);
     auto receive =
         tl.AddTask(send, &SwarmContainer::Receive, sc.get(), BoundaryCommSubset::all);
   }
 
-  TaskRegion &sync_region0 = tc.AddRegion(1);
+  TaskRegion &sync_region = tc.AddRegion(1);
   {
-    auto &tl = sync_region0[0];
-    auto stop_comm = tl.AddTask(none, StopCommunicationMesh, blocks);
+    auto &tl = sync_region[0];
+    auto check_completion = tl.AddTask(none, CountNumSent, blocks, t0 + dt, done);
   }
 
   return tc;
 }
 
+// TODO(BRR) to be replaced by iterative tasklist machinery
+TaskListStatus ParticleDriver::IterativeTransport() const {
+  TaskListStatus status;
+  bool transport_done = false;
+  int n_transport_iter = 0;
+  int n_transport_iter_max = 1000;
+  while (!transport_done) {
+    status = IterativeTransportTaskCollection(&transport_done).Execute();
+
+    n_transport_iter++;
+    PARTHENON_REQUIRE(n_transport_iter < n_transport_iter_max,
+                      "Too many transport iterations!");
+  }
+
+  return status;
+}
+
 TaskCollection ParticleDriver::MakeFinalizationTaskCollection() const {
   TaskCollection tc;
   TaskID none(0);
diff --git a/example/particles/particles.hpp b/example/particles/particles.hpp
index 6d74d59dbbb8..705572675247 100644
--- a/example/particles/particles.hpp
+++ b/example/particles/particles.hpp
@@ -33,7 +33,9 @@ class ParticleDriver : public EvolutionDriver {
   ParticleDriver(ParameterInput *pin, ApplicationInput *app_in, Mesh *pm)
       : EvolutionDriver(pin, app_in, pm), integrator(pin) {}
   TaskCollection MakeParticlesCreationTaskCollection() const;
-  TaskCollection MakeParticlesUpdateTaskCollection() const;
+  TaskCollection MakeParticlesTransportTaskCollection() const;
+  TaskListStatus IterativeTransport() const;
+  TaskCollection IterativeTransportTaskCollection(bool *done) const;
   TaskCollection MakeFinalizationTaskCollection() const;
   TaskListStatus Step();
 
diff --git a/src/interface/state_descriptor.cpp b/src/interface/state_descriptor.cpp
index aa3fa5e6af27..af73259ff6aa 100644
--- a/src/interface/state_descriptor.cpp
+++ b/src/interface/state_descriptor.cpp
@@ -452,6 +452,13 @@ StateDescriptor::CreateResolvedStateDescriptor(Packages_t &packages) {
     field_tracker.CategorizeCollection(name, field_dict, &field_provider);
     swarm_tracker.CategorizeCollection(name, package->AllSwarms(), &swarm_provider);
 
+    if (!package->AllSwarms().empty() && !std::is_same<Real, double>::value) {
+      PARTHENON_WARN(
+          "Swarms always use Real precision, even for ParticleVariables containing "
+          "time data, while Parthenon time variables are fixed to double precision. This "
+          "may cause inaccurate comparisons with cycle beginning and end times.")
+    }
+
     // Add package registered boundary conditions
     for (int i = 0; i < 6; ++i)
       state->UserBoundaryFunctions[i].insert(state->UserBoundaryFunctions[i].end(),
diff --git a/src/interface/swarm.cpp b/src/interface/swarm.cpp
index b6f1fb6fd991..8e715f1cb8db 100644
--- a/src/interface/swarm.cpp
+++ b/src/interface/swarm.cpp
@@ -65,38 +65,30 @@ SwarmDeviceContext Swarm::GetDeviceContext() const {
 
 Swarm::Swarm(const std::string &label, const Metadata &metadata, const int nmax_pool_in)
     : label_(label), m_(metadata), nmax_pool_(nmax_pool_in), mask_("mask", nmax_pool_),
-      marked_for_removal_("mfr", nmax_pool_), block_index_("block_index_", nmax_pool_),
+      marked_for_removal_("mfr", nmax_pool_),
+      empty_indices_("empty_indices_", nmax_pool_),
+      block_index_("block_index_", nmax_pool_),
       neighbor_indices_("neighbor_indices_", 4, 4, 4),
-      new_indices_("new_indices_", nmax_pool_),
-      from_to_indices_("from_to_indices_", nmax_pool_ + 1),
-      recv_neighbor_index_("recv_neighbor_index_", nmax_pool_),
-      recv_buffer_index_("recv_buffer_index_", nmax_pool_),
+      new_indices_("new_indices_", nmax_pool_), scratch_a_("scratch_a_", nmax_pool_),
+      scratch_b_("scratch_b_", nmax_pool_),
       num_particles_to_send_("num_particles_to_send_", NMAX_NEIGHBORS),
+      buffer_counters_("buffer_counters_", NMAX_NEIGHBORS),
+      neighbor_received_particles_("neighbor_received_particles_", NMAX_NEIGHBORS),
       cell_sorted_("cell_sorted_", nmax_pool_), mpiStatus(true) {
   PARTHENON_REQUIRE_THROWS(typeid(Coordinates_t) == typeid(UniformCartesian),
                            "SwarmDeviceContext only supports a uniform Cartesian mesh!");
 
   uid_ = get_uid_(label_);
 
+  // Add default swarm fields
   Add(swarm_position::x::name(), Metadata({Metadata::Real}));
   Add(swarm_position::y::name(), Metadata({Metadata::Real}));
   Add(swarm_position::z::name(), Metadata({Metadata::Real}));
+
+  // Initialize index metadata
   num_active_ = 0;
   max_active_index_ = inactive_max_active_index;
-
-  // TODO(BRR) Do this in a device kernel?
-  auto mask_h = Kokkos::create_mirror_view(HostMemSpace(), mask_);
-  auto marked_for_removal_h =
-      Kokkos::create_mirror_view(HostMemSpace(), marked_for_removal_);
-
-  for (int n = 0; n < nmax_pool_; n++) {
-    mask_h(n) = false;
-    marked_for_removal_h(n) = false;
-    free_indices_.push_back(n);
-  }
-
-  Kokkos::deep_copy(mask_, mask_h);
-  Kokkos::deep_copy(marked_for_removal_, marked_for_removal_h);
+  UpdateEmptyIndices();
 }
 
 void Swarm::Add(const std::vector<std::string> &label_array, const Metadata &metadata) {
@@ -196,25 +188,21 @@ void Swarm::Remove(const std::string &label) {
   }
 }
 
-void Swarm::setPoolMax(const std::int64_t nmax_pool) {
+void Swarm::SetPoolMax(const std::int64_t nmax_pool) {
   PARTHENON_REQUIRE(nmax_pool > nmax_pool_, "Must request larger pool size!");
-  std::int64_t n_new_begin = nmax_pool_;
   std::int64_t n_new = nmax_pool - nmax_pool_;
 
   auto pmb = GetBlockPointer();
   auto pm = pmb->pmy_mesh;
 
-  for (std::int64_t n = 0; n < n_new; n++) {
-    free_indices_.push_back(n + n_new_begin);
-  }
-
   // Rely on Kokkos setting the newly added values to false for these arrays
   Kokkos::resize(mask_, nmax_pool);
   Kokkos::resize(marked_for_removal_, nmax_pool);
+  Kokkos::resize(empty_indices_, nmax_pool);
   Kokkos::resize(new_indices_, nmax_pool);
-  Kokkos::resize(from_to_indices_, nmax_pool + 1);
-  Kokkos::resize(recv_neighbor_index_, nmax_pool);
-  Kokkos::resize(recv_buffer_index_, nmax_pool);
+  Kokkos::resize(scratch_a_, nmax_pool);
+  Kokkos::resize(scratch_b_, nmax_pool);
+
   pmb->LogMemUsage(2 * n_new * sizeof(bool));
 
   Kokkos::resize(cell_sorted_, nmax_pool);
@@ -240,7 +228,10 @@ void Swarm::setPoolMax(const std::int64_t nmax_pool) {
 
   nmax_pool_ = nmax_pool;
 
-  // Eliminate any cached SwarmPacks, as they will need to be rebuilt following setPoolMax
+  // Populate new empty indices
+  UpdateEmptyIndices();
+
+  // Eliminate any cached SwarmPacks, as they will need to be rebuilt following SetPoolMax
   pmb->meshblock_data.Get()->ClearSwarmCaches();
   pm->mesh_data.Get("base")->ClearSwarmCaches();
   for (auto &partition : pm->GetDefaultBlockPartitions()) {
@@ -251,129 +242,153 @@ void Swarm::setPoolMax(const std::int64_t nmax_pool) {
 NewParticlesContext Swarm::AddEmptyParticles(const int num_to_add) {
   PARTHENON_DEBUG_REQUIRE(num_to_add >= 0, "Cannot add negative numbers of particles!");
 
+  auto pmb = GetBlockPointer();
+
   if (num_to_add > 0) {
-    while (free_indices_.size() < num_to_add) {
-      increasePoolMax();
+    while (nmax_pool_ - num_active_ < num_to_add) {
+      IncreasePoolMax();
     }
 
-    // TODO(BRR) Use par_scan on device rather than do this on host
-    auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_);
-
-    auto block_index_h = block_index_.GetHostMirrorAndCopy();
+    auto &new_indices = new_indices_;
+    auto &empty_indices = empty_indices_;
+    auto &mask = mask_;
 
-    auto free_index = free_indices_.begin();
+    int max_new_active_index = 0;
+    parthenon::par_reduce(
+        PARTHENON_AUTO_LABEL, 0, num_to_add - 1,
+        KOKKOS_LAMBDA(const int n, int &max_ind) {
+          new_indices(n) = empty_indices(n);
+          mask(new_indices(n)) = true;
 
-    auto new_indices_h = new_indices_.GetHostMirror();
-
-    // Don't bother sanitizing the memory
-    for (int n = 0; n < num_to_add; n++) {
-      mask_h(*free_index) = true;
-      block_index_h(*free_index) = this_block_;
-      max_active_index_ = std::max<int>(max_active_index_, *free_index);
-      new_indices_h(n) = *free_index;
-
-      free_index = free_indices_.erase(free_index);
-    }
+          // Record vote for max active index
+          max_ind = new_indices(n);
+        },
+        Kokkos::Max<int>(max_new_active_index));
 
-    new_indices_.DeepCopy(new_indices_h);
+    // Update max active index if necessary
+    max_active_index_ = std::max(max_active_index_, max_new_active_index);
 
+    new_indices_max_idx_ = num_to_add - 1;
     num_active_ += num_to_add;
 
-    Kokkos::deep_copy(mask_, mask_h);
-    block_index_.DeepCopy(block_index_h);
-    new_indices_max_idx_ = num_to_add - 1;
+    UpdateEmptyIndices();
   } else {
     new_indices_max_idx_ = -1;
   }
 
+  // Create and return NewParticlesContext
   return NewParticlesContext(new_indices_max_idx_, new_indices_);
 }
 
+// Updates the empty_indices_ array so the first N elements contain an ascending list of
+// indices into empty elements of the swarm pool, where N is the number of empty indices
+void Swarm::UpdateEmptyIndices() {
+  auto &mask = mask_;
+  auto &empty_indices = empty_indices_;
+
+  // Associate scratch memory
+  auto &empty_indices_scan = scratch_a_;
+
+  // Calculate prefix sum of empty indices
+  parthenon::par_scan(
+      "Set empty indices prefix sum", 0, nmax_pool_ - 1,
+      KOKKOS_LAMBDA(const int n, int &update, const bool &final) {
+        const int val = !mask(n);
+        if (val) {
+          update += 1;
+        }
+
+        if (final) {
+          empty_indices_scan(n) = update;
+        }
+      });
+
+  // Update list of empty indices such that it is contiguous and in ascending order
+  parthenon::par_for(
+      PARTHENON_AUTO_LABEL, 0, nmax_pool_ - 1, KOKKOS_LAMBDA(const int n) {
+        if (!mask(n)) {
+          empty_indices(empty_indices_scan(n) - 1) = n;
+        }
+      });
+}
+
 // No active particles: nmax_active_index = inactive_max_active_index (= -1)
 // No particles removed: nmax_active_index unchanged
 // Particles removed: nmax_active_index is new max active index
 void Swarm::RemoveMarkedParticles() {
-  // TODO(BRR) Use par_scan to do this on device rather than host
-  auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_);
-  auto marked_for_removal_h =
-      Kokkos::create_mirror_view_and_copy(HostMemSpace(), marked_for_removal_);
-
-  // loop backwards to keep free_indices_ updated correctly
-  for (int n = max_active_index_; n >= 0; n--) {
-    if (mask_h(n)) {
-      if (marked_for_removal_h(n)) {
-        mask_h(n) = false;
-        free_indices_.push_front(n);
-        num_active_ -= 1;
-        if (n == max_active_index_) {
-          max_active_index_ -= 1;
+  int &max_active_index = max_active_index_;
+
+  auto &mask = mask_;
+  auto &marked_for_removal = marked_for_removal_;
+
+  // Update mask, count number of removed particles
+  int num_removed = 0;
+  parthenon::par_reduce(
+      PARTHENON_AUTO_LABEL, 0, max_active_index,
+      KOKKOS_LAMBDA(const int n, int &removed) {
+        if (mask(n)) {
+          if (marked_for_removal(n)) {
+            mask(n) = false;
+            marked_for_removal(n) = false;
+            removed += 1;
+          }
         }
-        marked_for_removal_h(n) = false;
-      }
-    }
-  }
+      },
+      Kokkos::Sum<int>(num_removed));
 
-  Kokkos::deep_copy(mask_, mask_h);
-  Kokkos::deep_copy(marked_for_removal_, marked_for_removal_h);
+  num_active_ -= num_removed;
+
+  UpdateEmptyIndices();
 }
 
 void Swarm::Defrag() {
   if (GetNumActive() == 0) {
     return;
   }
-  // TODO(BRR) Could this algorithm be more efficient? Does it matter?
-  // Add 1 to convert max index to max number
-  std::int64_t num_free = (max_active_index_ + 1) - num_active_;
-  auto pmb = GetBlockPointer();
 
-  auto from_to_indices_h = from_to_indices_.GetHostMirror();
+  // Associate scratch memory
+  auto &scan_scratch_toread = scratch_a_;
+  auto &map = scratch_b_;
 
-  auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_);
-
-  for (int n = 0; n <= max_active_index_; n++) {
-    from_to_indices_h(n) = unset_index_;
-  }
-
-  std::list<int> new_free_indices;
-
-  free_indices_.sort();
-
-  int index = max_active_index_;
-  int num_to_move = std::min<int>(num_free, num_active_);
-  for (int n = 0; n < num_to_move; n++) {
-    while (mask_h(index) == false) {
-      index--;
-    }
-    int index_to_move_from = index;
-    index--;
-
-    // Below this number "moved" particles should actually stay in place
-    if (index_to_move_from < num_active_) {
-      break;
-    }
-    int index_to_move_to = free_indices_.front();
-    free_indices_.pop_front();
-    new_free_indices.push_back(index_to_move_from);
-    from_to_indices_h(index_to_move_from) = index_to_move_to;
-  }
+  auto &mask = mask_;
 
-  // TODO(BRR) Not all these sorts may be necessary
-  new_free_indices.sort();
-  free_indices_.merge(new_free_indices);
+  const int &num_active = num_active_;
+  parthenon::par_scan(
+      "Set empty indices prefix sum", 0, nmax_pool_ - num_active_ - 1,
+      KOKKOS_LAMBDA(const int nn, int &update, const bool &final) {
+        const int n = nn + num_active;
+        const int val = mask(n);
+        if (val) {
+          update += 1;
+        }
+        if (final) scan_scratch_toread(n) = update;
+      });
 
-  from_to_indices_.DeepCopy(from_to_indices_h);
+  parthenon::par_for(
+      PARTHENON_AUTO_LABEL, 0, nmax_pool_ - 1, KOKKOS_LAMBDA(const int n) {
+        if (n >= num_active) {
+          if (mask(n)) {
+            map(scan_scratch_toread(n) - 1) = n;
+          }
+          mask(n) = false;
+        }
+      });
 
-  auto from_to_indices = from_to_indices_;
+  // Reuse scratch memory
+  auto &scan_scratch_towrite = scan_scratch_toread;
 
-  auto &mask = mask_;
-  pmb->par_for(
-      PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) {
-        if (from_to_indices(n) >= 0) {
-          mask(from_to_indices(n)) = mask(n);
-          mask(n) = false;
+  // Update list of empty indices
+  parthenon::par_scan(
+      "Set empty indices prefix sum", 0, num_active_ - 1,
+      KOKKOS_LAMBDA(const int n, int &update, const bool &final) {
+        const int val = !mask(n);
+        if (val) {
+          update += 1;
         }
+        if (final) scan_scratch_towrite(n) = update;
       });
 
+  // Get all dynamical variables in swarm
   auto &int_vector = std::get<getType<int>()>(vectors_);
   auto &real_vector = std::get<getType<Real>()>(vectors_);
   PackIndexMap real_imap;
@@ -387,15 +402,19 @@ void Swarm::Defrag() {
   const int realPackDim = vreal.GetDim(2);
   const int intPackDim = vint.GetDim(2);
 
-  pmb->par_for(
-      PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) {
-        if (from_to_indices(n) >= 0) {
+  // Loop over only the active number of particles, and if mask is empty, copy in particle
+  // using address from prefix sum
+  parthenon::par_for(
+      PARTHENON_AUTO_LABEL, 0, num_active_ - 1, KOKKOS_LAMBDA(const int n) {
+        if (!mask(n)) {
+          const int nread = map(scan_scratch_towrite(n) - 1);
           for (int vidx = 0; vidx < realPackDim; vidx++) {
-            vreal(vidx, from_to_indices(n)) = vreal(vidx, n);
+            vreal(vidx, n) = vreal(vidx, nread);
           }
           for (int vidx = 0; vidx < intPackDim; vidx++) {
-            vint(vidx, from_to_indices(n)) = vint(vidx, n);
+            vint(vidx, n) = vint(vidx, nread);
           }
+          mask(n) = true;
         }
       });
 
diff --git a/src/interface/swarm.hpp b/src/interface/swarm.hpp
index a3e366789986..8eb31a3d1cd7 100644
--- a/src/interface/swarm.hpp
+++ b/src/interface/swarm.hpp
@@ -161,10 +161,10 @@ class Swarm {
   std::string info() const { return info_; }
 
   /// Expand pool size geometrically as necessary
-  void increasePoolMax() { setPoolMax(2 * nmax_pool_); }
+  void IncreasePoolMax() { SetPoolMax(2 * nmax_pool_); }
 
   /// Set max pool size
-  void setPoolMax(const std::int64_t nmax_pool);
+  void SetPoolMax(const std::int64_t nmax_pool);
 
   /// Check whether metadata bit is set
   bool IsSet(const MetadataFlag bit) const { return m_.IsSet(bit); }
@@ -182,6 +182,9 @@ class Swarm {
   /// indicates gaps in the list.
   Real GetPackingEfficiency() const { return num_active_ / (max_active_index_ + 1); }
 
+  // Update sorted array of empty indices in the current memory pool
+  void UpdateEmptyIndices();
+
   /// Remove particles marked for removal and update internal indexing
   void RemoveMarkedParticles();
 
@@ -227,14 +230,10 @@ class Swarm {
   SwarmVariablePack<T> PackVariables(const std::vector<std::string> &name,
                                      PackIndexMap &vmap);
 
-  // Temporarily public
-  int num_particles_sent_;
-  bool finished_transport;
-
-  void LoadBuffers_(const int max_indices_size);
+  void LoadBuffers_();
   void UnloadBuffers_();
 
-  int CountParticlesToSend_(); // Must be public for launching kernel
+  void CountParticlesToSend_(); // Must be public for launching kernel
 
   template <typename T>
   const auto &GetVariableVector() const {
@@ -251,10 +250,6 @@ class Swarm {
 
   void SetNeighborIndices_();
 
-  void CountReceivedParticles_();
-  void UpdateNeighborBufferReceiveIndices_(ParArray1D<int> &neighbor_index,
-                                           ParArray1D<int> &buffer_index);
-
   template <class T>
   SwarmVariablePack<T> PackAllVariables_(PackIndexMap &vmap);
 
@@ -273,28 +268,25 @@ class Swarm {
 
   std::tuple<MapToParticle<int>, MapToParticle<Real>> maps_;
 
-  std::list<int> free_indices_;
   ParArray1D<bool> mask_;
   ParArray1D<bool> marked_for_removal_;
+  ParArray1D<int> empty_indices_; // Indices of empty slots in particle pool
   ParArrayND<int> block_index_; // Neighbor index for each particle. -1 for current block.
   ParArrayND<int> neighbor_indices_; // Indexing of vbvar's neighbor array. -1 for same.
                                      // k,j indices unused in 3D&2D, 2D, respectively
-  ParArray1D<int> new_indices_;     // Persistent array that provides the new indices when
-                                    // AddEmptyParticles is called. Always defragmented.
-  int new_indices_max_idx_;         // Maximum valid index of new_indices_ array.
-  ParArray1D<int> from_to_indices_; // Array used for sorting particles during defragment
-                                    // step (size nmax_pool + 1).
-  ParArray1D<int> recv_neighbor_index_; // Neighbor indices for received particles
-  ParArray1D<int> recv_buffer_index_;   // Buffer indices for received particles
+  ParArray1D<int> new_indices_; // Persistent array that provides the new indices when
+                                // AddEmptyParticles is called. Always defragmented.
+  int new_indices_max_idx_;     // Maximum valid index of new_indices_ array.
+  ParArray1D<int> scratch_a_;   // Scratch memory for index sorting
+  ParArray1D<int> scratch_b_;   // Scratch memory for index sorting
 
   constexpr static int no_block_ = -2;
   constexpr static int this_block_ = -1;
   constexpr static int unset_index_ = -1;
 
   ParArray1D<int> num_particles_to_send_;
-  ParArrayND<int> particle_indices_to_send_;
-
-  std::vector<int> neighbor_received_particles_;
+  ParArray1D<int> buffer_counters_;
+  ParArray1D<int> neighbor_received_particles_;
   int total_received_particles_;
 
   ParArrayND<int> neighbor_buffer_index_; // Map from neighbor index to neighbor bufid
diff --git a/src/interface/swarm_comms.cpp b/src/interface/swarm_comms.cpp
index 07be992ac727..650bcef09304 100644
--- a/src/interface/swarm_comms.cpp
+++ b/src/interface/swarm_comms.cpp
@@ -130,36 +130,6 @@ void Swarm::SetNeighborIndices_() {
     }
   }
 
-  // Draft alternative approach due to LFR utilizing mesh offset comparisons at the
-  // highest refinement level
-  // auto ll_block = pmb->loc.GetDaughter(0, 0, 0);
-  // int finest_level = pmb->loc.level() + 1;
-  // for (auto &n : pmb->neighbors) {
-  //  std::vector<LogicalLocation> dlocs;
-  //  auto &nloc =
-  //      n.loc; // Would need to use the location in the coordinates of the origin tree
-  //  if (nloc.level() == finest_level) {
-  //    dlocs.emplace_back(nloc);
-  //  } else if (nloc.level() == finest_level) {
-  //    dlocs = nloc.GetDaughters(ndim);
-  //  } else if (nloc.level() == finest_level - 2) {
-  //    auto tlocs = nloc.GetDaughters(ndim);
-  //    for (auto &t : tlocs) {
-  //      auto gdlocs = t.GetDaughters(ndim);
-  //      dlocs.insert(dlocs.end(), gdlocs.begin(), gdlocs.end());
-  //    }
-  //  } else {
-  //    PARTHENON_FAIL("Proper nesting is not being respected.");
-  //  }
-  //  for (auto &d : dlocs) {
-  //    const int k = d.lx3() - ll_block.lx3() + 1;
-  //    const int j = d.lx2() - ll_block.lx2() + 1;
-  //    const int i = d.lx1() - ll_block.lx1() + 1;
-  //    if (i >= 0 && i <= 3 && j >= 0 && j <= 3 && k >= 0 && k <= 3)
-  //      neighbor_indices_h(k, j, i) = n.gid;
-  //  }
-  //}
-
   neighbor_indices_.DeepCopy(neighbor_indices_h);
 }
 
@@ -174,8 +144,6 @@ void Swarm::SetupPersistentMPI() {
   // Build up convenience array of neighbor indices
   SetNeighborIndices_();
 
-  neighbor_received_particles_.resize(nbmax);
-
   // Build device array mapping neighbor index to neighbor bufid
   if (pmb->neighbors.size() > 0) {
     ParArrayND<int> neighbor_buffer_index("Neighbor buffer index", pmb->neighbors.size());
@@ -188,7 +156,7 @@ void Swarm::SetupPersistentMPI() {
   }
 }
 
-int Swarm::CountParticlesToSend_() {
+void Swarm::CountParticlesToSend_() {
   auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_);
   auto swarm_d = GetDeviceContext();
   auto pmb = GetBlockPointer();
@@ -197,10 +165,6 @@ int Swarm::CountParticlesToSend_() {
   // Fence to make sure particles aren't currently being transported locally
   // TODO(BRR) do this operation on device.
   pmb->exec_space.fence();
-  auto num_particles_to_send_h = num_particles_to_send_.GetHostMirror();
-  for (int n = 0; n < pmb->neighbors.size(); n++) {
-    num_particles_to_send_h(n) = 0;
-  }
   const int particle_size = GetParticleDataSize();
   vbswarm->particle_size = particle_size;
 
@@ -218,63 +182,31 @@ int Swarm::CountParticlesToSend_() {
         }
       });
 
-  int max_indices_size = 0;
-  int total_noblock_particles = 0;
-  auto block_index_h = block_index_.GetHostMirrorAndCopy();
-  for (int n = 0; n <= max_active_index_; n++) {
-    if (mask_h(n)) {
-      // This particle should be sent
-      if (block_index_h(n) >= 0) {
-        num_particles_to_send_h(block_index_h(n))++;
-        if (max_indices_size < num_particles_to_send_h(block_index_h(n))) {
-          max_indices_size = num_particles_to_send_h(block_index_h(n));
-        }
-      }
-      if (block_index_h(n) == no_block_) {
-        total_noblock_particles++;
-      }
-    }
-  }
-  // Size-0 arrays not permitted but we don't want to short-circuit subsequent logic
-  // that indicates completed communications
-  max_indices_size = std::max<int>(1, max_indices_size);
-
-  // Not a ragged-right array, just for convenience
-  if (total_noblock_particles > 0) {
-    auto noblock_indices =
-        ParArray1D<int>("Particles with no block", total_noblock_particles);
-    auto noblock_indices_h = noblock_indices.GetHostMirror();
-    int counter = 0;
-    for (int n = 0; n <= max_active_index_; n++) {
-      if (mask_h(n)) {
-        if (block_index_h(n) == no_block_) {
-          noblock_indices_h(counter) = n;
-          counter++;
+  // Facilitate lambda captures
+  auto &block_index = block_index_;
+  auto &num_particles_to_send = num_particles_to_send_;
+
+  // Zero out number of particles to send before accumulating
+  pmb->par_for(
+      PARTHENON_AUTO_LABEL, 0, NMAX_NEIGHBORS - 1,
+      KOKKOS_LAMBDA(const int n) { num_particles_to_send[n] = 0; });
+
+  parthenon::par_for(
+      PARTHENON_AUTO_LABEL, 0, max_active_index, KOKKOS_LAMBDA(const int n) {
+        if (swarm_d.IsActive(n)) {
+          bool on_current_mesh_block = true;
+          swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block);
+
+          if (block_index(n) >= 0) {
+            Kokkos::atomic_add(&num_particles_to_send(block_index(n)), 1);
+          }
         }
-      }
-    }
-    noblock_indices.DeepCopy(noblock_indices_h);
-  }
+      });
 
-  // TODO(BRR) don't allocate dynamically
-  particle_indices_to_send_ =
-      ParArrayND<int>("Particle indices to send", nbmax, max_indices_size);
-  auto particle_indices_to_send_h = particle_indices_to_send_.GetHostMirror();
-  std::vector<int> counter(nbmax, 0);
-  for (int n = 0; n <= max_active_index_; n++) {
-    if (mask_h(n)) {
-      if (block_index_h(n) >= 0) {
-        particle_indices_to_send_h(block_index_h(n), counter[block_index_h(n)]) = n;
-        counter[block_index_h(n)]++;
-      }
-    }
-  }
-  num_particles_to_send_.DeepCopy(num_particles_to_send_h);
-  particle_indices_to_send_.DeepCopy(particle_indices_to_send_h);
+  auto num_particles_to_send_h = num_particles_to_send_.GetHostMirrorAndCopy();
 
-  num_particles_sent_ = 0;
+  // Resize send buffers if too small
   for (int n = 0; n < pmb->neighbors.size(); n++) {
-    // Resize buffer if too small
     const int bufid = pmb->neighbors[n].bufid;
     auto sendbuf = vbswarm->bd_var_.send[bufid];
     if (sendbuf.extent(0) < num_particles_to_send_h(n) * particle_size) {
@@ -282,13 +214,10 @@ int Swarm::CountParticlesToSend_() {
       vbswarm->bd_var_.send[bufid] = sendbuf;
     }
     vbswarm->send_size[bufid] = num_particles_to_send_h(n) * particle_size;
-    num_particles_sent_ += num_particles_to_send_h(n);
   }
-
-  return max_indices_size;
 }
 
-void Swarm::LoadBuffers_(const int max_indices_size) {
+void Swarm::LoadBuffers_() {
   auto swarm_d = GetDeviceContext();
   auto pmb = GetBlockPointer();
   const int particle_size = GetParticleDataSize();
@@ -303,34 +232,45 @@ void Swarm::LoadBuffers_(const int max_indices_size) {
   const int realPackDim = vreal.GetDim(2);
   const int intPackDim = vint.GetDim(2);
 
-  // Pack index:
-  // [variable start] [swarm idx]
+  auto &x = Get<Real>(swarm_position::x::name()).Get();
+  auto &y = Get<Real>(swarm_position::y::name()).Get();
+  auto &z = Get<Real>(swarm_position::z::name()).Get();
+
+  // Zero buffer index counters
+  auto &buffer_counters = buffer_counters_;
+  pmb->par_for(
+      PARTHENON_AUTO_LABEL, 0, NMAX_NEIGHBORS - 1,
+      KOKKOS_LAMBDA(const int n) { buffer_counters[n] = 0; });
 
   auto &bdvar = vbswarm->bd_var_;
-  auto num_particles_to_send = num_particles_to_send_;
-  auto particle_indices_to_send = particle_indices_to_send_;
   auto neighbor_buffer_index = neighbor_buffer_index_;
+  // Loop over active particles and use atomic operations to find indices into buffers if
+  // this particle will be sent.
   pmb->par_for(
-      PARTHENON_AUTO_LABEL, 0, max_indices_size - 1,
-      KOKKOS_LAMBDA(const int n) {            // Max index
-        for (int m = 0; m < nneighbor; m++) { // Number of neighbors
+      PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) {
+        if (swarm_d.IsActive(n)) {
+          bool on_current_mesh_block = true;
+          const int m =
+              swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block);
           const int bufid = neighbor_buffer_index(m);
-          if (n < num_particles_to_send(m)) {
-            const int sidx = particle_indices_to_send(m, n);
-            int buffer_index = n * particle_size;
-            swarm_d.MarkParticleForRemoval(sidx);
+
+          if (m >= 0) {
+            const int bid = Kokkos::atomic_fetch_add(&buffer_counters(m), 1);
+            int buffer_index = bid * particle_size;
+            swarm_d.MarkParticleForRemoval(n);
             for (int i = 0; i < realPackDim; i++) {
-              bdvar.send[bufid](buffer_index) = vreal(i, sidx);
+              bdvar.send[bufid](buffer_index) = vreal(i, n);
               buffer_index++;
             }
             for (int i = 0; i < intPackDim; i++) {
-              bdvar.send[bufid](buffer_index) = static_cast<Real>(vint(i, sidx));
+              bdvar.send[bufid](buffer_index) = static_cast<Real>(vint(i, n));
               buffer_index++;
             }
           }
         }
       });
 
+  // Remove particles that were loaded to send to another block from this block
   RemoveMarkedParticles();
 }
 
@@ -339,99 +279,42 @@ void Swarm::Send(BoundaryCommSubset phase) {
   const int nneighbor = pmb->neighbors.size();
   auto swarm_d = GetDeviceContext();
 
-  if (nneighbor == 0) {
-    // TODO(BRR) Do we ever reach this branch?
-    // Process physical boundary conditions on "sent" particles
-    auto block_index_h = block_index_.GetHostMirrorAndCopy();
-    auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_);
-
-    int total_sent_particles = 0;
-    pmb->par_reduce(
-        PARTHENON_AUTO_LABEL, 0, max_active_index_,
-        KOKKOS_LAMBDA(int n, int &total_sent_particles) {
-          if (swarm_d.IsActive(n)) {
-            if (!swarm_d.IsOnCurrentMeshBlock(n)) {
-              total_sent_particles++;
-            }
-          }
-        },
-        Kokkos::Sum<int>(total_sent_particles));
-
-    if (total_sent_particles > 0) {
-      ParArray1D<int> new_indices("new indices", total_sent_particles);
-      auto new_indices_h = new_indices.GetHostMirrorAndCopy();
-      int sent_particle_index = 0;
-      for (int n = 0; n <= max_active_index_; n++) {
-        if (mask_h(n)) {
-          if (block_index_h(n) >= 0 || block_index_h(n) == no_block_) {
-            new_indices_h(sent_particle_index) = n;
-            sent_particle_index++;
-          }
-        }
-      }
-      new_indices.DeepCopy(new_indices_h);
-    }
-  } else {
-    // Query particles for those to be sent
-    int max_indices_size = CountParticlesToSend_();
+  // Query particles for those to be sent
+  CountParticlesToSend_();
 
-    // Prepare buffers for send operations
-    LoadBuffers_(max_indices_size);
+  // Prepare buffers for send operations
+  LoadBuffers_();
 
-    // Send buffer data
-    vbswarm->Send(phase);
-  }
+  // Send buffer data
+  vbswarm->Send(phase);
 }
 
-void Swarm::CountReceivedParticles_() {
+void Swarm::UnloadBuffers_() {
   auto pmb = GetBlockPointer();
+
+  // Count received particles
   total_received_particles_ = 0;
+  auto &neighbor_received_particles = neighbor_received_particles_;
+  auto neighbor_received_particles_h = neighbor_received_particles.GetHostMirror();
   for (int n = 0; n < pmb->neighbors.size(); n++) {
     const int bufid = pmb->neighbors[n].bufid;
     if (vbswarm->bd_var_.flag[bufid] == BoundaryStatus::arrived) {
       PARTHENON_DEBUG_REQUIRE(vbswarm->recv_size[bufid] % vbswarm->particle_size == 0,
                               "Receive buffer is not divisible by particle size!");
-      neighbor_received_particles_[n] =
+      neighbor_received_particles_h(n) =
           vbswarm->recv_size[bufid] / vbswarm->particle_size;
-      total_received_particles_ += neighbor_received_particles_[n];
+      total_received_particles_ += neighbor_received_particles_h(n);
     } else {
-      neighbor_received_particles_[n] = 0;
+      neighbor_received_particles_h(n) = 0;
     }
   }
-}
-
-void Swarm::UpdateNeighborBufferReceiveIndices_(ParArray1D<int> &neighbor_index,
-                                                ParArray1D<int> &buffer_index) {
-  auto pmb = GetBlockPointer();
-  auto neighbor_index_h = neighbor_index.GetHostMirror();
-  auto buffer_index_h =
-      buffer_index.GetHostMirror(); // Index of each particle in its received buffer
-
-  int id = 0;
-  for (int n = 0; n < pmb->neighbors.size(); n++) {
-    for (int m = 0; m < neighbor_received_particles_[n]; m++) {
-      neighbor_index_h(id) = n;
-      buffer_index_h(id) = m;
-      id++;
-    }
-  }
-  neighbor_index.DeepCopy(neighbor_index_h);
-  buffer_index.DeepCopy(buffer_index_h);
-}
-
-void Swarm::UnloadBuffers_() {
-  auto pmb = GetBlockPointer();
-
-  CountReceivedParticles_();
 
   auto &bdvar = vbswarm->bd_var_;
+  const int nbmax = vbswarm->bd_var_.nbmax;
 
   if (total_received_particles_ > 0) {
     auto newParticlesContext = AddEmptyParticles(total_received_particles_);
 
-    auto &recv_neighbor_index = recv_neighbor_index_;
-    auto &recv_buffer_index = recv_buffer_index_;
-    UpdateNeighborBufferReceiveIndices_(recv_neighbor_index, recv_buffer_index);
     auto neighbor_buffer_index = neighbor_buffer_index_;
 
     auto &int_vector = std::get<getType<int>()>(vectors_);
@@ -443,18 +326,37 @@ void Swarm::UnloadBuffers_() {
     int realPackDim = vreal.GetDim(2);
     int intPackDim = vint.GetDim(2);
 
-    // construct map from buffer index to swarm index (or just return vector of
-    // indices!)
     const int particle_size = GetParticleDataSize();
     auto swarm_d = GetDeviceContext();
 
+    // Change meaning of neighbor_received_particles from particles per neighbor to
+    // cumulative particles per neighbor
+    int val_prev = 0;
+    for (int n = 0; n < nbmax; n++) {
+      int val_curr = neighbor_received_particles_h(n);
+      neighbor_received_particles_h(n) += val_prev;
+      val_prev += val_curr;
+    }
+    neighbor_received_particles.DeepCopy(neighbor_received_particles_h);
+
+    auto &x = Get<Real>(swarm_position::x::name()).Get();
+    auto &y = Get<Real>(swarm_position::y::name()).Get();
+    auto &z = Get<Real>(swarm_position::z::name()).Get();
+
     pmb->par_for(
         PARTHENON_AUTO_LABEL, 0, newParticlesContext.GetNewParticlesMaxIndex(),
         // n is both new particle index and index over buffer values
         KOKKOS_LAMBDA(const int n) {
           const int sid = newParticlesContext.GetNewParticleIndex(n);
-          const int nid = recv_neighbor_index(n);
-          int bid = recv_buffer_index(n) * particle_size;
+          // Search for neighbor id over cumulative indices
+          int nid = 0;
+          while (n >= neighbor_received_particles(nid) && nid < nbmax - 1) {
+            nid++;
+          }
+
+          // Convert neighbor id to buffer id
+          int bid = nid == 0 ? n * particle_size
+                             : (n - neighbor_received_particles(nid - 1)) * particle_size;
           const int nbid = neighbor_buffer_index(nid);
           for (int i = 0; i < realPackDim; i++) {
             vreal(i, sid) = bdvar.recv[nbid](bid);
diff --git a/src/interface/swarm_container.cpp b/src/interface/swarm_container.cpp
index c449ef084c6a..9fc26a213e03 100644
--- a/src/interface/swarm_container.cpp
+++ b/src/interface/swarm_container.cpp
@@ -36,36 +36,6 @@ void SwarmContainer::Initialize(const std::shared_ptr<StateDescriptor> resolved_
   }
 }
 
-void SwarmContainer::InitializeBoundaries(const std::shared_ptr<MeshBlock> pmb) {
-  if (swarmVector_.empty()) {
-    // No Swarms in this container, so no need to initialize boundaries
-    // This allows default reflecting boundary conditions to be used when no
-    // swarms are present in a parthenon calculation.
-    // NOTE SwarmContainer::Initialize must have already been called.
-    return;
-  }
-
-  std::stringstream msg;
-  auto &bcs = pmb->pmy_mesh->mesh_bcs;
-  // Check that, if we are using user BCs, they are actually enrolled, and unsupported BCs
-  // are not being used
-  for (int iFace = 0; iFace < 6; iFace++) {
-    if (bcs[iFace] == BoundaryFlag::user) {
-      if (pmb->pmy_mesh->forest.GetTreePtr(pmb->loc.tree())->SwarmBndryFnctn[iFace] ==
-          nullptr) {
-        msg << (iFace % 2 == 0 ? "i" : "o") << "x" << iFace / 2 + 1
-            << " user boundary requested but provided function is null!";
-        PARTHENON_FAIL(msg);
-      }
-    } else if (bcs[iFace] != BoundaryFlag::outflow &&
-               bcs[iFace] != BoundaryFlag::periodic) {
-      msg << (iFace % 2 == 0 ? "i" : "o") << "x" << iFace / 2 + 1 << " boundary flag "
-          << static_cast<int>(bcs[iFace]) << " not supported!";
-      PARTHENON_FAIL(msg);
-    }
-  }
-}
-
 void SwarmContainer::Add(const std::vector<std::string> &labelArray,
                          const Metadata &metadata) {
   // generate the vector and call Add
diff --git a/src/interface/swarm_container.hpp b/src/interface/swarm_container.hpp
index 3e2072db3a64..1a1486b2a44c 100644
--- a/src/interface/swarm_container.hpp
+++ b/src/interface/swarm_container.hpp
@@ -90,8 +90,6 @@ class SwarmContainer {
   void Initialize(const std::shared_ptr<StateDescriptor> resolved_packages,
                   const std::shared_ptr<MeshBlock> pmb);
 
-  void InitializeBoundaries(const std::shared_ptr<MeshBlock> pmb);
-
   void Add(std::shared_ptr<Swarm> swarm) {
     swarmVector_.push_back(swarm);
     swarmMap_[swarm->label()] = swarm;
diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp
index d33fd757e353..fe1803882401 100644
--- a/src/mesh/forest/tree.cpp
+++ b/src/mesh/forest/tree.cpp
@@ -442,8 +442,6 @@ void Tree::EnrollBndryFncts(
       break;
     case BoundaryFlag::user:
       if (app_in->swarm_boundary_conditions[f] != nullptr) {
-        // This is checked to be non-null later in Swarm::AllocateBoundaries, in case user
-        // boundaries are requested but no swarms are used.
         SwarmBndryFnctn[f] = app_in->swarm_boundary_conditions[f];
       }
       break;
diff --git a/src/mesh/meshblock.cpp b/src/mesh/meshblock.cpp
index 3f84ce48293b..0bccf068576d 100644
--- a/src/mesh/meshblock.cpp
+++ b/src/mesh/meshblock.cpp
@@ -144,9 +144,6 @@ void MeshBlock::Initialize(int igid, int ilid, LogicalLocation iloc,
   auto &real_container = meshblock_data.Get();
   real_container->Initialize(shared_from_this());
 
-  // Initialize swarm boundary condition flags
-  real_container->GetSwarmData()->InitializeBoundaries(shared_from_this());
-
   // TODO(jdolence): Should these loops be moved to Variable creation
   // TODO(JMM): What variables should be in vars_cc_? They are used
   // for counting load-balance cost. Should it be different than the
diff --git a/src/tasks/tasks.hpp b/src/tasks/tasks.hpp
index 91b9e103d6aa..d432f0a7ea21 100644
--- a/src/tasks/tasks.hpp
+++ b/src/tasks/tasks.hpp
@@ -544,6 +544,8 @@ class TaskRegion {
 
  public:
   TaskRegion() = delete;
+  TaskRegion(const TaskRegion &) = delete; // Prevent copying TaskRegions during AddRegion
+                                           // calls which is a segfault
   explicit TaskRegion(const int num_lists) : task_lists(num_lists) {
     for (int i = 0; i < num_lists; i++)
       task_lists[i].SetID(i);
diff --git a/tst/unit/test_swarm.cpp b/tst/unit/test_swarm.cpp
index f365362ddcc6..bcff16106110 100644
--- a/tst/unit/test_swarm.cpp
+++ b/tst/unit/test_swarm.cpp
@@ -212,8 +212,8 @@ TEST_CASE("Swarm memory management", "[Swarm]") {
 
   // Check that data was moved during defrag
   x_h = swarm->Get<Real>(swarm_position::x::name()).Get().GetHostMirrorAndCopy();
-  REQUIRE(x_h(2) == 1.2);
-  REQUIRE(x_h(4) == 1.1);
+  REQUIRE(x_h(2) == 1.1);
+  REQUIRE(x_h(4) == 1.2);
   i_h = swarm->Get<int>("i").Get().GetHostMirrorAndCopy();
   REQUIRE(i_h(1) == 2);
 

From ec61c9cb102d6a67a4dd03e7d7fa4e10c82c1032 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Fri, 30 Aug 2024 09:42:12 +0200
Subject: [PATCH 40/40] Release 24.08 (#1164)

---
 CHANGELOG.md   | 20 +++++++++++++++++++-
 CMakeLists.txt |  2 +-
 2 files changed, 20 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index af70ab73376b..62b1ea9c5c78 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,8 +2,26 @@
 
 ## Current develop
 
+### Changed (changing behavior/API/variables/...)
+
+
+### Fixed (not changing behavior/API/variables/...)
+
+
+### Infrastructure (changes irrelevant to downstream codes)
+
+
+### Removed (removing behavior/API/varaibles/...)
+
+
+### Incompatibilities (i.e. breaking changes)
+
+
+## Release 24.08
+Date: 2024-08-30
+
 ### Added (new features/APIs/variables/...)
-- [[PR 1151]]((https://github.com/parthenon-hpc-lab/parthenon/pull/1151)) Add time offset `c` to LowStorageIntegrator
+- [[PR 1151]](https://github.com/parthenon-hpc-lab/parthenon/pull/1151) Add time offset `c` to LowStorageIntegrator
 - [[PR 1147]](https://github.com/parthenon-hpc-lab/parthenon/pull/1147) Add `par_reduce_inner` functions
 - [[PR 1159]](https://github.com/parthenon-hpc-lab/parthenon/pull/1159) Add additional timestep controllers in parthenon/time.
 - [[PR 1148]](https://github.com/parthenon-hpc-lab/parthenon/pull/1148) Add `GetPackDimension` to `StateDescriptor` for calculating pack sizes before `Mesh` initialization
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 768130ba1c0c..c398d34eae29 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -20,7 +20,7 @@ cmake_minimum_required(VERSION 3.16)
 # Imports machine-specific configuration
 include(cmake/MachineCfg.cmake)
 
-project(parthenon VERSION 24.03 LANGUAGES C CXX)
+project(parthenon VERSION 24.08 LANGUAGES C CXX)
 
 if (${CMAKE_VERSION} VERSION_GREATER_EQUAL 3.19.0)
   cmake_policy(SET CMP0110 NEW)