Skip to content

Fix field reflection across lower z #364

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
JamesMcClung merged 15 commits into from
Oct 24, 2025
Merged

Fix field reflection across lower z #364

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
ede9b01
test_reflective_bcs_integration: +IntegrationZ (fails)
JamesMcClung Oct 20, 2025
a59b46d
test_reflective_bcs_integration: specify Y
JamesMcClung Oct 20, 2025
089612c
test_reflective_bcs_integration: mv start pos
JamesMcClung Oct 20, 2025
736c7f8
test_reflective_bcs_integration: stricter detect
JamesMcClung Oct 20, 2025
a304166
test_reflective_bcs_integration: shorter nmax
JamesMcClung Oct 20, 2025
8a6e04a
test_reflective_bcs_integration: asserts
JamesMcClung Oct 20, 2025
226b556
test_reflective_bcs_integration: drop old assert
JamesMcClung Oct 20, 2025
823079e
test_reflective_bcs_integration: fix detect
JamesMcClung Oct 20, 2025
20be7b9
psc_bnd_fields_impl: change xbnds
JamesMcClung Oct 20, 2025
9640982
psc_bnd_fields_impl: fix conducting j
JamesMcClung Oct 20, 2025
671c967
psc_bnd_fields_impl: change a bound
JamesMcClung Oct 20, 2025
e0f1bb7
psc_bnd_fields_impl: combine/reorder loops
JamesMcClung Oct 20, 2025
97133ec
psc_bnd_fields_impl: consistently do 1 ghost cell
JamesMcClung Oct 20, 2025
f5765b8
psc_bnd_fields_impl: rm if debugs
JamesMcClung Oct 20, 2025
2027ac3
psc_bnd_fields_impl: restore if debugs
JamesMcClung Oct 24, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
114 changes: 43 additions & 71 deletions src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,18 +162,11 @@ struct BndFields_ : BndFieldsBase
ix++) {
F(EX, ix, 0, iz) = 0.;
F(EX, ix, -1, iz) = F(EX, ix, 1, iz);
F(EX, ix, -2, iz) = F(EX, ix, 2, iz);

F(EY, ix, -1, iz) = -F(EY, ix, 0, iz);
F(EY, ix, -2, iz) = -F(EY, ix, 1, iz);
}
}

for (int iz = -2; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(EZ, ix, 0, iz) = 0.;
F(EZ, ix, -1, iz) = F(EZ, ix, 1, iz);
F(EZ, ix, -2, iz) = F(EZ, ix, 2, iz);
}
}
} else if (d == 2) {
Expand All @@ -190,18 +183,15 @@ struct BndFields_ : BndFieldsBase
}
#endif
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(EX, ix, iy, 0) = 0.;
F(EX, ix, iy, -1) = F(EX, ix, iy, 1);
F(EZ, ix, iy, -1) = -F(EZ, ix, iy, 0);
F(EZ, ix, iy, -2) = -F(EZ, ix, iy, 1);
}
}

for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {
F(EY, ix, iy, 0) = 0.;
F(EY, ix, iy, -1) = F(EY, ix, iy, 1);

F(EZ, ix, iy, -1) = -F(EZ, ix, iy, 0);
}
}
} else {
Expand Down Expand Up @@ -235,13 +225,9 @@ struct BndFields_ : BndFieldsBase
ix++) {
F(EX, ix, my, iz) = 0.;
F(EX, ix, my + 1, iz) = F(EX, ix, my - 1, iz);

F(EY, ix, my, iz) = -F(EY, ix, my - 1, iz);
}
}

for (int iz = -2; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(EZ, ix, my, iz) = 0.;
F(EZ, ix, my + 1, iz) = F(EZ, ix, my - 1, iz);
}
Expand All @@ -261,18 +247,15 @@ struct BndFields_ : BndFieldsBase
}
#endif
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(EX, ix, iy, mz) = 0.;
F(EX, ix, iy, mz + 1) = F(EX, ix, iy, mz - 1);
F(EZ, ix, iy, mz) = -F(EZ, ix, iy, mz - 1);
F(EZ, ix, iy, mz + 1) = -F(EZ, ix, iy, mz - 2);
}
}

for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {
F(EY, ix, iy, mz) = 0.;
F(EY, ix, iy, mz + 1) = F(EY, ix, iy, mz - 1);

F(EZ, ix, iy, mz) = -F(EZ, ix, iy, mz - 1);
}
}
} else {
Expand Down Expand Up @@ -300,16 +283,13 @@ struct BndFields_ : BndFieldsBase
}
}
#endif
for (int iz = -1; iz < ldims[2] + 1; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(HY, ix, -1, iz) = F(HY, ix, 1, iz);
F(HX, ix, -1, iz) = -F(HX, ix, 0, iz);
}
}
for (int iz = -1; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(HX, ix, -1, iz) = -F(HX, ix, 0, iz);

F(HY, ix, -1, iz) = F(HY, ix, 1, iz);

F(HZ, ix, -1, iz) = -F(HZ, ix, 0, iz);
}
}
Expand All @@ -327,16 +307,13 @@ struct BndFields_ : BndFieldsBase
}
#endif
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {
F(HZ, ix, iy, -1) = F(HZ, ix, iy, 1);
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(HX, ix, iy, -1) = -F(HX, ix, iy, 0);
F(HX, ix, iy, -2) = -F(HX, ix, iy, 1);
}
}
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {

F(HY, ix, iy, -1) = -F(HY, ix, iy, 0);
F(HY, ix, iy, -2) = -F(HY, ix, iy, 1);

F(HZ, ix, iy, -1) = F(HZ, ix, iy, 1);
}
}
} else {
Expand Down Expand Up @@ -368,16 +345,11 @@ struct BndFields_ : BndFieldsBase
for (int iz = -2; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(HY, ix, my + 1, iz) = F(HY, ix, my - 1, iz);
F(HX, ix, my, iz) = -F(HX, ix, my - 1, iz);
F(HX, ix, my + 1, iz) = -F(HX, ix, my - 2, iz);
}
}
for (int iz = -2; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {

F(HY, ix, my + 1, iz) = F(HY, ix, my - 1, iz);

F(HZ, ix, my, iz) = -F(HZ, ix, my - 1, iz);
F(HZ, ix, my + 1, iz) = -F(HZ, ix, my - 2, iz);
}
}
} else if (d == 2) {
Expand All @@ -394,16 +366,13 @@ struct BndFields_ : BndFieldsBase
}
#endif
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {
F(HZ, ix, iy, mz + 1) = F(HZ, ix, iy, mz - 1);
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(HX, ix, iy, mz) = -F(HX, ix, iy, mz - 1);
F(HX, ix, iy, mz + 1) = -F(HX, ix, iy, mz - 2);
}
}
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = -2; ix < ldims[0] + 2; ix++) {

F(HY, ix, iy, mz) = -F(HY, ix, iy, mz - 1);
F(HY, ix, iy, mz + 1) = -F(HY, ix, iy, mz - 2);

F(HZ, ix, iy, mz + 1) = F(HZ, ix, iy, mz - 1);
}
}
} else {
Expand All @@ -421,12 +390,12 @@ struct BndFields_ : BndFieldsBase
for (int iz = -2; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(JYI, ix, 1, iz) -= F(JYI, ix, -2, iz);
F(JYI, ix, 0, iz) -= F(JYI, ix, -1, iz);
F(JYI, ix, -1, iz) = 0.;
F(JYI, ix, -2, iz) = 0.;
F(JXI, ix, 1, iz) += F(JXI, ix, -1, iz);
F(JXI, ix, -1, iz) = 0.;

F(JYI, ix, 0, iz) -= F(JYI, ix, -1, iz);
F(JYI, ix, -1, iz) = 0.;

F(JZI, ix, 1, iz) += F(JZI, ix, -1, iz);
F(JZI, ix, -1, iz) = 0.;
}
Expand All @@ -435,13 +404,14 @@ struct BndFields_ : BndFieldsBase
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(JZI, ix, iy, 0) -= F(JZI, ix, iy, -1);
F(JZI, ix, iy, 0) -= F(JZI, ix, iy, -1);
F(JZI, ix, iy, -1) = 0.;
F(JXI, ix, iy, 1) += F(JXI, ix, iy, -1);
F(JXI, ix, iy, -1) = 0.;

F(JYI, ix, iy, 1) += F(JYI, ix, iy, -1);
F(JYI, ix, iy, -1) = 0.;

F(JZI, ix, iy, 0) -= F(JZI, ix, iy, -1);
F(JZI, ix, iy, -1) = 0.;
}
}
} else {
Expand All @@ -460,12 +430,12 @@ struct BndFields_ : BndFieldsBase
for (int iz = -2; iz < ldims[2] + 2; iz++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(JYI, ix, my - 2, iz) -= F(JYI, ix, my + 1, iz);
F(JYI, ix, my - 1, iz) -= F(JYI, ix, my, iz);
F(JYI, ix, my, iz) = 0.;
F(JYI, ix, my + 1, iz) = 0.;
F(JXI, ix, my - 1, iz) += F(JXI, ix, my + 1, iz);
F(JXI, ix, my + 1, iz) = 0.;

F(JYI, ix, my - 1, iz) -= F(JYI, ix, my, iz);
F(JYI, ix, my, iz) = 0.;

F(JZI, ix, my - 1, iz) += F(JZI, ix, my + 1, iz);
F(JZI, ix, my + 1, iz) = 0.;
}
Expand All @@ -475,12 +445,14 @@ struct BndFields_ : BndFieldsBase
for (int iy = -2; iy < ldims[1] + 2; iy++) {
for (int ix = MAX(-2, ib[0]); ix < MIN(ldims[0] + 2, ib[0] + im[0]);
ix++) {
F(JZI, ix, iy, mz - 1) -= F(JZI, ix, iy, mz);
F(JZI, ix, iy, mz) = 0.;
F(JXI, ix, iy, mz - 1) += F(JXI, ix, iy, mz + 1);
F(JXI, ix, iy, mz + 1) = 0.;

F(JYI, ix, iy, mz - 1) += F(JYI, ix, iy, mz + 1);
F(JYI, ix, iy, mz + 1) = 0.;

F(JZI, ix, iy, mz - 1) -= F(JZI, ix, iy, mz);
F(JZI, ix, iy, mz) = 0.;
}
}
} else {
Expand Down
118 changes: 102 additions & 16 deletions src/libpsc/tests/test_reflective_bcs_integration.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,11 +49,11 @@ Grid_t* setupGrid()
{0, 0, 0}, // location of lower corner
{1, 1, 1}}; // n patches

auto bc =
psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL, BND_FLD_PERIODIC},
{BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL, BND_FLD_PERIODIC},
{BND_PRT_PERIODIC, BND_PRT_REFLECTING, BND_PRT_PERIODIC},
{BND_PRT_PERIODIC, BND_PRT_REFLECTING, BND_PRT_PERIODIC}};
auto bc = psc::grid::BC{
{BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL, BND_FLD_CONDUCTING_WALL},
{BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL, BND_FLD_CONDUCTING_WALL},
{BND_PRT_PERIODIC, BND_PRT_REFLECTING, BND_PRT_REFLECTING},
{BND_PRT_PERIODIC, BND_PRT_REFLECTING, BND_PRT_REFLECTING}};

auto kinds = Grid_t::Kinds(NR_KINDS);
kinds[KIND_ELECTRON] = {-1.0, 1.0, "e"};
Expand Down Expand Up @@ -83,12 +83,12 @@ Grid_t* setupGrid()
// ======================================================================
// Integration test (pun not intended)

TEST(ReflectiveBcsTest, Integration)
TEST(ReflectiveBcsTest, IntegrationY)
{
// ----------------------------------------------------------------------
// setup

psc_params.nmax = 100;
psc_params.nmax = 20;
psc_params.stats_every = 1;
psc_params.cfl = .75;

Expand Down Expand Up @@ -125,10 +125,10 @@ TEST(ReflectiveBcsTest, Integration)
auto injector = mprts.injector();
auto inj = injector[p];
double vy = 100; // fast enough to escape electric attraction
double y_center = grid.domain.length[1] / 2.0;
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_center, 5}, {0, -vy, 0}, 1, 0}); // electron
inj({{0, y_center, 5}, {0, vy, 0}, 1, 1}); // positron
inj({{0, y_quarter, 5}, {0, -vy, 0}, 1, 0}); // electron
inj({{0, y_quarter, 5}, {0, vy, 0}, 1, 1}); // positron
}

// ----------------------------------------------------------------------
Expand All @@ -151,24 +151,110 @@ TEST(ReflectiveBcsTest, Integration)
bool reflected = false;

for (; grid.timestep_ < psc_params.nmax; grid.timestep_++) {
about_to_reflect =
prts[0].x()[1] < grid.domain.dx[1] && prts[0].u()[1] < 0.0;

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);

if (prts[0].u()[1] > 0.0) {
if (prts[0].x()[1] < grid.domain.dx[1] && prts[0].u()[1] > 0.0) {
reflected = about_to_reflect;
break;
}
}

ASSERT_TRUE(reflected) << "timestep " << grid.timestep_;

ASSERT_EQ(prts.size(), 2);
ASSERT_GT(prts[0].u()[1], 0.0);
}

TEST(ReflectiveBcsTest, IntegrationZ)
{
// ----------------------------------------------------------------------
// 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<MfieldsState, Mparticles, Dim> outf{grid, {}};
OutputParticles outp{grid, {}};
DiagEnergies oute{grid.comm(), 0};
auto diagnostics = makeDiagnosticsDefault(outf, outp, oute);

auto psc =
makePscIntegrator<PscConfig>(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 vz = 100; // fast enough to escape electric attraction
double z_quarter = grid.domain.length[2] / 4.0;
// inject 2 particles at same location to satisfy Gauss' law at t=0
inj({{0, 5, z_quarter}, {0, 0, -vz}, 1, 0}); // electron
inj({{0, 5, z_quarter}, {0, 0, vz}, 1, 1}); // positron
}

// ----------------------------------------------------------------------
// run the simulation

auto accessor = mprts.accessor();
auto prts = accessor[p];

ASSERT_EQ(prts.size(), 2);
ASSERT_LT(prts[0].u()[2], 0.0);
ASSERT_GT(prts[1].u()[2], 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);

bool about_to_reflect = false;
bool reflected = false;

for (; grid.timestep_ < psc_params.nmax; grid.timestep_++) {
about_to_reflect =
prts[0].x()[1] < grid.domain.dx[1] / 2.0 && prts[0].u()[1] < 0.0;
prts[0].x()[2] < grid.domain.dx[2] && prts[0].u()[2] < 0.0;

psc.step();
ASSERT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold);
ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold);

if (prts[0].x()[2] < grid.domain.dx[2] && prts[0].u()[2] > 0.0) {
reflected = about_to_reflect;
break;
}
}

ASSERT_TRUE(reflected) << "timestep " << grid.timestep_;

ASSERT_EQ(prts.size(), 2);
ASSERT_GT(prts[0].u()[1], 0.0);
ASSERT_LT(prts[1].u()[1], 0.0);
ASSERT_GT(prts[0].u()[2], 0.0);
}

// ======================================================================
Expand Down