From 74dd2914ca7b866f98d22512b0d9c1c52eaa1486 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 6 Nov 2025 09:39:56 -0500 Subject: [PATCH 01/36] +test_open_bcs_integration it crashes when the particles are dropped, it seems --- src/libpsc/tests/CMakeLists.txt | 1 + .../tests/test_open_bcs_integration.cxx | 170 ++++++++++++++++++ 2 files changed, 171 insertions(+) create mode 100644 src/libpsc/tests/test_open_bcs_integration.cxx diff --git a/src/libpsc/tests/CMakeLists.txt b/src/libpsc/tests/CMakeLists.txt index e46e74a4b..6ed4f2882 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_open_bcs_integration.cxx b/src/libpsc/tests/test_open_bcs_integration.cxx new file mode 100644 index 000000000..92e36fcde --- /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; +} From 247d5e1bbf66f040ca4b21e2c4cdfc1fd84ba055 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 6 Nov 2025 09:50:51 -0500 Subject: [PATCH 02/36] bnd_particles_impl: drop for open BCs --- src/include/bnd_particles_impl.hxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/include/bnd_particles_impl.hxx b/src/include/bnd_particles_impl.hxx index 653b3be54..cbe30b2ec 100644 --- a/src/include/bnd_particles_impl.hxx +++ b/src/include/bnd_particles_impl.hxx @@ -153,6 +153,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, pxi[d] = -pxi[d]; dir[d] = 0; break; + case BND_PRT_OPEN: case BND_PRT_ABSORBING: drop = true; break; default: assert(0); } @@ -177,6 +178,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, } break; } + case BND_PRT_OPEN: case BND_PRT_ABSORBING: drop = true; break; default: assert(0); } From a23b3edd7ffcdbdab1bc9be76efb7129bc714b14 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 12:26:31 -0500 Subject: [PATCH 03/36] bnd_particles_impl: no fallthrough --- src/include/bnd_particles_impl.hxx | 23 ++++++++++++++++++----- 1 file changed, 18 insertions(+), 5 deletions(-) diff --git a/src/include/bnd_particles_impl.hxx b/src/include/bnd_particles_impl.hxx index cbe30b2ec..4173e0e8a 100644 --- a/src/include/bnd_particles_impl.hxx +++ b/src/include/bnd_particles_impl.hxx @@ -148,13 +148,20 @@ 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_OPEN: - case BND_PRT_ABSORBING: drop = true; break; + } + case BND_PRT_OPEN: { + drop = true; + break; + } + case BND_PRT_ABSORBING: { + drop = true; + break; + } default: assert(0); } } @@ -178,8 +185,14 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, } break; } - case BND_PRT_OPEN: - case BND_PRT_ABSORBING: drop = true; break; + case BND_PRT_OPEN: { + drop = true; + break; + } + case BND_PRT_ABSORBING: { + drop = true; + break; + } default: assert(0); } } From ffb670ad18a65fdbdad90d195f194ac481cc4498 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 29 Dec 2025 09:24:20 -0500 Subject: [PATCH 04/36] psc_config: -vpic stuff --- src/psc_config.hxx | 85 ---------------------------------------------- 1 file changed, 85 deletions(-) diff --git a/src/psc_config.hxx b/src/psc_config.hxx index 6401a481e..9f1e40ac2 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 -}; From 9a12ddb5fd34e6010d3bc598bcecd752ba4fc340 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 13:29:03 -0500 Subject: [PATCH 05/36] -bnd_particles_vpic; * --- src/libpsc/tests/testing.hxx | 20 -------------------- src/libpsc/vpic/bnd_particles_vpic.hxx | 12 ------------ 2 files changed, 32 deletions(-) delete mode 100644 src/libpsc/vpic/bnd_particles_vpic.hxx diff --git a/src/libpsc/tests/testing.hxx b/src/libpsc/tests/testing.hxx index 6dc2584bc..1922b34cd 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) {} -}; From ff84046e968066596ba283d9909e02d77bc3c1ce Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 13:29:44 -0500 Subject: [PATCH 06/36] -bnd_particles; *: -BndParticlesBase --- src/include/bnd_particles.hxx | 10 ---------- src/include/bnd_particles_impl.hxx | 4 ++-- src/libpsc/integrate.cxx | 1 - src/libpsc/psc_balance/psc_balance_impl.hxx | 1 - 4 files changed, 2 insertions(+), 14 deletions(-) delete mode 100644 src/include/bnd_particles.hxx diff --git a/src/include/bnd_particles.hxx b/src/include/bnd_particles.hxx deleted file mode 100644 index a578670e2..000000000 --- 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 4173e0e8a..93dc3d09c 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,7 +17,7 @@ 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; diff --git a/src/libpsc/integrate.cxx b/src/libpsc/integrate.cxx index 6f67ba3b9..92430a1b3 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 78164a445..2799464ed 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" From 7f73dbe03601e1e20bddc1503084f1744cd7946f Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 14:36:56 -0500 Subject: [PATCH 07/36] bnd_particles_impl: xm -> patch_size --- src/include/bnd_particles_impl.hxx | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/include/bnd_particles_impl.hxx b/src/include/bnd_particles_impl.hxx index 93dc3d09c..668534706 100644 --- a/src/include/bnd_particles_impl.hxx +++ b/src/include/bnd_particles_impl.hxx @@ -21,6 +21,7 @@ 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; @@ -98,10 +99,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, 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++) { @@ -139,7 +137,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]) { @@ -167,7 +165,7 @@ void BndParticlesCommon::process_patch(const Grid_t& grid, } } 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) { @@ -176,7 +174,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); @@ -206,7 +204,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) { From 0af790c0cf4d7ff6594ff76a82fb10fdd8ffbf6d Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 14:39:24 -0500 Subject: [PATCH 08/36] bnd_particles_impl: use grid.ldims --- src/include/bnd_particles_impl.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/bnd_particles_impl.hxx b/src/include/bnd_particles_impl.hxx index 668534706..81abd3312 100644 --- a/src/include/bnd_particles_impl.hxx +++ b/src/include/bnd_particles_impl.hxx @@ -97,8 +97,8 @@ 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(); Real3 patch_size = Real3(gpatch.xe - gpatch.xb); auto* dpatch = &ddcp->patches_[p]; From 9c6db4c15ed80852e49eeef00915a8e77a0537b4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 14:43:28 -0500 Subject: [PATCH 09/36] domain: +dx_inv --- src/include/grid/domain.hxx | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/include/grid/domain.hxx b/src/include/grid/domain.hxx index ce7f20f6a..2e4661c63 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 From 5a47c6a6fa78ef2d4a5c776554fb87987ca9ede1 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 14:56:03 -0500 Subject: [PATCH 10/36] domain; *: +dx_inv --- src/include/boundary_injector.hxx | 2 +- src/include/particle_indexer.hxx | 2 +- src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx | 2 +- src/libpsc/psc_push_particles/inc_curr_1vb_var1.cxx | 3 +-- src/libpsc/psc_push_particles/inc_curr_zigzag.cxx | 2 +- src/libpsc/psc_push_particles/push_particles_1vb.hxx | 2 +- src/libpsc/psc_push_particles/push_particles_esirkepov.hxx | 2 +- src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx | 2 +- 8 files changed, 8 insertions(+), 9 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index c7b0258e6..8b6cfd402 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -97,7 +97,7 @@ 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++) { diff --git a/src/include/particle_indexer.hxx b/src/include/particle_indexer.hxx index e8e7b8453..1263fc225 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]; } 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 c6af98748..00ef3cc58 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 78d5034cc..0c13effb7 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 51e07eb0e..f6cd6b433 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 d2c080300..cc48976e8 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -89,7 +89,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 58285a867..085668ab1 100644 --- a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx +++ b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx @@ -36,7 +36,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 = 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/pushp_current_esirkepov.hxx b/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx index 8e611ebc0..76f01174a 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; From e7ebf89820bdf81b30847f2fa857b9472d326c15 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 16:07:53 -0500 Subject: [PATCH 11/36] particle_indexer: +isValidCellPosition --- src/include/particle_indexer.hxx | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/src/include/particle_indexer.hxx b/src/include/particle_indexer.hxx index 1263fc225..4efd41eaa 100644 --- a/src/include/particle_indexer.hxx +++ b/src/include/particle_indexer.hxx @@ -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", From 1192b0b324ab2e033942bed159129e258ac1144e Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 20 Jan 2026 16:09:18 -0500 Subject: [PATCH 12/36] bnd_particles_impl: use isValidCellPosition --- src/include/bnd_particles_impl.hxx | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/include/bnd_particles_impl.hxx b/src/include/bnd_particles_impl.hxx index 81abd3312..abf18020a 100644 --- a/src/include/bnd_particles_impl.hxx +++ b/src/include/bnd_particles_impl.hxx @@ -119,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; } From 3776316e994a8cff0fbf016e6c8e26f45febaa04 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 21 Jan 2026 09:23:57 -0500 Subject: [PATCH 13/36] push_particles_1vb: one more Int3 --- src/libpsc/psc_push_particles/push_particles_1vb.hxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index cc48976e8..4f5c005b6 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -67,18 +67,18 @@ struct PushParticlesVb Int3 final_index = final_pos_normalized.fint(); // 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); } } } From 7675046947d99908ca09425c6d7c6b8a782a5176 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 21 Jan 2026 15:04:11 -0500 Subject: [PATCH 14/36] push_particles_1vb: +exit_to_edge --- .../psc_push_particles/push_particles_1vb.hxx | 41 +++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index 4f5c005b6..a63177f87 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -2,6 +2,45 @@ // 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_& grid, int p) +{ + 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,6 +105,8 @@ 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 Int3 initial_index; if (!Dim::InvarX::value) { From 0c89b61bdac53be52e62e7406cdf42ce4b20ed57 Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 21 Jan 2026 15:20:27 -0500 Subject: [PATCH 15/36] push_particles_esirkepov: +handle_exiting (not really) --- .../push_particles_esirkepov.hxx | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx index 085668ab1..34605d747 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 @@ -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); } } } From 5d63c8f0aed2dcd2ed3176de00e8963a7bb67b8d Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 21 Jan 2026 16:49:11 -0500 Subject: [PATCH 16/36] boundary_injector: rename --- src/include/boundary_injector.hxx | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 8b6cfd402..4a443b801 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -113,16 +113,17 @@ public: 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]++) { + for (Int3 initial_idx = ilo; initial_idx[0] < ihi[0]; initial_idx[0]++) { + for (initial_idx[2] = ilo[2]; initial_idx[2] < ihi[2]; + initial_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); + flds.storage()(initial_idx[0] + grid.ibn[0], -1 + grid.ibn[1], + initial_idx[2] + grid.ibn[2], JXI + 1); - cell_idx[INJECT_DIM_IDX_] = -1; + initial_idx[INJECT_DIM_IDX_] = -1; Real3 cell_corner = - grid.domain.corner + Double3(cell_idx) * grid.domain.dx; + grid.domain.corner + Double3(initial_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; @@ -156,7 +157,7 @@ public: // 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); + final_idx, initial_idx, qni_wni, v); charge_injected += qni_wni; } @@ -164,8 +165,8 @@ public: // 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) = + flds.storage()(initial_idx[0] + grid.ibn[0], -1 + grid.ibn[1], + initial_idx[2] + grid.ibn[2], JXI + 1) = boundary_current_before + charge_injected * grid.domain.dx[1] / grid.dt / prts_per_unit_density; } From 83a5d15b6a3fed4b694aa1061ea7c5a58cec80a0 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 08:08:17 -0500 Subject: [PATCH 17/36] test_boundary_injector: integrate 2 steps (fails) --- src/libpsc/tests/test_boundary_injector.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 44e540fc1..4b899fef3 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -112,7 +112,7 @@ TEST(BoundaryInjectorTest, Integration1Particle) PscParams psc_params; - psc_params.nmax = 1; + psc_params.nmax = 2; psc_params.stats_every = 1; psc_params.cfl = .75; @@ -175,7 +175,7 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) PscParams psc_params; - psc_params.nmax = 1; + psc_params.nmax = 2; psc_params.stats_every = 1; psc_params.cfl = .75; @@ -238,7 +238,7 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) PscParams psc_params; - psc_params.nmax = 1; + psc_params.nmax = 2; psc_params.stats_every = 1; psc_params.cfl = .75; From 33f2539177c717f4c417c4d227c3168eea93c3be Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:02:55 -0500 Subject: [PATCH 18/36] boundary_injector: simplify deposition tests show the same magnitude of failure as before --- src/include/boundary_injector.hxx | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 4a443b801..4924335be 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -116,17 +116,12 @@ public: for (Int3 initial_idx = ilo; initial_idx[0] < ihi[0]; initial_idx[0]++) { for (initial_idx[2] = ilo[2]; initial_idx[2] < ihi[2]; initial_idx[2]++) { - auto boundary_current_before = - flds.storage()(initial_idx[0] + grid.ibn[0], -1 + grid.ibn[1], - initial_idx[2] + grid.ibn[2], JXI + 1); - initial_idx[INJECT_DIM_IDX_] = -1; Real3 cell_corner = grid.domain.corner + Double3(initial_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++) { @@ -148,6 +143,8 @@ public: // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts() 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; Int3 final_idx; Real3 final_normalized_pos; @@ -158,17 +155,7 @@ public: 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); - - charge_injected += qni_wni; } - - // override whatever current was deposited in the boundary to be - // consistent with the charge having started completely out of the - // domain - flds.storage()(initial_idx[0] + grid.ibn[0], -1 + grid.ibn[1], - initial_idx[2] + grid.ibn[2], JXI + 1) = - boundary_current_before + charge_injected * grid.domain.dx[1] / - grid.dt / prts_per_unit_density; } } } From 4a22298801bcda515cae3934fb7d4f0849661489 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:19:01 -0500 Subject: [PATCH 19/36] scripts/run_test: mpirun --- scripts/run_test.sh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/scripts/run_test.sh b/scripts/run_test.sh index c6a5dbf3e..4d464d30f 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}" From 6a6be6911424aebc3a745352b786346a3f86f5f2 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:22:22 -0500 Subject: [PATCH 20/36] boundary_injector: Real3=Vec3 --- src/include/boundary_injector.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 4924335be..bc50c5517 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -74,7 +74,7 @@ public: using Current = typename PushParticles::Current; using Dim = typename PushParticles::Dim; using real_t = typename PushParticles::real_t; - using Real3 = typename PushParticles::Real3; + using Real3 = Vec3; using checks_order = typename PushParticles::checks_order; BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) From 53055b1deaf8ff3c5ddf4dbf023c5ada12e920bd Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:26:56 -0500 Subject: [PATCH 21/36] boundary_injector: cleanup types --- src/include/boundary_injector.hxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index bc50c5517..08a015d89 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -72,10 +72,8 @@ 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 = Vec3; - using checks_order = typename PushParticles::checks_order; BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) : partice_generator{particle_generator}, @@ -163,6 +161,6 @@ public: private: ParticleGenerator partice_generator; - AdvanceParticle advance; + AdvanceParticle_t advance; real_t prts_per_unit_density; }; From 55c9cd45931fb35680a10ed9761b9d557b532bd4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:44:05 -0500 Subject: [PATCH 22/36] boundary_injector: -find_idx_pos_1st_rel --- src/include/boundary_injector.hxx | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 08a015d89..42071a45b 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -144,10 +144,8 @@ public: // pretend it came from the edge to inject proper current initial_normalized_pos[INJECT_DIM_IDX_] = -1.0; - Int3 final_idx; - Real3 final_normalized_pos; - find_idx_pos_1st_rel(prt.x, dxi, final_idx, final_normalized_pos, - real_t(0.)); + Real3 final_normalized_pos = prt.x * dxi; + Int3 final_idx = final_normalized_pos.fint(); // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt real_t qni_wni = grid.kinds[prt.kind].q * prt.w; From 2db0ab4c7a1b585e29964fb647a2e51ae516098f Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:48:40 -0500 Subject: [PATCH 23/36] boundary_injector: use VecRange and maybe fix bounds --- src/include/boundary_injector.hxx | 71 +++++++++++++++---------------- 1 file changed, 34 insertions(+), 37 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 42071a45b..6ec1aefd1 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 @@ -99,59 +100,55 @@ public: 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 initial_idx = ilo; initial_idx[0] < ihi[0]; initial_idx[0]++) { - for (initial_idx[2] = ilo[2]; initial_idx[2] < ihi[2]; - initial_idx[2]++) { - initial_idx[INJECT_DIM_IDX_] = -1; - - Real3 cell_corner = - grid.domain.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 (Int3 initial_idx : VecRange(ilo, ihi)) { + Real3 cell_corner = + grid.domain.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++) { - psc::particle::Inject prt = - partice_generator.get(cell_corner, grid.domain.dx); + 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); + 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; - } + 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); + injector(prt); - // Update currents - // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts() + // Update currents + // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts() - 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; + 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; - Real3 final_normalized_pos = prt.x * dxi; - Int3 final_idx = final_normalized_pos.fint(); + Real3 final_normalized_pos = prt.x * dxi; + Int3 final_idx = final_normalized_pos.fint(); - // 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); - } + // 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); } } } From 6ceead6b66e3f0b1904078cffd134fe54e2af4e8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:52:07 -0500 Subject: [PATCH 24/36] boundary_injector: cell_corner is patch-frame --- src/include/boundary_injector.hxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 6ec1aefd1..974ca85aa 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -115,8 +115,7 @@ public: typename Current::fields_t J(flds); for (Int3 initial_idx : VecRange(ilo, ihi)) { - Real3 cell_corner = - grid.domain.corner + Double3(initial_idx) * grid.domain.dx; + 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); @@ -128,8 +127,8 @@ public: 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 + if (prt.x[INJECT_DIM_IDX_] < 0.0) { + // don't inject a particle that fails to enter the patch continue; } From 3fed928794baf77d8770c7aa268d4cfbc746a279 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 09:54:54 -0500 Subject: [PATCH 25/36] boundary_injector: typo --- src/include/boundary_injector.hxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 974ca85aa..4f732fb96 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -77,7 +77,7 @@ public: using Real3 = Vec3; BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) - : partice_generator{particle_generator}, + : particle_generator{particle_generator}, advance{grid.dt}, prts_per_unit_density{grid.norm.prts_per_unit_density} {} @@ -121,7 +121,7 @@ public: 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); + particle_generator.get(cell_corner, grid.domain.dx); Real3 v = advance.calc_v(prt.u); Real3 initial_x = prt.x; @@ -154,7 +154,7 @@ public: } private: - ParticleGenerator partice_generator; + ParticleGenerator particle_generator; AdvanceParticle_t advance; real_t prts_per_unit_density; }; From 504bbd478589e0620f495ad03692827fd02a898e Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 12:57:59 -0500 Subject: [PATCH 26/36] test_boundary_injector: unit values --- src/libpsc/tests/test_boundary_injector.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 4b899fef3..e26cddac1 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; From d6599af4684cb1111e663a8824c3d05dcd4f17c1 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 12:59:34 -0500 Subject: [PATCH 27/36] psc_bnd_fields_impl: revive open_H --- .../psc_bnd_fields/psc_bnd_fields_impl.hxx | 69 ++++++++++--------- 1 file changed, 36 insertions(+), 33 deletions(-) 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 1c670e5aa..3a4751687 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) + @@ -644,7 +647,7 @@ struct BndFieldsNone : BndFieldsBase { using Mfields = MF; - void fill_ghosts_E(Mfields& mflds){}; - void fill_ghosts_H(Mfields& mflds){}; - void add_ghosts_J(Mfields& mflds){}; + void fill_ghosts_E(Mfields& mflds) {}; + void fill_ghosts_H(Mfields& mflds) {}; + void add_ghosts_J(Mfields& mflds) {}; }; From d9affd0a2f35d4f0ef6be4d00a303449ec78cec1 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 13:30:51 -0500 Subject: [PATCH 28/36] boundary_injector: don't specify dt_frac --- src/include/boundary_injector.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 4f732fb96..5d22b43cb 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -125,7 +125,7 @@ public: Real3 v = advance.calc_v(prt.u); Real3 initial_x = prt.x; - advance.push_x(prt.x, v, 1.0); + 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 From 9f60fa708573f42c91da523465b631462d6c2515 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 13:34:31 -0500 Subject: [PATCH 29/36] test_boundary_injector: assert checks --- src/libpsc/tests/test_boundary_injector.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index e26cddac1..5af54f888 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -161,8 +161,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); @@ -224,8 +224,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); @@ -292,8 +292,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; From 44230a6d975a87ea68870c81628f50bbee9b0ff5 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 14:13:09 -0500 Subject: [PATCH 30/36] *: explain fake gauss error from ghost corners --- src/include/boundary_injector.hxx | 6 ++++++ src/libpsc/psc_push_particles/push_particles_1vb.hxx | 5 +++++ src/libpsc/tests/test_boundary_injector.cxx | 6 ++++++ 3 files changed, 17 insertions(+) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 5d22b43cb..93c3885da 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -120,6 +120,12 @@ public: 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); diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index a63177f87..95690fa13 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -25,6 +25,11 @@ template void exit_to_edge(Vec3& final_x_normed, const Int3& final_i3, const Grid_& 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) { diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index 5af54f888..ef51fd70f 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -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; From 646370d93de53e82e6443eddd42eca7054dfb3ce Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 16:04:45 -0500 Subject: [PATCH 31/36] push_particles_1vb: fix grid t --- src/libpsc/psc_push_particles/push_particles_1vb.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index 95690fa13..f1e96b211 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -23,7 +23,7 @@ namespace */ template void exit_to_edge(Vec3& final_x_normed, const Int3& final_i3, - const Grid_& grid, int p) + 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 From 880c71927f7415bfcca43e0a8c16570d4f6c4912 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 16:12:04 -0500 Subject: [PATCH 32/36] push_particles_esirkepov: fix dxi --- src/libpsc/psc_push_particles/push_particles_esirkepov.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx index 34605d747..4b292f082 100644 --- a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx +++ b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx @@ -60,7 +60,7 @@ struct PushParticlesEsirkepov static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi = grid.domain.dx_inv; + Real3 dxi(grid.domain.dx_inv); real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; assert(kinds.size() <= MAX_NR_KINDS); From d47f8972a6134dc36c980bc35731da3da5dbcde9 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 22 Jan 2026 16:23:12 -0500 Subject: [PATCH 33/36] psc_bnd_fields_impl: i win --- src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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 3a4751687..2e8c18516 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -640,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; + // clang-format off void fill_ghosts_E(Mfields& mflds) {}; void fill_ghosts_H(Mfields& mflds) {}; void add_ghosts_J(Mfields& mflds) {}; + // clang-format on }; From d6f85a477d21fca67f749a2324e90fad683c9fd2 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 30 Jan 2026 13:24:46 -0500 Subject: [PATCH 34/36] push_particles_esirkepov: un-"fix" dxi This reverts commit 880c71927f7415bfcca43e0a8c16570d4f6c4912. --- src/libpsc/psc_push_particles/push_particles_esirkepov.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx index 4b292f082..34605d747 100644 --- a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx +++ b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx @@ -60,7 +60,7 @@ struct PushParticlesEsirkepov static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi(grid.domain.dx_inv); + Real3 dxi = grid.domain.dx_inv; real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; assert(kinds.size() <= MAX_NR_KINDS); From 6f4fda89d12f458fa3d2bf6190dabf6bca01935d Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 30 Jan 2026 13:26:27 -0500 Subject: [PATCH 35/36] boundary_injector: append _ to private fields --- src/include/boundary_injector.hxx | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx index 93c3885da..72f9dc4eb 100644 --- a/src/include/boundary_injector.hxx +++ b/src/include/boundary_injector.hxx @@ -77,9 +77,9 @@ public: using Real3 = Vec3; BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) - : particle_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 @@ -117,7 +117,7 @@ public: 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); + 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) @@ -127,11 +127,11 @@ public: // 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); + particle_generator_.get(cell_corner, grid.domain.dx); - Real3 v = advance.calc_v(prt.u); + Real3 v = advance_.calc_v(prt.u); Real3 initial_x = prt.x; - advance.push_x(prt.x, v); + 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 @@ -160,7 +160,7 @@ public: } private: - ParticleGenerator particle_generator; - AdvanceParticle_t advance; - real_t prts_per_unit_density; + ParticleGenerator particle_generator_; + AdvanceParticle_t advance_; + real_t prts_per_unit_density_; }; From 83865e982be9e0100740642fbd826c82e5dda8cc Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 30 Jan 2026 13:33:25 -0500 Subject: [PATCH 36/36] push_particles_esirkepov: un-unfix dxi --- src/libpsc/psc_push_particles/push_particles_esirkepov.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx index 34605d747..22617b280 100644 --- a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx +++ b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx @@ -60,7 +60,7 @@ struct PushParticlesEsirkepov static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { const auto& grid = mprts.grid(); - Real3 dxi = grid.domain.dx_inv; + 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);