diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index d4831de71..c7b0258e6 100644 --- a/src/include/boundary_injector.hxx +++ b/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 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 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; } } } diff --git a/src/include/composite_injector.hxx b/src/include/composite_injector.hxx new file mode 100644 index 000000000..105c605bf --- /dev/null +++ b/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 +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 +CompositeInjector make_composite(INJECTOR_1 injector_1, + INJECTOR_2 injector_2) +{ + return CompositeInjector(injector_1, injector_2); +} diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 36508d96f..d424c8354 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/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 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, PscConfig::PushParticles>{{}, grid}; + BoundaryInjector{ + ParticleGenerator(1, 1), grid}; auto psc = makePscIntegrator(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, PscConfig::PushParticles>{{}, grid}; + BoundaryInjector{ + ParticleGenerator(-1, 1), grid}; auto psc = makePscIntegrator(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 outf{grid, {}}; + OutputParticles outp{grid, {}}; + DiagEnergies oute{grid.comm(), 0}; + auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + + auto inject_electrons = + BoundaryInjector{ + ParticleGenerator(-1, 0), grid}; + auto inject_ions = + BoundaryInjector{ + ParticleGenerator(-1, 1), grid}; + + auto inject = make_composite(inject_electrons, inject_ions); + + auto psc = makePscIntegrator(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