Composite injector

src/include/boundary_injector.hxx

@@ -11,6 +11,8 @@
 #include "setup_particles.hxx"
 #include "../libpsc/psc_push_particles/inc_push.cxx"
 
+/// @brief A particle generator for use with @ref BoundaryInjector. Samples
+/// particles from a (possibly shifted) Maxwellian distribution.
 class ParticleGeneratorMaxwellian
 {
 public:
@@ -49,6 +51,14 @@ private:
   rng::Uniform<Real> uniform_dist{0.0, 1.0};
 };
 
+/// @brief Injects particles on a given boundary, sampling from a given particle
+/// generator. For precise control over multiple particle species, use one
+/// BoundaryInjector per species, combined with @ref CompositeInjector.
+/// @tparam PARTICLE_GENERATOR a type that defines `get(min_pos, pos_range)` and
+/// returns an injectable particle within that range of positions (usually a
+/// grid cell); see @ref ParticleGeneratorMaxwellian
+/// @tparam PUSH_PARTICLES type that provides the types `Mparticles`,
+/// `MfieldsState`, `Current`, `real_t`, etc.
 template <typename PARTICLE_GENERATOR, typename PUSH_PARTICLES>
 class BoundaryInjector
 {
@@ -105,6 +115,10 @@ public:
 
       for (Int3 cell_idx = ilo; cell_idx[0] < ihi[0]; cell_idx[0]++) {
         for (cell_idx[2] = ilo[2]; cell_idx[2] < ihi[2]; cell_idx[2]++) {
+          auto boundary_current_before =
+            flds.storage()(cell_idx[0] + grid.ibn[0], -1 + grid.ibn[1],
+                           cell_idx[2] + grid.ibn[2], JXI + 1);
+
           cell_idx[INJECT_DIM_IDX_] = -1;
 
           Real3 cell_corner =
@@ -147,11 +161,13 @@ public:
             charge_injected += qni_wni;
           }
 
-          // set current in boundary to account for injected charge
+          // override whatever current was deposited in the boundary to be
+          // consistent with the charge having started completely out of the
+          // domain
           flds.storage()(cell_idx[0] + grid.ibn[0], -1 + grid.ibn[1],
                          cell_idx[2] + grid.ibn[2], JXI + 1) =
-            charge_injected * grid.domain.dx[1] / grid.dt /
-            prts_per_unit_density;
+            boundary_current_before + charge_injected * grid.domain.dx[1] /
+                                        grid.dt / prts_per_unit_density;
         }
       }
     }

src/include/composite_injector.hxx

@@ -0,0 +1,48 @@
+#pragma once
+
+/// @brief An injector comprising two other injectors. Enables multiple
+/// injectors in a single psc case. Useful for injecting multiple species with
+/// @ref BoundaryInjector.
+/// @tparam INJECTOR_1 first type of injector
+/// @tparam INJECTOR_2 second type of injector
+template <typename INJECTOR_1, typename INJECTOR_2>
+class CompositeInjector
+{
+public:
+  using Mparticles = typename INJECTOR_1::Mparticles;
+  using MfieldsState = typename INJECTOR_1::MfieldsState;
+
+  CompositeInjector(INJECTOR_1 injector_1, INJECTOR_2 injector_2)
+    : injector_1{injector_1}, injector_2{injector_2}
+  {}
+
+  /// @brief Calls both member injectors in order.
+  /// @param mprts particles
+  /// @param mflds fields
+  void operator()(Mparticles& mprts, MfieldsState& mflds)
+  {
+    injector_1(mprts, mflds);
+    injector_2(mprts, mflds);
+  }
+
+private:
+  INJECTOR_1 injector_1;
+  INJECTOR_2 injector_2;
+};
+
+/// @brief Helper method to deduce the type of a composite injector, enabling
+/// e.g. `auto composite_injector = make_composite(injector_1, injector_2);`
+///
+/// Wouldn't be necessary in C++17; see
+/// https://isocpp.github.io/CppCoreGuidelines/CppCoreGuidelines#Rt-deduce.
+/// @tparam INJECTOR_1 first type of injector
+/// @tparam INJECTOR_2 second type of injector
+/// @param injector_1 first injector
+/// @param injector_2 second injector
+/// @return
+template <typename INJECTOR_1, typename INJECTOR_2>
+CompositeInjector<INJECTOR_1, INJECTOR_2> make_composite(INJECTOR_1 injector_1,
+                                                         INJECTOR_2 injector_2)
+{
+  return CompositeInjector<INJECTOR_1, INJECTOR_2>(injector_1, injector_2);
+}

src/libpsc/tests/test_boundary_injector.cxx

@@ -3,6 +3,7 @@
 #include "test_common.hxx"
 
 #include "boundary_injector.hxx"
+#include "composite_injector.hxx"
 
 #include "psc.hxx"
 #include "DiagnosticsDefault.h"
@@ -81,16 +82,19 @@ Grid_t* setupGrid()
   return new Grid_t{domain, bc, kinds, norm, dt, -1, ibn};
 }
 
-template <int PRT_COUNT>
 struct ParticleGenerator
 {
   using Real = psc::particle::Inject::Real;
   using Real3 = psc::particle::Inject::Real3;
 
+  ParticleGenerator(int max_n_injected, int kind_idx)
+    : max_n_injected(max_n_injected), kind_idx(kind_idx)
+  {}
+
   psc::particle::Inject get(Real3 min_pos, Real3 pos_range)
   {
     Real uy = 2.0;
-    if (PRT_COUNT > 0 && n_injected++ >= PRT_COUNT) {
+    if (max_n_injected > 0 && n_injected++ >= max_n_injected) {
       // Setting uy=0 means the particle won't enter the domain, and thus won't
       // be injected
       uy = 0.0;
@@ -99,13 +103,14 @@ struct ParticleGenerator
     Real3 x = min_pos + pos_range * Real3{0, .999, 0.};
     Real3 u{0.0, uy, 0.0};
     Real w = 1.0;
-    int kind_idx = 1;
     psc::particle::Tag tag = 0;
 
     return {x, u, w, kind_idx, tag};
   }
 
   int n_injected = 0;
+  int max_n_injected;
+  int kind_idx;
 };
 
 TEST(BoundaryInjectorTest, Integration1Particle)
@@ -140,7 +145,8 @@ TEST(BoundaryInjectorTest, Integration1Particle)
   auto diagnostics = makeDiagnosticsDefault(outf, outp, oute);
 
   auto inject_particles =
-    BoundaryInjector<ParticleGenerator<1>, PscConfig::PushParticles>{{}, grid};
+    BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>{
+      ParticleGenerator(1, 1), grid};
 
   auto psc = makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts,
                                           balance, collision, checks, marder,
@@ -202,7 +208,8 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles)
   auto diagnostics = makeDiagnosticsDefault(outf, outp, oute);
 
   auto inject_particles =
-    BoundaryInjector<ParticleGenerator<-1>, PscConfig::PushParticles>{{}, grid};
+    BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>{
+      ParticleGenerator(-1, 1), grid};
 
   auto psc = makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts,
                                           balance, collision, checks, marder,
@@ -232,6 +239,82 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles)
   ASSERT_GT(prts.size(), 1);
 }
 
+TEST(BoundaryInjectorTest, IntegrationManySpecies)
+{
+  // ----------------------------------------------------------------------
+  // setup
+
+  PscParams psc_params;
+
+  psc_params.nmax = 1;
+  psc_params.stats_every = 1;
+  psc_params.cfl = .75;
+
+  auto grid_ptr = setupGrid();
+  auto& grid = *grid_ptr;
+
+  MfieldsState mflds{grid};
+  Mparticles mprts{grid};
+
+  ChecksParams checks_params{};
+  checks_params.continuity.check_interval = 1;
+  checks_params.gauss.check_interval = 1;
+  Checks checks{grid, MPI_COMM_WORLD, checks_params};
+
+  Balance balance{.1};
+  Collision collision{grid, 0, 0.1};
+  Marder marder(grid, 0.9, 3, false);
+
+  OutputFields<MfieldsState, Mparticles, Dim> outf{grid, {}};
+  OutputParticles outp{grid, {}};
+  DiagEnergies oute{grid.comm(), 0};
+  auto diagnostics = makeDiagnosticsDefault(outf, outp, oute);
+
+  auto inject_electrons =
+    BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>{
+      ParticleGenerator(-1, 0), grid};
+  auto inject_ions =
+    BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>{
+      ParticleGenerator(-1, 1), grid};
+
+  auto inject = make_composite(inject_electrons, inject_ions);
+
+  auto psc = makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts,
+                                          balance, collision, checks, marder,
+                                          diagnostics, inject);
+
+  // ----------------------------------------------------------------------
+  // set up initial conditions
+
+  ASSERT_EQ(grid.n_patches(), 1);
+  int p = 0;
+
+  // ----------------------------------------------------------------------
+  // run the simulation
+
+  auto accessor = mprts.accessor();
+  auto prts = accessor[p];
+
+  ASSERT_EQ(prts.size(), 0);
+
+  for (; grid.timestep_ < psc_params.nmax; grid.timestep_++) {
+    psc.step();
+
+    EXPECT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold);
+    EXPECT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold);
+  }
+
+  bool found_electrons = false;
+  bool found_ions = false;
+  for (auto prt : prts) {
+    found_electrons |= prt.kind() == 0;
+    found_ions |= prt.kind() == 1;
+  }
+
+  ASSERT_TRUE(found_electrons);
+  ASSERT_TRUE(found_ions);
+}
+
 // ======================================================================
 // main