diff --git a/src/include/particle.h b/src/include/particle.h index 225addd9a..d07b1e918 100644 --- a/src/include/particle.h +++ b/src/include/particle.h @@ -37,3 +37,11 @@ struct Inject } // namespace particle } // namespace psc + +inline std::ostream& operator<<(std::ostream& os, + const psc::particle::Inject& inj) +{ + os << "Inject{x=" << inj.x << ", u=" << inj.u << ", w=" << inj.w + << ", kind=" << inj.kind << ", tag=" << inj.tag << "}"; + return os; +} \ No newline at end of file diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx index bad72d12c..c522f3454 100644 --- a/src/include/setup_particles.hxx +++ b/src/include/setup_particles.hxx @@ -9,10 +9,10 @@ struct psc_particle_npt { - int kind; ///< particle kind - double n; ///< density - double p[3]; ///< momentum - double T[3]; ///< temperature + int kind; ///< particle kind + double n; ///< density + Double3 p; ///< momentum + Double3 T; ///< temperature psc::particle::Tag tag; }; @@ -241,6 +241,18 @@ struct SetupParticles }); } + // ---------------------------------------------------------------------- + // getWeight + + real_t getWeight(const psc_particle_np& np, int n_in_cell) + { + if (fractional_n_particles_per_cell) { + return 1.; + } else { + return np.n / (n_in_cell * norm_.cori); + } + } + // ---------------------------------------------------------------------- // setupParticles @@ -269,12 +281,7 @@ struct SetupParticles op_cellwise(grid, p, init_np, [&](int n_in_cell, psc_particle_np& np, Double3& pos) { for (int cnt = 0; cnt < n_in_cell; cnt++) { - real_t wni; - if (fractional_n_particles_per_cell) { - wni = 1.; - } else { - wni = np.n / (n_in_cell * norm_.cori); - } + real_t wni = getWeight(np, n_in_cell); auto prt = setupParticle(np, pos, wni); injector(prt); } diff --git a/src/libpsc/tests/CMakeLists.txt b/src/libpsc/tests/CMakeLists.txt index 913d29da9..1c29b8f49 100644 --- a/src/libpsc/tests/CMakeLists.txt +++ b/src/libpsc/tests/CMakeLists.txt @@ -24,6 +24,7 @@ add_psc_test(test_output_particles) add_psc_test(test_mfields) add_psc_test(test_mfields_cuda) add_psc_test(test_bnd) +add_psc_test(test_bnd_prt_inflow) add_psc_test(test_deposit) add_psc_test(test_current_deposition) add_psc_test(test_push_particles) diff --git a/src/libpsc/tests/test_bnd_prt_inflow.cxx b/src/libpsc/tests/test_bnd_prt_inflow.cxx new file mode 100644 index 000000000..ff1629340 --- /dev/null +++ b/src/libpsc/tests/test_bnd_prt_inflow.cxx @@ -0,0 +1,424 @@ + +#include + +#include "test_common.hxx" + +#include "psc_particles_double.h" +#include "psc_particles_single.h" +#include "particle_with_id.h" +#include "setup_particles.hxx" +#include "pushp.hxx" +#include "dim.hxx" +#ifdef USE_CUDA +#include "../libpsc/cuda/mparticles_cuda.hxx" +#include "../libpsc/cuda/mparticles_cuda.inl" +#endif +#include "particles_simple.inl" +#include + +template +struct Config +{ + using Mparticles = _Mparticles; + using MakeGrid = _MakeGrid; +}; + +using MparticlesInflowTestTypes = + ::testing::Types, + Config +#ifdef USE_CUDA + , + Config, MakeTestGridYZ1>, + Config, MakeTestGridYZ> +#endif + >; + +TYPED_TEST_SUITE(MparticlesInflowTest, MparticlesInflowTestTypes); + +TEST(TestSetupParticlesInflow, SetupSimple) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 1, 1}}; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, .1}; + Mparticles mprts{grid}; + + SetupParticles setup_particles(grid); + + Double3 pos = {5., 5., 5.}; + Double3 u = Double3{0.0, 0.0, 0.0}; + int tag = 0; + psc_particle_np np = {0, 1.0, [u]() { return u; }, tag}; + double wni = 1.0; + auto prt = setup_particles.setupParticle(np, pos, wni); + + EXPECT_EQ(prt.x, pos); + EXPECT_EQ(prt.u, u); + EXPECT_EQ(prt.w, wni); + EXPECT_EQ(prt.kind, 0); + EXPECT_EQ(prt.tag, tag); +} + +TEST(TestSetupParticlesInflow, SetupMaxwellian) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 1, 1}}; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, .1}; + Mparticles mprts{grid}; + + SetupParticles setup_particles(grid); + + Double3 pos = {5., 5., 5.}; + Double3 u = {0.0, 0.5, 0.0}; + Double3 T = {0.0, 0.0, 0.0}; + int tag = 0; + + psc_particle_npt npt = {0, 1.0, u, T, tag}; + auto p = setup_particles.createMaxwellian(npt); + psc_particle_np np = {0, 1.0, p, tag}; + double wni = 1.0; + auto prt = setup_particles.setupParticle(np, pos, wni); + + EXPECT_EQ(prt.x, pos); + EXPECT_EQ(prt.u, u); + EXPECT_EQ(prt.w, wni); + EXPECT_EQ(prt.kind, 0); + EXPECT_EQ(prt.tag, tag); +} + +template +class Inflow +{ +public: + using Mparticles = MP; + using real_t = typename Mparticles::real_t; + using Dim = DIM; + + // offset_in_cell: () -> double + Inflow(const Grid_t& grid, psc_particle_npt npt, + double (*offset_in_cell_dist)()) + : grid_(grid), + setup_particles_(grid), + advance_(grid.dt), + // FIXME np.p is a std::function and is called many times; better to use a + // lambda + np_{npt.kind, npt.n, setup_particles_.createMaxwellian(npt)}, + offset_in_cell_dist_(offset_in_cell_dist) + { + assert(INJECT_DIM_IDX_ >= 0); + } + + auto get_advanced_prt(Double3 pos, real_t wni) + { + auto prt = setup_particles_.setupParticle(np_, pos, wni); + + auto v = advance_.calc_v(prt.u); + advance_.push_x(prt.x, v, 1.0); + + return prt; + } + + template + void inject_into_boundary_cell(Injector& injector, + Int3 boundary_cell_global_idx) + { + assert(boundary_cell_global_idx[INJECT_DIM_IDX_] == 0); + boundary_cell_global_idx[INJECT_DIM_IDX_] = -1; + + int n_in_cell = setup_particles_.get_n_in_cell(np_); + double wni = setup_particles_.getWeight(np_, n_in_cell); + + for (int cnt = 0; cnt < n_in_cell; cnt++) { + Double3 offset = {offset_in_cell_dist_(), offset_in_cell_dist_(), + offset_in_cell_dist_()}; + auto pos = + (Double3(boundary_cell_global_idx) + offset) * grid_.domain.dx + + grid_.domain.corner; + auto prt = get_advanced_prt(pos, wni); + + if (prt.x[INJECT_DIM_IDX_] < grid_.domain.corner[INJECT_DIM_IDX_]) { + continue; + } + + injector(prt); + } + } + + template + void inject_into_boundary_patch(Injector& injector, + const Grid_t::Patch& boundary_patch) + { + Int3 ilo = boundary_patch.off; + Int3 ihi = ilo + grid_.ldims; + + assert(ilo[INJECT_DIM_IDX_] == 0); + + int dim1 = std::min((INJECT_DIM_IDX_ + 1) % 3, (INJECT_DIM_IDX_ + 2) % 3); + int dim2 = std::max((INJECT_DIM_IDX_ + 1) % 3, (INJECT_DIM_IDX_ + 2) % 3); + + for (Int3 cell_idx = ilo; cell_idx[dim1] < ihi[dim1]; cell_idx[dim1]++) { + for (cell_idx[dim2] = ilo[dim2]; cell_idx[dim2] < ihi[dim2]; + cell_idx[dim2]++) { + inject_into_boundary_cell(injector, cell_idx); + } + } + } + + template + void inject_into_boundary(Injectors& injectors_by_patch) + { + for (int patch_idx = 0; patch_idx < grid_.n_patches(); patch_idx++) { + const auto& patch = grid_.patches[patch_idx]; + + if (patch.off[INJECT_DIM_IDX_] != 0) { + continue; + } + + auto& injector = injectors_by_patch[patch_idx]; + inject_into_boundary_patch(injector, patch); + } + } + + const Grid_t& grid_; + AdvanceParticle advance_; + SetupParticles setup_particles_; + psc_particle_np np_; + double (*offset_in_cell_dist_)(); + const static int INJECT_DIM_IDX_ = (!Dim::InvarX::value ? 0 + : !Dim::InvarY::value ? 1 + : !Dim::InvarZ::value ? 2 + : -1); +}; + +class TestInjector +{ +public: + void operator()(psc::particle::Inject prt) { prts.push_back(prt); } + + std::vector prts; +}; + +static double half() { return 0.5; } + +TEST(TestSetupParticlesInflow, Advance) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 1, 1}}; + double dt = 10.; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, dt}; + Mparticles mprts{grid}; + + Double3 pos = {5., -2., 5.}; + Double3 u = {0.0, 0.5, 0.0}; + Double3 T = {0.0, 0.0, 0.0}; + psc_particle_npt npt = {0, 1.0, u, T}; + + Inflow inflow(grid, npt, *half); + + auto prt = inflow.get_advanced_prt(pos, 1.0); + + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 2.47214, 1e-5); + EXPECT_NEAR(prt.x[2], 5., 1e-5); +} + +TEST(TestSetupParticlesInflow, InjectIntoCell) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 1, 1}}; + double dt = 10.; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, dt}; + Mparticles mprts{grid}; + + Double3 u = {0.0, 0.8, 0.0}; + Double3 T = {0.0, 0.0, 0.0}; + psc_particle_npt npt = {0, 1.0, u, T}; + + Inflow inflow(grid, npt, *half); + + TestInjector injector; + inflow.inject_into_boundary_cell(injector, {0, 0, 0}); + + EXPECT_EQ(injector.prts.size(), prm.nicell); + + for (auto prt : injector.prts) { + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 1.246950, 1e-5); + EXPECT_NEAR(prt.x[2], 5., 1e-5); + } +} + +TEST(TestSetupParticlesInflow, InjectIntoCellFilter) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 1, 1}}; + double dt = 10.; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, dt}; + Mparticles mprts{grid}; + + // too slow to enter domain + Double3 u = {0.0, 0.5, 0.0}; + Double3 T = {0.0, 0.0, 0.0}; + psc_particle_npt npt = {0, 1.0, u, T}; + + Inflow inflow(grid, npt, *half); + + TestInjector injector; + inflow.inject_into_boundary_cell(injector, {0, 0, 0}); + + EXPECT_EQ(injector.prts.size(), 0); +} + +TEST(TestSetupParticlesInflow, InjectIntoPatch) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 1, 1}}; + double dt = 10.; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, dt}; + Mparticles mprts{grid}; + + Double3 u = {0.0, 0.8, 0.0}; + Double3 T = {0.0, 0.0, 0.0}; + psc_particle_npt npt = {0, 1.0, u, T}; + + Inflow inflow(grid, npt, *half); + + TestInjector injector; + inflow.inject_into_boundary_patch(injector, grid.patches[0]); + + EXPECT_EQ(injector.prts.size(), + domain.ldims[0] * domain.ldims[2] * prm.nicell); + + for (int i = 0; i < injector.prts.size(); i++) { + auto prt = injector.prts[i]; + + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 1.246950, 1e-5); + + if (i < prm.nicell) { + EXPECT_NEAR(prt.x[2], 5., 1e-5); + } else { + EXPECT_NEAR(prt.x[2], 15., 1e-5); + } + } +} + +TEST(TestSetupParticlesInflow, InjectIntoBoundaryY) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 2, 2}}; + double dt = 10.; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, dt}; + Mparticles mprts{grid}; + + Double3 u = {0.0, 0.8, 0.0}; + Double3 T = {0.0, 0.0, 0.0}; + psc_particle_npt npt = {0, 1.0, u, T}; + + Inflow inflow(grid, npt, *half); + + std::vector injectors = {TestInjector(), TestInjector(), + TestInjector(), TestInjector()}; + inflow.inject_into_boundary(injectors); + + EXPECT_EQ(injectors[0].prts.size(), + domain.ldims[0] * domain.ldims[2] * prm.nicell); + for (auto prt : injectors[0].prts) { + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 1.246950, 1e-5); + EXPECT_NEAR(prt.x[2], 5., 1e-5); + } + + EXPECT_EQ(injectors[1].prts.size(), 0); + + EXPECT_EQ(injectors[2].prts.size(), + domain.ldims[0] * domain.ldims[2] * prm.nicell); + for (auto prt : injectors[2].prts) { + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 1.246950, 1e-5); + EXPECT_NEAR(prt.x[2], 15., 1e-5); + } + + EXPECT_EQ(injectors[3].prts.size(), 0); +} + +TEST(TestSetupParticlesInflow, InjectIntoBoundaryZ) +{ + using Mparticles = MparticlesDouble; + + auto domain = Grid_t::Domain{{1, 2, 2}, {10., 20., 20.}, {}, {1, 2, 2}}; + double dt = 10.; + auto kinds = Grid_t::Kinds{{1., 100., "i"}}; + auto prm = Grid_t::NormalizationParams::dimensionless(); + prm.nicell = 2; + Grid_t grid{domain, {}, kinds, {prm}, dt}; + Mparticles mprts{grid}; + + Double3 u = {0.0, 0.0, 0.8}; + Double3 T = {0.0, 0.0, 0.0}; + psc_particle_npt npt = {0, 1.0, u, T}; + + Inflow inflow(grid, npt, *half); + + std::vector injectors = {TestInjector(), TestInjector(), + TestInjector(), TestInjector()}; + inflow.inject_into_boundary(injectors); + + EXPECT_EQ(injectors[0].prts.size(), + domain.ldims[0] * domain.ldims[1] * prm.nicell); + for (auto prt : injectors[0].prts) { + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 5., 1e-5); + EXPECT_NEAR(prt.x[2], 1.246950, 1e-5); + } + + EXPECT_EQ(injectors[1].prts.size(), + domain.ldims[0] * domain.ldims[1] * prm.nicell); + for (auto prt : injectors[1].prts) { + EXPECT_NEAR(prt.x[0], 5., 1e-5); + EXPECT_NEAR(prt.x[1], 15., 1e-5); + EXPECT_NEAR(prt.x[2], 1.246950, 1e-5); + } + + EXPECT_EQ(injectors[2].prts.size(), 0); + + EXPECT_EQ(injectors[3].prts.size(), 0); +} + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + + MPI_Finalize(); + return rc; +}