Enable user boundary conditions for particles

example/particles/main.cpp

@@ -35,16 +35,24 @@ int main(int argc, char *argv[]) {
   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;
+        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;
+        particles_example::SetSwarmOX1UserBC;
   }
   pman.ParthenonInitPackagesAndMesh();
 

example/particles/particles.cpp

@@ -26,12 +26,12 @@
 namespace particles_example {
 
 std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>>
-SetSwarmIx1UserBC() {
+SetSwarmIX1UserBC() {
   return DeviceAllocate<ParticleBoundIX1User>();
 }
 
 std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>>
-SetSwarmOx1UserBC() {
+SetSwarmOX1UserBC() {
   return DeviceAllocate<ParticleBoundOX1User>();
 }
 

example/particles/particles.hpp

@@ -66,9 +66,9 @@ class ParticleDriver : public EvolutionDriver {
 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>> SetSwarmIX1UserBC();
 
-std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>> SetSwarmOx1UserBC();
+std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>> SetSwarmOX1UserBC();
 
 namespace Particles {
 

src/interface/swarm.cpp

@@ -92,107 +92,47 @@ Swarm::Swarm(const std::string &label, const Metadata &metadata, const int nmax_
   marked_for_removal_.data.DeepCopy(marked_for_removal_h);
 }
 
-void Swarm::AllocateBoundaries() {
-  auto pmb = GetBlockPointer();
+template <class BOutflow, class BPeriodic, int iFace>
+void Swarm::AllocateBoundariesImpl_(MeshBlock *pmb) {
   std::stringstream msg;
-
   auto &bcs = pmb->pmy_mesh->mesh_bcs;
-
-  if (bcs[0] == BoundaryFlag::outflow) {
-    bounds_uptrs[0] = DeviceAllocate<ParticleBoundIX1Outflow>();
-  } else if (bcs[0] == BoundaryFlag::periodic) {
-    bounds_uptrs[0] = DeviceAllocate<ParticleBoundIX1Periodic>();
-  } else if (bcs[0] == BoundaryFlag::user) {
-    if (pmb->pmy_mesh->SwarmBndryFnctn[0] != nullptr) {
-      bounds_uptrs[0] = pmb->pmy_mesh->SwarmBndryFnctn[0]();
-    } else {
-      msg << "ix1 user boundary requested but provided function is null!";
-      PARTHENON_THROW(msg);
-    }
-  } else {
-    msg << "ix1 boundary flag " << static_cast<int>(bcs[0]) << " not supported!";
-    PARTHENON_THROW(msg);
-  }
-
-  if (bcs[1] == BoundaryFlag::outflow) {
-    bounds_uptrs[1] = DeviceAllocate<ParticleBoundOX1Outflow>();
-  } else if (bcs[1] == BoundaryFlag::periodic) {
-    bounds_uptrs[1] = DeviceAllocate<ParticleBoundOX1Periodic>();
-  } else if (bcs[1] == BoundaryFlag::user) {
-    if (pmb->pmy_mesh->SwarmBndryFnctn[1] != nullptr) {
-      bounds_uptrs[1] = pmb->pmy_mesh->SwarmBndryFnctn[1]();
+  if (bcs[iFace] == BoundaryFlag::outflow) {
+    bounds_uptrs[iFace] = DeviceAllocate<BOutflow>();
+  } else if (bcs[iFace] == BoundaryFlag::periodic) {
+    bounds_uptrs[iFace] = DeviceAllocate<BPeriodic>();
+  } else if (bcs[iFace] == BoundaryFlag::user) {
+    if (pmb->pmy_mesh->SwarmBndryFnctn[iFace] != nullptr) {
+      bounds_uptrs[iFace] = pmb->pmy_mesh->SwarmBndryFnctn[iFace]();
     } else {
-      msg << "ox1 user boundary requested but provided function is null!";
+      msg << "ix" << iFace + 1
+          << " user boundary requested but provided function is null!";
       PARTHENON_THROW(msg);
     }
   } else {
-    msg << "ox1 boundary flag " << static_cast<int>(bcs[1]) << " not supported!";
-    PARTHENON_THROW(msg);
-  }
-
-  if (bcs[2] == BoundaryFlag::outflow) {
-    bounds_uptrs[2] = DeviceAllocate<ParticleBoundIX2Outflow>();
-  } else if (bcs[2] == BoundaryFlag::periodic) {
-    bounds_uptrs[2] = DeviceAllocate<ParticleBoundIX2Periodic>();
-  } else if (bcs[2] == BoundaryFlag::user) {
-    if (pmb->pmy_mesh->SwarmBndryFnctn[2] != nullptr) {
-      bounds_uptrs[2] = pmb->pmy_mesh->SwarmBndryFnctn[2]();
-    } else {
-      msg << "ix2 user boundary requested but provided function is null!";
-      PARTHENON_THROW(msg);
-    }
-  } else {
-    msg << "ix2 boundary flag " << static_cast<int>(bcs[2]) << " not supported!";
+    msg << "ix" << iFace + 1 << " boundary flag " << static_cast<int>(bcs[iFace])
+        << " not supported!";
     PARTHENON_THROW(msg);
   }
+}
 
-  if (bcs[3] == BoundaryFlag::outflow) {
-    bounds_uptrs[3] = DeviceAllocate<ParticleBoundOX2Outflow>();
-  } else if (bcs[3] == BoundaryFlag::periodic) {
-    bounds_uptrs[3] = DeviceAllocate<ParticleBoundOX2Periodic>();
-  } else if (bcs[3] == BoundaryFlag::user) {
-    if (pmb->pmy_mesh->SwarmBndryFnctn[3] != nullptr) {
-      bounds_uptrs[3] = pmb->pmy_mesh->SwarmBndryFnctn[3]();
-    } else {
-      msg << "ox2 user boundary requested but provided function is null!";
-      PARTHENON_THROW(msg);
-    }
-  } else {
-    msg << "ox2 boundary flag " << static_cast<int>(bcs[3]) << " not supported!";
-    PARTHENON_THROW(msg);
-  }
+void Swarm::AllocateBoundaries() {
+  auto pmb = GetBlockPointer();
+  std::stringstream msg;
 
-  if (bcs[4] == BoundaryFlag::outflow) {
-    bounds_uptrs[4] = DeviceAllocate<ParticleBoundIX3Outflow>();
-  } else if (bcs[4] == BoundaryFlag::periodic) {
-    bounds_uptrs[4] = DeviceAllocate<ParticleBoundIX3Periodic>();
-  } else if (bcs[4] == BoundaryFlag::user) {
-    if (pmb->pmy_mesh->SwarmBndryFnctn[4] != nullptr) {
-      bounds_uptrs[4] = pmb->pmy_mesh->SwarmBndryFnctn[4]();
-    } else {
-      msg << "ix3 user boundary requested but provided function is null!";
-      PARTHENON_THROW(msg);
-    }
-  } else {
-    msg << "ix3 boundary flag " << static_cast<int>(bcs[4]) << " not supported!";
-    PARTHENON_THROW(msg);
-  }
+  auto &bcs = pmb->pmy_mesh->mesh_bcs;
 
-  if (bcs[5] == BoundaryFlag::outflow) {
-    bounds_uptrs[5] = DeviceAllocate<ParticleBoundOX3Outflow>();
-  } else if (bcs[5] == BoundaryFlag::periodic) {
-    bounds_uptrs[5] = DeviceAllocate<ParticleBoundOX3Periodic>();
-  } else if (bcs[5] == BoundaryFlag::user) {
-    if (pmb->pmy_mesh->SwarmBndryFnctn[5] != nullptr) {
-      bounds_uptrs[5] = pmb->pmy_mesh->SwarmBndryFnctn[5]();
-    } else {
-      msg << "ox3 user boundary requested but provided function is null!";
-      PARTHENON_THROW(msg);
-    }
-  } else {
-    msg << "ox3 boundary flag " << static_cast<int>(bcs[5]) << " not supported!";
-    PARTHENON_THROW(msg);
-  }
+  AllocateBoundariesImpl_<ParticleBoundIX1Outflow, ParticleBoundIX1Periodic, 0>(
+      pmb.get());
+  AllocateBoundariesImpl_<ParticleBoundOX1Outflow, ParticleBoundOX1Periodic, 1>(
+      pmb.get());
+  AllocateBoundariesImpl_<ParticleBoundIX2Outflow, ParticleBoundIX2Periodic, 2>(
+      pmb.get());
+  AllocateBoundariesImpl_<ParticleBoundOX2Outflow, ParticleBoundOX2Periodic, 3>(
+      pmb.get());
+  AllocateBoundariesImpl_<ParticleBoundIX3Outflow, ParticleBoundIX3Periodic, 4>(
+      pmb.get());
+  AllocateBoundariesImpl_<ParticleBoundOX3Outflow, ParticleBoundOX3Periodic, 5>(
+      pmb.get());
 
   for (int n = 0; n < 6; n++) {
     bounds_d.bounds[n] = bounds_uptrs[n].get();

src/interface/swarm.hpp

@@ -210,6 +210,9 @@ class Swarm {
   template <class T>
   vpack_types::SwarmVarList<T> MakeVarListAll_();
 
+  template <class BOutflow, class BPeriodic, int iFace>
+  void AllocateBoundariesImpl_(MeshBlock *pmb);
+
   void SetNeighborIndices1D_();
   void SetNeighborIndices2D_();
   void SetNeighborIndices3D_();

src/mesh/mesh.cpp

@@ -968,7 +968,7 @@ void Mesh::EnrollBndryFncts_(ApplicationInput *app_in) {
       } else {
         std::stringstream msg;
         msg << "A user boundary condition for face " << f
-            << " was requested, but not swarm condition was enrolled." << std::endl;
+            << " was requested, but no swarm condition was enrolled." << std::endl;
         PARTHENON_THROW(msg);
       }
       break;

tst/unit/test_swarm.cpp

@@ -24,6 +24,7 @@
 
 #include <catch2/catch.hpp>
 
+#include "bvals/bvals_interfaces.hpp"
 #include "interface/swarm.hpp"
 #include "mesh/mesh.hpp"
 
@@ -32,17 +33,37 @@
 
 using Real = double;
 using parthenon::ApplicationInput;
+using parthenon::BoundaryFlag;
+using parthenon::DeviceAllocate;
+using parthenon::DeviceDeleter;
 using parthenon::Mesh;
 using parthenon::MeshBlock;
 using parthenon::Metadata;
 using parthenon::Packages_t;
 using parthenon::ParameterInput;
 using parthenon::ParArrayND;
+using parthenon::ParticleBound;
 using parthenon::Swarm;
+using parthenon::SwarmDeviceContext;
 using std::endl;
 
 constexpr int NUMINIT = 10;
 
+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);
+    }
+  }
+};
+
+std::unique_ptr<ParticleBound, DeviceDeleter<parthenon::DevMemSpace>>
+SetSwarmIX1UserBC() {
+  return DeviceAllocate<ParticleBoundIX1User>();
+}
+
 TEST_CASE("Swarm memory management", "[Swarm]") {
   std::stringstream is;
   is << "<parthenon/mesh>" << endl;
@@ -61,10 +82,16 @@ TEST_CASE("Swarm memory management", "[Swarm]") {
   Packages_t packages;
   auto meshblock = std::make_shared<MeshBlock>(1, 1);
   auto mesh = std::make_shared<Mesh>(pin.get(), app_in.get(), packages, 1);
+  mesh->mesh_bcs[0] = BoundaryFlag::user;
+  mesh->SwarmBndryFnctn[0] = SetSwarmIX1UserBC;
+  for (int i = 1; i < 6; i++) {
+    mesh->mesh_bcs[i] = BoundaryFlag::outflow;
+  }
   meshblock->pmy_mesh = mesh.get();
   Metadata m;
   auto swarm = std::make_shared<Swarm>("test swarm", m, NUMINIT);
   swarm->SetBlockPointer(meshblock);
+  swarm->AllocateBoundaries();
   auto swarm_d = swarm->GetDeviceContext();
   REQUIRE(swarm->GetNumActive() == 0);
   REQUIRE(swarm->GetMaxActiveIndex() == 0);