Particle BCs without relocatable device code

CHANGELOG.md

@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
 - [[PR 1037]](https://github.com/parthenon-hpc-lab/parthenon/pull/1037) Add SwarmPacks
 - [[PR 1068]](https://github.com/parthenon-hpc-lab/parthenon/pull/1068) Add ability to dump sparse pack contents as a string
 - [[PR 1062]](https://github.com/parthenon-hpc-lab/parthenon/pull/1062) UserWorkBeforeRestartOutput #1062
@@ -37,6 +38,7 @@
 
 
 ### Incompatibilities (i.e. breaking changes)
+- [[PR 1026]](https://github.com/parthenon-hpc-lab/parthenon/pull/1026) Particle BCs without relocatable device code
 - [[PR 1037]](https://github.com/parthenon-hpc-lab/parthenon/pull/1037) Add SwarmPacks
 - [[PR 1042]](https://github.com/parthenon-hpc-lab/parthenon/pull/1042) Use Offset class and clean up of NeighborBlock
 - [[PR 1019]](https://github.com/parthenon-hpc-lab/parthenon/pull/1019) Remove support for file formats < 3

doc/sphinx/src/particles.rst

@@ -190,17 +190,23 @@ For example ``SwarmPack`` usage, see the ``particle_leapfrog`` example.
 Boundary conditions
 -------------------
 
-Particle boundary conditions are not applied in separate kernel calls;
-instead, inherited classes containing boundary condition functions for
-updating particles or removing them when they are in boundary regions
-are allocated depending on the boundary flags specified in the input
-file. Currently, outflow and periodic boundaries are supported natively.
-User-specified boundary conditions must be set by specifying the “user”
-flag in the input parameter file and then updating the appropriate
-Swarm::bounds array entries to factory functions that allocate
-device-side boundary condition objects. An example is given in the
-``particles`` example when ix1 and ox1 are set to ``user`` in the input
-parameter file.
+Particle boundary conditions are applied in per-block per-boundary kernel
+launches analogous to grid-based variables. Outflow and periodic boundaries
+are supported natively, but other boundary conditions (including reflecting)
+must be provided by the downstream application. Particle boundary conditions are
+enrolled by setting entries in ``ApplicationInput::swarm_boundary_conditions``
+to per-boundary (inner ``x1``, outer ``x2``, etc.) custom boundary functions
+with signature
+
+.. code:: cpp
+
+   void SwarmUserInnerX1(std::shared_ptr<Swarm> &swarm);
+
+The ``particles`` example demonstrates how to create and enroll custom particle
+boundary conditions.
+
+Note that periodic boundary conditions cannot be enrolled by the user; the
+default ``periodic`` option for Parthenon must be requested in the input file.
 
 Outputs
 --------

example/particle_leapfrog/particle_leapfrog.cpp

@@ -262,11 +262,6 @@ TaskStatus TransportParticles(MeshData<Real> *md, const StagedIntegrator *integr
           pack_pos(b, swarm_position::x(), n) += pack_v(b, iv + 0, n) * 0.5 * dt;
           pack_pos(b, swarm_position::y(), n) += pack_v(b, iv + 1, n) * 0.5 * dt;
           pack_pos(b, swarm_position::z(), n) += pack_v(b, iv + 2, n) * 0.5 * dt;
-
-          bool on_current_mesh_block;
-          swarm_d.GetNeighborBlockIndex(
-              n, pack_pos(b, swarm_position::x(), n), pack_pos(b, swarm_position::y(), n),
-              pack_pos(b, swarm_position::z(), n), on_current_mesh_block);
         }
       });
 
@@ -276,10 +271,9 @@ TaskStatus TransportParticles(MeshData<Real> *md, const StagedIntegrator *integr
 // Custom step function to allow for looping over MPI-related tasks until complete
 TaskListStatus ParticleDriver::Step() {
   TaskListStatus status;
-  integrator.dt = tm.dt;
 
-  BlockList_t &blocks = pmesh->block_list;
-  auto num_task_lists_executed_independently = blocks.size();
+  // Enforce fixed timestep
+  integrator.dt = tm.dt;
 
   status = MakeParticlesUpdateTaskCollection().Execute();
 

example/particle_tracers/particle_tracers.cpp

@@ -144,12 +144,6 @@ std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin) {
   int num_tracers = pin->GetOrAddReal("Tracers", "num_tracers", 100);
   pkg->AddParam<>("num_tracers", num_tracers);
 
-  // Initialize random number generator pool
-  int rng_seed = pin->GetOrAddInteger("Tracers", "rng_seed", 1273);
-  pkg->AddParam<>("rng_seed", rng_seed);
-  RNGPool rng_pool(rng_seed);
-  pkg->AddParam<>("rng_pool", rng_pool);
-
   // Add swarm of tracer particles
   std::string swarm_name = "tracers";
   Metadata swarm_metadata({Metadata::Provides, Metadata::None});
@@ -187,9 +181,6 @@ TaskStatus AdvectTracers(MeshBlock *pmb, const StagedIntegrator *integrator) {
           x(n) += vx * dt;
           y(n) += vy * dt;
           z(n) += vz * dt;
-
-          bool on_current_mesh_block = true;
-          swarm_d.GetNeighborBlockIndex(n, x(n), y(n), z(n), on_current_mesh_block);
         }
       });
 
@@ -237,7 +228,8 @@ TaskStatus DepositTracers(MeshBlock *pmb) {
             k = static_cast<int>(std::floor((z(n) - minx_k) / dx_k) + kb.s);
           }
 
-          // For testing in this example we make sure the indices are correct
+          // For testing in this example we make sure the indices are correct; these could
+          // be demoted to Debug-only calls
           if (i >= ib.s && i <= ib.e && j >= jb.s && j <= jb.e && k >= kb.s &&
               k <= kb.e) {
             Kokkos::atomic_add(&tracer_dep(k, j, i), 1.0);
@@ -322,7 +314,8 @@ void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
   auto &advected = mbd->Get("advected").data;
   auto &swarm = pmb->meshblock_data.Get()->GetSwarmData()->Get("tracers");
   const auto num_tracers = tr_pkg->Param<int>("num_tracers");
-  auto rng_pool = tr_pkg->Param<RNGPool>("rng_pool");
+  auto rng_pool =
+      RNGPool(pmb->gid); // Seed is meshblock gid for consistency across MPI decomposition
 
   const int ndim = pmb->pmy_mesh->ndim;
   PARTHENON_REQUIRE(ndim <= 2, "Tracer particles example only supports <= 2D!");
@@ -515,7 +508,7 @@ TaskCollection ParticleDriver::MakeTaskCollection(BlockList_t &blocks, int stage
       auto deposit = tl.AddTask(receive, tracers_example::DepositTracers, pmb.get());
 
       // Defragment if swarm memory pool occupancy is 90%
-      auto defrag = tl.AddTask(none, &SwarmContainer::Defrag, sc.get(), 0.9);
+      auto defrag = tl.AddTask(deposit, &SwarmContainer::Defrag, sc.get(), 0.9);
     }
   }
 

example/particles/main.cpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2022. 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
@@ -13,8 +13,24 @@
 
 #include "parthenon_manager.hpp"
 
+#include "bvals/boundary_conditions.hpp"
+#include "bvals/boundary_conditions_generic.hpp"
 #include "particles.hpp"
 
+using namespace parthenon::BoundaryFunction;
+
+// Example inner boundary condition (this just reuses existing features) to show how to
+// create and enroll a user swarm boundary condition. Note that currently both Swarm and
+// field boundary conditions must be provided when "user" is specified.
+// Note that BCType::Periodic cannot be enrolled as a user boundary condition.
+void SwarmUserInnerX1(std::shared_ptr<Swarm> &swarm) {
+  GenericSwarmBC<X1DIR, BCSide::Inner, BCType::Outflow>(swarm);
+}
+
+void SwarmUserOuterX1(std::shared_ptr<Swarm> &swarm) {
+  GenericSwarmBC<X1DIR, BCSide::Outer, BCType::Outflow>(swarm);
+}
+
 int main(int argc, char *argv[]) {
   using parthenon::ParthenonManager;
   using parthenon::ParthenonStatus;
@@ -33,27 +49,13 @@ int main(int argc, char *argv[]) {
 
   // Redefine parthenon defaults
   pman.app_input->ProcessPackages = particles_example::ProcessPackages;
-  pman.app_input->ProblemGenerator = particles_example::ProblemGenerator;
-  if (pman.pinput->GetString("parthenon/mesh", "ix1_bc") == "user") {
-    // In this case, we are setting a custom swarm boundary condition while still using
-    // a default parthenon boundary condition for cell variables. In general, one can
-    // provide both custom cell variable and swarm boundary conditions. However, to use
-    // custom boundary conditions for either cell variables or swarms, the parthenon
-    // boundary must be set to "user" and both cell variable and swarm boundaries provided
-    // as here.
-    pman.app_input->boundary_conditions[parthenon::BoundaryFace::inner_x1] =
-        parthenon::BoundaryFunction::OutflowInnerX1;
-    pman.app_input->swarm_boundary_conditions[parthenon::BoundaryFace::inner_x1] =
-        particles_example::SetSwarmIX1UserBC;
-  }
-  if (pman.pinput->GetString("parthenon/mesh", "ox1_bc") == "user") {
-    // Again, we use a default parthenon boundary condition for cell variables but a
-    // custom swarm boundary condition.
-    pman.app_input->boundary_conditions[parthenon::BoundaryFace::outer_x1] =
-        parthenon::BoundaryFunction::OutflowOuterX1;
-    pman.app_input->swarm_boundary_conditions[parthenon::BoundaryFace::outer_x1] =
-        particles_example::SetSwarmOX1UserBC;
-  }
+  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;
+  // Note that this example does not use a ProblemGenerator
   pman.ParthenonInitPackagesAndMesh();
 
   // This needs to be scoped so that the driver object is destructed before Finalize

example/particles/parthinput.particles

@@ -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
@@ -24,8 +24,10 @@ refinement = none
 nx1 = 16
 x1min = -0.5
 x1max = 0.5
-ix1_bc = periodic
-ox1_bc = periodic
+ix1_bc = user
+ox1_bc = user
+# ix1_bc = periodic # Optionally use periodic boundary conditions everywhere
+# ox1_bc = periodic
 
 nx2 = 16
 x2min = -0.5
@@ -40,8 +42,8 @@ ix3_bc = periodic
 ox3_bc = periodic
 
 <parthenon/meshblock>
-nx1 = 16
-nx2 = 16
+nx1 = 8
+nx2 = 8
 
 <parthenon/output0>
 file_type = hdf5

example/particles/particles.cpp

@@ -11,40 +11,29 @@
 // the public, perform publicly and display publicly, and to permit others to do so.
 //========================================================================================
 
-#include "particles.hpp"
-
 #include <algorithm>
 #include <limits>
 #include <memory>
 #include <string>
 #include <utility>
 #include <vector>
 
+#include "particles.hpp"
+
 // *************************************************//
 // redefine some internal parthenon functions      *//
 // *************************************************//
 namespace particles_example {
 
-std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>>
-SetSwarmIX1UserBC() {
-  return DeviceAllocate<ParticleBoundIX1User>();
-}
-
-std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>>
-SetSwarmOX1UserBC() {
-  return DeviceAllocate<ParticleBoundOX1User>();
-}
+using namespace parthenon;
+using namespace parthenon::BoundaryFunction;
 
 Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin) {
   Packages_t packages;
   packages.Add(particles_example::Particles::Initialize(pin.get()));
   return packages;
 }
 
-void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin) {
-  // Don't do anything for now
-}
-
 enum class DepositionMethod { per_particle, per_cell };
 
 // *************************************************//
@@ -534,8 +523,7 @@ TaskStatus StopCommunicationMesh(const BlockList_t &blocks) {
       // TODO(BRR) May want logic like this if we have non-blocking TaskRegions
       // if (nb.snb.rank != Globals::my_rank) {
       //  if (swarm->vbswarm->bd_var_.flag[nb.bufid] != BoundaryStatus::completed) {
-      //    printf("[%i] Neighbor %i not complete!\n", Globals::my_rank, n);
-      //    //return TaskStatus::incomplete;
+      //    return TaskStatus::incomplete;
       //  }
       //}
 

example/particles/particles.hpp

@@ -1,5 +1,5 @@
 //========================================================================================
-// (C) (or copyright) 2020-2022. 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
@@ -26,28 +26,6 @@ using namespace parthenon;
 
 namespace particles_example {
 
-// Duplicates of ParticleBoundIX1Outflow and ParticleBoundOX1Outflow to illustrate custom
-// particle boundary conditions
-class ParticleBoundIX1User : public ParticleBound {
- public:
-  KOKKOS_INLINE_FUNCTION void Apply(const int n, double &x, double &y, double &z,
-                                    const SwarmDeviceContext &swarm_d) const override {
-    if (x < swarm_d.x_min_global_) {
-      swarm_d.MarkParticleForRemoval(n);
-    }
-  }
-};
-
-class ParticleBoundOX1User : public ParticleBound {
- public:
-  KOKKOS_INLINE_FUNCTION void Apply(const int n, double &x, double &y, double &z,
-                                    const SwarmDeviceContext &swarm_d) const override {
-    if (x > swarm_d.x_max_global_) {
-      swarm_d.MarkParticleForRemoval(n);
-    }
-  }
-};
-
 typedef Kokkos::Random_XorShift64_Pool<> RNGPool;
 
 class ParticleDriver : public EvolutionDriver {
@@ -63,13 +41,8 @@ class ParticleDriver : public EvolutionDriver {
   LowStorageIntegrator integrator;
 };
 
-void ProblemGenerator(MeshBlock *pmb, ParameterInput *pin);
 Packages_t ProcessPackages(std::unique_ptr<ParameterInput> &pin);
 
-std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>> SetSwarmIX1UserBC();
-
-std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>> SetSwarmOX1UserBC();
-
 namespace Particles {
 
 std::shared_ptr<StateDescriptor> Initialize(ParameterInput *pin);

src/CMakeLists.txt

@@ -144,7 +144,6 @@ add_library(parthenon
   interface/state_descriptor.cpp
   interface/swarm.cpp
   interface/swarm.hpp
-  interface/swarm_boundaries.hpp
   interface/swarm_comms.cpp
   interface/swarm_container.cpp
   interface/swarm_default_names.hpp