Polymorphic injectors, external currents, and diagnostics

src/include/boundary_injector.hxx

@@ -11,6 +11,7 @@
 #include "setup_particles.hxx"
 #include "kg/VecRange.hxx"
 #include "../libpsc/psc_push_particles/inc_push.cxx"
+#include "injector_base.hxx"
 
 /// @brief A particle generator for use with @ref BoundaryInjector. Samples
 /// particles from a (possibly shifted) Maxwellian distribution.
@@ -54,14 +55,16 @@ private:
 
 /// @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.
+/// BoundaryInjector per species.
 /// @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
+  : public InjectorBase<typename PUSH_PARTICLES::Mparticles,
+                        typename PUSH_PARTICLES::MfieldsState>
 {
   static const int INJECT_DIM_IDX_ = 1;
 
@@ -82,6 +85,11 @@ public:
       prts_per_unit_density_{grid.norm.prts_per_unit_density}
   {}
 
+  void inject(Mparticles& mprts, MfieldsState& mflds) override
+  {
+    (*this)(mprts, mflds);
+  }
+
   /// Injects particles at the lower y-bound as if there were a population of
   /// particles just beyond the edge. The imaginary particle population has unit
   /// density, and individual particles from that population are sampled using

src/include/external_current_base.hxx

@@ -0,0 +1,24 @@
+#pragma once
+
+template <typename MFIELDS_STATE>
+struct ExternalCurrentBase
+{
+  using MfieldsState = MFIELDS_STATE;
+
+  virtual void inject_current(MfieldsState& mflds) = 0;
+};
+
+template <typename MFIELDS_STATE>
+struct ExternalCurrentFromLambda : ExternalCurrentBase<MFIELDS_STATE>
+{
+  using MfieldsState = MFIELDS_STATE;
+
+  ExternalCurrentFromLambda(std::function<void(MfieldsState&)> lambda)
+    : lambda{lambda}
+  {}
+
+  void inject_current(MfieldsState& mflds) override { return lambda(mflds); }
+
+private:
+  std::function<void(MfieldsState&)> lambda;
+};

src/include/injector_base.hxx

@@ -0,0 +1,29 @@
+#pragma once
+
+template <typename MPARTICLES, typename MFIELDS_STATE>
+struct InjectorBase
+{
+  using Mparticles = MPARTICLES;
+  using MfieldsState = MFIELDS_STATE;
+
+  virtual void inject(Mparticles& mprts, MfieldsState& mflds) = 0;
+};
+
+template <typename MPARTICLES, typename MFIELDS_STATE>
+struct InjectFromLambda : InjectorBase<MPARTICLES, MFIELDS_STATE>
+{
+  using Mparticles = MPARTICLES;
+  using MfieldsState = MFIELDS_STATE;
+
+  InjectFromLambda(std::function<void(Mparticles&, MfieldsState&)> lambda)
+    : lambda{lambda}
+  {}
+
+  void inject(Mparticles& mprts, MfieldsState& mflds) override
+  {
+    return lambda(mprts, mflds);
+  }
+
+private:
+  std::function<void(Mparticles&, MfieldsState&)> lambda;
+};

src/include/psc.hxx

@@ -8,6 +8,8 @@
 #include <setup_particles.hxx>
 
 #include "../libpsc/vpic/fields_item_vpic.hxx"
+#include "injector_base.hxx"
+#include "external_current_base.hxx"
 #include <checks_params.hxx>
 #include <output_particles.hxx>
 #include <push_particles.hxx>
@@ -99,8 +101,7 @@ inline double courant_length(const Grid_t::Domain& domain)
 // ======================================================================
 // Psc
 
-template <typename PscConfig, typename Diagnostics, typename InjectParticles,
-          typename ExtCurrent>
+template <typename PscConfig, typename Diagnostics>
 struct Psc
 {
   using Mparticles = typename PscConfig::Mparticles;
@@ -116,14 +117,15 @@ struct Psc
   using BndFields = typename PscConfig::BndFields;
   using BndParticles = typename PscConfig::BndParticles;
   using Dim = typename PscConfig::Dim;
+  using InjectorBaseT = InjectorBase<Mparticles, MfieldsState>;
+  using ExternalCurrentBaseT = ExternalCurrentBase<MfieldsState>;
 
   // ----------------------------------------------------------------------
   // ctor
 
   Psc(const PscParams& params, Grid_t& grid, MfieldsState& mflds,
       Mparticles& mprts, Balance& balance, Collision& collision, Checks& checks,
-      Marder& marder, Diagnostics& diagnostics,
-      InjectParticles& inject_particles, ExtCurrent& ext_current)
+      Marder& marder, Diagnostics& diagnostics)
     : p_{params},
       grid_{&grid},
       mflds_{mflds},
@@ -134,8 +136,6 @@ struct Psc
       marder_{marder},
       bndp_{grid},
       diagnostics_{diagnostics},
-      inject_particles_{inject_particles},
-      ext_current_{ext_current},
       checkpointing_{params.write_checkpoint_every_step}
   {
     time_start_ = MPI_Wtime();
@@ -159,6 +159,23 @@ struct Psc
     initialize();
   }
 
+  // ----------------------------------------------------------------------
+  // API for modifying various internal components
+
+  void add_injector(InjectorBaseT* injector)
+  {
+    if (injector) {
+      injectors_.push_back(injector);
+    }
+  }
+
+  void add_external_current(ExternalCurrentBaseT* external_current)
+  {
+    if (external_current) {
+      external_currents_.push_back(external_current);
+    }
+  }
+
   // ----------------------------------------------------------------------
   // initialize_stats
 
@@ -343,12 +360,16 @@ struct Psc
 
     // === particle injection
     prof_start(pr_inject_prts);
-    inject_particles_(mprts_, mflds_);
+    for (auto injector : injectors_) {
+      injector->inject(mprts_, mflds_);
+    }
     prof_stop(pr_inject_prts);
 
     // === external current
     prof_start(pr_external_current);
-    this->ext_current_(grid(), mflds_);
+    for (auto external_current : external_currents_) {
+      external_current->inject_current(mflds_);
+    }
     prof_stop(pr_external_current);
 
     mpi_printf(comm, "***** Bnd particles...\n");
@@ -485,8 +506,8 @@ protected:
   Checks& checks_;
   Marder& marder_;
   Diagnostics& diagnostics_;
-  InjectParticles& inject_particles_;
-  ExtCurrent& ext_current_;
+  std::vector<InjectorBaseT*> injectors_;
+  std::vector<ExternalCurrentBaseT*> external_currents_;
 
   Sort sort_;
   PushParticles pushp_;
@@ -507,59 +528,19 @@ protected:
   int st_time_step;
 };
 
-// ======================================================================
-// InjectParticlesNone
-
-class InjectParticlesNone
-{
-public:
-  template <typename Mparticles, typename MfieldsState>
-  void operator()(Mparticles& mprts, MfieldsState& mflds)
-  {}
-};
-
-namespace
-{
-
-InjectParticlesNone injectParticlesNone;
-
-} // namespace
-
-// ======================================================================
-// ExtCurrentNone
-
-class ExtCurrentNone
-{
-public:
-  template <typename MfieldsState>
-  void operator()(const Grid_t& grid, MfieldsState& mflds)
-  {}
-};
-
-namespace
-{
-
-ExtCurrentNone extCurrentNone;
-
-} // namespace
 // ======================================================================
 // makePscIntegrator
 
 template <typename PscConfig, typename MfieldsState, typename Mparticles,
           typename Balance, typename Collision, typename Checks,
-          typename Marder, typename Diagnostics,
-          typename InjectParticles = InjectParticlesNone,
-          typename ExtCurrent = ExtCurrentNone>
-Psc<PscConfig, Diagnostics, InjectParticles, ExtCurrent> makePscIntegrator(
+          typename Marder, typename Diagnostics>
+Psc<PscConfig, Diagnostics> makePscIntegrator(
   const PscParams& params, Grid_t& grid, MfieldsState& mflds, Mparticles& mprts,
   Balance& balance, Collision& collision, Checks& checks, Marder& marder,
-  Diagnostics& diagnostics,
-  InjectParticles& inject_particles = injectParticlesNone,
-  ExtCurrent& ext_current = extCurrentNone)
+  Diagnostics& diagnostics)
 {
-  return {params,     grid,   mflds,  mprts,       balance,
-          collision,  checks, marder, diagnostics, inject_particles,
-          ext_current};
+  return {params,    grid,   mflds,  mprts,      balance,
+          collision, checks, marder, diagnostics};
 }
 
 // ======================================================================

src/libpsc/tests/test_boundary_injector.cxx

@@ -3,7 +3,6 @@
 #include "test_common.hxx"
 
 #include "boundary_injector.hxx"
-#include "composite_injector.hxx"
 
 #include "psc.hxx"
 #include "DiagnosticsDefault.h"
@@ -142,13 +141,13 @@ TEST(BoundaryInjectorTest, Integration1Particle)
   DiagEnergies oute{grid.comm(), 0};
   auto diagnostics = makeDiagnosticsDefault(outf, outp, oute);
 
-  auto inject_particles =
-    BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>{
-      ParticleGenerator(1, 1), grid};
+  auto psc =
+    makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts, balance,
+                                 collision, checks, marder, diagnostics);
 
-  auto psc = makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts,
-                                          balance, collision, checks, marder,
-                                          diagnostics, inject_particles);
+  psc.add_injector(
+    new BoundaryInjector<ParticleGenerator, typename PscConfig::PushParticles>(
+      ParticleGenerator(1, 1), grid));
 
   // ----------------------------------------------------------------------
   // set up initial conditions
@@ -205,13 +204,13 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles)
   DiagEnergies oute{grid.comm(), 0};
   auto diagnostics = makeDiagnosticsDefault(outf, outp, oute);
 
-  auto inject_particles =
-    BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>{
-      ParticleGenerator(-1, 1), grid};
+  auto psc =
+    makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts, balance,
+                                 collision, checks, marder, diagnostics);
 
-  auto psc = makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts,
-                                          balance, collision, checks, marder,
-                                          diagnostics, inject_particles);
+  psc.add_injector(
+    new BoundaryInjector<ParticleGenerator, PscConfig::PushParticles>(
+      ParticleGenerator(-1, 1), grid));
 
   // ----------------------------------------------------------------------
   // set up initial conditions
@@ -275,11 +274,12 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies)
     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);
 
-  auto psc = makePscIntegrator<PscConfig>(psc_params, grid, mflds, mprts,
-                                          balance, collision, checks, marder,
-                                          diagnostics, inject);
+  psc.add_injector(&inject_ions);
+  psc.add_injector(&inject_electrons);
 
   // ----------------------------------------------------------------------
   // set up initial conditions

src/psc_2d_shock.cxx

@@ -608,9 +608,12 @@ void run()
   // ----------------------------------------------------------------------
   // hand off to PscIntegrator to run the simulation
 
-  auto psc = makePscIntegrator<PscConfig>(psc_params, *grid_ptr, mflds, mprts,
-                                          balance, collision, checks, marder,
-                                          diagnostics, lf_inject_heat);
+  auto psc =
+    makePscIntegrator<PscConfig>(psc_params, *grid_ptr, mflds, mprts, balance,
+                                 collision, checks, marder, diagnostics);
+
+  psc.add_injector(
+    new InjectFromLambda<Mparticles, MfieldsState>(lf_inject_heat));
 
   psc.integrate();
 }