diff --git a/scripts/run_test.sh b/scripts/run_test.sh index c6a5dbf3ee..4d464d30f0 100755 --- a/scripts/run_test.sh +++ b/scripts/run_test.sh @@ -1,6 +1,6 @@ #!/bin/bash -# Usage: scripts/run_test.sh name_of_test +# Usage: scripts/run_test.sh name_of_test [gtest args] # Example: scripts/run_test.sh test_injector_boundary_inflow # Assumes you are in the root psc directory and have a build directory. @@ -10,4 +10,4 @@ make $1 mkdir -p runs cp ../bits/adios2cfg.xml runs/ cd runs -../src/libpsc/tests/$1 "${@:2}" +mpirun -np 1 ../src/libpsc/tests/$1 "${@:2}" diff --git a/src/include/bnd_particles.hxx b/src/include/bnd_particles.hxx deleted file mode 100644 index a578670e2a..0000000000 --- a/src/include/bnd_particles.hxx +++ /dev/null @@ -1,10 +0,0 @@ - -#pragma once - -#include "particles.hxx" - -// ====================================================================== -// BndParticlesBase - -struct BndParticlesBase -{}; diff --git a/src/include/bnd_particles_impl.hxx b/src/include/bnd_particles_impl.hxx index 653b3be541..abf18020ad 100644 --- a/src/include/bnd_particles_impl.hxx +++ b/src/include/bnd_particles_impl.hxx @@ -7,7 +7,7 @@ #include "balance.hxx" #include "ddc_particles.hxx" -#include "bnd_particles.hxx" +#include "particles.hxx" #include "particles_simple.hxx" extern int pr_time_step_no_comm; @@ -17,10 +17,11 @@ extern double* psc_balance_comp_time_by_patch; // BndParticlesCommon template -struct BndParticlesCommon : BndParticlesBase +struct BndParticlesCommon { using Mparticles = MP; using real_t = typename Mparticles::real_t; + using Real3 = Vec3; using ddcp_t = ddc_particles; using BndBuffers = typename Mparticles::BndBuffers; @@ -96,12 +97,9 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, // New-style boundary requirements. // These will need revisiting when it comes to non-periodic domains. + const Int3& ldims = grid.ldims; const auto& gpatch = grid.patches[p]; - const Int3& ldims = pi.ldims(); - real_t xm[3]; - for (int d = 0; d < 3; d++) { - xm[d] = gpatch.xe[d] - gpatch.xb[d]; - } + Real3 patch_size = Real3(gpatch.xe - gpatch.xb); auto* dpatch = &ddcp->patches_[p]; for (int dir1 = 0; dir1 < N_DIR; dir1++) { @@ -121,12 +119,9 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, Int3 pos = pi.cellPosition(xi); - if (pos[0] >= 0 && - pos[0] < ldims[0] && // OPT, could be optimized with casts to unsigned - pos[1] >= 0 && pos[1] < ldims[1] && pos[2] >= 0 && pos[2] < ldims[2]) { + if (pi.isValidCellPosition(pos)) { // fast path // particle is still inside patch: move into right position - pi.validCellIndex(xi); buf[head++] = *prt; continue; } @@ -139,7 +134,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, for (int d = 0; d < 3; d++) { if (pos[d] < 0) { if (!grid.atBoundaryLo(p, d) || grid.bc.prt_lo[d] == BND_PRT_PERIODIC) { - xi[d] += xm[d]; + xi[d] += patch_size[d]; dir[d] = -1; int ci = pi.cellPosition(xi[d], d); if (ci >= ldims[d]) { @@ -148,18 +143,26 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, } } else { switch (grid.bc.prt_lo[d]) { - case BND_PRT_REFLECTING: + case BND_PRT_REFLECTING: { xi[d] = -xi[d]; pxi[d] = -pxi[d]; dir[d] = 0; break; - case BND_PRT_ABSORBING: drop = true; break; + } + case BND_PRT_OPEN: { + drop = true; + break; + } + case BND_PRT_ABSORBING: { + drop = true; + break; + } default: assert(0); } } } else if (pos[d] >= ldims[d]) { if (!grid.atBoundaryHi(p, d) || grid.bc.prt_hi[d] == BND_PRT_PERIODIC) { - xi[d] -= xm[d]; + xi[d] -= patch_size[d]; dir[d] = +1; int ci = pi.cellPosition(xi[d], d); if (ci < 0) { @@ -168,7 +171,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, } else { switch (grid.bc.prt_hi[d]) { case BND_PRT_REFLECTING: { - xi[d] = 2.f * xm[d] - xi[d]; + xi[d] = 2.f * patch_size[d] - xi[d]; pxi[d] = -pxi[d]; dir[d] = 0; int ci = pi.cellPosition(xi[d], d); @@ -177,7 +180,14 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, } break; } - case BND_PRT_ABSORBING: drop = true; break; + case BND_PRT_OPEN: { + drop = true; + break; + } + case BND_PRT_ABSORBING: { + drop = true; + break; + } default: assert(0); } } @@ -191,7 +201,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, xi[d] = 0.f; } assert(xi[d] >= 0.f); - assert(xi[d] <= xm[d]); + assert(xi[d] <= patch_size[d]); } } if (!drop) { diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index c7b0258e68..72f9dc4eb6 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -9,6 +9,7 @@ #include "pushp.hxx" #include "dim.hxx" #include "setup_particles.hxx" +#include "kg/VecRange.hxx" #include "../libpsc/psc_push_particles/inc_push.cxx" /// @brief A particle generator for use with @ref BoundaryInjector. Samples @@ -72,15 +73,13 @@ public: 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; + using Real3 = Vec3; BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) - : partice_generator{particle_generator}, - advance{grid.dt}, - prts_per_unit_density{grid.norm.prts_per_unit_density} + : particle_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 @@ -97,84 +96,71 @@ public: const Grid_t& grid = mprts.grid(); auto injectors_by_patch = mprts.injector(); - Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx); + Real3 dxi = grid.domain.dx_inv; 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) { + if (!grid.atBoundaryLo(patch_idx, INJECT_DIM_IDX_)) { continue; } - Int3 ihi = ilo + grid.ldims; + Int3 ilo = {0, 0, 0}; + Int3 ihi = grid.ldims; + + ilo[INJECT_DIM_IDX_] = -1; + ihi[INJECT_DIM_IDX_] = 0; 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]++) { - 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 = - 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() + for (Int3 initial_idx : VecRange(ilo, ihi)) { + Real3 cell_corner = Double3(initial_idx) * grid.domain.dx; + int n_prts_to_try_inject = + get_n_in_cell(1.0, prts_per_unit_density_, true); + + for (int prt_count = 0; prt_count < n_prts_to_try_inject; prt_count++) { + // FIXME #948112531345 (also see the other FIXMEs with this id) + // Depositing current in ghost corners leads to false-positive gauss + // errors. This could be avoided here by artificially constraining the + // trajectories of injected particles. However, assuming the particle + // boundary condition is "open", outflowing particles cause the same + // problem, and there's nothing to be done about that. + psc::particle::Inject prt = + particle_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); + + if (prt.x[INJECT_DIM_IDX_] < 0.0) { + // don't inject a particle that fails to enter the patch + continue; + } - Real3 initial_normalized_pos = initial_x * dxi; + injector(prt); - Int3 final_idx; - Real3 final_normalized_pos; - find_idx_pos_1st_rel(prt.x, dxi, final_idx, final_normalized_pos, - real_t(0.)); + // Update currents + // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts() - // 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); + Real3 initial_normalized_pos = initial_x * dxi; + // pretend it came from the edge to inject proper current + initial_normalized_pos[INJECT_DIM_IDX_] = -1.0; - charge_injected += qni_wni; - } + Real3 final_normalized_pos = prt.x * dxi; + Int3 final_idx = final_normalized_pos.fint(); - // 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) = - boundary_current_before + charge_injected * grid.domain.dx[1] / - grid.dt / prts_per_unit_density; + // 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, initial_idx, qni_wni, v); } } } } private: - ParticleGenerator partice_generator; - AdvanceParticle advance; - real_t prts_per_unit_density; + ParticleGenerator particle_generator_; + AdvanceParticle_t advance_; + real_t prts_per_unit_density_; }; diff --git a/src/include/grid/domain.hxx b/src/include/grid/domain.hxx index ce7f20f6a2..2e4661c634 100644 --- a/src/include/grid/domain.hxx +++ b/src/include/grid/domain.hxx @@ -43,6 +43,7 @@ struct Domain ldims = gdims / np; dx = length / Real3(gdims); + dx_inv = Real3(gdims) / length; } void view() const @@ -65,6 +66,8 @@ struct Domain Int3 ldims; /// @brief Side lengths of each grid cell, in physical units. Real3 dx; + /// @brief Inverses of side lengths of each grid cell, in physical units. + Real3 dx_inv; }; template diff --git a/src/include/particle_indexer.hxx b/src/include/particle_indexer.hxx index e8e7b84537..4efd41eaaf 100644 --- a/src/include/particle_indexer.hxx +++ b/src/include/particle_indexer.hxx @@ -66,7 +66,7 @@ struct ParticleIndexer using Real3 = Vec3; ParticleIndexer(const Grid_t& grid) - : dxi_(Real3{1., 1., 1.} / Real3(grid.domain.dx)), ldims_(grid.ldims) + : dxi_(grid.domain.dx_inv), ldims_(grid.ldims) { n_cells_ = ldims_[0] * ldims_[1] * ldims_[2]; } @@ -98,10 +98,25 @@ struct ParticleIndexer return cellIndex(cpos); } + bool isValidCellPosition(const Int3& cpos) const + { + for (int d = 0; d < 3; d++) { + // optimization trick: if cpos[d]<0, it becomes larger than any positive + // int after casting both to uint + // TODO (cursed?) store ldims as uints, so this happens implicitly + if (uint(cpos[d]) >= ldims_[d]) { + return false; + } + } + return true; + } + int validCellIndex(const real_t* pos) const { Int3 cpos = cellPosition(pos); for (int d = 0; d < 3; d++) { + // optimization trick: if cpos[d]<0, it becomes larger than any positive + // int after casting both to uint if (uint(cpos[d]) >= ldims_[d]) { printf("validCellIndex: cpos[%d] = %d ldims_[%d] = %d // pos[%d] = %g " "pos[%d]*dxi_[%d] = %g\n", diff --git a/src/libpsc/integrate.cxx b/src/libpsc/integrate.cxx index 6f67ba3b91..92430a1b38 100644 --- a/src/libpsc/integrate.cxx +++ b/src/libpsc/integrate.cxx @@ -5,7 +5,6 @@ #include "push_particles.hxx" #include "sort.hxx" #include "collision.hxx" -#include "bnd_particles.hxx" #include "checks_params.hxx" #include diff --git a/src/libpsc/psc_balance/psc_balance_impl.hxx b/src/libpsc/psc_balance/psc_balance_impl.hxx index 78164a4457..2799464ed8 100644 --- a/src/libpsc/psc_balance/psc_balance_impl.hxx +++ b/src/libpsc/psc_balance/psc_balance_impl.hxx @@ -4,7 +4,6 @@ #include "fields.hxx" #include "fields_traits.hxx" #include "balance.hxx" -#include "bnd_particles.hxx" #include "bnd.hxx" #include "mpi_dtype_traits.hxx" 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 1c670e5aaa..2e8c18516c 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -469,18 +469,22 @@ struct BndFields_ : BndFieldsBase void open_H_lo(Mfields& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); + const Grid_t& grid = mflds.grid(); + Int3 ldims = grid.ldims; + Int3 ib = mflds.ib(); + Int3 im = mflds.im(); + real_t dt = grid.dt; + real_t dx = grid.domain.dx[0]; + real_t dy = grid.domain.dx[1]; + real_t dz = grid.domain.dx[2]; - // TODO - // assert(0); -#if 0 +#if 1 if (d == 1) { #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { fields_t_set_nan(&F(HX, ix, -1, iz)); fields_t_set_nan(&F(HX, ix, -2, iz)); fields_t_set_nan(&F(HY, ix, -1, iz)); @@ -490,11 +494,10 @@ struct BndFields_ : BndFieldsBase } } #endif - real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], - dz = F.grid().domain.dx[2]; + for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HX, ix, -1, iz) = (/* + 4.f * C_s_pulse_y1(x,y,z+0.5*dz,t) */ -2.f * F(EZ, ix, 0, iz) @@ -512,8 +515,8 @@ struct BndFields_ : BndFieldsBase } else if (d == 2) { #ifdef DEBUG for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { fields_t_set_nan(&F(HX, ix, iy, -1)); fields_t_set_nan(&F(HX, ix, iy, -2)); fields_t_set_nan(&F(HY, ix, iy, -1)); @@ -523,11 +526,9 @@ struct BndFields_ : BndFieldsBase } } #endif - real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], - dz = F.grid().domain.dx[2]; for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HY, ix, iy, -1) = (/* + 4.f * C_s_pulse_z1(x+0.5*dx,y,z,t) */ -2.f * F(EX, ix, iy, 0) - @@ -552,18 +553,24 @@ struct BndFields_ : BndFieldsBase void open_H_hi(Mfields& mflds, int p, int d) { auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); + const Grid_t& grid = mflds.grid(); + Int3 ldims = grid.ldims; + Int3 ib = mflds.ib(); + Int3 im = mflds.im(); + real_t dt = grid.dt; + real_t dx = grid.domain.dx[0]; + real_t dy = grid.domain.dx[1]; + real_t dz = grid.domain.dx[2]; // TODO // assert(0); -#if 0 +#if 1 if (d == 1) { int my _mrc_unused = ldims[1]; #ifdef DEBUG for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { fields_t_set_nan(&F(HX, ix, my, iz)); fields_t_set_nan(&F(HX, ix, my + 1, iz)); fields_t_set_nan(&F(HY, ix, my, iz)); @@ -572,11 +579,9 @@ struct BndFields_ : BndFieldsBase } } #endif - real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], - dz = F.grid().domain.dx[2]; for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HX, ix, my, iz) = (/* + 4.f * C_s_pulse_y2(x,y,z+0.5*dz,t) */ +2.f * F(EZ, ix, my, iz) @@ -606,11 +611,9 @@ struct BndFields_ : BndFieldsBase } } #endif - real_t dt = F.grid().dt, dy = F.grid().domain.dx[1], - dz = F.grid().domain.dx[2]; for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, F.ib_[0]); - ix < std::min(ldims[0] + 2, F.ib_[0] + F.im_[0]); ix++) { + for (int ix = std::max(-2, ib[0]); + ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { F(HY, ix, iy, mz) = (/* - 4.f * C_s_pulse_z2(x+0.5*dx,y,z,t) */ +2.f * F(EX, ix, iy, mz) + @@ -637,14 +640,16 @@ struct BndFields_ : BndFieldsBase // ====================================================================== // BndFieldsNone -// FIXME, this is mostly useful at most for testing and maybe should go away +// used by CUDA template struct BndFieldsNone : BndFieldsBase { using Mfields = MF; - void fill_ghosts_E(Mfields& mflds){}; - void fill_ghosts_H(Mfields& mflds){}; - void add_ghosts_J(Mfields& mflds){}; + // clang-format off + void fill_ghosts_E(Mfields& mflds) {}; + void fill_ghosts_H(Mfields& mflds) {}; + void add_ghosts_J(Mfields& mflds) {}; + // clang-format on }; diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx index c6af987484..00ef3cc583 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx @@ -18,7 +18,7 @@ struct Current1vbSplit Current1vbSplit(const Grid_t& grid) : dt_(grid.dt), - dxi_{grid.domain.dx.inv()}, + dxi_{grid.domain.dx_inv}, deposition_(real_t(grid.norm.fnqs / grid.dt) * Real3(grid.domain.dx)) {} diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx index 78d5034cc9..0c13effb7b 100644 --- a/src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx @@ -19,8 +19,7 @@ struct Current1vbVar1 using real_t = typename fields_t::real_t; using Real3 = Vec3; - Current1vbVar1(const Grid_t& grid) - : dt_(grid.dt), dxi_{Real3{1., 1., 1.} / Real3(grid.domain.dx)} + Current1vbVar1(const Grid_t& grid) : dt_(grid.dt), dxi_{grid.domain.dx_inv} { fnqxs_ = grid.domain.dx[0] * grid.norm.fnqs / grid.dt; fnqys_ = grid.domain.dx[1] * grid.norm.fnqs / grid.dt; diff --git a/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx b/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx index 51e07eb0e4..f6cd6b4334 100644 --- a/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx +++ b/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx @@ -19,7 +19,7 @@ struct CurrentZigzag CurrentZigzag(const Grid_t& grid) : dt_(grid.dt), - dxi_{Real3{1., 1., 1.} / Real3(grid.domain.dx)}, + dxi_{grid.domain.dx_inv}, deposition_(real_t(grid.norm.fnqs / grid.dt) * Real3(grid.domain.dx)) {} diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index d2c0803003..f1e96b2110 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -2,6 +2,50 @@ // J. Villasenor and O. Buneman, "Rigorous charge conservation for local // electromagnetic field solvers", Computer Physics Communications 69 (1992) 306 +namespace +{ +/** + * @brief "Pushes" an exiting particle to the outer edge of the first layer of + * ghost cells. This should happen before current deposition, so the full + * exiting current is deposited. Note that higher order particles would need to + * be pushed deeper into the ghost region. + * + * Actually, this function works with the grid-normalized particle position, not + * the particle itself (or even its actual position). The grid-normalized + * position is what's actually used to deposit current, and the particle is + * dropped later, so its position doesn't need to be updated. + * @tparam Real `float` or `double` + * @param final_x_normed the particle's grid-normalized position, which is + * possibly mutated + * @param final_i3 the particle's 3d cell index + * @param grid the grid + * @param p the patch index + */ +template +void exit_to_edge(Vec3& final_x_normed, const Int3& final_i3, + const Grid_t& grid, int p) +{ + // FIXME #948112531345 (also see the other FIXMEs with this id) + // Current deposited in ghost corners isn't sent to other patches, and + // thus leads to false-positive gauss errors. + // This is a general problem with non-periodic boundaries that have currents + // in ghost cells (i.e., just open BCs as of now), not just this function. + for (int d = 0; d < 3; d++) { + if (grid.bc.prt_lo[d] == BND_PRT_OPEN && grid.atBoundaryLo(p, d)) { + if (final_i3[d] < 0) { + final_x_normed[d] = -1.0; + } + } + + if (grid.bc.prt_hi[d] == BND_PRT_OPEN && grid.atBoundaryHi(p, d)) { + if (final_i3[d] >= grid.ldims[d]) { + final_x_normed[d] = grid.ldims[d] + 1.0; + } + } + } +} +} // namespace + // ====================================================================== // PushParticlesVb @@ -66,19 +110,21 @@ struct PushParticlesVb Real3 final_pos_normalized = prt.x() * dxi; Int3 final_index = final_pos_normalized.fint(); + exit_to_edge(final_pos_normalized, final_index, grid, p); + // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt - int lg[3]; + Int3 initial_index; if (!Dim::InvarX::value) { - lg[0] = ip.cx.g.l; + initial_index[0] = ip.cx.g.l; } if (!Dim::InvarY::value) { - lg[1] = ip.cy.g.l; + initial_index[1] = ip.cy.g.l; } if (!Dim::InvarZ::value) { - lg[2] = ip.cz.g.l; + initial_index[2] = ip.cz.g.l; } current.calc_j(J, initial_pos_normalized, final_pos_normalized, - final_index, lg, prt.qni_wni(), v); + final_index, initial_index, prt.qni_wni(), v); } } } @@ -89,7 +135,7 @@ struct PushParticlesVb static void stagger_mprts_patch(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx); + Real3 dxi = grid.domain.dx_inv; real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; assert(kinds.size() <= MAX_NR_KINDS); diff --git a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx index 58285a867e..22617b2806 100644 --- a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx +++ b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx @@ -8,6 +8,30 @@ #include "pushp_current_esirkepov.hxx" #include "../libpsc/psc_checks/checks_impl.hxx" +namespace +{ +/** + * @brief Handle exiting particles. When a particle crosses an open boundary, it + * is dropped. Before that happens, an exiting current needs to be deposited, or + * else virtual charge will accumulate at the boundary. + * + * This function currently asserts that there are no open boundaries (and is + * thus a noop in release builds). See `push_particles_1vb.hxx` for an actual + * implementation for 1st-order particles. + * @tparam Real 'float' or 'double' + * @param grid the grid + * @param p the patch index + */ +template +void handle_exiting(const Grid_& grid, int p) +{ + for (int d = 0; d < 3; d++) { + assert(grid.bc.prt_lo[d] != BND_PRT_OPEN); + assert(grid.bc.prt_hi[d] != BND_PRT_OPEN); + } +} +} // namespace + // ====================================================================== // PushParticlesEsirkepov @@ -36,7 +60,7 @@ struct PushParticlesEsirkepov static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi = Real3{1., 1., 1.} / Real3(grid.domain.dx); + Real3 dxi = Real3(grid.domain.dx_inv); // need to explicitly convert dtype real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; assert(kinds.size() <= MAX_NR_KINDS); @@ -86,6 +110,8 @@ struct PushParticlesEsirkepov // CURRENT DENSITY AT (n+1.0)*dt current.prep(prt.qni_wni(), v); current.calc(J); + + handle_exiting(grid, p); } } } diff --git a/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx b/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx index 8e611ebc07..76f01174a5 100644 --- a/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx +++ b/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx @@ -177,7 +177,7 @@ struct CurrentEsirkepov using fields_t = Fields; CurrentEsirkepov(const Grid_t& grid) - : dxi_{Real3{1., 1., 1.} / Real3(grid.domain.dx)}, fnqs_(grid.norm.fnqs) + : dxi_{grid.domain.dx_inv}, fnqs_(grid.norm.fnqs) { fnqxs_ = grid.domain.dx[0] * fnqs_ / grid.dt; fnqys_ = grid.domain.dx[1] * fnqs_ / grid.dt; diff --git a/src/libpsc/tests/CMakeLists.txt b/src/libpsc/tests/CMakeLists.txt index e46e74a4b9..6ed4f28823 100644 --- a/src/libpsc/tests/CMakeLists.txt +++ b/src/libpsc/tests/CMakeLists.txt @@ -25,6 +25,7 @@ add_psc_test(test_mparticles) add_psc_test(test_output_particles) add_psc_test(test_reflective_bcs) add_psc_test(test_reflective_bcs_integration) +add_psc_test(test_open_bcs_integration) add_psc_test(test_mfields) add_psc_test(test_mfields_cuda) add_psc_test(test_bnd) diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 44e540fc19..ef51fd70f1 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -45,10 +45,10 @@ 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 domain = Grid_t::Domain{{1, 8, 2}, // n grid points + {1, 8, 2}, // physical lengths + {0, 0, 0}, // location of lower corner + {1, 1, 1}}; // n patches auto bc = // FIXME wrong BCs @@ -63,9 +63,9 @@ Grid_t* setupGrid() // --- generic setup auto norm_params = Grid_t::NormalizationParams::dimensionless(); - norm_params.nicell = 2; + norm_params.nicell = 1; - double dt = .1; + double dt = 1; Grid_t::Normalization norm{norm_params}; int n_ghosts = 2; @@ -92,6 +92,12 @@ struct ParticleGenerator uy = 0.0; } + // FIXME #948112531345 (also see the other FIXMEs with this id) + // If x[2] != 0, the injected particle deposits some + // current in the cells above its path. If one of those cells is in a ghost + // corner, the current deposited in that cell IS NOT sent across patches. + // That current is necessary to compute div E correctly. The run itself + // should be fine; it only throws off the gauss check. Real3 x = min_pos + pos_range * Real3{0, .999, 0.}; Real3 u{0.0, uy, 0.0}; Real w = 1.0; @@ -112,7 +118,7 @@ TEST(BoundaryInjectorTest, Integration1Particle) PscParams psc_params; - psc_params.nmax = 1; + psc_params.nmax = 2; psc_params.stats_every = 1; psc_params.cfl = .75; @@ -161,8 +167,8 @@ TEST(BoundaryInjectorTest, Integration1Particle) for (; grid.timestep_ < psc_params.nmax;) { 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_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); } ASSERT_EQ(prts.size(), 1); @@ -175,7 +181,7 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) PscParams psc_params; - psc_params.nmax = 1; + psc_params.nmax = 2; psc_params.stats_every = 1; psc_params.cfl = .75; @@ -224,8 +230,8 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) for (; grid.timestep_ < psc_params.nmax;) { 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_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); } ASSERT_GT(prts.size(), 1); @@ -238,7 +244,7 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) PscParams psc_params; - psc_params.nmax = 1; + psc_params.nmax = 2; psc_params.stats_every = 1; psc_params.cfl = .75; @@ -292,8 +298,8 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) for (; grid.timestep_ < psc_params.nmax;) { 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_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); } bool found_electrons = false; diff --git a/src/libpsc/tests/test_open_bcs_integration.cxx b/src/libpsc/tests/test_open_bcs_integration.cxx new file mode 100644 index 0000000000..92e36fcde3 --- /dev/null +++ b/src/libpsc/tests/test_open_bcs_integration.cxx @@ -0,0 +1,170 @@ +#include "gtest/gtest.h" + +#include +#include +#include + +#include "DiagnosticsDefault.h" +#include "OutputFieldsDefault.h" +#include "add_ghosts_reflecting.hxx" +#include "../psc_config.hxx" + +// ====================================================================== +// PSC configuration +// +// This sets up compile-time configuration for the code, in particular +// what data structures and algorithms to use +// +// EDIT to change order / floating point type / cuda / 2d/3d + +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; + +// ====================================================================== +// Global parameters + +namespace +{ +PscParams psc_params; +} // namespace + +// ====================================================================== +// setupGrid + +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 = psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_OPEN, BND_FLD_OPEN}, + {BND_FLD_PERIODIC, BND_FLD_OPEN, BND_FLD_OPEN}, + {BND_PRT_PERIODIC, BND_PRT_OPEN, BND_PRT_OPEN}, + {BND_PRT_PERIODIC, BND_PRT_OPEN, BND_PRT_OPEN}}; + + 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 = 1; + + double dt = psc_params.cfl * courant_length(domain); + 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}; +} + +// ====================================================================== +// Start with 2 particles, send them flying off towards the boundary, and demand +// that: +// a) they get dropped +// b) there are no continuity or gauss errors in the fields + +TEST(OpenBcsTest, IntegrationY) +{ + // ---------------------------------------------------------------------- + // setup + + psc_params.nmax = 20; + 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 psc = + makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, + collision, checks, marder, diagnostics); + + // ---------------------------------------------------------------------- + // set up initial conditions + + EXPECT_EQ(grid.n_patches(), 1); + int p = 0; + + { + auto injector = mprts.injector(); + auto inj = injector[p]; + double vy = 100; // fast enough to escape electric attraction + double y_quarter = grid.domain.length[1] / 4.0; + // inject 2 particles at same location to satisfy Gauss' law at t=0 + inj({{0, y_quarter, 5}, {0, -vy, 0}, 1, 0}); // electron + inj({{0, y_quarter, 5}, {0, vy, 0}, 1, 1}); // positron + } + + // ---------------------------------------------------------------------- + // run the simulation + + auto accessor = mprts.accessor(); + auto prts = accessor[p]; + + ASSERT_EQ(prts.size(), 2); + ASSERT_LT(prts[0].u()[1], 0.0); + ASSERT_GT(prts[1].u()[1], 0.0); + + ASSERT_EQ(prts[0].m(), 1.0); + ASSERT_EQ(prts[1].m(), 1.0); + + ASSERT_EQ(prts[0].q(), -1.0); + ASSERT_EQ(prts[1].q(), 1.0); + + for (; grid.timestep_ < psc_params.nmax;) { + psc.step(); + ASSERT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); + } + + ASSERT_TRUE(prts.size() < 2); +} + +// ====================================================================== +// main + +int main(int argc, char** argv) +{ + psc_init(argc, argv); + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + psc_finalize(); + return rc; +} diff --git a/src/libpsc/tests/testing.hxx b/src/libpsc/tests/testing.hxx index 6dc2584bc1..1922b34cd7 100644 --- a/src/libpsc/tests/testing.hxx +++ b/src/libpsc/tests/testing.hxx @@ -14,10 +14,6 @@ #include "setup_fields.hxx" #include "setup_particles.hxx" -#include "../libpsc/vpic/mfields_state_psc.hxx" -#include "../libpsc/vpic/push_particles_vpic.hxx" -#include "../libpsc/vpic/bnd_particles_vpic.hxx" - #ifdef USE_CUDA #include "../libpsc/cuda/push_particles_cuda_impl.hxx" #include "../libpsc/cuda/setup_fields_cuda.hxx" @@ -91,22 +87,6 @@ using TestConfig1vbec3dSingleXZ = Config1vbecSplit>, checks_order_1st>; -using VpicConfig = VpicConfigPsc; - -using _MfieldsStateVpic = typename VpicConfig::MfieldsState; -using _MparticlesVpic = typename VpicConfig::Mparticles; -using _PushParticlesVpic = PushParticlesVpic< - _MparticlesVpic, _MfieldsStateVpic, typename VpicConfig::ParticlesOps, - typename VpicConfig::AccumulatorOps, typename VpicConfig::AccumulateOps, - typename VpicConfig::InterpolatorOps>; - -using TestConfigVpic = - TestConfig, - BndParticlesVpic<_MparticlesVpic>>; - #ifdef USE_CUDA using TestConfig1vbec3dCuda = TestConfig -struct BndParticlesVpic : BndParticlesBase -{ - BndParticlesVpic(const Grid_t& grid) {} - - void operator()(Mparticles& mprts) {} -}; diff --git a/src/psc_config.hxx b/src/psc_config.hxx index 6401a481e0..9f1e40ac2e 100644 --- a/src/psc_config.hxx +++ b/src/psc_config.hxx @@ -181,88 +181,3 @@ struct PscConfig1vbecCuda {}; #endif - -#include "../libpsc/vpic/bnd_fields_vpic.hxx" -#include "../libpsc/vpic/bnd_particles_vpic.hxx" -#include "../libpsc/vpic/bnd_vpic.hxx" -#include "../libpsc/vpic/checks_vpic.hxx" -#include "../libpsc/vpic/collision_vpic.hxx" -#include "../libpsc/vpic/marder_vpic.hxx" -#include "../libpsc/vpic/push_fields_vpic.hxx" -#include "../libpsc/vpic/push_particles_vpic.hxx" -#include "../libpsc/vpic/sort_vpic.hxx" - -#ifdef USE_VPIC -struct PscConfigVpicWrap -{ - using VpicConfig = VpicConfigWrap; - - using MfieldsState = typename VpicConfig::MfieldsState; - using Mparticles = typename VpicConfig::Mparticles; - using MfieldsHydro = typename VpicConfig::MfieldsHydro; - - using Balance = Balance_; - using Sort = SortVpicWrap; - using Collision = PscCollisionVpic; - using PushParticles = PushParticlesVpic< - Mparticles, MfieldsState, typename VpicConfig::ParticlesOps, - typename VpicConfig::AccumulatorOps, typename VpicConfig::AccumulateOps, - typename VpicConfig::InterpolatorOps>; - using PushFields = PushFieldsVpicWrap; - using Bnd = BndVpic; - using BndFields = BndFieldsVpic; - using BndParticles = BndParticlesVpic; - using Checks = ChecksVpic; - using Marder = MarderVpicWrap; - using OutputParticles = OutputParticlesHdf5; - using OutputHydro = - OutputHydroVpicWrap; - using Dim = dim_xyz; - -#if 0 - using DiagMixin = VpicDiagMixin; -#else - using DiagMixin = NoneDiagMixin; -#endif -}; -#endif - -struct PscConfigVpicPsc -{ - using VpicConfig = VpicConfigPsc; - - using MfieldsState = typename VpicConfig::MfieldsState; - using Mparticles = typename VpicConfig::Mparticles; - using MfieldsHydro = MfieldsHydroQ; - - using Balance = Balance_; - using Sort = SortVpic; - using Collision = PscCollisionVpic; - using PushParticles = PushParticlesVpic< - Mparticles, MfieldsState, typename VpicConfig::ParticlesOps, - typename VpicConfig::AccumulatorOps, typename VpicConfig::AccumulateOps, - typename VpicConfig::InterpolatorOps>; - using PushFields = PushFieldsVpic; - using Bnd = BndVpic; - using BndFields = BndFieldsVpic; - using BndParticles = BndParticlesVpic; - using Checks = ChecksVpic; - using Marder = MarderVpic; - using OutputParticles = OutputParticlesHdf5; - using OutputHydro = OutputHydroQ; - using Dim = dim_xz; - -#if 0 - using DiagMixin = VpicDiagMixin; -#else - using DiagMixin = - NoneDiagMixin; -#endif -};