From 86e2f0aa535f4b1b1f3b22f9a7b0f7ff01915e32 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 16 Oct 2025 12:53:11 -0400 Subject: [PATCH 1/7] +composite_injector --- src/include/composite_injector.hxx | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) create mode 100644 src/include/composite_injector.hxx diff --git a/src/include/composite_injector.hxx b/src/include/composite_injector.hxx new file mode 100644 index 000000000..e0d3a5466 --- /dev/null +++ b/src/include/composite_injector.hxx @@ -0,0 +1,30 @@ +#pragma once + +/// @brief An injector comprising two other injectors. Enables multiple +/// injectors in a single psc case. +/// @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; +}; From b98144e4eda1e7b728a6f0e1af51d0e9f4d4da1c Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 17 Oct 2025 17:02:57 -0400 Subject: [PATCH 2/7] composite_injector: +make_composite --- src/include/composite_injector.hxx | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/src/include/composite_injector.hxx b/src/include/composite_injector.hxx index e0d3a5466..7a3fe3648 100644 --- a/src/include/composite_injector.hxx +++ b/src/include/composite_injector.hxx @@ -28,3 +28,20 @@ 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); +} From 7ca11ea11f41676199cf47f65489e446d60963ef Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 20 Oct 2025 09:57:02 -0400 Subject: [PATCH 3/7] test_boundary_injector: +IntegrationManySpecies --- src/libpsc/tests/test_boundary_injector.cxx | 89 +++++++++++++++++++-- 1 file changed, 84 insertions(+), 5 deletions(-) diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 36508d96f..52532226f 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,7 +82,8 @@ Grid_t* setupGrid() return new Grid_t{domain, bc, kinds, norm, dt, -1, ibn}; } -template +// TODO this template is dumb, just make a proper ctor +template struct ParticleGenerator { using Real = psc::particle::Inject::Real; @@ -99,10 +101,9 @@ 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}; + return {x, u, w, KIND_IDX, tag}; } int n_injected = 0; @@ -140,7 +141,8 @@ TEST(BoundaryInjectorTest, Integration1Particle) auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); auto inject_particles = - BoundaryInjector, PscConfig::PushParticles>{{}, grid}; + BoundaryInjector, PscConfig::PushParticles>{{}, + grid}; auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, balance, collision, checks, marder, @@ -202,7 +204,8 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); auto inject_particles = - BoundaryInjector, PscConfig::PushParticles>{{}, grid}; + BoundaryInjector, PscConfig::PushParticles>{{}, + grid}; auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, balance, collision, checks, marder, @@ -232,6 +235,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, PscConfig::PushParticles>{{}, + grid}; + auto inject_ions = + BoundaryInjector, PscConfig::PushParticles>{{}, + 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 From 3301567dfbe3f9817146a89bc09956ee2db219e3 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 24 Oct 2025 14:47:09 -0400 Subject: [PATCH 4/7] boundary_injector: don't override total boundary current --- src/include/boundary_injector.hxx | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index d4831de71..cc6ef07bf 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -105,6 +105,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 +151,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; } } } From 087842d0e6b97055bc37deb70122a63fb930a3c3 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 24 Oct 2025 14:52:20 -0400 Subject: [PATCH 5/7] test_boundary_injector: don't template ParticleGenerator --- src/libpsc/tests/test_boundary_injector.cxx | 28 ++++++++++++--------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 52532226f..d424c8354 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -82,17 +82,19 @@ Grid_t* setupGrid() return new Grid_t{domain, bc, kinds, norm, dt, -1, ibn}; } -// TODO this template is dumb, just make a proper ctor -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; @@ -103,10 +105,12 @@ struct ParticleGenerator Real w = 1.0; psc::particle::Tag tag = 0; - return {x, u, w, KIND_IDX, tag}; + return {x, u, w, kind_idx, tag}; } int n_injected = 0; + int max_n_injected; + int kind_idx; }; TEST(BoundaryInjectorTest, Integration1Particle) @@ -141,8 +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, @@ -204,8 +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, @@ -267,11 +271,11 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); auto inject_electrons = - BoundaryInjector, PscConfig::PushParticles>{{}, - grid}; + BoundaryInjector{ + ParticleGenerator(-1, 0), grid}; auto inject_ions = - BoundaryInjector, PscConfig::PushParticles>{{}, - grid}; + BoundaryInjector{ + ParticleGenerator(-1, 1), grid}; auto inject = make_composite(inject_electrons, inject_ions); From 55d3182ca8ad38456aa2fd791b9c157cb06be40e Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 24 Oct 2025 15:15:54 -0400 Subject: [PATCH 6/7] boundary_injector: more comments --- src/include/boundary_injector.hxx | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index cc6ef07bf..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 { From 2f7dcf82a4a8ffbecb3cff03e4329a78ca7781a8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 24 Oct 2025 15:16:06 -0400 Subject: [PATCH 7/7] composite_injector: update comments --- src/include/composite_injector.hxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/include/composite_injector.hxx b/src/include/composite_injector.hxx index 7a3fe3648..105c605bf 100644 --- a/src/include/composite_injector.hxx +++ b/src/include/composite_injector.hxx @@ -1,7 +1,8 @@ #pragma once /// @brief An injector comprising two other injectors. Enables multiple -/// injectors in a single psc case. +/// 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