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 85e49486ea..b8ee7318a7 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -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) { @@ -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 { @@ -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); } @@ -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 { @@ -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); } } @@ -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 { @@ -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) { @@ -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 { @@ -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.; } @@ -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 { @@ -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.; } @@ -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 { diff --git a/src/libpsc/tests/test_reflective_bcs_integration.cxx b/src/libpsc/tests/test_reflective_bcs_integration.cxx index 1dd70e6d21..44d4901d03 100644 --- a/src/libpsc/tests/test_reflective_bcs_integration.cxx +++ b/src/libpsc/tests/test_reflective_bcs_integration.cxx @@ -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"}; @@ -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; @@ -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 } // ---------------------------------------------------------------------- @@ -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 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 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); } // ======================================================================