diff --git a/scripts/run_test.sh b/scripts/run_test.sh new file mode 100755 index 0000000000..c6a5dbf3ee --- /dev/null +++ b/scripts/run_test.sh @@ -0,0 +1,13 @@ +#!/bin/bash + +# Usage: scripts/run_test.sh name_of_test +# Example: scripts/run_test.sh test_injector_boundary_inflow +# Assumes you are in the root psc directory and have a build directory. + +set -e +cd build +make $1 +mkdir -p runs +cp ../bits/adios2cfg.xml runs/ +cd runs +../src/libpsc/tests/$1 "${@:2}" diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx new file mode 100644 index 0000000000..d4831de714 --- /dev/null +++ b/src/include/boundary_injector.hxx @@ -0,0 +1,164 @@ +#pragma once + +#include + +#include "grid.hxx" +#include "rng.hxx" +#include "particle.h" +#include +#include "pushp.hxx" +#include "dim.hxx" +#include "setup_particles.hxx" +#include "../libpsc/psc_push_particles/inc_push.cxx" + +class ParticleGeneratorMaxwellian +{ +public: + using Real = psc::particle::Inject::Real; + using Real3 = psc::particle::Inject::Real3; + + // FIXME would be nice to just pass 1 thing for kind-related info + ParticleGeneratorMaxwellian(int kind_idx, Grid_t::Kind kind, Real3 mean_u, + Real3 temperature) + : kind_idx{kind_idx} + { + for (int d = 0; d < 3; d++) { + Real stdev_u = sqrt(temperature[d] / kind.m); + vdfs[d] = VelocityDistributionFunction{mean_u[d], stdev_u}; + } + } + + psc::particle::Inject get(Real3 min_pos, Real3 pos_range) + { + Real3 x; + for (int d = 0; d < 3; d++) { + x[d] = min_pos[d] + uniform_dist.get() * pos_range[d]; + } + + Real3 u{vdfs[0].get(), vdfs[1].get(), vdfs[2].get()}; + Real w = 1.0; + psc::particle::Tag tag = 0; + + return {x, u, w, kind_idx, tag}; + } + +private: + using VelocityDistributionFunction = rng::Normal; + Vec3 vdfs; + int kind_idx; + rng::Uniform uniform_dist{0.0, 1.0}; +}; + +template +class BoundaryInjector +{ + static const int INJECT_DIM_IDX_ = 1; + +public: + using ParticleGenerator = PARTICLE_GENERATOR; + using PushParticles = PUSH_PARTICLES; + + using Mparticles = typename PushParticles::Mparticles; + using MfieldsState = typename PushParticles::MfieldsState; + using AdvanceParticle_t = typename PushParticles::AdvanceParticle_t; + using Current = typename PushParticles::Current; + using Dim = typename PushParticles::Dim; + using real_t = typename PushParticles::real_t; + using Real3 = typename PushParticles::Real3; + using checks_order = typename PushParticles::checks_order; + + BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) + : partice_generator{particle_generator}, + advance{grid.dt}, + prts_per_unit_density{grid.norm.prts_per_unit_density} + {} + + /// Injects particles at the lower y-bound as if there were a population of + /// particles just beyond the edge. The imaginary particle population has unit + /// density, and individual particles from that population are sampled using + /// the given ParticleGenerator. + /// + /// Some of these limitations may be removed in the future. + void operator()(Mparticles& mprts, MfieldsState& mflds) + { + static_assert(INJECT_DIM_IDX_ == 1, + "only injection at lower bound of y is supported"); + + const Grid_t& grid = mprts.grid(); + auto injectors_by_patch = mprts.injector(); + + Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx); + Current current(grid); + + for (int patch_idx = 0; patch_idx < grid.n_patches(); patch_idx++) { + Int3 ilo = grid.patches[patch_idx].off; + + if (ilo[INJECT_DIM_IDX_] != 0) { + continue; + } + + Int3 ihi = ilo + grid.ldims; + + auto&& injector = injectors_by_patch[patch_idx]; + auto flds = mflds[patch_idx]; + typename Current::fields_t J(flds); + + 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]++) { + cell_idx[INJECT_DIM_IDX_] = -1; + + Real3 cell_corner = + grid.domain.corner + Double3(cell_idx) * grid.domain.dx; + int n_prts_to_try_inject = + get_n_in_cell(1.0, prts_per_unit_density, true); + real_t charge_injected = 0; + + for (int prt_count = 0; prt_count < n_prts_to_try_inject; + prt_count++) { + psc::particle::Inject prt = + partice_generator.get(cell_corner, grid.domain.dx); + + Real3 v = advance.calc_v(prt.u); + Real3 initial_x = prt.x; + advance.push_x(prt.x, v, 1.0); + + if (prt.x[INJECT_DIM_IDX_] < grid.domain.corner[INJECT_DIM_IDX_]) { + // don't inject a particle that fails to enter the domain + continue; + } + + injector(prt); + + // Update currents + // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts() + + Real3 initial_normalized_pos = initial_x * dxi; + + Int3 final_idx; + Real3 final_normalized_pos; + find_idx_pos_1st_rel(prt.x, dxi, final_idx, final_normalized_pos, + real_t(0.)); + + // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt + real_t qni_wni = grid.kinds[prt.kind].q * prt.w; + current.calc_j(J, initial_normalized_pos, final_normalized_pos, + final_idx, cell_idx, qni_wni, v); + + charge_injected += qni_wni; + } + + // set current in boundary to account for injected charge + 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; + } + } + } + } + +private: + ParticleGenerator partice_generator; + AdvanceParticle advance; + real_t prts_per_unit_density; +}; diff --git a/src/include/grid.hxx b/src/include/grid.hxx index a5b8fb0da6..e11cc9ff4c 100644 --- a/src/include/grid.hxx +++ b/src/include/grid.hxx @@ -286,6 +286,7 @@ struct Grid_::Normalization double vt = sqrt(prm.tt / prm.mm); double wp = sqrt(sqr(prm.qq) * prm.n0 / prm.eps0 / prm.mm); + prts_per_unit_density = prm.nicell; cori = 1. / prm.nicell; double alpha_ = wp / wl; beta = vt / cc; @@ -299,6 +300,7 @@ struct Grid_::Normalization real_t eta = {1.}; real_t beta = {1.}; real_t cori = {1.}; + real_t prts_per_unit_density; real_t e0; real_t b0; diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx index d8e95c280b..174d40bc02 100644 --- a/src/include/setup_particles.hxx +++ b/src/include/setup_particles.hxx @@ -4,7 +4,11 @@ #include #include #include + +#include + #include "centering.hxx" +#include "particle.h" #include "rng.hxx" struct psc_particle_npt @@ -24,9 +28,6 @@ struct psc_particle_np psc::particle::Tag tag; }; -// ====================================================================== -// SetupParticles - namespace { const centering::Centerer defaultCenterer(centering::CC); @@ -92,6 +93,28 @@ private: std::function func; }; +// ====================================================================== +// get_n_in_cell +// +// helper function for partition / particle setup + +template +int get_n_in_cell(real_t density, real_t prts_per_unit_density, + bool fractional_n_particles_per_cell) +{ + static rng::Uniform dist{0, 1}; + if (density == 0.0) { + return 0; + } + if (fractional_n_particles_per_cell) { + return density * prts_per_unit_density + dist.get(); + } + return std::max(1, int(density * prts_per_unit_density + .5)); +} + +// ====================================================================== +// SetupParticles + template struct SetupParticles { @@ -111,19 +134,11 @@ struct SetupParticles // ---------------------------------------------------------------------- // get_n_in_cell - // - // helper function for partition / particle setup - int get_n_in_cell(const psc_particle_np& np) + int get_n_in_cell(real_t density) { - static rng::Uniform dist{0, 1}; - if (np.n == 0) { - return 0; - } - if (fractional_n_particles_per_cell) { - return np.n / norm_.cori + dist.get(); - } - return std::max(1, int(np.n / norm_.cori + .5)); + return ::get_n_in_cell(density, real_t(norm_.prts_per_unit_density), + fractional_n_particles_per_cell); } // ---------------------------------------------------------------------- @@ -161,7 +176,7 @@ struct SetupParticles int n_in_cell; if (pop != neutralizing_population) { - n_in_cell = get_n_in_cell(np); + n_in_cell = get_n_in_cell(np.n); n_q_in_cell += kinds_[np.kind].q * n_in_cell; } else { // FIXME, should handle the case where not the last population @@ -180,9 +195,9 @@ struct SetupParticles // setupParticle psc::particle::Inject setupParticle(const psc_particle_np& np, Double3 pos, - double wni) + double weight) { - return psc::particle::Inject{pos, np.p(), wni, np.kind, np.tag}; + return psc::particle::Inject{pos, np.p(), weight, np.kind, np.tag}; } // ---------------------------------------------------------------------- @@ -244,12 +259,12 @@ struct SetupParticles // ---------------------------------------------------------------------- // getWeight - real_t getWeight(const psc_particle_np& np, int n_in_cell) + real_t getWeight(real_t density, int n_in_cell) { if (fractional_n_particles_per_cell) { return 1.; } else { - return np.n / (n_in_cell * norm_.cori); + return density * norm_.prts_per_unit_density / n_in_cell; } } @@ -281,8 +296,8 @@ 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 = getWeight(np, n_in_cell); - auto prt = setupParticle(np, pos, wni); + real_t weight = getWeight(np.n, n_in_cell); + auto prt = setupParticle(np, pos, weight); injector(prt); } }); diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index b8ee7318a7..5bec697cf6 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -472,7 +472,8 @@ struct BndFields_ : BndFieldsBase const int* ldims = mflds.grid().ldims; Int3 ib = mflds.ib(), im = mflds.im(); - assert(0); + // TODO + // assert(0); #if 0 if (d == 1) { @@ -554,7 +555,8 @@ struct BndFields_ : BndFieldsBase const int* ldims = mflds.grid().ldims; Int3 ib = mflds.ib(), im = mflds.im(); - assert(0); + // TODO + // assert(0); #if 0 if (d == 1) { int my _mrc_unused = ldims[1]; diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index 85dc83ad32..2eb70759ac 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/src/libpsc/psc_checks/checks_impl.hxx @@ -82,6 +82,7 @@ public: MPI_Allreduce(&local_err, &max_err, 1, MPI_DOUBLE, MPI_MAX, grid.comm()); if (should_print_diffs(max_err)) { + mpi_printf(grid.comm(), "continuity: drho -- -dt*div j\n"); psc::helper::print_diff(d_rho, -dt_divj, err_threshold); } @@ -160,6 +161,7 @@ public: max_err = std::max(max_err, patch_err); if (should_print_diffs(patch_err)) { + mpi_printf(grid.comm(), "gauss: rho -- div E\n"); psc::helper::print_diff(patch_rho, patch_dive, err_threshold); } } diff --git a/src/libpsc/psc_push_particles/inc_curr.c b/src/libpsc/psc_push_particles/inc_curr.cxx similarity index 87% rename from src/libpsc/psc_push_particles/inc_curr.c rename to src/libpsc/psc_push_particles/inc_curr.cxx index 10ae07edf7..781335e05e 100644 --- a/src/libpsc/psc_push_particles/inc_curr.c +++ b/src/libpsc/psc_push_particles/inc_curr.cxx @@ -25,10 +25,10 @@ GT_INLINE void deposit(Curr& curr, const int _i[3], const real_t fnqs[3], deposition(curr, i, qni_wni, dx, xa, dim_xyz{}); } -#include "inc_curr_1vb_split.c" -#include "inc_curr_1vb_var1.c" -#include "inc_curr_1vb_2d.c" -#include "inc_curr_zigzag.c" +#include "inc_curr_1vb_split.cxx" +#include "inc_curr_1vb_var1.cxx" +#include "inc_curr_1vb_2d.cxx" +#include "inc_curr_zigzag.cxx" template struct Current1vb; diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_2d.c b/src/libpsc/psc_push_particles/inc_curr_1vb_2d.cxx similarity index 100% rename from src/libpsc/psc_push_particles/inc_curr_1vb_2d.c rename to src/libpsc/psc_push_particles/inc_curr_1vb_2d.cxx diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.c b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx similarity index 100% rename from src/libpsc/psc_push_particles/inc_curr_1vb_split.c rename to src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_var1.c b/src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx similarity index 100% rename from src/libpsc/psc_push_particles/inc_curr_1vb_var1.c rename to src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx diff --git a/src/libpsc/psc_push_particles/inc_curr_zigzag.c b/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx similarity index 100% rename from src/libpsc/psc_push_particles/inc_curr_zigzag.c rename to src/libpsc/psc_push_particles/inc_curr_zigzag.cxx diff --git a/src/libpsc/psc_push_particles/inc_push.c b/src/libpsc/psc_push_particles/inc_push.cxx similarity index 100% rename from src/libpsc/psc_push_particles/inc_push.c rename to src/libpsc/psc_push_particles/inc_push.cxx diff --git a/src/libpsc/psc_push_particles/push_config.hxx b/src/libpsc/psc_push_particles/push_config.hxx index 916fb469ff..32ff205716 100644 --- a/src/libpsc/psc_push_particles/push_config.hxx +++ b/src/libpsc/psc_push_particles/push_config.hxx @@ -6,8 +6,8 @@ #include "inc_defs.h" #include "interpolate.hxx" -#include "inc_curr.c" -#include "inc_push.c" +#include "inc_curr.cxx" +#include "inc_push.cxx" #include "push_particles.hxx" #include "push_particles_esirkepov.hxx" #include "push_particles_1vb.hxx" diff --git a/src/libpsc/tests/CMakeLists.txt b/src/libpsc/tests/CMakeLists.txt index 3b6b895339..e46e74a4b9 100644 --- a/src/libpsc/tests/CMakeLists.txt +++ b/src/libpsc/tests/CMakeLists.txt @@ -28,7 +28,7 @@ add_psc_test(test_reflective_bcs_integration) 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_boundary_injector) 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 deleted file mode 100644 index ff16293403..0000000000 --- a/src/libpsc/tests/test_bnd_prt_inflow.cxx +++ /dev/null @@ -1,424 +0,0 @@ - -#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; -} diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx new file mode 100644 index 0000000000..36508d96fd --- /dev/null +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -0,0 +1,245 @@ +#include + +#include "test_common.hxx" + +#include "boundary_injector.hxx" + +#include "psc.hxx" +#include "DiagnosticsDefault.h" +#include "OutputFieldsDefault.h" +#include "../psc_config.hxx" + +TEST(BoundaryInjectorTest, ParticleGeneratorMaxwellianTest) +{ + int kind_idx = 15; + Grid_t::Kind kind{1.0, 1836.0, "ion"}; + ParticleGeneratorMaxwellian::Real w = 1.0; + ParticleGeneratorMaxwellian::Real3 mean_u{0.0, 5.0, 15.0}; + ParticleGeneratorMaxwellian::Real3 temperature{0.0, 0.0, 0.0}; + ParticleGeneratorMaxwellian::Real3 pos{1.0, 2.0, 5.0}; + + ParticleGeneratorMaxwellian gen{kind_idx, kind, mean_u, temperature}; + + auto prt = gen.get(pos, {0.0, 0.0, 0.0}); + + ASSERT_EQ(prt.kind, kind_idx); + ASSERT_EQ(prt.w, w); + ASSERT_EQ(prt.tag, 0); + ASSERT_EQ(prt.u, mean_u); // zero temperature => exact velocity + ASSERT_EQ(prt.x, pos); +} + +// ====================================================================== + +using Dim = dim_yz; +using PscConfig = PscConfig1vbecDouble; + +using MfieldsState = PscConfig::MfieldsState; +using Mparticles = PscConfig::Mparticles; +using Balance = PscConfig::Balance; +using Collision = PscConfig::Collision; +using Checks = PscConfig::Checks; +using Marder = PscConfig::Marder; +using OutputParticles = PscConfig::OutputParticles; + +Grid_t* setupGrid() +{ + auto domain = Grid_t::Domain{{1, 8, 2}, // n grid points + {10.0, 80.0, 20.0}, // physical lengths + {0, 0, 0}, // location of lower corner + {1, 1, 1}}; // n patches + + auto bc = + // FIXME wrong BCs + psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_OPEN, BND_FLD_PERIODIC}, + {BND_FLD_PERIODIC, BND_FLD_OPEN, BND_FLD_PERIODIC}, + {BND_PRT_PERIODIC, BND_PRT_OPEN, BND_PRT_PERIODIC}, + {BND_PRT_PERIODIC, BND_PRT_OPEN, BND_PRT_PERIODIC}}; + + auto kinds = Grid_t::Kinds(NR_KINDS); + kinds[KIND_ELECTRON] = {-1.0, 1.0, "e"}; + kinds[KIND_ION] = {1.0, 1.0, "i"}; + + // --- generic setup + auto norm_params = Grid_t::NormalizationParams::dimensionless(); + norm_params.nicell = 2; + + double dt = .1; + Grid_t::Normalization norm{norm_params}; + + Int3 ibn = {2, 2, 2}; + if (Dim::InvarX::value) { + ibn[0] = 0; + } + if (Dim::InvarY::value) { + ibn[1] = 0; + } + if (Dim::InvarZ::value) { + ibn[2] = 0; + } + + 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; + + psc::particle::Inject get(Real3 min_pos, Real3 pos_range) + { + Real uy = 2.0; + if (PRT_COUNT > 0 && n_injected++ >= PRT_COUNT) { + // Setting uy=0 means the particle won't enter the domain, and thus won't + // be injected + uy = 0.0; + } + + 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; +}; + +TEST(BoundaryInjectorTest, Integration1Particle) +{ + // ---------------------------------------------------------------------- + // 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_particles = + BoundaryInjector, PscConfig::PushParticles>{{}, grid}; + + auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, + balance, collision, checks, marder, + diagnostics, inject_particles); + + // ---------------------------------------------------------------------- + // 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); + } + + ASSERT_EQ(prts.size(), 1); +} + +TEST(BoundaryInjectorTest, IntegrationManyParticles) +{ + // ---------------------------------------------------------------------- + // 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_particles = + BoundaryInjector, PscConfig::PushParticles>{{}, grid}; + + auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, + balance, collision, checks, marder, + diagnostics, inject_particles); + + // ---------------------------------------------------------------------- + // 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); + } + + ASSERT_GT(prts.size(), 1); +} + +// ====================================================================== +// main + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + MPI_Finalize(); + return rc; +} diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 2a321e0666..18cb4206b1 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -22,7 +22,6 @@ using PscConfig = PscConfig1vbecDouble; // ---------------------------------------------------------------------- -using BgkMfields = PscConfig::Mfields; using MfieldsState = PscConfig::MfieldsState; using Mparticles = PscConfig::Mparticles; using Balance = PscConfig::Balance; diff --git a/src/libpsc/tests/test_reflective_bcs_integration.cxx b/src/libpsc/tests/test_reflective_bcs_integration.cxx index 44d4901d03..8519c5690a 100644 --- a/src/libpsc/tests/test_reflective_bcs_integration.cxx +++ b/src/libpsc/tests/test_reflective_bcs_integration.cxx @@ -22,7 +22,6 @@ using PscConfig = PscConfig1vbecDouble; // ---------------------------------------------------------------------- -using BgkMfields = PscConfig::Mfields; using MfieldsState = PscConfig::MfieldsState; using Mparticles = PscConfig::Mparticles; using Balance = PscConfig::Balance;