Pr/inject species ratio

src/include/inject.hxx

@@ -11,13 +11,12 @@
 
 struct InjectBase
 {
-  InjectBase(int interval, int tau, int kind_n)
-    : interval(interval), tau(tau), kind_n(kind_n)
+
+  InjectBase(int interval)
+    : interval(interval)
   {}
 
   // param
   const int interval; // inject every so many steps
-  const int tau; // in steps
-  const int kind_n; // the number of particles to inject are based on this kind's density
 };
 

src/include/inject_impl.hxx

@@ -1,14 +1,14 @@
 
 
 #include "../libpsc/psc_output_fields/fields_item_moments_1st.hxx"
+#include "../libpsc/psc_output_fields/fields_item_fields.hxx"
 #include <bnd.hxx>
 #include <fields.hxx>
 #include <fields_item.hxx>
 #include <inject.hxx>
 
 #include <stdlib.h>
 #include <string>
-
 // ======================================================================
 // Inject_
 
@@ -21,16 +21,23 @@ struct Inject_ : InjectBase
   using real_t = typename Mparticles::real_t;
   using ItemMoment_t = _ItemMoment;
   using SetupParticles = ::SetupParticles<Mparticles>;
+#ifdef USE_CUDA
+  using MFields = HMFields;  
+#else
+  using MFields = MfieldsC;
+#endif
 
   // ----------------------------------------------------------------------
   // ctor
-
-  Inject_(const Grid_t& grid, int interval, int tau, int kind_n,
-          Target_t target, SetupParticles& setup_particles)
-    : InjectBase{interval, tau, kind_n},
-      target_{target},
+  template <typename FUNC>
+  Inject_( Grid_t& grid, int interval,
+          SetupParticles& setup_particles,
+          FUNC init_npt)
+          //std::function<void (int, int, psc_particle_npt&, MFields)> init_npt)
+    : InjectBase{interval},
       moment_n_{grid},
-      setup_particles_{setup_particles}
+      setup_particles_{setup_particles},
+      init_npt_{init_npt}
   {}
 
   // ----------------------------------------------------------------------
@@ -60,36 +67,22 @@ struct Inject_ : InjectBase
     auto mf_n = evalMfields(moment_n_);
     prof_stop(pr_2);
 
-    prof_start(pr_3);
-    //auto& mf_n = mres.template get_as<Mfields>(kind_n, kind_n + 1);
-    prof_stop(pr_3);
-
-    real_t fac = (interval * grid.dt / tau) / (1. + interval * grid.dt / tau);
-
     auto lf_init_npt = [&](int kind, Double3 pos, int p, Int3 idx,
                            psc_particle_npt& npt) {
-      if (target_.is_inside(pos)) {
-        target_.init_npt(kind, pos, npt);
-        npt.n -= mf_n[p](kind_n, idx[0], idx[1], idx[2]);
-        if (npt.n < 0) {
-          npt.n = 0;
-        }
-        npt.n *= fac;
-      }
-    };
 
+        init_npt_(kind, pos, p, idx, npt, mf_n);
+    };
     prof_start(pr_4);
     setup_particles_.setupParticles(mprts, lf_init_npt);
     prof_stop(pr_4);
 
-    //mres.put_as(mf_n, 0, 0);
     prof_stop(pr);
   }
 
 private:
-  Target_t target_;
   ItemMoment_t moment_n_;
   SetupParticles setup_particles_;
+  std::function<void (int, Double3, int, Int3,  psc_particle_npt&, MFields&)> init_npt_;
 };
 
 // ======================================================================

src/libpsc/tests/test_inject.cxx

@@ -99,7 +99,7 @@ TYPED_TEST(InjectTest, Test1)
   const int inject_tau = 10;
   auto target = InjectTestTarget{};
   auto setup_particles = SetupParticles<Mparticles>{grid};
-  Inject inject{grid, inject_interval, inject_tau, 0, target, setup_particles}; // FIXME, can't use "auto inject = Inject{...}", though I want to
+  Inject inject{grid, inject_interval, inject_tau, target, setup_particles, 1, -1, 0.0}; // FIXME, can't use "auto inject = Inject{...}", though I want to
 
   // let's start with no particles
   Mparticles mprts{grid};

src/psc_flatfoil_yz.cxx

@@ -11,7 +11,9 @@
 #include "heating_spot_foil.hxx"
 #include "inject_impl.hxx"
 
-#define DIM_3D
+#include "../libpsc/psc_output_fields/fields_item_fields.hxx"
+
+//#define DIM_3D
 
 // ======================================================================
 // Particle kinds
@@ -199,7 +201,7 @@ void setupParameters()
   // -- set some generic PSC parameters
   psc_params.nmax = 10000001; // 5001;
   psc_params.cfl = 0.75;
-  psc_params.write_checkpoint_every_step = 1000;
+  psc_params.write_checkpoint_every_step = 100;
   psc_params.stats_every = 1;
 
   // -- start from checkpoint:
@@ -249,9 +251,12 @@ Grid_t* setupGrid()
   Int3 gdims = {160, 160, 3 * 160};        // global number of grid points
   Int3 np = {5, 5, 3 * 5};                 // division into patches
 #else
-  Grid_t::Real3 LL = {1., 800., 3. * 800.}; // domain size (in d_e)
-  Int3 gdims = {1, 1600, 3 * 1600};         // global number of grid points
-  Int3 np = {1, 50, 3 * 50};                // division into patches
+  //Grid_t::Real3 LL = {1., 800., 3. * 800.}; // domain size (in d_e)
+  //Int3 gdims = {1, 1600, 3 * 1600};         // global number of grid points
+  //Int3 np = {1, 50, 3 * 50};                // division into patches
+    Grid_t::Real3 LL = {1., 32., 3200.}; // domain size (in d_e)
+    Int3 gdims = {1, 32, 2*3200};        // global number of grid points
+    Int3 np = {1, 1, 2*100};                 // division into patches
 #endif
 
   Grid_t::Domain domain{gdims, LL, -.5 * LL, np};
@@ -478,8 +483,35 @@ void run()
   setup_particles.fractional_n_particles_per_cell = true;
   setup_particles.neutralizing_population = MY_ION;
 
-  Inject inject{grid,        g.inject_interval, inject_tau,
-                MY_ELECTRON, inject_target,     setup_particles};
+#ifdef USE_CUDA
+  using MFields = HMFields;  
+#else
+  using MFields = MfieldsC;
+#endif
+  double fac = (g.inject_interval * grid.dt / inject_tau) / (1. + g.inject_interval * grid.dt / inject_tau);
+  auto lf_init_npt = [&fac, &inject_target](int kind, Double3 pos, int p, Int3 idx,
+                         psc_particle_npt& npt, MFields& mf_n) {
+
+    if (inject_target.is_inside(pos)) {
+
+      if(kind == MY_ELECTRON_HE){
+        inject_target.init_npt(MY_ELECTRON, pos, npt);
+        npt.n -= mf_n[p](MY_ELECTRON, idx[0], idx[1], idx[2]);
+        npt.n *= g.electron_HE_ratio;
+      }else{
+        inject_target.init_npt(kind, pos, npt);
+        npt.n -= mf_n[p](kind, idx[0], idx[1], idx[2]);
+        if(kind == MY_ELECTRON)
+          npt.n *= (1 - g.electron_HE_ratio);
+      }
+      if (npt.n < 0) {
+        npt.n = 0;
+      }
+      npt.n *= fac;
+    }
+  };
+
+  Inject inject{grid, g.inject_interval, setup_particles, lf_init_npt};
 
   auto lf_inject = [&](const Grid_t& grid, Mparticles& mprts) {
     static int pr_inject, pr_heating;