diff --git a/src/include/inject.hxx b/src/include/inject.hxx index b61c8f6862..b7ac7b5e92 100644 --- a/src/include/inject.hxx +++ b/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 }; diff --git a/src/include/inject_impl.hxx b/src/include/inject_impl.hxx index 910f28dbc3..e293d2b944 100644 --- a/src/include/inject_impl.hxx +++ b/src/include/inject_impl.hxx @@ -1,6 +1,7 @@ #include "../libpsc/psc_output_fields/fields_item_moments_1st.hxx" +#include "../libpsc/psc_output_fields/fields_item_fields.hxx" #include #include #include @@ -8,7 +9,6 @@ #include #include - // ====================================================================== // Inject_ @@ -21,16 +21,23 @@ struct Inject_ : InjectBase using real_t = typename Mparticles::real_t; using ItemMoment_t = _ItemMoment; using SetupParticles = ::SetupParticles; +#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 + Inject_( Grid_t& grid, int interval, + SetupParticles& setup_particles, + FUNC init_npt) + //std::function 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(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 init_npt_; }; // ====================================================================== diff --git a/src/libpsc/tests/test_inject.cxx b/src/libpsc/tests/test_inject.cxx index b63ab59fc9..59231b3335 100644 --- a/src/libpsc/tests/test_inject.cxx +++ b/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{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}; diff --git a/src/psc_flatfoil_yz.cxx b/src/psc_flatfoil_yz.cxx index 9a6ddf4b1d..511b1c0aff 100644 --- a/src/psc_flatfoil_yz.cxx +++ b/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;