From 7a9b2f855f559f38c89649f641421e55c85f3d60 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 12 Sep 2024 19:15:06 -0600
Subject: [PATCH 01/23] it compiles at least

---
 .../boundary_exchange_driver.cpp              | 16 +--
 example/particles/main.cpp                    | 14 +--
 src/CMakeLists.txt                            |  3 +
 src/application_input.cpp                     | 97 +++++++++++++++++++
 src/application_input.hpp                     | 39 ++++++--
 src/bvals/boundary_conditions.cpp             |  4 +-
 src/bvals/boundary_conditions.hpp             |  6 +-
 src/bvals/boundary_flag.cpp                   | 48 +--------
 src/bvals/bvals.hpp                           |  1 -
 src/defs.hpp                                  |  6 +-
 src/mesh/forest/forest.hpp                    | 13 +--
 src/mesh/forest/tree.cpp                      | 73 +++-----------
 src/mesh/forest/tree.hpp                      | 12 ++-
 src/mesh/mesh.cpp                             | 50 +++++++---
 src/mesh/mesh.hpp                             | 13 +++
 src/outputs/parthenon_hdf5.cpp                |  8 +-
 src/outputs/parthenon_hdf5.hpp                |  4 +
 src/outputs/parthenon_hdf5_attributes.cpp     | 11 +++
 src/outputs/restart.hpp                       |  1 -
 src/outputs/restart_hdf5.cpp                  |  2 -
 tst/unit/test_swarm.cpp                       |  2 +-
 21 files changed, 256 insertions(+), 167 deletions(-)
 create mode 100644 src/application_input.cpp

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/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..87721ef36072
--- /dev/null
+++ b/src/application_input.cpp
@@ -0,0 +1,97 @@
+//========================================================================================
+// (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::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..7abd6990e5be 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,27 @@ 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 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/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 f88c6fc645e6..db199245bbc1 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -149,13 +149,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>>{
@@ -173,7 +168,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:
@@ -193,13 +189,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);
@@ -217,7 +209,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);
 }
@@ -1180,4 +1173,29 @@ 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")};
+  mesh_swarm_bc_names = {
+      pin->GetOrAddString("parthenon/mesh", "ix1_bc", mesh_bc_names[0]),
+      pin->GetOrAddString("parthenon/mesh", "ox1_bc", mesh_bc_names[1]),
+      pin->GetOrAddString("parthenon/mesh", "ix2_bc", mesh_bc_names[2]),
+      pin->GetOrAddString("parthenon/mesh", "ox2_bc", mesh_bc_names[3]),
+      pin->GetOrAddString("parthenon/mesh", "ix3_bc", mesh_bc_names[4]),
+      pin->GetOrAddString("parthenon/mesh", "ox3_bc", mesh_bc_names[5])};
+}
+
+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 684a897aad56..918d50c175cb 100644
--- a/src/mesh/mesh.hpp
+++ b/src/mesh/mesh.hpp
@@ -116,6 +116,15 @@ 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
+  // 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.
   std::array<BoundaryFlag, BOUNDARY_NFACES> mesh_bcs;
   int ndim; // number of dimensions
   const bool adaptive, multilevel, multigrid;
@@ -297,6 +306,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/parthenon_hdf5.cpp b/src/outputs/parthenon_hdf5.cpp
index 3f7c88835c69..f0a4d876514e 100644
--- a/src/outputs/parthenon_hdf5.cpp
+++ b/src/outputs/parthenon_hdf5.cpp
@@ -185,12 +185,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 f745dd30d0ba..5248ddf5392a 100644
--- a/src/outputs/parthenon_hdf5_attributes.cpp
+++ b/src/outputs/parthenon_hdf5_attributes.cpp
@@ -32,6 +32,7 @@
 #include <tuple>
 #include <vector>
 
+#include "defs.hpp"
 #include "kokkos_abstraction.hpp"
 #include "utils/concepts_lite.hpp"
 #include "utils/error_checking.hpp"
@@ -93,6 +94,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 f16bfccc52de..3938b97ae436 100644
--- a/src/outputs/restart_hdf5.cpp
+++ b/src/outputs/restart_hdf5.cpp
@@ -125,8 +125,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");
diff --git a/tst/unit/test_swarm.cpp b/tst/unit/test_swarm.cpp
index bcff16106110..8221b7040bed 100644
--- a/tst/unit/test_swarm.cpp
+++ b/tst/unit/test_swarm.cpp
@@ -92,7 +92,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();

From 401e0c2022285bd9017b78a96edb14457031f23f Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 12 Sep 2024 19:29:42 -0600
Subject: [PATCH 02/23] add easy machinery to register reflecting BCs

---
 doc/sphinx/src/boundary_conditions.rst | 31 +++++++++++++++++---------
 src/application_input.cpp              | 10 +++++++++
 src/application_input.hpp              |  1 +
 3 files changed, 32 insertions(+), 10 deletions(-)

diff --git a/doc/sphinx/src/boundary_conditions.rst b/doc/sphinx/src/boundary_conditions.rst
index cb0e08532492..2f23bfe1f744 100644
--- a/doc/sphinx/src/boundary_conditions.rst
+++ b/doc/sphinx/src/boundary_conditions.rst
@@ -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 int he 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/src/application_input.cpp b/src/application_input.cpp
index 87721ef36072..42f7eb116272 100644
--- a/src/application_input.cpp
+++ b/src/application_input.cpp
@@ -47,6 +47,16 @@ ApplicationInput::ApplicationInput() {
   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) {
diff --git a/src/application_input.hpp b/src/application_input.hpp
index 7abd6990e5be..9c91585bfb3e 100644
--- a/src/application_input.hpp
+++ b/src/application_input.hpp
@@ -81,6 +81,7 @@ class ApplicationInput {
   void RegisterBoundaryCondition(BoundaryFace face, BValFunc condition) {
     RegisterBoundaryCondition(face, "user", condition);
   }
+  void RegisterDefaultReflectingBoundaryConditions();
 
   void RegisterSwarmBoundaryCondition(BoundaryFace face, const std::string &name,
                                       SBValFunc condition);

From 1d2ae51e7405ce11cab63396d2f16ea9fa62b599 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Thu, 12 Sep 2024 19:33:21 -0600
Subject: [PATCH 03/23] changelog

---
 CHANGELOG.md | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 510e3a6def57..041b071a4149 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -6,6 +6,7 @@
 - [[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/...)
+- [[PR1177]](https://github.com/parthenon-hpc-lab/parthenon/pull/1177) Make mesh-level boundary conditions usable without the "user" flag
 - [[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
 
@@ -19,7 +20,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

From c656f13455000b870acd60aa7d5f1b09808c8153 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 13 Sep 2024 08:33:34 -0600
Subject: [PATCH 04/23] swarm bcs differetn from mesh bcs

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

diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index db199245bbc1..af20881b2d09 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -1181,12 +1181,12 @@ void Mesh::SetBCNames_(ParameterInput *pin) {
                    pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"),
                    pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")};
   mesh_swarm_bc_names = {
-      pin->GetOrAddString("parthenon/mesh", "ix1_bc", mesh_bc_names[0]),
-      pin->GetOrAddString("parthenon/mesh", "ox1_bc", mesh_bc_names[1]),
-      pin->GetOrAddString("parthenon/mesh", "ix2_bc", mesh_bc_names[2]),
-      pin->GetOrAddString("parthenon/mesh", "ox2_bc", mesh_bc_names[3]),
-      pin->GetOrAddString("parthenon/mesh", "ix3_bc", mesh_bc_names[4]),
-      pin->GetOrAddString("parthenon/mesh", "ox3_bc", mesh_bc_names[5])};
+      pin->GetOrAddString("parthenon/mesh", "ix1_swarm_bc", mesh_bc_names[0]),
+      pin->GetOrAddString("parthenon/mesh", "ox1_swarm_bc", mesh_bc_names[1]),
+      pin->GetOrAddString("parthenon/mesh", "ix2_swarm_bc", mesh_bc_names[2]),
+      pin->GetOrAddString("parthenon/mesh", "ox2_swarm_bc", mesh_bc_names[3]),
+      pin->GetOrAddString("parthenon/mesh", "ix3_swarm_bc", mesh_bc_names[4]),
+      pin->GetOrAddString("parthenon/mesh", "ox3_swarm_bc", mesh_bc_names[5])};
 }
 
 std::array<BoundaryFlag, BOUNDARY_NFACES>

From 30f0ecc8962cff2cbca0ce640267b183bc885131 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 13 Sep 2024 08:33:58 -0600
Subject: [PATCH 05/23] use new input block

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

diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index af20881b2d09..a822f9b03797 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -1181,12 +1181,12 @@ void Mesh::SetBCNames_(ParameterInput *pin) {
                    pin->GetOrAddString("parthenon/mesh", "ix3_bc", "outflow"),
                    pin->GetOrAddString("parthenon/mesh", "ox3_bc", "outflow")};
   mesh_swarm_bc_names = {
-      pin->GetOrAddString("parthenon/mesh", "ix1_swarm_bc", mesh_bc_names[0]),
-      pin->GetOrAddString("parthenon/mesh", "ox1_swarm_bc", mesh_bc_names[1]),
-      pin->GetOrAddString("parthenon/mesh", "ix2_swarm_bc", mesh_bc_names[2]),
-      pin->GetOrAddString("parthenon/mesh", "ox2_swarm_bc", mesh_bc_names[3]),
-      pin->GetOrAddString("parthenon/mesh", "ix3_swarm_bc", mesh_bc_names[4]),
-      pin->GetOrAddString("parthenon/mesh", "ox3_swarm_bc", mesh_bc_names[5])};
+      pin->GetOrAddString("parthenon/swarm", "ix1_bc", mesh_bc_names[0]),
+      pin->GetOrAddString("parthenon/swarm", "ox1_bc", mesh_bc_names[1]),
+      pin->GetOrAddString("parthenon/swarm", "ix2_bc", mesh_bc_names[2]),
+      pin->GetOrAddString("parthenon/swarm", "ox2_bc", mesh_bc_names[3]),
+      pin->GetOrAddString("parthenon/swarm", "ix3_bc", mesh_bc_names[4]),
+      pin->GetOrAddString("parthenon/swarm", "ox3_bc", mesh_bc_names[5])};
 }
 
 std::array<BoundaryFlag, BOUNDARY_NFACES>

From 1df88056482bda344e3748f787a615ba16b3a218 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 27 Sep 2024 17:02:24 -0600
Subject: [PATCH 06/23] typo

---
 doc/sphinx/src/boundary_conditions.rst | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/doc/sphinx/src/boundary_conditions.rst b/doc/sphinx/src/boundary_conditions.rst
index 2f23bfe1f744..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
 ===================
@@ -108,7 +108,7 @@ are available `here <https://github.com/parthenon-hpc-lab/parthenon/blob/develop
   ``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 int he case where a single variable
+  reflecting boundary conditions in the case where a single variable
   is used as the state vector.
 
 

From 21c0c5d0a403b7d80caecef4ab40e0687a9231b0 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 27 Sep 2024 17:18:37 -0600
Subject: [PATCH 07/23] add error checking for swarm/mesh BC consistency

---
 src/mesh/mesh.cpp | 10 ++++++++++
 src/mesh/mesh.hpp |  4 ----
 2 files changed, 10 insertions(+), 4 deletions(-)

diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index 49b6d5a520ef..789c1723b25e 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -1189,6 +1189,16 @@ void Mesh::SetBCNames_(ParameterInput *pin) {
       pin->GetOrAddString("parthenon/swarm", "ox2_bc", mesh_bc_names[3]),
       pin->GetOrAddString("parthenon/swarm", "ix3_bc", mesh_bc_names[4]),
       pin->GetOrAddString("parthenon/swarm", "ox3_bc", 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 == "periodic",
+                        "If the mesh is periodic, swarms must be also.");
+    }
+  }
 }
 
 std::array<BoundaryFlag, BOUNDARY_NFACES>
diff --git a/src/mesh/mesh.hpp b/src/mesh/mesh.hpp
index 918d50c175cb..9f50e9a12b28 100644
--- a/src/mesh/mesh.hpp
+++ b/src/mesh/mesh.hpp
@@ -121,10 +121,6 @@ class Mesh {
   BValNames_t mesh_swarm_bc_names;
 
   // these are flags not boundary functions
-  // 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.
   std::array<BoundaryFlag, BOUNDARY_NFACES> mesh_bcs;
   int ndim; // number of dimensions
   const bool adaptive, multilevel, multigrid;

From 69b767e775cdb12d6a3510c912bf945255bda044 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 27 Sep 2024 17:19:11 -0600
Subject: [PATCH 08/23] typo

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

diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index 789c1723b25e..3f490702323d 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -1195,7 +1195,7 @@ void Mesh::SetBCNames_(ParameterInput *pin) {
   // boundaries are decoupled.
   for (int i = 0; i < BOUNDARY_NFACES; ++i) {
     if (mesh_bc_names[i] == "periodic") {
-      PARTHENON_REQUIRE(mesh_swarm_bc_names == "periodic",
+      PARTHENON_REQUIRE(mesh_swarm_bc_names[i] == "periodic",
                         "If the mesh is periodic, swarms must be also.");
     }
   }

From 5925d1d10168ef5d8c513052d08e13e81fd5d572 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 27 Sep 2024 18:01:56 -0600
Subject: [PATCH 09/23] phdf diff

---
 .../packages/parthenon_tools/parthenon_tools/phdf_diff.py  | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

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")

From 2bb344f06d35174688c33852b6fec7d65913a514 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Sat, 28 Sep 2024 10:55:57 -0600
Subject: [PATCH 10/23] Register reflecting BCs for advection examples

---
 example/advection/main.cpp        | 1 +
 example/sparse_advection/main.cpp | 1 +
 2 files changed, 2 insertions(+)

diff --git a/example/advection/main.cpp b/example/advection/main.cpp
index 146093bc22eb..43764c35f774 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;
+  mpan.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/sparse_advection/main.cpp b/example/sparse_advection/main.cpp
index 585d1a35ab1e..4e19a6ac8034 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;
+  mpan.app_input->RegisterDefaultReflectingBoundaryConditions();
 
   // call ParthenonInit to initialize MPI and Kokkos, parse the input deck, and set up
   auto manager_status = pman.ParthenonInitEnv(argc, argv);

From 9ca8a2e6c336652246aefe5e6c7d7c1d703c2650 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Sat, 28 Sep 2024 13:03:18 -0600
Subject: [PATCH 11/23] typo

---
 example/advection/main.cpp        | 2 +-
 example/sparse_advection/main.cpp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/example/advection/main.cpp b/example/advection/main.cpp
index 43764c35f774..5ed943b07474 100644
--- a/example/advection/main.cpp
+++ b/example/advection/main.cpp
@@ -25,7 +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;
-  mpan.app_input->RegisterDefaultReflectingBoundaryConditions();
+  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/sparse_advection/main.cpp b/example/sparse_advection/main.cpp
index 4e19a6ac8034..ce6b197e335a 100644
--- a/example/sparse_advection/main.cpp
+++ b/example/sparse_advection/main.cpp
@@ -27,7 +27,7 @@ int main(int argc, char *argv[]) {
   pman.app_input->ProblemGenerator = sparse_advection_example::ProblemGenerator;
   pman.app_input->PostStepDiagnosticsInLoop =
       sparse_advection_example::PostStepDiagnosticsInLoop;
-  mpan.app_input->RegisterDefaultReflectingBoundaryConditions();
+  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);

From 33f63c1c557c4a61252f1994605b247cbd512f3c Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Sat, 28 Sep 2024 15:03:01 -0600
Subject: [PATCH 12/23] silly backwards compatibility thing to make it so you
 don't have to specify swarm BCs if you're not using swarms

---
 src/mesh/mesh.cpp | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/mesh/mesh.cpp b/src/mesh/mesh.cpp
index 3f490702323d..d71c1b50dde8 100644
--- a/src/mesh/mesh.cpp
+++ b/src/mesh/mesh.cpp
@@ -1182,13 +1182,18 @@ void Mesh::SetBCNames_(ParameterInput *pin) {
                    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", mesh_bc_names[0]),
-      pin->GetOrAddString("parthenon/swarm", "ox1_bc", mesh_bc_names[1]),
-      pin->GetOrAddString("parthenon/swarm", "ix2_bc", mesh_bc_names[2]),
-      pin->GetOrAddString("parthenon/swarm", "ox2_bc", mesh_bc_names[3]),
-      pin->GetOrAddString("parthenon/swarm", "ix3_bc", mesh_bc_names[4]),
-      pin->GetOrAddString("parthenon/swarm", "ox3_bc", mesh_bc_names[5])};
+      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

From b10d7b103f57451ace50f405f648da7c530d6074 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Wed, 16 Oct 2024 01:30:55 +0200
Subject: [PATCH 13/23] Address CUDA MPI/ICP issue with Kokkos <=4.4.1 (#1189)

---
 CHANGELOG.md   | 3 ++-
 CMakeLists.txt | 7 +++++++
 2 files changed, 9 insertions(+), 1 deletion(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 326c1e6c55e6..1f046faabe75 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -12,11 +12,12 @@
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR 1187]](https://github.com/parthenon-hpc-lab/parthenon/pull/1187) Make DataCollection::Add safer and generalize MeshBlockData::Initialize
-- [[PR 1186]](https://github.com/parthenon-hpc-lab/parthenon/pull/1186) Bump Kokkos submodule to 4.4.1
+- [[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 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.
 
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 24a7a3b2ebfb..61fd4fd3cfff 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -291,6 +291,13 @@ if (Kokkos_ENABLE_CUDA)
   if(CHECK_REGISTRY_PRESSURE)
     add_compile_options(-Xptxas=-v)
   endif()
+
+  # Async malloc are enabled by default until Kokkos <= 4.4.1 but
+  # cause issues with IPCs and MPI,
+  # see https://github.com/parthenon-hpc-lab/athenapk/pull/114#issuecomment-2358981857
+  # Given that it's also not expected to be a significant performance gain
+  # it's hard disabled now.
+  set(Kokkos_ENABLE_IMPL_CUDA_MALLOC_ASYNC OFF CACHE BOOL "Disable Cuda async malloc (to address MPI/IPC issues)")
 endif()
 # Note that these options may not play nice with nvcc wrapper
 if (TEST_INTEL_OPTIMIZATION)

From 6a0097a519c235d9487a1af2f0050a79fd8a804f Mon Sep 17 00:00:00 2001
From: Adam Dempsey <adempsey@lanl.gov>
Date: Thu, 31 Oct 2024 07:38:34 -0600
Subject: [PATCH 14/23] Remove else from if constexpr when there are returns

---
 src/utils/hash.hpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

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>

From f397d14c3f80481c6cd6452748148f52c9bf1de6 Mon Sep 17 00:00:00 2001
From: Ben Ryan <brryan@lanl.gov>
Date: Thu, 31 Oct 2024 10:54:09 -0600
Subject: [PATCH 15/23] Jonah's fix for this CI issue

---
 .github/workflows/docs.yml | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

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}"

From 9ec0cf8fc6c485da6bb404a7a7d6879333530a84 Mon Sep 17 00:00:00 2001
From: Ben Ryan <brryan@lanl.gov>
Date: Thu, 31 Oct 2024 10:55:43 -0600
Subject: [PATCH 16/23] CHANGELOG

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

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 7405fd170490..bbd37785d197 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@
 - [[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/...)
+- [[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

From b5364b7e2777fff9850b1c2823ccd136d0ad4c0b Mon Sep 17 00:00:00 2001
From: Alex Long <along@lanl.gov>
Date: Thu, 31 Oct 2024 22:35:00 -0600
Subject: [PATCH 17/23] Consolidate buffer packing functions with less atomics
 (#1199)

* Consolidate buffer packing functions with less atomics

+ Fold in buffer sizing to the load buffer function
+ Use a sort and a discontinuity check instead of an
atomic_fetch_add inside of LoadBuffer_ to get particle
index into buffer
+ Reduce transcendental functions call in particle sourcing in
particle example

* Address PR comments

+ Call the member variable in SwarmKey the sort_key
+ Remove CountParticlesInBuffer function
+ Add buffer_start and buffer_sorted as swarm member variables

* Update example/particles/particles.cpp

Co-authored-by: Ben Ryan <bryan10@illinois.edu>

* Update src/interface/swarm_comms.cpp

Co-authored-by: Ben Ryan <bryan10@illinois.edu>

---------

Co-authored-by: Ben Ryan <bryan10@illinois.edu>
Co-authored-by: Ben Ryan <brryan@lanl.gov>
---
 example/particles/particles.cpp        |  30 ++--
 src/interface/swarm.cpp                |  23 +--
 src/interface/swarm.hpp                |   5 +-
 src/interface/swarm_comms.cpp          | 190 ++++++++++++-------------
 src/interface/swarm_device_context.hpp |   7 +-
 5 files changed, 131 insertions(+), 124 deletions(-)

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/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..7b236b6e4e1f 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,98 @@ 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();
+  }
 }
 
 void Swarm::Send(BoundaryCommSubset phase) {
@@ -279,10 +277,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_;

From ddc6500d438fa8adbd522cdb2c3a2f0b4f0c2525 Mon Sep 17 00:00:00 2001
From: Philipp Grete <pgrete@hs.uni-hamburg.de>
Date: Mon, 4 Nov 2024 18:55:01 +0100
Subject: [PATCH 18/23] [Trivial] Fix type used for array init (#1170)

* Fix type used for array init

* init array with constexpr expression

* CC

---------

Co-authored-by: Jonah Miller <jonahm@lanl.gov>
---
 CHANGELOG.md                 | 1 +
 src/outputs/output_utils.hpp | 2 +-
 src/outputs/restart_hdf5.cpp | 4 ++--
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index bbd37785d197..17f8a009242d 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -19,6 +19,7 @@
 - [[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.
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/restart_hdf5.cpp b/src/outputs/restart_hdf5.cpp
index cababddbfafe..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
@@ -213,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;

From f2171ff36d9f3c627d9702b794a965963d23e790 Mon Sep 17 00:00:00 2001
From: Ben Ryan <brryan@lanl.gov>
Date: Thu, 7 Nov 2024 20:43:27 -0700
Subject: [PATCH 19/23] Leapfrog fix (#1206)

* Missing send size init

* cleanup, CHANGELOG

* verbose CI

* further CI debugging

* This should be working...

* This should be fixed... but I get a segfault on GPU

* Is it my AMD GPU thats wrong?

* Missing a return statement

* retest

* Oops missing statement

* Revert test

* revert workflow
---
 CHANGELOG.md                  |  1 +
 src/interface/swarm_comms.cpp |  5 +++++
 src/utils/sort.hpp            | 17 +++++++++++++++--
 3 files changed, 21 insertions(+), 2 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 17f8a009242d..c29941247e43 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -11,6 +11,7 @@
 - [[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
diff --git a/src/interface/swarm_comms.cpp b/src/interface/swarm_comms.cpp
index 7b236b6e4e1f..054ded34fabb 100644
--- a/src/interface/swarm_comms.cpp
+++ b/src/interface/swarm_comms.cpp
@@ -269,6 +269,11 @@ void Swarm::LoadBuffers_() {
 
     // 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;
+    }
   }
 }
 
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);

From 64473d70df5461fcf10ab52fd603aad310cb46b1 Mon Sep 17 00:00:00 2001
From: Ben Ryan <brryan@lanl.gov>
Date: Fri, 8 Nov 2024 09:55:59 -0700
Subject: [PATCH 20/23] CHANGELOG

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

diff --git a/CHANGELOG.md b/CHANGELOG.md
index c29941247e43..554debf3c568 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1103]](https://github.com/parthenon-hpc-lab/parthenon/pull/1103) Add sparsity to vector wave equation test
 - [[PR 1185]](https://github.com/parthenon-hpc-lab/parthenon/pull/1185/files) Bugfix to particle defragmentation
 - [[PR 1184]](https://github.com/parthenon-hpc-lab/parthenon/pull/1184) Fix swarm block neighbor indexing in 1D, 2D
 - [[PR 1183]](https://github.com/parthenon-hpc-lab/parthenon/pull/1183) Fix particle leapfrog example initialization data

From 0f6cec4e113cc0c8241497521e633fadf3734791 Mon Sep 17 00:00:00 2001
From: Ben Ryan <brryan@lanl.gov>
Date: Fri, 8 Nov 2024 09:57:08 -0700
Subject: [PATCH 21/23] copyright

---
 example/fine_advection/parthenon_app_inputs.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/example/fine_advection/parthenon_app_inputs.cpp b/example/fine_advection/parthenon_app_inputs.cpp
index 9f8c7f80e5c6..35236015e391 100644
--- a/example/fine_advection/parthenon_app_inputs.cpp
+++ b/example/fine_advection/parthenon_app_inputs.cpp
@@ -1,5 +1,5 @@
 //========================================================================================
-// (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

From 5272708ee9ea031455c4e742bd2962d956e924ce Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 12 Nov 2024 13:19:16 -0700
Subject: [PATCH 22/23] Philipp comments

---
 example/fine_advection/advection_package.cpp    |  4 ++--
 example/fine_advection/advection_package.hpp    |  9 +++++----
 example/fine_advection/parthenon_app_inputs.cpp | 12 ++++++------
 3 files changed, 13 insertions(+), 12 deletions(-)

diff --git a/example/fine_advection/advection_package.cpp b/example/fine_advection/advection_package.cpp
index 648dfe4ac83f..428d91b9245f 100644
--- a/example/fine_advection/advection_package.cpp
+++ b/example/fine_advection/advection_package.cpp
@@ -225,9 +225,9 @@ TaskStatus FillDerived(MeshData<Real> *md) {
           const int fi = (ndim > 0) ? (i - nghost) * 2 + nghost : i;
           pack(b, Conserved::phi_fine_restricted(), k, j, i) = 0.0;
           Real ntot = 0.0;
-          for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+          for (int koff = 0; koff <= (ndim > 2); ++koff)
             for (int joff = 0; joff <= (ndim > 1); ++joff)
-              for (int koff = 0; koff <= (ndim > 2); ++koff) {
+              for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                 ntot += 1.0;
                 pack(b, Conserved::phi_fine_restricted(), k, j, i) +=
                     pack(b, Conserved::phi_fine(), kf + koff, jf + joff, fi + ioff);
diff --git a/example/fine_advection/advection_package.hpp b/example/fine_advection/advection_package.hpp
index 6e3f8e03f232..12b6b9b65aad 100644
--- a/example/fine_advection/advection_package.hpp
+++ b/example/fine_advection/advection_package.hpp
@@ -80,11 +80,12 @@ TaskStatus CalculateFluxes(pack_desc_t &desc, parthenon::TopologicalElement FACE
   IndexRange ib = md->GetBoundsI(cl, IndexDomain::interior, FACE);
   IndexRange jb = md->GetBoundsJ(cl, IndexDomain::interior, FACE);
   IndexRange kb = md->GetBoundsK(cl, IndexDomain::interior, FACE);
-  constexpr int scratch_size = 0;
-  constexpr int scratch_level = 1;
+  constexpr int unused_scratch_size = 0;
+  constexpr int unused_scratch_level = 1;
   parthenon::par_for_outer(
-      PARTHENON_AUTO_LABEL, scratch_size, scratch_level, 0, pack.GetNBlocks() - 1, kb.s,
-      kb.e, KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
+      PARTHENON_AUTO_LABEL, unused_scratch_size, unused_scratch_level, 0,
+      pack.GetNBlocks() - 1, kb.s, kb.e,
+      KOKKOS_LAMBDA(parthenon::team_mbr_t member, const int b, const int k) {
         parthenon::Indexer2D idxer({jb.s, jb.e}, {ib.s, ib.e});
         for (int l = pack.GetLowerBound(b); l <= pack.GetUpperBound(b); ++l) {
           parthenon::par_for_inner(member, 0, idxer.size() - 1, [&](const int idx) {
diff --git a/example/fine_advection/parthenon_app_inputs.cpp b/example/fine_advection/parthenon_app_inputs.cpp
index 35236015e391..8025970f35ca 100644
--- a/example/fine_advection/parthenon_app_inputs.cpp
+++ b/example/fine_advection/parthenon_app_inputs.cpp
@@ -119,9 +119,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
             Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) +
                        coords.Xc<2>(j) * coords.Xc<2>(j) +
                        coords.Xc<3>(k) * coords.Xc<3>(k);
-            for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+            for (int koff = 0; koff <= (ndim > 2); ++koff)
               for (int joff = 0; joff <= (ndim > 1); ++joff)
-                for (int koff = 0; koff <= (ndim > 2); ++koff) {
+                for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                   pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) =
                       1. + amp * exp(-100.0 * rsq);
                 }
@@ -129,16 +129,16 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
             Real rsq = coords.Xc<1>(i) * coords.Xc<1>(i) +
                        coords.Xc<2>(j) * coords.Xc<2>(j) +
                        coords.Xc<3>(k) * coords.Xc<3>(k);
-            for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+            for (int koff = 0; koff <= (ndim > 2); ++koff)
               for (int joff = 0; joff <= (ndim > 1); ++joff)
-                for (int koff = 0; koff <= (ndim > 2); ++koff) {
+                for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                   pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) =
                       (rsq < 0.15 * 0.15 ? 1.0 : 0.0);
                 }
           } else {
-            for (int ioff = 0; ioff <= (ndim > 0); ++ioff)
+            for (int koff = 0; koff <= (ndim > 2); ++koff)
               for (int joff = 0; joff <= (ndim > 1); ++joff)
-                for (int koff = 0; koff <= (ndim > 2); ++koff) {
+                for (int ioff = 0; ioff <= (ndim > 0); ++ioff) {
                   pack(b, phi_fine(), kf + koff, jf + joff, fi + ioff) = 0.0;
                 }
           }

From 84fc520f7360ff41402d10363dbd3d951da35a57 Mon Sep 17 00:00:00 2001
From: Luke Roberts <lfroberts@lanl.gov>
Date: Tue, 12 Nov 2024 13:22:47 -0700
Subject: [PATCH 23/23] more Philipp comments

---
 example/fine_advection/parthinput.advection | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/example/fine_advection/parthinput.advection b/example/fine_advection/parthinput.advection
index 31ee255c5053..b3bf7ca10d0c 100644
--- a/example/fine_advection/parthinput.advection
+++ b/example/fine_advection/parthinput.advection
@@ -57,6 +57,13 @@ profile = hard_sphere
 refine_tol = 0.3    # control the package specific refinement tagging function
 derefine_tol = 0.03
 
+# Include advection equation for regular resolution CC fields
+do_regular_advection = true
+# Include advection equation for fine resolution CC fields
+do_fine_advection = true
+# Include vector wave equation using CT
+do_CT_advection = true
+
 <parthenon/output1>
 file_type = rst
 dt = 0.05