diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml
index 4b95511d82f5..9e1cbcfad76f 100644
--- a/.github/workflows/docs.yml
+++ b/.github/workflows/docs.yml
@@ -12,9 +12,12 @@ on:
 jobs:
     docs:
       name: build and deploy docs
-      runs-on: ubuntu-latest
+      # not using latest due to issues with pip user packages, see
+      # https://github.com/actions/runner-images/issues/10781 and
+      # https://github.com/actions/runner-images/issues/10636
+      runs-on: ubuntu-22.04
 
-      steps:      
+      steps:
         - name: Checkout code
           uses: actions/checkout@v2
           with:
@@ -23,9 +26,9 @@ jobs:
           run: export DEBIAN_FRONTEND=noninteractive
         - name: install dependencies
           run: |
-            pip install --break-system-packages sphinx
-            pip install --break-system-packages sphinx-rtd-theme
-            pip install --break-system-packages sphinx-multiversion
+            pip install sphinx
+            pip install sphinx-rtd-theme
+            pip install sphinx-multiversion
         - name: build docs
           run: |
             echo "Repo = ${GITHUB_REPOSITORY}"
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1f046faabe75..c29941247e43 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,12 +11,16 @@
 - [[PR 1161]](https://github.com/parthenon-hpc-lab/parthenon/pull/1161) Make flux field Metadata accessible, add Metadata::CellMemAligned flag, small perfomance upgrades
 
 ### Changed (changing behavior/API/variables/...)
+- [[PR 1206]](https://github.com/parthenon-hpc-lab/parthenon/pull/1206) Leapfrog fix
+- [[PR1203]](https://github.com/parthenon-hpc-lab/parthenon/pull/1203) Pin Ubuntu CI image
+- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag
 - [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize
 - [[Issue 1165]](https://github.com/parthenon-hpc-lab/parthenon/issues/1165) Bump Kokkos submodule to 4.4.1
 - [[PR 1171]](https://github.com/parthenon-hpc-lab/parthenon/pull/1171) Add PARTHENON_USE_SYSTEM_PACKAGES build option
 - [[PR 1172]](https://github.com/parthenon-hpc-lab/parthenon/pull/1172) Make parthenon manager robust against external MPI init and finalize calls
 
 ### Fixed (not changing behavior/API/variables/...)
+- [[PR 1170]](https://github.com/parthenon-hpc-lab/parthenon/pull/1170) Fixed incorrect initialization of array by a const not constexpr
 - [[PR 1189]](https://github.com/parthenon-hpc-lab/parthenon/pull/1189) Address CUDA MPI/ICP issue with Kokkos <=4.4.1
 - [[PR 1178]](https://github.com/parthenon-hpc-lab/parthenon/pull/1178) Fix issue with mesh pointer when using relative residual tolerance in BiCGSTAB solver.
 - [[PR 1173]](https://github.com/parthenon-hpc-lab/parthenon/pull/1173) Make debugging easier by making parthenon throw an error if ParameterInput is different on multiple MPI ranks.
@@ -28,7 +32,7 @@
 
 
 ### Incompatibilities (i.e. breaking changes)
-
+- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag
 
 ## Release 24.08
 Date: 2024-08-30
diff --git a/doc/sphinx/src/boundary_conditions.rst b/doc/sphinx/src/boundary_conditions.rst
index cb0e08532492..9537f0c65728 100644
--- a/doc/sphinx/src/boundary_conditions.rst
+++ b/doc/sphinx/src/boundary_conditions.rst
@@ -1,4 +1,4 @@
-.. _sphinx-doc:
+.. _boundary-conds:
 
 Boundary Conditions
 ===================
@@ -10,7 +10,6 @@ Natively, Parthenon supports three kinds of boundary conditions:
 
 - ``periodic``
 - ``outflow``
-- ``reflecting``
 
 which are all imposed on variables with the ``Metadata::FillGhost``
 metadata flag. To set the boundaries in each direction, set the
@@ -22,8 +21,8 @@ metadata flag. To set the boundaries in each direction, set the
    ix1_bc = outflow
    ox1_bc = outflow
 
-   ix2_bc = reflecting
-   ox2_bc = reflecting
+   ix2_bc = outflow
+   ox2_bc = outflow
 
    ix3_bc = periodic
    ox3_bc = periodic
@@ -40,7 +39,9 @@ for your ``parthenon_manager``. e.g.,
 
 .. code:: c++
 
-   pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x1] = MyBoundaryInnerX1;
+   pman.app_input->RegisterBoundaryCondition(
+	  parthenon::BoundaryFace::inner_x1,
+	  "my_bc_name", MyBoundaryInnerX1);
 
 where ``BoundaryFace`` is an enum defined in ``defs.hpp`` as
 
@@ -58,13 +59,13 @@ where ``BoundaryFace`` is an enum defined in ``defs.hpp`` as
      outer_x3 = 5
    };
 
-You can then set this boundary condition via the ``user`` flag in the
-input file:
+You can then set this boundary condition by using the name you
+registered in the input file:
 
 ::
 
    <parthenon/mesh>
-   ix1_bc = user
+   ix1_bc = my_bc_name
 
 Boundary conditions so defined should look roughly like
 
@@ -100,9 +101,19 @@ Other than these requirements, the ``Boundary`` object can do whatever
 you like. Reference implementations of the standard boundary conditions
 are available `here <https://github.com/parthenon-hpc-lab/parthenon/blob/develop/src/bvals/boundary_conditions.cpp>`__.
 
+.. note::
 
-Per package user-defined boundary conditions
---------------------------------------------
+  A per-variable reflecting boundary condition is available, but you
+  must register it manually. To do so, simply call
+  ``app_in->RegisterDefaultReflectingBoundaryConditions()`` and it
+  will be available as a mesh boundary with the name ``reflecting``.
+  The reason manual registration is required is to support custom
+  reflecting boundary conditions in the case where a single variable
+  is used as the state vector.
+
+
+Per package user-defined boundary conditions.
+---------------------------------
 
 In addition to user defined *global* boundary conditions, Parthenon also supports 
 registration of boundary conditions at the *per package* level. These per package 
@@ -133,4 +144,4 @@ for a more complete example):
     pkg->UserBoundaryFunctions[BF::inner_x1].push_back(GetMyBC<X1DIR, BCSide::Inner>());
     pkg->UserBoundaryFunctions[BF::inner_x2].push_back(GetMyBC<X2DIR, BCSide::Inner>());
     ...
-  }
\ No newline at end of file
+  }
diff --git a/example/advection/main.cpp b/example/advection/main.cpp
index 146093bc22eb..5ed943b07474 100644
--- a/example/advection/main.cpp
+++ b/example/advection/main.cpp
@@ -25,6 +25,7 @@ int main(int argc, char *argv[]) {
   pman.app_input->ProblemGenerator = advection_example::ProblemGenerator;
   pman.app_input->UserWorkAfterLoop = advection_example::UserWorkAfterLoop;
   pman.app_input->UserMeshWorkBeforeOutput = advection_example::UserMeshWorkBeforeOutput;
+  pman.app_input->RegisterDefaultReflectingBoundaryConditions();
 
   // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
   auto manager_status = pman.ParthenonInitEnv(argc, argv);
diff --git a/example/boundary_exchange/boundary_exchange_driver.cpp b/example/boundary_exchange/boundary_exchange_driver.cpp
index bdb43d7e66ca..640705ca77ef 100644
--- a/example/boundary_exchange/boundary_exchange_driver.cpp
+++ b/example/boundary_exchange/boundary_exchange_driver.cpp
@@ -76,23 +76,23 @@ int main(int argc, char *argv[]) {
                      ar3_t{1.0, 1.0, 1.0});
   // forest_def.AddFace(0, {n[0], n[1], n[3], n[2]}, ar3_t{0.0, 0.0, 0.0},
   // ar3_t{1.0, 1.0, 1.0});
-  forest_def.AddBC(edge_t({n[0], n[1]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[0], n[3]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[0], n[1]}));
+  forest_def.AddBC(edge_t({n[0], n[3]}));
 
   forest_def.AddFace(1, {n[1], n[4], n[2], n[5]}, ar3_t{2.0, 0.0, 0.0},
                      ar3_t{3.0, 1.0, 1.0});
-  forest_def.AddBC(edge_t({n[1], n[4]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[4], n[5]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[1], n[4]}));
+  forest_def.AddBC(edge_t({n[4], n[5]}));
 
   forest_def.AddFace(3, {n[3], n[2], n[6], n[7]}, ar3_t{0.0, 2.0, 0.0},
                      ar3_t{1.0, 3.0, 1.0});
-  forest_def.AddBC(edge_t({n[6], n[7]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[3], n[6]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[6], n[7]}));
+  forest_def.AddBC(edge_t({n[3], n[6]}));
 
   forest_def.AddFace(4, {n[2], n[5], n[7], n[8]}, ar3_t{2.0, 2.0, 0.0},
                      ar3_t{3.0, 3.0, 1.0});
-  forest_def.AddBC(edge_t({n[5], n[8]}), parthenon::BoundaryFlag::outflow);
-  forest_def.AddBC(edge_t({n[7], n[8]}), parthenon::BoundaryFlag::outflow);
+  forest_def.AddBC(edge_t({n[5], n[8]}));
+  forest_def.AddBC(edge_t({n[7], n[8]}));
 
   forest_def.AddInitialRefinement(parthenon::LogicalLocation(0, 1, 0, 0, 0));
   pman.ParthenonInitPackagesAndMesh(forest_def);
diff --git a/example/particles/main.cpp b/example/particles/main.cpp
index 45b84fcfe823..e733033608c3 100644
--- a/example/particles/main.cpp
+++ b/example/particles/main.cpp
@@ -49,12 +49,14 @@ int main(int argc, char *argv[]) {
 
   // Redefine parthenon defaults
   pman.app_input->ProcessPackages = particles_example::ProcessPackages;
-  pman.app_input->boundary_conditions[BoundaryFace::inner_x1] =
-      parthenon::BoundaryFunction::OutflowInnerX1;
-  pman.app_input->boundary_conditions[BoundaryFace::outer_x1] =
-      parthenon::BoundaryFunction::OutflowOuterX1;
-  pman.app_input->swarm_boundary_conditions[BoundaryFace::inner_x1] = SwarmUserInnerX1;
-  pman.app_input->swarm_boundary_conditions[BoundaryFace::outer_x1] = SwarmUserOuterX1;
+  pman.app_input->RegisterBoundaryCondition(BoundaryFace::inner_x1,
+                                            parthenon::BoundaryFunction::OutflowInnerX1);
+  pman.app_input->RegisterBoundaryCondition(BoundaryFace::outer_x1,
+                                            parthenon::BoundaryFunction::OutflowOuterX1);
+  pman.app_input->RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1,
+                                                 SwarmUserInnerX1);
+  pman.app_input->RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1,
+                                                 SwarmUserOuterX1);
   // Note that this example does not use a ProblemGenerator
   pman.ParthenonInitPackagesAndMesh();
 
diff --git a/example/particles/particles.cpp b/example/particles/particles.cpp
index a01a121f75fa..0f250e9ab02f 100644
--- a/example/particles/particles.cpp
+++ b/example/particles/particles.cpp
@@ -288,17 +288,18 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
           } while (r > 0.5);
 
           // Randomly sample direction perpendicular to origin
-          Real theta = acos(2. * rng_gen.drand() - 1.);
-          Real phi = 2. * M_PI * rng_gen.drand();
-          v(0, n) = sin(theta) * cos(phi);
-          v(1, n) = sin(theta) * sin(phi);
-          v(2, n) = cos(theta);
+          const Real mu = 2.0 * rng_gen.drand() - 1.0;
+          const Real phi = 2. * M_PI * rng_gen.drand();
+          const Real stheta = std::sqrt(1.0 - mu * mu);
+          v(0, n) = stheta * cos(phi);
+          v(1, n) = stheta * sin(phi);
+          v(2, n) = mu;
           // Project v onto plane normal to sphere
           Real vdN = v(0, n) * x(n) + v(1, n) * y(n) + v(2, n) * z(n);
-          Real NdN = r * r;
-          v(0, n) = v(0, n) - vdN / NdN * x(n);
-          v(1, n) = v(1, n) - vdN / NdN * y(n);
-          v(2, n) = v(2, n) - vdN / NdN * z(n);
+          Real inverse_NdN = 1. / (r * r);
+          v(0, n) = v(0, n) - vdN * inverse_NdN * x(n);
+          v(1, n) = v(1, n) - vdN * inverse_NdN * y(n);
+          v(2, n) = v(2, n) - vdN * inverse_NdN * z(n);
 
           // Normalize
           Real v_tmp = sqrt(v(0, n) * v(0, n) + v(1, n) * v(1, n) + v(2, n) * v(2, n));
@@ -335,11 +336,12 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
           z(n) = minx_k + nx_k * dx_k * rng_gen.drand();
 
           // Randomly sample direction on the unit sphere, fixing speed
-          Real theta = acos(2. * rng_gen.drand() - 1.);
-          Real phi = 2. * M_PI * rng_gen.drand();
-          v(0, n) = vel * sin(theta) * cos(phi);
-          v(1, n) = vel * sin(theta) * sin(phi);
-          v(2, n) = vel * cos(theta);
+          const Real mu = 2.0 * rng_gen.drand() - 1.0;
+          const Real phi = 2. * M_PI * rng_gen.drand();
+          const Real stheta = std::sqrt(1.0 - mu * mu);
+          v(0, n) = vel * stheta * cos(phi);
+          v(1, n) = vel * stheta * sin(phi);
+          v(2, n) = vel * mu;
 
           // Create particles at the beginning of the timestep
           t(n) = t0;
diff --git a/example/sparse_advection/main.cpp b/example/sparse_advection/main.cpp
index 585d1a35ab1e..ce6b197e335a 100644
--- a/example/sparse_advection/main.cpp
+++ b/example/sparse_advection/main.cpp
@@ -27,6 +27,7 @@ int main(int argc, char *argv[]) {
   pman.app_input->ProblemGenerator = sparse_advection_example::ProblemGenerator;
   pman.app_input->PostStepDiagnosticsInLoop =
       sparse_advection_example::PostStepDiagnosticsInLoop;
+  pman.app_input->RegisterDefaultReflectingBoundaryConditions();
 
   // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
   auto manager_status = pman.ParthenonInitEnv(argc, argv);
diff --git a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py
index b7cb8b1d74c4..dd1ff9051d81 100644
--- a/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py
+++ b/scripts/python/packages/parthenon_tools/parthenon_tools/phdf_diff.py
@@ -211,7 +211,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e
         if one:
             return ret_code
 
-    # Compare the names of attributes in /Info, except "Time"
+    # Compare the names of attributes in /Info, except those we know
+    # may vary safely
     f0_Info = {
         key: value
         for key, value in f0.Info.items()
@@ -219,6 +220,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e
         and key != "BlocksPerPE"
         and key != "WallTime"
         and key != "OutputFormatVersion"
+        and key != "BoundaryConditions"
+        and key != "SwarmBoundaryConditions"
     }
     f1_Info = {
         key: value
@@ -227,6 +230,8 @@ def compare_metadata(f0, f1, quiet=False, one=False, check_input=False, tol=1.0e
         and key != "BlocksPerPE"
         and key != "WallTime"
         and key != "OutputFormatVersion"
+        and key != "BoundaryConditions"
+        and key != "SwarmBoundaryConditions"
     }
     if sorted(f0_Info.keys()) != sorted(f1_Info.keys()):
         print("Names of attributes in '/Info' of differ")
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 133279563126..37883ddcb7af 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -95,6 +95,9 @@ set(COORDINATE_TYPE UniformCartesian) # TODO: Make this an option when more are
 configure_file(config.hpp.in generated/config.hpp @ONLY)
 
 add_library(parthenon
+  application_input.cpp
+  application_input.hpp
+
   bvals/comms/bvals_in_one.hpp
   bvals/comms/bvals_utils.hpp
   bvals/comms/build_boundary_buffers.cpp
diff --git a/src/application_input.cpp b/src/application_input.cpp
new file mode 100644
index 000000000000..42f7eb116272
--- /dev/null
+++ b/src/application_input.cpp
@@ -0,0 +1,107 @@
+//========================================================================================
+// (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.
+//========================================================================================
+
+#include <string>
+
+#include "application_input.hpp"
+#include "basic_types.hpp"
+#include "bvals/boundary_conditions.hpp"
+#include "utils/error_checking.hpp"
+
+namespace parthenon {
+ApplicationInput::ApplicationInput() {
+  using namespace BoundaryFunction;
+  const std::string OUTFLOW = "outflow";
+  const std::string PERIODIC = "periodic";
+
+  // Periodic handled separately for mesh
+  RegisterBoundaryCondition(BoundaryFace::inner_x1, OUTFLOW, &OutflowInnerX1);
+  RegisterBoundaryCondition(BoundaryFace::outer_x1, OUTFLOW, &OutflowOuterX1);
+  RegisterBoundaryCondition(BoundaryFace::inner_x2, OUTFLOW, &OutflowInnerX2);
+  RegisterBoundaryCondition(BoundaryFace::outer_x2, OUTFLOW, &OutflowOuterX2);
+  RegisterBoundaryCondition(BoundaryFace::inner_x3, OUTFLOW, &OutflowInnerX3);
+  RegisterBoundaryCondition(BoundaryFace::outer_x3, OUTFLOW, &OutflowOuterX3);
+
+  // Periodic is explicit function for swarms
+  RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1, OUTFLOW, &SwarmOutflowInnerX1);
+  RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1, OUTFLOW, &SwarmOutflowOuterX1);
+  RegisterSwarmBoundaryCondition(BoundaryFace::inner_x2, OUTFLOW, &SwarmOutflowInnerX2);
+  RegisterSwarmBoundaryCondition(BoundaryFace::outer_x2, OUTFLOW, &SwarmOutflowOuterX2);
+  RegisterSwarmBoundaryCondition(BoundaryFace::inner_x3, OUTFLOW, &SwarmOutflowInnerX3);
+  RegisterSwarmBoundaryCondition(BoundaryFace::outer_x3, OUTFLOW, &SwarmOutflowOuterX3);
+  RegisterSwarmBoundaryCondition(BoundaryFace::inner_x1, PERIODIC, &SwarmPeriodicInnerX1);
+  RegisterSwarmBoundaryCondition(BoundaryFace::outer_x1, PERIODIC, &SwarmPeriodicOuterX1);
+  RegisterSwarmBoundaryCondition(BoundaryFace::inner_x2, PERIODIC, &SwarmPeriodicInnerX2);
+  RegisterSwarmBoundaryCondition(BoundaryFace::outer_x2, PERIODIC, &SwarmPeriodicOuterX2);
+  RegisterSwarmBoundaryCondition(BoundaryFace::inner_x3, PERIODIC, &SwarmPeriodicInnerX3);
+  RegisterSwarmBoundaryCondition(BoundaryFace::outer_x3, PERIODIC, &SwarmPeriodicOuterX3);
+}
+
+void ApplicationInput::RegisterDefaultReflectingBoundaryConditions() {
+  using namespace BoundaryFunction;
+  const std::string REFLECTING = "reflecting";
+  RegisterBoundaryCondition(BoundaryFace::inner_x1, REFLECTING, &ReflectInnerX1);
+  RegisterBoundaryCondition(BoundaryFace::outer_x1, REFLECTING, &ReflectOuterX1);
+  RegisterBoundaryCondition(BoundaryFace::inner_x2, REFLECTING, &ReflectInnerX2);
+  RegisterBoundaryCondition(BoundaryFace::outer_x2, REFLECTING, &ReflectOuterX2);
+  RegisterBoundaryCondition(BoundaryFace::inner_x3, REFLECTING, &ReflectInnerX3);
+  RegisterBoundaryCondition(BoundaryFace::outer_x3, REFLECTING, &ReflectOuterX3);
+}
+void ApplicationInput::RegisterBoundaryCondition(BoundaryFace face,
+                                                 const std::string &name,
+                                                 BValFunc condition) {
+  if (boundary_conditions_[face].count(name) > 0) {
+    PARTHENON_THROW("Boundary condition " + name + " at face " + std::to_string(face) +
+                    "already registered.");
+  }
+  boundary_conditions_[face][name] = condition;
+}
+void ApplicationInput::RegisterSwarmBoundaryCondition(BoundaryFace face,
+                                                      const std::string &name,
+                                                      SBValFunc condition) {
+  if (swarm_boundary_conditions_[face].count(name) > 0) {
+    PARTHENON_THROW("Swarm boundary condition " + name + " at face " +
+                    std::to_string(face) + "already registered.");
+  }
+  swarm_boundary_conditions_[face][name] = condition;
+}
+
+BValFunc ApplicationInput::GetBoundaryCondition(BoundaryFace face,
+                                                const std::string &name) const {
+  if (boundary_conditions_[face].count(name) == 0) {
+    std::stringstream msg;
+    msg << "Boundary condition " << name << " at face " << face << "not registered!\n"
+        << "Available conditions for this face are:\n";
+    for (const auto &[name, func] : boundary_conditions_[face]) {
+      msg << name << "\n";
+    }
+    PARTHENON_THROW(msg);
+  }
+  return boundary_conditions_[face].at(name);
+}
+SBValFunc ApplicationInput::GetSwarmBoundaryCondition(BoundaryFace face,
+                                                      const std::string &name) const {
+  if (swarm_boundary_conditions_[face].count(name) == 0) {
+    std::stringstream msg;
+    msg << "Swarm boundary condition " << name << " at face " << face
+        << "not registered!\n"
+        << "Available conditions for this face are:\n";
+    for (const auto &[name, func] : swarm_boundary_conditions_[face]) {
+      msg << name << "\n";
+    }
+    PARTHENON_THROW(msg);
+  }
+  return swarm_boundary_conditions_[face].at(name);
+}
+
+} // namespace parthenon
diff --git a/src/application_input.hpp b/src/application_input.hpp
index 3a8c9bbbee06..9c91585bfb3e 100644
--- a/src/application_input.hpp
+++ b/src/application_input.hpp
@@ -19,17 +19,25 @@
 #include <string>
 #include <vector>
 
+#include "basic_types.hpp"
 #include "bvals/boundary_conditions.hpp"
 #include "defs.hpp"
-#include "interface/state_descriptor.hpp"
-#include "outputs/output_parameters.hpp"
-#include "parameter_input.hpp"
 #include "parthenon_arrays.hpp"
 
 namespace parthenon {
+class Mesh;
+class MeshBlock;
+class MeshBlockApplicationData;
+class OutputParameters;
+class Packages_t;
+class ParameterInput;
+template <typename T>
+class MeshData;
 
-struct ApplicationInput {
+class ApplicationInput {
  public:
+  ApplicationInput();
+
   // ParthenonManager functions
   std::function<Packages_t(std::unique_ptr<ParameterInput> &)> ProcessPackages = nullptr;
 
@@ -57,8 +65,6 @@ struct ApplicationInput {
 
   std::function<void(Mesh *, ParameterInput *, SimTime &)> UserWorkAfterLoop = nullptr;
   std::function<void(Mesh *, ParameterInput *, SimTime &)> UserWorkBeforeLoop = nullptr;
-  BValFunc boundary_conditions[BOUNDARY_NFACES] = {nullptr};
-  SBValFunc swarm_boundary_conditions[BOUNDARY_NFACES] = {nullptr};
 
   // MeshBlock functions
   std::function<std::unique_ptr<MeshBlockApplicationData>(MeshBlock *, ParameterInput *)>
@@ -68,6 +74,28 @@ struct ApplicationInput {
   std::function<void(MeshBlock *, ParameterInput *)> PostInitialization = nullptr;
   std::function<void(MeshBlock *, ParameterInput *, const SimTime &)>
       MeshBlockUserWorkBeforeOutput = nullptr;
+
+  // Physical boundaries
+  void RegisterBoundaryCondition(BoundaryFace face, const std::string &name,
+                                 BValFunc condition);
+  void RegisterBoundaryCondition(BoundaryFace face, BValFunc condition) {
+    RegisterBoundaryCondition(face, "user", condition);
+  }
+  void RegisterDefaultReflectingBoundaryConditions();
+
+  void RegisterSwarmBoundaryCondition(BoundaryFace face, const std::string &name,
+                                      SBValFunc condition);
+  void RegisterSwarmBoundaryCondition(BoundaryFace face, SBValFunc condition) {
+    RegisterSwarmBoundaryCondition(face, "user", condition);
+  }
+
+  // Getters
+  BValFunc GetBoundaryCondition(BoundaryFace face, const std::string &name) const;
+  SBValFunc GetSwarmBoundaryCondition(BoundaryFace face, const std::string &name) const;
+
+ private:
+  Dictionary<BValFunc> boundary_conditions_[BOUNDARY_NFACES];
+  Dictionary<SBValFunc> swarm_boundary_conditions_[BOUNDARY_NFACES];
 };
 
 } // namespace parthenon
diff --git a/src/bvals/boundary_conditions.cpp b/src/bvals/boundary_conditions.cpp
index 8f6777f5e705..648f04ac7478 100644
--- a/src/bvals/boundary_conditions.cpp
+++ b/src/bvals/boundary_conditions.cpp
@@ -207,7 +207,7 @@ bool DoPhysicalBoundary_(const BoundaryFlag flag, const BoundaryFace face,
     return false;
   } // ndim always at least 1
 
-  return true; // reflect, outflow, user, dims correct
+  return true; // user, dims correct
 }
 
 bool DoPhysicalSwarmBoundary_(const BoundaryFlag flag, const BoundaryFace face,
@@ -215,7 +215,7 @@ bool DoPhysicalSwarmBoundary_(const BoundaryFlag flag, const BoundaryFace face,
   if (flag == BoundaryFlag::undef) return false;
   if (flag == BoundaryFlag::block) return false;
 
-  return true; // outflow, periodic, user, dims (particles always 3D) correct
+  return true; // periodic, user, dims (particles always 3D) correct
 }
 
 } // namespace boundary_cond_impl
diff --git a/src/bvals/boundary_conditions.hpp b/src/bvals/boundary_conditions.hpp
index 020c20bcb720..b14ad1955e7f 100644
--- a/src/bvals/boundary_conditions.hpp
+++ b/src/bvals/boundary_conditions.hpp
@@ -14,11 +14,13 @@
 #ifndef BVALS_BOUNDARY_CONDITIONS_HPP_
 #define BVALS_BOUNDARY_CONDITIONS_HPP_
 
+#include <array>
 #include <functional>
 #include <memory>
-#include <string>
+#include <vector>
 
 #include "basic_types.hpp"
+#include "defs.hpp"
 
 namespace parthenon {
 
@@ -32,6 +34,8 @@ class Swarm;
 // Physical boundary conditions
 using BValFunc = std::function<void(std::shared_ptr<MeshBlockData<Real>> &, bool)>;
 using SBValFunc = std::function<void(std::shared_ptr<Swarm> &)>;
+using BValFuncArray_t = std::array<std::vector<BValFunc>, BOUNDARY_NFACES>;
+using SBValFuncArray_t = std::array<std::vector<SBValFunc>, BOUNDARY_NFACES>;
 
 TaskStatus ApplyBoundaryConditionsOnCoarseOrFine(std::shared_ptr<MeshBlockData<Real>> &rc,
                                                  bool coarse);
diff --git a/src/bvals/boundary_flag.cpp b/src/bvals/boundary_flag.cpp
index 649ba7bd7c5b..37dfa229939a 100644
--- a/src/bvals/boundary_flag.cpp
+++ b/src/bvals/boundary_flag.cpp
@@ -3,7 +3,7 @@
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (C) (or copyright) 2020. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -35,54 +35,14 @@ namespace parthenon {
 //  condition. Typically called in Mesh() ctor and in pgen/*.cpp files.
 
 BoundaryFlag GetBoundaryFlag(const std::string &input_string) {
-  if (input_string == "reflecting") {
-    return BoundaryFlag::reflect;
-  } else if (input_string == "outflow") {
-    return BoundaryFlag::outflow;
-  } else if (input_string == "periodic") {
+  if (input_string == "periodic") {
     return BoundaryFlag::periodic;
   } else if (input_string == "none") {
     return BoundaryFlag::undef;
   } else if (input_string == "block") {
     return BoundaryFlag::block;
-  } else if (input_string == "user") {
-    return BoundaryFlag::user;
   } else {
-    std::stringstream msg;
-    msg << "### FATAL ERROR in GetBoundaryFlag" << std::endl
-        << "Input string=" << input_string << "\n"
-        << "is an invalid boundary type" << std::endl;
-    PARTHENON_FAIL(msg);
-  }
-}
-
-//----------------------------------------------------------------------------------------
-//! \fn GetBoundaryString(BoundaryFlag input_flag)
-//  \brief Parses enumerated type BoundaryFlag internal integer representation to return
-//  string describing the boundary condition. Typicall used to format descriptive errors
-//  or diagnostics. Inverse of GetBoundaryFlag().
-
-std::string GetBoundaryString(BoundaryFlag input_flag) {
-  switch (input_flag) {
-  case BoundaryFlag::block: // -1
-    return "block";
-  case BoundaryFlag::undef: // 0
-    return "none";
-  case BoundaryFlag::reflect:
-    return "reflecting";
-  case BoundaryFlag::outflow:
-    return "outflow";
-  case BoundaryFlag::periodic:
-    return "periodic";
-  case BoundaryFlag::user:
-    return "user";
-  default:
-    std::stringstream msg;
-    msg << "### FATAL ERROR in GetBoundaryString" << std::endl
-        << "Input enum class BoundaryFlag=" << static_cast<int>(input_flag) << "\n"
-        << "is an invalid boundary type" << std::endl;
-    PARTHENON_FAIL(msg);
-    break;
+    return BoundaryFlag::user;
   }
 }
 
@@ -96,7 +56,7 @@ std::string GetBoundaryString(BoundaryFlag input_flag) {
 void CheckBoundaryFlag(BoundaryFlag block_flag, CoordinateDirection dir) {
   std::stringstream msg;
   msg << "### FATAL ERROR in CheckBoundaryFlag" << std::endl
-      << "Attempting to set invalid MeshBlock boundary= " << GetBoundaryString(block_flag)
+      << "Attempting to set invalid MeshBlock boundary"
       << "\nin x" << dir << " direction" << std::endl;
   switch (dir) {
   case CoordinateDirection::X1DIR:
diff --git a/src/bvals/bvals.hpp b/src/bvals/bvals.hpp
index 0c93ae858f5b..bc6ed0ebed56 100644
--- a/src/bvals/bvals.hpp
+++ b/src/bvals/bvals.hpp
@@ -50,7 +50,6 @@ struct RegionSize;
 
 // free functions to return boundary flag given input string, and vice versa
 BoundaryFlag GetBoundaryFlag(const std::string &input_string);
-std::string GetBoundaryString(BoundaryFlag input_flag);
 // + confirming that the MeshBlock's boundaries are all valid selections
 void CheckBoundaryFlag(BoundaryFlag block_flag, CoordinateDirection dir);
 
diff --git a/src/defs.hpp b/src/defs.hpp
index d058b14abfc6..cd656d05e037 100644
--- a/src/defs.hpp
+++ b/src/defs.hpp
@@ -3,7 +3,7 @@
 // Copyright(C) 2014 James M. Stone <jmstone@princeton.edu> and other code contributors
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (C) (or copyright) 2020-2023. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -19,6 +19,7 @@
 //! \file defs.hpp
 //  \brief contains Athena++ general purpose types, structures, enums, etc.
 
+#include <array>
 #include <cmath>
 #include <cstdint>
 #include <memory>
@@ -121,7 +122,7 @@ struct RegionSize {
 // tasks.hpp, ???
 
 // identifiers for boundary conditions
-enum class BoundaryFlag { block = -1, undef, reflect, outflow, periodic, user };
+enum class BoundaryFlag { block = -1, undef, periodic, user };
 
 // identifiers for all 6 faces of a MeshBlock
 constexpr int BOUNDARY_NFACES = 6;
@@ -134,6 +135,7 @@ enum BoundaryFace {
   inner_x3 = 4,
   outer_x3 = 5
 };
+using BValNames_t = std::array<std::string, BOUNDARY_NFACES>;
 
 inline BoundaryFace GetInnerBoundaryFace(CoordinateDirection dir) {
   if (dir == X1DIR) {
diff --git a/src/interface/swarm.cpp b/src/interface/swarm.cpp
index 294338dee0d2..39f1e8a28123 100644
--- a/src/interface/swarm.cpp
+++ b/src/interface/swarm.cpp
@@ -33,6 +33,7 @@ SwarmDeviceContext Swarm::GetDeviceContext() const {
   context.block_index_ = block_index_;
   context.neighbor_indices_ = neighbor_indices_;
   context.cell_sorted_ = cell_sorted_;
+  context.buffer_sorted_ = buffer_sorted_;
   context.cell_sorted_begin_ = cell_sorted_begin_;
   context.cell_sorted_number_ = cell_sorted_number_;
 
@@ -73,9 +74,10 @@ Swarm::Swarm(const std::string &label, const Metadata &metadata, const int nmax_
       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),
+      buffer_start_("buffer_start_", NMAX_NEIGHBORS),
       neighbor_received_particles_("neighbor_received_particles_", NMAX_NEIGHBORS),
-      cell_sorted_("cell_sorted_", nmax_pool_), mpiStatus(true) {
+      cell_sorted_("cell_sorted_", nmax_pool_),
+      buffer_sorted_("buffer_sorted_", nmax_pool_), mpiStatus(true) {
   PARTHENON_REQUIRE_THROWS(typeid(Coordinates_t) == typeid(UniformCartesian),
                            "SwarmDeviceContext only supports a uniform Cartesian mesh!");
 
@@ -209,6 +211,9 @@ void Swarm::SetPoolMax(const std::int64_t nmax_pool) {
   Kokkos::resize(cell_sorted_, nmax_pool);
   pmb->LogMemUsage(n_new * sizeof(SwarmKey));
 
+  Kokkos::resize(buffer_sorted_, nmax_pool);
+  pmb->LogMemUsage(n_new * sizeof(SwarmKey));
+
   block_index_.Resize(nmax_pool);
   pmb->LogMemUsage(n_new * sizeof(int));
 
@@ -490,35 +495,35 @@ void Swarm::SortParticlesByCell() {
             break;
           }
 
-          if (cell_sorted(start_index).cell_idx_1d_ == cell_idx_1d) {
+          if (cell_sorted(start_index).sort_idx_ == cell_idx_1d) {
             if (start_index == 0) {
               break;
-            } else if (cell_sorted(start_index - 1).cell_idx_1d_ != cell_idx_1d) {
+            } else if (cell_sorted(start_index - 1).sort_idx_ != cell_idx_1d) {
               break;
             } else {
               start_index--;
               continue;
             }
           }
-          if (cell_sorted(start_index).cell_idx_1d_ >= cell_idx_1d) {
+          if (cell_sorted(start_index).sort_idx_ >= cell_idx_1d) {
             start_index--;
             if (start_index < 0) {
               start_index = -1;
               break;
             }
-            if (cell_sorted(start_index).cell_idx_1d_ < cell_idx_1d) {
+            if (cell_sorted(start_index).sort_idx_ < cell_idx_1d) {
               start_index = -1;
               break;
             }
             continue;
           }
-          if (cell_sorted(start_index).cell_idx_1d_ < cell_idx_1d) {
+          if (cell_sorted(start_index).sort_idx_ < cell_idx_1d) {
             start_index++;
             if (start_index > max_active_index) {
               start_index = -1;
               break;
             }
-            if (cell_sorted(start_index).cell_idx_1d_ > cell_idx_1d) {
+            if (cell_sorted(start_index).sort_idx_ > cell_idx_1d) {
               start_index = -1;
               break;
             }
@@ -532,7 +537,7 @@ void Swarm::SortParticlesByCell() {
           int number = 0;
           int current_index = start_index;
           while (current_index <= max_active_index &&
-                 cell_sorted(current_index).cell_idx_1d_ == cell_idx_1d) {
+                 cell_sorted(current_index).sort_idx_ == cell_idx_1d) {
             current_index++;
             number++;
             cell_sorted_number(k, j, i) = number;
diff --git a/src/interface/swarm.hpp b/src/interface/swarm.hpp
index 2058dc58f894..49cd906c2a15 100644
--- a/src/interface/swarm.hpp
+++ b/src/interface/swarm.hpp
@@ -289,7 +289,7 @@ class Swarm {
   constexpr static int unset_index_ = -1;
 
   ParArray1D<int> num_particles_to_send_;
-  ParArray1D<int> buffer_counters_;
+  ParArray1D<int> buffer_start_;
   ParArray1D<int> neighbor_received_particles_;
   int total_received_particles_;
 
@@ -298,6 +298,9 @@ class Swarm {
   ParArray1D<SwarmKey>
       cell_sorted_; // 1D per-cell sorted array of key-value swarm memory indices
 
+  ParArray1D<SwarmKey>
+      buffer_sorted_; // 1D per-buffer sorted array of key-value swarm memory indices
+
   ParArrayND<int>
       cell_sorted_begin_; // Per-cell array of starting indices in cell_sorted_
 
diff --git a/src/interface/swarm_comms.cpp b/src/interface/swarm_comms.cpp
index 3ee4e2af7512..054ded34fabb 100644
--- a/src/interface/swarm_comms.cpp
+++ b/src/interface/swarm_comms.cpp
@@ -156,72 +156,14 @@ void Swarm::SetupPersistentMPI() {
   }
 }
 
-void Swarm::CountParticlesToSend_() {
-  auto mask_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), mask_);
-  auto swarm_d = GetDeviceContext();
-  auto pmb = GetBlockPointer();
-  const int nbmax = vbswarm->bd_var_.nbmax;
-
-  // Fence to make sure particles aren't currently being transported locally
-  // TODO(BRR) do this operation on device.
-  pmb->exec_space.fence();
-  const int particle_size = GetParticleDataSize();
-  vbswarm->particle_size = particle_size;
-
-  // TODO(BRR) This kernel launch should be folded into the subsequent logic once we
-  // convert that to kernel-based reductions
-  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();
-  const int max_active_index = GetMaxActiveIndex();
-  pmb->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);
-        }
-      });
-
-  // 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);
-          }
-        }
-      });
-
-  auto num_particles_to_send_h = num_particles_to_send_.GetHostMirrorAndCopy();
-
-  // Resize send buffers if too small
-  for (int n = 0; n < pmb->neighbors.size(); n++) {
-    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) {
-      sendbuf = BufArray1D<Real>("Buffer", num_particles_to_send_h(n) * particle_size);
-      vbswarm->bd_var_.send[bufid] = sendbuf;
-    }
-    vbswarm->send_size[bufid] = num_particles_to_send_h(n) * particle_size;
-  }
-}
-
 void Swarm::LoadBuffers_() {
   auto swarm_d = GetDeviceContext();
   auto pmb = GetBlockPointer();
   const int particle_size = GetParticleDataSize();
+  vbswarm->particle_size = particle_size;
   const int nneighbor = pmb->neighbors.size();
+  // Fence to make sure particles aren't currently being transported locally
+  pmb->exec_space.fence();
 
   auto &int_vector = std::get<getType<int>()>(vectors_);
   auto &real_vector = std::get<getType<Real>()>(vectors_);
@@ -236,42 +178,103 @@ void Swarm::LoadBuffers_() {
   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; });
+  if (max_active_index_ >= 0) {
+    auto &buffer_sorted = buffer_sorted_;
+    auto &buffer_start = buffer_start_;
 
-  auto &bdvar = vbswarm->bd_var_;
-  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_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 (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, n);
-              buffer_index++;
-            }
-            for (int i = 0; i < intPackDim; i++) {
-              bdvar.send[bufid](buffer_index) = static_cast<Real>(vint(i, n));
-              buffer_index++;
+    pmb->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;
+            const int m =
+                swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block);
+            buffer_sorted(n) = SwarmKey(m, n);
+          } else {
+            buffer_sorted(n) = SwarmKey(this_block_, n);
+          }
+        });
+
+    // sort by buffer index
+    sort(buffer_sorted, SwarmKeyComparator(), 0, max_active_index_);
+
+    // use discontinuity check to determine start of a buffer in buffer_sorted array
+    auto &num_particles_to_send = num_particles_to_send_;
+    auto max_active_index = max_active_index_;
+
+    // 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;
+          buffer_start[n] = 0;
+        });
+
+    pmb->par_for(
+        PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) {
+          auto m = buffer_sorted(n).sort_idx_;
+          // start checks (used for index of particle in buffer)
+          if (m >= 0 && n == 0) {
+            buffer_start(m) = 0;
+          } else if (m >= 0 && m != buffer_sorted(n - 1).sort_idx_) {
+            buffer_start(m) = n;
+          }
+          // end checks (used to to size particle buffers)
+          if (m >= 0 && n == max_active_index) {
+            num_particles_to_send(m) = n + 1;
+          } else if (m >= 0 && m != buffer_sorted(n + 1).sort_idx_) {
+            num_particles_to_send(m) = n + 1;
+          }
+        });
+
+    // copy values back to host for buffer sizing
+    auto num_particles_to_send_h = num_particles_to_send_.GetHostMirrorAndCopy();
+    auto buffer_start_h = buffer_start.GetHostMirrorAndCopy();
+
+    // Resize send buffers if too small
+    for (int n = 0; n < pmb->neighbors.size(); n++) {
+      num_particles_to_send_h(n) -= buffer_start_h(n);
+      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) {
+        sendbuf = BufArray1D<Real>("Buffer", num_particles_to_send_h(n) * particle_size);
+        vbswarm->bd_var_.send[bufid] = sendbuf;
+      }
+      vbswarm->send_size[bufid] = num_particles_to_send_h(n) * particle_size;
+    }
+
+    auto &bdvar = vbswarm->bd_var_;
+    auto neighbor_buffer_index = neighbor_buffer_index_;
+    // Loop over active particles buffer_sorted, use index n and buffer_start to
+    /// get index in buffer m, pack that particle in buffer
+    pmb->par_for(
+        PARTHENON_AUTO_LABEL, 0, max_active_index_, KOKKOS_LAMBDA(const int n) {
+          auto p_index = buffer_sorted(n).swarm_idx_;
+          if (swarm_d.IsActive(p_index)) {
+            const int m = buffer_sorted(n).sort_idx_;
+            const int bufid = neighbor_buffer_index(m);
+            if (m >= 0) {
+              const int bid = n - buffer_start[m];
+              int buffer_index = bid * particle_size;
+              swarm_d.MarkParticleForRemoval(p_index);
+              for (int i = 0; i < realPackDim; i++) {
+                bdvar.send[bufid](buffer_index) = vreal(i, p_index);
+                buffer_index++;
+              }
+              for (int i = 0; i < intPackDim; i++) {
+                bdvar.send[bufid](buffer_index) = static_cast<Real>(vint(i, p_index));
+                buffer_index++;
+              }
             }
           }
-        }
-      });
+        });
 
-  // Remove particles that were loaded to send to another block from this block
-  RemoveMarkedParticles();
+    // Remove particles that were loaded to send to another block from this block
+    RemoveMarkedParticles();
+  } else {
+    for (int n = 0; n < pmb->neighbors.size(); n++) {
+      const int bufid = pmb->neighbors[n].bufid;
+      vbswarm->send_size[bufid] = 0;
+    }
+  }
 }
 
 void Swarm::Send(BoundaryCommSubset phase) {
@@ -279,10 +282,8 @@ void Swarm::Send(BoundaryCommSubset phase) {
   const int nneighbor = pmb->neighbors.size();
   auto swarm_d = GetDeviceContext();
 
-  // Query particles for those to be sent
-  CountParticlesToSend_();
-
-  // Prepare buffers for send operations
+  // Potentially resize buffer, get consistent index from particle array, get ready to
+  // send
   LoadBuffers_();
 
   // Send buffer data
diff --git a/src/interface/swarm_device_context.hpp b/src/interface/swarm_device_context.hpp
index 936d2d56ad35..b16daf0d0cdf 100644
--- a/src/interface/swarm_device_context.hpp
+++ b/src/interface/swarm_device_context.hpp
@@ -25,16 +25,16 @@ struct SwarmKey {
   SwarmKey() {}
   KOKKOS_INLINE_FUNCTION
   SwarmKey(const int cell_idx_1d, const int swarm_idx_1d)
-      : cell_idx_1d_(cell_idx_1d), swarm_idx_(swarm_idx_1d) {}
+      : sort_idx_(cell_idx_1d), swarm_idx_(swarm_idx_1d) {}
 
-  int cell_idx_1d_;
+  int sort_idx_;
   int swarm_idx_;
 };
 
 struct SwarmKeyComparator {
   KOKKOS_INLINE_FUNCTION
   bool operator()(const SwarmKey &s1, const SwarmKey &s2) {
-    return s1.cell_idx_1d_ < s2.cell_idx_1d_;
+    return s1.sort_idx_ < s2.sort_idx_;
   }
 };
 
@@ -139,6 +139,7 @@ class SwarmDeviceContext {
   ParArrayND<int> block_index_;
   ParArrayND<int> neighbor_indices_; // 4x4x4 array of possible block AMR regions
   ParArray1D<SwarmKey> cell_sorted_;
+  ParArray1D<SwarmKey> buffer_sorted_;
   ParArrayND<int> cell_sorted_begin_;
   ParArrayND<int> cell_sorted_number_;
   int ndim_;
diff --git a/src/mesh/forest/forest.hpp b/src/mesh/forest/forest.hpp
index cf5f2780cbc1..db9de2e0025b 100644
--- a/src/mesh/forest/forest.hpp
+++ b/src/mesh/forest/forest.hpp
@@ -62,7 +62,8 @@ class ForestDefinition {
     face_sizes.emplace_back(xmin, xmax, ar3_t{1.0, 1.0, 1.0}, ai3_t{1, 1, 1});
   }
 
-  void AddBC(Edge edge, BoundaryFlag bf, std::optional<Edge> periodic_connection = {}) {
+  void AddBC(Edge edge, BoundaryFlag bf = BoundaryFlag::user,
+             std::optional<Edge> periodic_connection = {}) {
     if (bf == BoundaryFlag::periodic)
       PARTHENON_REQUIRE(periodic_connection,
                         "Must specify another edge for periodic boundary conditions.");
@@ -134,12 +135,12 @@ class Forest {
     return trees.at(loc.tree())->GetBlockBCs(loc);
   }
 
-  void EnrollBndryFncts(
-      ApplicationInput *app_in,
-      std::array<std::vector<BValFunc>, BOUNDARY_NFACES> UserBoundaryFunctions_in,
-      std::array<std::vector<SBValFunc>, BOUNDARY_NFACES> UserSwarmBoundaryFunctions_in) {
+  void EnrollBndryFncts(ApplicationInput *app_in, const BValNames_t &names,
+                        const BValNames_t &swarm_names,
+                        const BValFuncArray_t &UserBoundaryFunctions_in,
+                        const SBValFuncArray_t &UserSwarmBoundaryFunctions_in) {
     for (auto &[id, ptree] : trees)
-      ptree->EnrollBndryFncts(app_in, UserBoundaryFunctions_in,
+      ptree->EnrollBndryFncts(app_in, names, swarm_names, UserBoundaryFunctions_in,
                               UserSwarmBoundaryFunctions_in);
   }
 
diff --git a/src/mesh/forest/tree.cpp b/src/mesh/forest/tree.cpp
index fe1803882401..d1d019ee5bd4 100644
--- a/src/mesh/forest/tree.cpp
+++ b/src/mesh/forest/tree.cpp
@@ -24,12 +24,14 @@
 #include <utility>
 #include <vector>
 
+#include "application_input.hpp"
 #include "basic_types.hpp"
 #include "defs.hpp"
 #include "mesh/forest/forest_topology.hpp"
 #include "mesh/forest/logical_coordinate_transformation.hpp"
 #include "mesh/forest/logical_location.hpp"
 #include "mesh/forest/tree.hpp"
+#include "parameter_input.hpp"
 #include "utils/bit_hacks.hpp"
 #include "utils/indexer.hpp"
 
@@ -383,70 +385,21 @@ std::int64_t Tree::GetOldGid(const LogicalLocation &loc) const {
   return -1;
 }
 
-void Tree::EnrollBndryFncts(
-    ApplicationInput *app_in,
-    std::array<std::vector<BValFunc>, BOUNDARY_NFACES> UserBoundaryFunctions_in,
-    std::array<std::vector<SBValFunc>, BOUNDARY_NFACES> UserSwarmBoundaryFunctions_in) {
+void Tree::EnrollBndryFncts(ApplicationInput *app_in, const BValNames_t &names,
+                            const BValNames_t &swarm_names,
+                            const BValFuncArray_t &UserBoundaryFunctions_in,
+                            const SBValFuncArray_t &UserSwarmBoundaryFunctions_in) {
   UserBoundaryFunctions = UserBoundaryFunctions_in;
   UserSwarmBoundaryFunctions = UserSwarmBoundaryFunctions_in;
-  static const BValFunc outflow[6] = {
-      BoundaryFunction::OutflowInnerX1, BoundaryFunction::OutflowOuterX1,
-      BoundaryFunction::OutflowInnerX2, BoundaryFunction::OutflowOuterX2,
-      BoundaryFunction::OutflowInnerX3, BoundaryFunction::OutflowOuterX3};
-  static const BValFunc reflect[6] = {
-      BoundaryFunction::ReflectInnerX1, BoundaryFunction::ReflectOuterX1,
-      BoundaryFunction::ReflectInnerX2, BoundaryFunction::ReflectOuterX2,
-      BoundaryFunction::ReflectInnerX3, BoundaryFunction::ReflectOuterX3};
-  static const SBValFunc soutflow[6] = {
-      BoundaryFunction::SwarmOutflowInnerX1, BoundaryFunction::SwarmOutflowOuterX1,
-      BoundaryFunction::SwarmOutflowInnerX2, BoundaryFunction::SwarmOutflowOuterX2,
-      BoundaryFunction::SwarmOutflowInnerX3, BoundaryFunction::SwarmOutflowOuterX3};
-  static const SBValFunc speriodic[6] = {
-      BoundaryFunction::SwarmPeriodicInnerX1, BoundaryFunction::SwarmPeriodicOuterX1,
-      BoundaryFunction::SwarmPeriodicInnerX2, BoundaryFunction::SwarmPeriodicOuterX2,
-      BoundaryFunction::SwarmPeriodicInnerX3, BoundaryFunction::SwarmPeriodicOuterX3};
 
   for (int f = 0; f < BOUNDARY_NFACES; f++) {
-    switch (boundary_conditions[f]) {
-    case BoundaryFlag::reflect:
-      MeshBndryFnctn[f] = reflect[f];
-      break;
-    case BoundaryFlag::outflow:
-      MeshBndryFnctn[f] = outflow[f];
-      SwarmBndryFnctn[f] = soutflow[f];
-      break;
-    case BoundaryFlag::user:
-      if (app_in->boundary_conditions[f] != nullptr) {
-        MeshBndryFnctn[f] = app_in->boundary_conditions[f];
-      } else {
-        std::stringstream msg;
-        msg << "A user boundary condition for face " << f
-            << " was requested. but no condition was enrolled." << std::endl;
-        PARTHENON_THROW(msg);
-      }
-      break;
-    default: // periodic/block BCs handled elsewhere.
-      break;
-    }
-
-    switch (boundary_conditions[f]) {
-    case BoundaryFlag::outflow:
-      SwarmBndryFnctn[f] = soutflow[f];
-      break;
-    case BoundaryFlag::periodic:
-      SwarmBndryFnctn[f] = speriodic[f];
-      break;
-    case BoundaryFlag::reflect:
-      // Default "reflect" boundaries not available for swarms; catch later on if swarms
-      // are present
-      break;
-    case BoundaryFlag::user:
-      if (app_in->swarm_boundary_conditions[f] != nullptr) {
-        SwarmBndryFnctn[f] = app_in->swarm_boundary_conditions[f];
-      }
-      break;
-    default: // Default BCs handled elsewhere
-      break;
+    auto flag = boundary_conditions[f];
+    auto face = static_cast<BoundaryFace>(f);
+    if (flag == BoundaryFlag::user) {
+      MeshBndryFnctn[f] = app_in->GetBoundaryCondition(face, names[f]);
+      SwarmBndryFnctn[f] = app_in->GetSwarmBoundaryCondition(face, swarm_names[f]);
+    } else if (flag == BoundaryFlag::periodic) {
+      SwarmBndryFnctn[f] = app_in->GetSwarmBoundaryCondition(face, swarm_names[f]);
     }
   }
 }
diff --git a/src/mesh/forest/tree.hpp b/src/mesh/forest/tree.hpp
index 7b53073ac92b..f423b0180f78 100644
--- a/src/mesh/forest/tree.hpp
+++ b/src/mesh/forest/tree.hpp
@@ -23,7 +23,6 @@
 #include <utility>
 #include <vector>
 
-#include "application_input.hpp"
 #include "basic_types.hpp"
 #include "bvals/boundary_conditions.hpp"
 #include "defs.hpp"
@@ -35,6 +34,9 @@
 #include "utils/indexer.hpp"
 
 namespace parthenon {
+class ApplicationInput;
+class ParameterInput;
+
 namespace forest {
 
 // We don't explicitly allow for periodic boundaries, since we can encode periodicity
@@ -117,10 +119,10 @@ class Tree : public std::enable_shared_from_this<Tree> {
   std::vector<std::shared_ptr<Node>> forest_nodes;
 
   // Boundary Functions
-  void EnrollBndryFncts(
-      ApplicationInput *app_in,
-      std::array<std::vector<BValFunc>, BOUNDARY_NFACES> UserBoundaryFunctions_in,
-      std::array<std::vector<SBValFunc>, BOUNDARY_NFACES> UserSwarmBoundaryFunctions_in);
+  void EnrollBndryFncts(ApplicationInput *app_in, const BValNames_t &names,
+                        const BValNames_t &swarm_names,
+                        const BValFuncArray_t &UserBoundaryFunctions_in,
+                        const SBValFuncArray_t &UserSwarmBoundaryFunctions_in);
   BValFunc MeshBndryFnctn[BOUNDARY_NFACES];
   SBValFunc SwarmBndryFnctn[BOUNDARY_NFACES];
   std::array<std::vector<BValFunc>, BOUNDARY_NFACES> UserBoundaryFunctions;
diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index c45f4f9679f2..86d727bae46b 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -151,13 +151,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
        pin->GetInteger("parthenon/mesh", "nx3")},
       {false, pin->GetInteger("parthenon/mesh", "nx2") == 1,
        pin->GetInteger("parthenon/mesh", "nx3") == 1});
-  mesh_bcs = {
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))};
+  SetBCNames_(pin);
+  mesh_bcs = GetBCsFromNames_(mesh_bc_names);
   ndim = (mesh_size.nx(X3DIR) > 1) ? 3 : ((mesh_size.nx(X2DIR) > 1) ? 2 : 1);
 
   for (auto &[dir, label] : std::vector<std::tuple<CoordinateDirection, std::string>>{
@@ -175,7 +170,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
   // Load balancing flag and parameters
   forest = forest::Forest::HyperRectangular(mesh_size, base_block_size, mesh_bcs);
   root_level = forest.root_level;
-  forest.EnrollBndryFncts(app_in, resolved_packages->UserBoundaryFunctions,
+  forest.EnrollBndryFncts(app_in, mesh_bc_names, mesh_swarm_bc_names,
+                          resolved_packages->UserBoundaryFunctions,
                           resolved_packages->UserSwarmBoundaryFunctions);
 
   // SMR / AMR:
@@ -195,13 +191,9 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
     : Mesh(pin, app_in, packages, base_constructor_selector_t()) {
   mesh_size =
       RegionSize({0, 0, 0}, {1, 1, 0}, {1, 1, 1}, {1, 1, 1}, {false, false, true});
-  mesh_bcs = {
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix1_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox1_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix2_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox2_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ix3_bc", "reflecting")),
-      GetBoundaryFlag(pin->GetOrAddString("parthenon/mesh", "ox3_bc", "reflecting"))};
+
+  SetBCNames_(pin);
+  mesh_bcs = GetBCsFromNames_(mesh_bc_names);
   for (auto &[dir, label] : std::vector<std::tuple<CoordinateDirection, std::string>>{
            {X1DIR, "nx1"}, {X2DIR, "nx2"}, {X3DIR, "nx3"}}) {
     base_block_size.xrat(dir) = mesh_size.xrat(dir);
@@ -219,7 +211,8 @@ Mesh::Mesh(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
   // Load balancing flag and parameters
   forest = forest::Forest::Make2D(forest_def);
   root_level = forest.root_level;
-  forest.EnrollBndryFncts(app_in, resolved_packages->UserBoundaryFunctions,
+  forest.EnrollBndryFncts(app_in, mesh_bc_names, mesh_swarm_bc_names,
+                          resolved_packages->UserBoundaryFunctions,
                           resolved_packages->UserSwarmBoundaryFunctions);
   BuildBlockList(pin, app_in, packages, 0);
 }
@@ -1183,4 +1176,44 @@ Mesh::GetLevelsAndLogicalLocationsFlat() const noexcept {
   return std::make_pair(levels, logicalLocations);
 }
 
+void Mesh::SetBCNames_(ParameterInput *pin) {
+  mesh_bc_names = {pin->GetOrAddString("parthenon/mesh", "ix1_bc", "outflow"),
+                   pin->GetOrAddString("parthenon/mesh", "ox1_bc", "outflow"),
+                   pin->GetOrAddString("parthenon/mesh", "ix2_bc", "outflow"),
+                   pin->GetOrAddString("parthenon/mesh", "ox2_bc", "outflow"),
+                   pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"),
+                   pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")};
+  // JMM: This is needed because not all BCs are necessarily
+  // implemented for swarms
+  auto maybe = [](const std::string &s) {
+    return ((s == "outflow") || (s == "periodic")) ? s : "outflow";
+  };
+  mesh_swarm_bc_names = {
+      pin->GetOrAddString("parthenon/swarm", "ix1_bc", maybe(mesh_bc_names[0])),
+      pin->GetOrAddString("parthenon/swarm", "ox1_bc", maybe(mesh_bc_names[1])),
+      pin->GetOrAddString("parthenon/swarm", "ix2_bc", maybe(mesh_bc_names[2])),
+      pin->GetOrAddString("parthenon/swarm", "ox2_bc", maybe(mesh_bc_names[3])),
+      pin->GetOrAddString("parthenon/swarm", "ix3_bc", maybe(mesh_bc_names[4])),
+      pin->GetOrAddString("parthenon/swarm", "ox3_bc", maybe(mesh_bc_names[5]))};
+  // JMM: A consequence of having only one boundary flag array but
+  // multiple boundary function arrays is that swarms *must* be
+  // periodic if the mesh is periodic but otherwise mesh and swarm
+  // boundaries are decoupled.
+  for (int i = 0; i < BOUNDARY_NFACES; ++i) {
+    if (mesh_bc_names[i] == "periodic") {
+      PARTHENON_REQUIRE(mesh_swarm_bc_names[i] == "periodic",
+                        "If the mesh is periodic, swarms must be also.");
+    }
+  }
+}
+
+std::array<BoundaryFlag, BOUNDARY_NFACES>
+Mesh::GetBCsFromNames_(const std::array<std::string, BOUNDARY_NFACES> &names) const {
+  std::array<BoundaryFlag, BOUNDARY_NFACES> out;
+  for (int f = 0; f < BOUNDARY_NFACES; ++f) {
+    out[f] = GetBoundaryFlag(names[f]);
+  }
+  return out;
+}
+
 } // namespace parthenon
diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp
index 88b99c333ad4..ed26be1c6a1b 100644
--- a/src/mesh/mesh.hpp
+++ b/src/mesh/mesh.hpp
@@ -116,6 +116,11 @@ class Mesh {
   bool is_restart;
   RegionSize mesh_size;
   RegionSize base_block_size;
+
+  BValNames_t mesh_bc_names;
+  BValNames_t mesh_swarm_bc_names;
+
+  // these are flags not boundary functions
   std::array<BoundaryFlag, BOUNDARY_NFACES> mesh_bcs;
   int ndim; // number of dimensions
   const bool adaptive, multilevel, multigrid;
@@ -297,6 +302,10 @@ class Mesh {
   std::unordered_map<std::string, MPI_Comm> mpi_comm_map_;
 #endif
 
+  void SetBCNames_(ParameterInput *pin);
+  std::array<BoundaryFlag, BOUNDARY_NFACES>
+  GetBCsFromNames_(const BValNames_t &names) const;
+
   // functions
   void CheckMeshValidity() const;
   void BuildBlockList(ParameterInput *pin, ApplicationInput *app_in, Packages_t &packages,
diff --git a/src/outputs/output_utils.hpp b/src/outputs/output_utils.hpp
index 5f95fbbaedc2..5c094ef688d7 100644
--- a/src/outputs/output_utils.hpp
+++ b/src/outputs/output_utils.hpp
@@ -50,7 +50,7 @@ namespace OutputUtils {
 // Helper struct containing some information about a variable
 struct VarInfo {
  public:
-  static constexpr int VNDIM = MAX_VARIABLE_DIMENSION;
+  static constexpr const int VNDIM = MAX_VARIABLE_DIMENSION;
   std::string label;
   int num_components;
   int tensor_rank; // 0- to 3-D for cell-centered variables, 0- to 6-D for arbitrary shape
diff --git a/src/outputs/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp
index 7d65d366cf77..cd5009d490d9 100644
--- a/src/outputs/parthenon_hdf5.cpp
+++ b/src/outputs/parthenon_hdf5.cpp
@@ -190,12 +190,8 @@ void PHDF5Output::WriteOutputFileImpl(Mesh *pm, ParameterInput *pin, SimTime *tm
                        info_group);
 
     // Boundary conditions
-    std::vector<std::string> boundary_condition_str(BOUNDARY_NFACES);
-    for (size_t i = 0; i < boundary_condition_str.size(); i++) {
-      boundary_condition_str[i] = GetBoundaryString(pm->mesh_bcs[i]);
-    }
-
-    HDF5WriteAttribute("BoundaryConditions", boundary_condition_str, info_group);
+    HDF5WriteAttribute("BoundaryConditions", pm->mesh_bc_names, info_group);
+    HDF5WriteAttribute("SwarmBoundaryConditions", pm->mesh_swarm_bc_names, info_group);
     Kokkos::Profiling::popRegion(); // write Info
   }                                 // Info section
 
diff --git a/src/outputs/parthenon_hdf5.hpp b/src/outputs/parthenon_hdf5.hpp
index 72deaab90681..e99c58b14815 100644
--- a/src/outputs/parthenon_hdf5.hpp
+++ b/src/outputs/parthenon_hdf5.hpp
@@ -18,6 +18,7 @@
 #define OUTPUTS_PARTHENON_HDF5_HPP_
 
 #include "config.hpp"
+#include "defs.hpp"
 
 #include "kokkos_abstraction.hpp"
 #include "parthenon_arrays.hpp"
@@ -122,6 +123,9 @@ void HDF5WriteAttribute(const std::string &name, size_t num_values, const T *dat
 // In CPP file
 void HDF5WriteAttribute(const std::string &name, const std::string &value,
                         hid_t location);
+void HDF5WriteAttribute(const std::string &name,
+                        const std::array<std::string, BOUNDARY_NFACES> &values,
+                        hid_t location);
 
 template <typename T>
 void HDF5WriteAttribute(const std::string &name, const std::vector<T> &values,
diff --git a/src/outputs/parthenon_hdf5_attributes.cpp b/src/outputs/parthenon_hdf5_attributes.cpp
index 45cece79881c..ef1d3c74c462 100644
--- a/src/outputs/parthenon_hdf5_attributes.cpp
+++ b/src/outputs/parthenon_hdf5_attributes.cpp
@@ -33,6 +33,7 @@
 #include <tuple>
 #include <vector>
 
+#include "defs.hpp"
 #include "kokkos_abstraction.hpp"
 #include "utils/concepts_lite.hpp"
 #include "utils/error_checking.hpp"
@@ -94,6 +95,16 @@ void HDF5WriteAttribute(const std::string &name, const std::vector<std::string>
   HDF5WriteAttribute(name, char_ptrs, location);
 }
 
+void HDF5WriteAttribute(const std::string &name,
+                        const std::array<std::string, BOUNDARY_NFACES> &values,
+                        hid_t location) {
+  std::vector<const char *> char_ptrs(values.size());
+  for (size_t i = 0; i < values.size(); ++i) {
+    char_ptrs[i] = values[i].c_str();
+  }
+  HDF5WriteAttribute(name, char_ptrs, location);
+}
+
 template <>
 std::vector<std::string> HDF5ReadAttributeVec(hid_t location, const std::string &name) {
   // get strings as char pointers, HDF5 will allocate the memory and we need to free it
diff --git a/src/outputs/restart.hpp b/src/outputs/restart.hpp
index 22b059fc3c5e..138949dd1667 100644
--- a/src/outputs/restart.hpp
+++ b/src/outputs/restart.hpp
@@ -86,7 +86,6 @@ class RestartReader {
 
   struct MeshInfo {
     int nbnew, nbdel, nbtotal, root_level, includes_ghost, n_ghost;
-    std::vector<std::string> bound_cond;
     std::vector<int> block_size;
     std::vector<Real> grid_dim;
     std::vector<int64_t> lx123;
diff --git a/src/outputs/restart_hdf5.cpp b/src/outputs/restart_hdf5.cpp
index c7584f954b8d..8078fbd5fa9a 100644
--- a/src/outputs/restart_hdf5.cpp
+++ b/src/outputs/restart_hdf5.cpp
@@ -3,7 +3,7 @@
 // Copyright(C) 2020-2024 The Parthenon collaboration
 // Licensed under the 3-clause BSD License, see LICENSE file for details
 //========================================================================================
-// (C) (or copyright) 2020-2021. Triad National Security, LLC. All rights reserved.
+// (C) (or copyright) 2020-2024. Triad National Security, LLC. All rights reserved.
 //
 // This program was produced under U.S. Government contract 89233218CNA000001 for Los
 // Alamos National Laboratory (LANL), which is operated by Triad National Security, LLC
@@ -126,8 +126,6 @@ RestartReaderHDF5::MeshInfo RestartReaderHDF5::GetMeshInfo() const {
   mesh_info.nbtotal = GetAttr<int>("Info", "NumMeshBlocks");
   mesh_info.root_level = GetAttr<int>("Info", "RootLevel");
 
-  mesh_info.bound_cond = GetAttrVec<std::string>("Info", "BoundaryConditions");
-
   mesh_info.block_size = GetAttrVec<int>("Info", "MeshBlockSize");
   mesh_info.includes_ghost = GetAttr<int>("Info", "IncludesGhost");
   mesh_info.n_ghost = GetAttr<int>("Info", "NGhost");
@@ -215,7 +213,7 @@ void RestartReaderHDF5::ReadBlocks(const std::string &name, IndexRange range,
 #else  // HDF5 enabled
   auto hdl = OpenDataset<Real>(name);
 
-  const int VNDIM = info.VNDIM;
+  constexpr int VNDIM = OutputUtils::VarInfo::VNDIM;
 
   /** Select hyperslab in dataset **/
   int total_dim = 0;
diff --git a/src/utils/hash.hpp b/src/utils/hash.hpp
index 77642a64dba3..4f1ff6557c74 100644
--- a/src/utils/hash.hpp
+++ b/src/utils/hash.hpp
@@ -31,9 +31,8 @@ std::size_t hash_combine(std::size_t lhs, const T &v, Rest &&...rest) {
   lhs ^= rhs + 0x9e3779b9 + (lhs << 6) + (lhs >> 2);
   if constexpr (sizeof...(Rest) > 0) {
     return hash_combine(lhs, std::forward<Rest>(rest)...);
-  } else {
-    return lhs;
   }
+  return lhs;
 }
 
 template <class Tup, std::size_t I = std::tuple_size<Tup>::value - 1>
diff --git a/src/utils/sort.hpp b/src/utils/sort.hpp
index 97e9c77a88e4..a4ab8139dab3 100644
--- a/src/utils/sort.hpp
+++ b/src/utils/sort.hpp
@@ -61,7 +61,7 @@ void sort(ParArray1D<Key> data, KeyComparator comparator, size_t min_idx,
           size_t max_idx) {
   PARTHENON_DEBUG_REQUIRE(min_idx < data.extent(0), "Invalid minimum sort index!");
   PARTHENON_DEBUG_REQUIRE(max_idx < data.extent(0), "Invalid maximum sort index!");
-#ifdef KOKKOS_ENABLE_CUDA
+#if defined(KOKKOS_ENABLE_CUDA)
 #ifdef __clang__
   PARTHENON_FAIL("sort is using thrust and there exists an incompatibility with clang, "
                  "see https://github.com/lanl/parthenon/issues/647 for more details. We "
@@ -74,6 +74,13 @@ void sort(ParArray1D<Key> data, KeyComparator comparator, size_t min_idx,
   thrust::device_ptr<Key> last_d = thrust::device_pointer_cast(data.data()) + max_idx + 1;
   thrust::sort(first_d, last_d, comparator);
 #endif
+#elif defined(KOKKOS_ENABLE_HIP)
+  auto data_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), data);
+  std::sort(data_h.data() + min_idx, data_h.data() + max_idx + 1, comparator);
+  Kokkos::deep_copy(data, data_h);
+  // TODO(BRR) With Kokkos 4.4, switch to Kokkos::sort
+  // auto sub_data = Kokkos::subview(data, std::make_pair(min_idx, max_idx + 1));
+  // Kokkos::sort(sub_data, comparator);
 #else
   if (std::is_same<DevExecSpace, HostExecSpace>::value) {
     std::sort(data.data() + min_idx, data.data() + max_idx + 1, comparator);
@@ -89,7 +96,7 @@ template <class Key>
 void sort(ParArray1D<Key> data, size_t min_idx, size_t max_idx) {
   PARTHENON_DEBUG_REQUIRE(min_idx < data.extent(0), "Invalid minimum sort index!");
   PARTHENON_DEBUG_REQUIRE(max_idx < data.extent(0), "Invalid maximum sort index!");
-#ifdef KOKKOS_ENABLE_CUDA
+#if defined(KOKKOS_ENABLE_CUDA) || defined(KOKKOS_ENABLE_HIP)
 #ifdef __clang__
   PARTHENON_FAIL("sort is using thrust and there exists an incompatibility with clang, "
                  "see https://github.com/lanl/parthenon/issues/647 for more details. We "
@@ -102,6 +109,12 @@ void sort(ParArray1D<Key> data, size_t min_idx, size_t max_idx) {
   thrust::device_ptr<Key> last_d = thrust::device_pointer_cast(data.data()) + max_idx + 1;
   thrust::sort(first_d, last_d);
 #endif
+  auto data_h = Kokkos::create_mirror_view_and_copy(HostMemSpace(), data);
+  std::sort(data_h.data() + min_idx, data_h.data() + max_idx + 1);
+  Kokkos::deep_copy(data, data_h);
+  // TODO(BRR) With Kokkos 4.4, switch to Kokkos::sort
+  // auto sub_data = Kokkos::subview(data, std::make_pair(min_idx, max_idx + 1));
+  // Kokkos::sort(sub_data);
 #else
   if (std::is_same<DevExecSpace, HostExecSpace>::value) {
     std::sort(data.data() + min_idx, data.data() + max_idx + 1);
diff --git a/tst/unit/test_swarm.cpp b/tst/unit/test_swarm.cpp
index 22f6b1f4d735..ada71b112537 100644
--- a/tst/unit/test_swarm.cpp
+++ b/tst/unit/test_swarm.cpp
@@ -94,7 +94,7 @@ TEST_CASE("Swarm memory management", "[Swarm]") {
   mesh->mesh_bcs[0] = BoundaryFlag::user;
   meshblock->boundary_flag[0] = BoundaryFlag::user;
   for (int i = 1; i < 6; i++) {
-    mesh->mesh_bcs[i] = BoundaryFlag::outflow;
+    mesh->mesh_bcs[i] = BoundaryFlag::user;
     meshblock->boundary_flag[i] = BoundaryFlag::user;
   }
   meshblock->pmy_mesh = mesh.get();