Add SwarmPacks

example/particle_leapfrog/particle_leapfrog.cpp

@@ -79,7 +79,7 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 
 Real EstimateTimestepBlock(MeshBlockData<Real> *rc) {
   auto pmb = rc->GetBlockPointer();
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
   auto pkg = pmb->packages.Get("particles_package");
   const auto &cfl = pkg->Param<Real>("cfl");
 
@@ -135,7 +135,7 @@ const Kokkos::Array<Kokkos::Array<Real, 6>, num_test_particles> particles_ic = {
 
 void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto pkg = pmb->packages.Get("particles_package");
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
 
   const IndexRange &ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
   const IndexRange &jb = pmb->cellbounds.GetBoundsJ(IndexDomain::interior);
@@ -178,9 +178,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto new_particles_context = swarm->AddEmptyParticles(num_particles_this_block);
 
   auto &id = swarm->Get<int>("id").Get();
-  auto &x = swarm->Get<Real>("x").Get();
-  auto &y = swarm->Get<Real>("y").Get();
-  auto &z = swarm->Get<Real>("z").Get();
+  auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   auto &v = swarm->Get<Real>("v").Get();
   auto &vv = swarm->Get<Real>("vv").Get();
 
@@ -208,17 +208,17 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
 }
 
 TaskStatus TransportParticles(MeshBlock *pmb, const StagedIntegrator *integrator) {
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
   auto pkg = pmb->packages.Get("particles_package");
 
   int max_active_index = swarm->GetMaxActiveIndex();
 
   Real dt = integrator->dt;
 
   auto &id = swarm->Get<int>("id").Get();
-  auto &x = swarm->Get<Real>("x").Get();
-  auto &y = swarm->Get<Real>("y").Get();
-  auto &z = swarm->Get<Real>("z").Get();
+  auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   auto &v = swarm->Get<Real>("v").Get();
 
   auto swarm_d = swarm->GetDeviceContext();
@@ -280,7 +280,7 @@ TaskCollection ParticleDriver::MakeParticlesUpdateTaskCollection() const {
     for (int i = 0; i < blocks.size(); i++) {
       auto &tl = sync_region0[0];
       auto &pmb = blocks[i];
-      auto &sc = pmb->swarm_data.Get();
+      auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
       auto reset_comms = tl.AddTask(none, &SwarmContainer::ResetCommunication, sc.get());
     }
   }
@@ -289,7 +289,7 @@ TaskCollection ParticleDriver::MakeParticlesUpdateTaskCollection() const {
   for (int i = 0; i < blocks.size(); i++) {
     auto &pmb = blocks[i];
 
-    auto &sc = pmb->swarm_data.Get();
+    auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
 
     auto &tl = async_region0[i];
 
@@ -315,7 +315,7 @@ TaskCollection ParticleDriver::MakeFinalizationTaskCollection() const {
 
   for (int i = 0; i < blocks.size(); i++) {
     auto &pmb = blocks[i];
-    auto &sc = pmb->swarm_data.Get();
+    auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
     auto &tl = async_region1[i];
 
     // Defragment if swarm memory pool occupancy is 90%

example/particle_tracers/particle_tracers.cpp

@@ -165,16 +165,16 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
 } // namespace particles_package
 
 TaskStatus AdvectTracers(MeshBlock *pmb, const StagedIntegrator *integrator) {
-  auto swarm = pmb->swarm_data.Get()->Get("tracers");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("tracers");
   auto adv_pkg = pmb->packages.Get("advection_package");
 
   int max_active_index = swarm->GetMaxActiveIndex();
 
   Real dt = integrator->dt;
 
-  auto &x = swarm->Get<Real>("x").Get();
-  auto &y = swarm->Get<Real>("y").Get();
-  auto &z = swarm->Get<Real>("z").Get();
+  auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
 
   const auto &vx = adv_pkg->Param<Real>("vx");
   const auto &vy = adv_pkg->Param<Real>("vy");
@@ -197,7 +197,7 @@ TaskStatus AdvectTracers(MeshBlock *pmb, const StagedIntegrator *integrator) {
 }
 
 TaskStatus DepositTracers(MeshBlock *pmb) {
-  auto swarm = pmb->swarm_data.Get()->Get("tracers");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("tracers");
 
   // Meshblock geometry
   const IndexRange &ib = pmb->cellbounds.GetBoundsI(IndexDomain::interior);
@@ -211,9 +211,9 @@ TaskStatus DepositTracers(MeshBlock *pmb) {
   const Real &minx_j = pmb->coords.Xf<2>(jb.s);
   const Real &minx_k = pmb->coords.Xf<3>(kb.s);
 
-  const auto &x = swarm->Get<Real>("x").Get();
-  const auto &y = swarm->Get<Real>("y").Get();
-  const auto &z = swarm->Get<Real>("z").Get();
+  const auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  const auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  const auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   auto swarm_d = swarm->GetDeviceContext();
 
   auto &tracer_dep = pmb->meshblock_data.Get()->Get("tracer_deposition").data;
@@ -321,7 +321,7 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto &tr_pkg = pmb->packages.Get("particles_package");
   auto &mbd = pmb->meshblock_data.Get();
   auto &advected = mbd->Get("advected").data;
-  auto &swarm = pmb->swarm_data.Get()->Get("tracers");
+  auto &swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("tracers");
   const auto num_tracers = tr_pkg->Param<int>("num_tracers");
   auto rng_pool = tr_pkg->Param<RNGPool>("rng_pool");
 
@@ -377,9 +377,9 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
 
   auto new_particles_context = swarm->AddEmptyParticles(num_tracers_meshblock);
 
-  auto &x = swarm->Get<Real>("x").Get();
-  auto &y = swarm->Get<Real>("y").Get();
-  auto &z = swarm->Get<Real>("z").Get();
+  auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   auto &id = swarm->Get<int>("id").Get();
 
   auto swarm_d = swarm->GetDeviceContext();
@@ -493,7 +493,7 @@ TaskCollection ParticleDriver::MakeTaskCollection(BlockList_t &blocks, int stage
       for (int i = 0; i < blocks.size(); i++) {
         auto &tl = sync_region0[0];
         auto &pmb = blocks[i];
-        auto &sc = pmb->swarm_data.Get();
+        auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
         auto reset_comms =
             tl.AddTask(none, &SwarmContainer::ResetCommunication, sc.get());
       }
@@ -503,7 +503,7 @@ TaskCollection ParticleDriver::MakeTaskCollection(BlockList_t &blocks, int stage
     for (int n = 0; n < nblocks; n++) {
       auto &tl = async_region1[n];
       auto &pmb = blocks[n];
-      auto &sc = pmb->swarm_data.Get();
+      auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
       auto tracerAdvect =
           tl.AddTask(none, tracers_example::AdvectTracers, pmb.get(), integrator.get());
 

example/particles/particles.cpp

@@ -133,7 +133,7 @@ TaskStatus DestroySomeParticles(MeshBlock *pmb) {
   PARTHENON_INSTRUMENT
 
   auto pkg = pmb->packages.Get("particles_package");
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
   auto rng_pool = pkg->Param<RNGPool>("rng_pool");
   const auto destroy_particles_frac = pkg->Param<Real>("destroy_particles_frac");
 
@@ -163,7 +163,7 @@ TaskStatus SortParticlesIfUsingPerCellDeposition(MeshBlock *pmb) {
   auto pkg = pmb->packages.Get("particles_package");
   const auto deposition_method = pkg->Param<DepositionMethod>("deposition_method");
   if (deposition_method == DepositionMethod::per_cell) {
-    auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+    auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
     swarm->SortParticlesByCell();
   }
 
@@ -173,7 +173,7 @@ TaskStatus SortParticlesIfUsingPerCellDeposition(MeshBlock *pmb) {
 TaskStatus DepositParticles(MeshBlock *pmb) {
   PARTHENON_INSTRUMENT
 
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
 
   auto pkg = pmb->packages.Get("particles_package");
   const auto deposition_method = pkg->Param<DepositionMethod>("deposition_method");
@@ -189,9 +189,9 @@ TaskStatus DepositParticles(MeshBlock *pmb) {
   const Real &minx_j = pmb->coords.Xf<2>(jb.s);
   const Real &minx_k = pmb->coords.Xf<3>(kb.s);
 
-  const auto &x = swarm->Get<Real>("x").Get();
-  const auto &y = swarm->Get<Real>("y").Get();
-  const auto &z = swarm->Get<Real>("z").Get();
+  const auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  const auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  const auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   const auto &weight = swarm->Get<Real>("weight").Get();
   auto swarm_d = swarm->GetDeviceContext();
 
@@ -244,7 +244,7 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
   PARTHENON_INSTRUMENT
 
   auto pkg = pmb->packages.Get("particles_package");
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
   auto rng_pool = pkg->Param<RNGPool>("rng_pool");
   auto num_particles = pkg->Param<int>("num_particles");
   auto vel = pkg->Param<Real>("particle_speed");
@@ -268,9 +268,9 @@ TaskStatus CreateSomeParticles(MeshBlock *pmb, const double t0) {
   const Real &minx_k = pmb->coords.Xf<3>(kb.s);
 
   auto &t = swarm->Get<Real>("t").Get();
-  auto &x = swarm->Get<Real>("x").Get();
-  auto &y = swarm->Get<Real>("y").Get();
-  auto &z = swarm->Get<Real>("z").Get();
+  auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   auto &v = swarm->Get<Real>("v").Get();
   auto &weight = swarm->Get<Real>("weight").Get();
 
@@ -355,7 +355,7 @@ TaskStatus TransportParticles(MeshBlock *pmb, const StagedIntegrator *integrator
                               const double t0) {
   PARTHENON_INSTRUMENT
 
-  auto swarm = pmb->swarm_data.Get()->Get("my_particles");
+  auto swarm = pmb->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
   auto pkg = pmb->packages.Get("particles_package");
   const auto orbiting_particles = pkg->Param<bool>("orbiting_particles");
 
@@ -364,9 +364,9 @@ TaskStatus TransportParticles(MeshBlock *pmb, const StagedIntegrator *integrator
   Real dt = integrator->dt;
 
   auto &t = swarm->Get<Real>("t").Get();
-  auto &x = swarm->Get<Real>("x").Get();
-  auto &y = swarm->Get<Real>("y").Get();
-  auto &z = swarm->Get<Real>("z").Get();
+  auto &x = swarm->Get<Real>(swarm_position::x::name()).Get();
+  auto &y = swarm->Get<Real>(swarm_position::y::name()).Get();
+  auto &z = swarm->Get<Real>(swarm_position::z::name()).Get();
   auto &v = swarm->Get<Real>("v").Get();
 
   const Real &dx_i = pmb->coords.Dxf<1>(pmb->cellbounds.is(IndexDomain::interior));
@@ -499,7 +499,7 @@ TaskListStatus ParticleDriver::Step() {
     for (auto &block : blocks) {
       // TODO(BRR) Despite this "my_particles"-specific call, this function feels like it
       // should be generalized
-      auto swarm = block->swarm_data.Get()->Get("my_particles");
+      auto swarm = block->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
       if (!swarm->finished_transport) {
         particles_update_done = false;
       }
@@ -519,7 +519,7 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
 
   int num_sent_local = 0;
   for (auto &block : blocks) {
-    auto sc = block->swarm_data.Get();
+    auto sc = block->meshblock_data.Get()->swarm_data.Get();
     auto swarm = sc->Get("my_particles");
     swarm->finished_transport = false;
     num_sent_local += swarm->num_particles_sent_;
@@ -528,7 +528,7 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
   int num_sent_global = num_sent_local; // potentially overwritten by following Allreduce
 #ifdef MPI_PARALLEL
   for (auto &block : blocks) {
-    auto swarm = block->swarm_data.Get()->Get("my_particles");
+    auto swarm = block->meshblock_data.Get()->swarm_data.Get()->Get("my_particles");
     for (int n = 0; n < block->pbval->nneighbor; n++) {
       NeighborBlock &nb = block->pbval->neighbor[n];
       // TODO(BRR) May want logic like this if we have non-blocking TaskRegions
@@ -552,7 +552,7 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
   if (num_sent_global == 0) {
     for (auto &block : blocks) {
       auto &pmb = block;
-      auto sc = pmb->swarm_data.Get();
+      auto sc = pmb->meshblock_data.Get()->swarm_data.Get();
       auto swarm = sc->Get("my_particles");
       swarm->finished_transport = true;
     }
@@ -561,7 +561,7 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
   // Reset boundary statuses
   for (auto &block : blocks) {
     auto &pmb = block;
-    auto sc = pmb->swarm_data.Get();
+    auto sc = pmb->meshblock_data.Get()->swarm_data.Get();
     auto swarm = sc->Get("my_particles");
     for (int n = 0; n < swarm->vbswarm->bd_var_.nbmax; n++) {
       auto &nb = pmb->pbval->neighbor[n];
@@ -601,7 +601,7 @@ TaskCollection ParticleDriver::MakeParticlesUpdateTaskCollection() const {
   for (int i = 0; i < blocks.size(); i++) {
     auto &pmb = blocks[i];
 
-    auto &sc = pmb->swarm_data.Get();
+    auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
 
     auto &tl = async_region0[i];
 
@@ -633,7 +633,7 @@ TaskCollection ParticleDriver::MakeFinalizationTaskCollection() const {
 
   for (int i = 0; i < blocks.size(); i++) {
     auto &pmb = blocks[i];
-    auto &sc = pmb->swarm_data.Get();
+    auto &sc = pmb->meshblock_data.Get()->swarm_data.Get();
     auto &sc1 = pmb->meshblock_data.Get();
     auto &tl = async_region1[i];
 

scripts/python/packages/parthenon_tools/parthenon_tools/phdf.py

@@ -59,15 +59,15 @@ def Get(self, variable):
 
     @property
     def x(self):
-        return self.Get("x")
+        return self.Get("swarm.x")
 
     @property
     def y(self):
-        return self.Get("y")
+        return self.Get("swarm.y")
 
     @property
     def z(self):
-        return self.Get("z")
+        return self.Get("swarm.z")
 
     def __getitem__(self, key):
         return self.Get(key)

src/CMakeLists.txt

@@ -101,10 +101,10 @@ add_library(parthenon
   bvals/comms/bnd_info.cpp
   bvals/comms/bnd_info.hpp
   bvals/comms/boundary_communication.cpp
-  bvals/comms/flux_correction.cpp 
-  bvals/comms/tag_map.cpp 
-  bvals/comms/tag_map.hpp 
-  
+  bvals/comms/flux_correction.cpp
+  bvals/comms/tag_map.cpp
+  bvals/comms/tag_map.hpp
+
   bvals/boundary_conditions_generic.hpp
   bvals/boundary_conditions.cpp
   bvals/boundary_conditions.hpp
@@ -131,13 +131,15 @@ add_library(parthenon
   interface/mesh_data.hpp
   interface/meshblock_data.cpp
   interface/meshblock_data.hpp
-  interface/swarm_comms.cpp
   interface/swarm_container.cpp
   interface/swarm.cpp
   interface/swarm.hpp
   interface/swarm_boundaries.hpp
+  interface/swarm_comms.cpp
   interface/swarm_device_context.hpp
+  interface/swarm_variable_types.hpp
   interface/make_pack_descriptor.hpp
+  interface/make_swarm_pack_descriptor.hpp
   interface/metadata.cpp
   interface/metadata.hpp
   interface/packages.hpp
@@ -150,6 +152,8 @@ add_library(parthenon
   interface/sparse_pool.hpp
   interface/state_descriptor.hpp
   interface/state_descriptor.cpp
+  interface/swarm_pack_base.hpp
+  interface/swarm_pack.hpp
   interface/update.cpp
   interface/update.hpp
   interface/var_id.hpp