diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx new file mode 100644 index 000000000..c318ba33b --- /dev/null +++ b/src/include/add_ghosts_reflecting.hxx @@ -0,0 +1,71 @@ +#pragma once + +#include "../kg/include/kg/Vec3.h" +#include "grid.hxx" + +template +void add_ghosts_reflecting_lo(const Int3& ldims, FE& mres_gt, const Int3& ib, + int p, int d, int mb, int me) +{ + int bx = ldims[0] == 1 ? 0 : 1; + if (d == 1) { + for (int iz = -1; iz < ldims[2] + 1; iz++) { + for (int ix = -bx; ix < ldims[0] + bx; ix++) { + int iy = 0; + { + for (int m = mb; m < me; m++) { + mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += + mres_gt(ix - ib[0], iy - 1 - ib[1], iz - ib[2], m, p); + } + } + } + } + } else if (d == 2) { + for (int iy = 0 * -1; iy < ldims[1] + 0 * 1; iy++) { + for (int ix = -bx; ix < ldims[0] + bx; ix++) { + int iz = 0; + { + for (int m = mb; m < me; m++) { + mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += + mres_gt(ix - ib[0], iy - ib[1], iz - 1 - ib[2], m, p); + } + } + } + } + } else { + assert(0); + } +} + +template +void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, + int p, int d, int mb, int me) +{ + // FIXME only need to scan 1 cell into the ghost region for 1st-order moments + + Int3 idx_end = ldims - ib; + Int3 idx_begin = ib; + + idx_begin[d] = ldims[d] + ib[d]; + idx_end[d] = ldims[d]; + + Int3 idx; + + for (int m = mb; m < me; m++) { + for (idx[2] = idx_begin[2]; idx[2] < idx_end[2]; idx[2]++) { + for (idx[1] = idx_begin[1]; idx[1] < idx_end[1]; idx[1]++) { + for (idx[0] = idx_begin[0]; idx[0] < idx_end[0]; idx[0]++) { + Int3 idx_reflected = idx; + idx_reflected[d] = 2 * ldims[d] - idx[d] - 1; + + auto val_reflected = + mres_gt(idx_reflected[0] - ib[0], idx_reflected[1] - ib[1], + idx_reflected[2] - ib[2], m, p); + + mres_gt(idx[0] - ib[0], idx[1] - ib[1], idx[2] - ib[2], m, p) += + val_reflected; + } + } + } + } +} \ No newline at end of file diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index cbcbab3ba..854ce4c7a 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -7,6 +7,7 @@ #include "fields3d.hxx" #include "particles.hxx" #include "psc_fields_c.h" +#include "add_ghosts_reflecting.hxx" #include @@ -53,74 +54,6 @@ private: // ---------------------------------------------------------------------- // boundary stuff FIXME, should go elsewhere... - template - void add_ghosts_reflecting_lo(const Int3& ldims, FE& mres_gt, const Int3& ib, - int p, int d, int mb, int me) - { - int bx = ldims[0] == 1 ? 0 : 1; - if (d == 1) { - for (int iz = -1; iz < ldims[2] + 1; iz++) { - for (int ix = -bx; ix < ldims[0] + bx; ix++) { - int iy = 0; - { - for (int m = mb; m < me; m++) { - mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += - mres_gt(ix - ib[0], iy - 1 - ib[1], iz - ib[2], m, p); - } - } - } - } - } else if (d == 2) { - for (int iy = 0 * -1; iy < ldims[1] + 0 * 1; iy++) { - for (int ix = -bx; ix < ldims[0] + bx; ix++) { - int iz = 0; - { - for (int m = mb; m < me; m++) { - mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += - mres_gt(ix - ib[0], iy - ib[1], iz - 1 - ib[2], m, p); - } - } - } - } - } else { - assert(0); - } - } - - template - void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, - int p, int d, int mb, int me) - { - int bx = ldims[0] == 1 ? 0 : 1; - if (d == 1) { - for (int iz = -1; iz < ldims[2] + 1; iz++) { - for (int ix = -bx; ix < ldims[0] + bx; ix++) { - int iy = ldims[1] - 1; - { - for (int m = mb; m < me; m++) { - mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += - mres_gt(ix - ib[0], iy - 1 + ib[1], iz - ib[2], m, p); - } - } - } - } - } else if (d == 2) { - for (int iy = 0 * -1; iy < ldims[1] + 0 * 1; iy++) { - for (int ix = -bx; ix < ldims[0] + bx; ix++) { - int iz = ldims[2] - 1; - { - for (int m = mb; m < me; m++) { - mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += - mres_gt(ix - ib[0], iy - ib[1], iz + 1 - ib[2], m, p); - } - } - } - } - } else { - assert(0); - } - } - template void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, int p, int mb, int me) diff --git a/src/libpsc/psc_checks/checks_impl.hxx b/src/libpsc/psc_checks/checks_impl.hxx index 33ddbe55e..85dc83ad3 100644 --- a/src/libpsc/psc_checks/checks_impl.hxx +++ b/src/libpsc/psc_checks/checks_impl.hxx @@ -209,11 +209,10 @@ struct checks_order_2nd using Moment_rho_nc = Moment_rho_2nd_nc; }; -template +template class ChecksCommon { public: - using Mparticles = MP; using storage_type = S; using item_rho_type = ITEM_RHO; @@ -228,5 +227,5 @@ public: template using Checks_ = - ChecksCommon>; diff --git a/src/libpsc/tests/CMakeLists.txt b/src/libpsc/tests/CMakeLists.txt index 1c29b8f49..35f990dc5 100644 --- a/src/libpsc/tests/CMakeLists.txt +++ b/src/libpsc/tests/CMakeLists.txt @@ -21,6 +21,7 @@ add_psc_test(test_rng) add_psc_test(test_mparticles_cuda) add_psc_test(test_mparticles) add_psc_test(test_output_particles) +add_psc_test(test_reflective_bcs) add_psc_test(test_mfields) add_psc_test(test_mfields_cuda) add_psc_test(test_bnd) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx new file mode 100644 index 000000000..5f839c6d1 --- /dev/null +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -0,0 +1,307 @@ +#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 BgkMfields = PscConfig::Mfields; +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_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 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}; +} + +// ====================================================================== +// run + +static void run(int argc, char** argv) +{ + mpi_printf(MPI_COMM_WORLD, "*** Setting up...\n"); + + // ---------------------------------------------------------------------- + // setup various parameters first + + psc_params.nmax = 100; + psc_params.stats_every = 1; + psc_params.cfl = .75; + + // ---------------------------------------------------------------------- + // Set up grid, state fields, particles + + auto grid_ptr = setupGrid(); + auto& grid = *grid_ptr; + MfieldsState mflds{grid}; + Mparticles mprts{grid}; + + // ---------------------------------------------------------------------- + // Set up various objects needed to run this case + + // -- Balance + psc_params.balance_interval = 0; + Balance balance{.1}; + + // -- Sort + psc_params.sort_interval = 100; + + // -- Collision + int collision_interval = 0; + double collision_nu = .1; + Collision collision{grid, collision_interval, collision_nu}; + + // -- Checks + ChecksParams checks_params{}; + checks_params.gauss.check_interval = 1; + checks_params.continuity.check_interval = 1; + + Checks checks{grid, MPI_COMM_WORLD, checks_params}; + + // -- Marder correction + double marder_diffusion = 0.9; + int marder_loop = 3; + bool marder_dump = false; + psc_params.marder_interval = -1; + Marder marder(grid, marder_diffusion, marder_loop, marder_dump); + + // ---------------------------------------------------------------------- + // Set up output + + // -- output fields + OutputFieldsParams outf_params{}; + outf_params.fields.pfield.out_interval = -1; + outf_params.moments.pfield.out_interval = -1; + OutputFields outf{grid, outf_params}; + + // -- output particles + OutputParticlesParams outp_params{}; + outp_params.every_step = -1; + OutputParticles outp{grid, outp_params}; + + int oute_interval = -100; + DiagEnergies oute{grid.comm(), oute_interval}; + + auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + + // ---------------------------------------------------------------------- + // set up initial conditions + + { + auto injector = mprts.injector(); + auto inj = injector[0]; + inj({{0, 15, 5}, {0, -.001, 0}, 1, 0}); + } + + // ---------------------------------------------------------------------- + // run the simulation + + auto psc = + makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, + collision, checks, marder, diagnostics); + + for (; grid.timestep_ < psc_params.nmax; grid.timestep_++) { + psc.step(); + ASSERT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); + } +} + +// ====================================================================== +// Integration test (pun not intended) + +TEST(ReflectiveBcsTest, Integration) +{ + int argc = 1; + char** argv = nullptr; + + run(argc, argv); +} + +// ====================================================================== +// Unit tests + +int ravel_idx(Int3 idx, Int3 shape) +{ + return idx[0] + shape[0] * (idx[1] + shape[1] * (idx[2])); +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) +{ + Grid_t* grid_ptr = setupGrid(); + auto& grid = *grid_ptr; + + Int3 ldims = grid.ldims; + Int3 shape = ldims + 2 * grid.ibn; + Int3 ib = -grid.ibn; + + auto mres = psc::mflds::zeros(grid, 1, ib); + + EXPECT_EQ(grid.n_patches(), 1); + int p = 0; + + for (int i = 0; i < shape[0]; i++) { + for (int j = 0; j < shape[1]; j++) { + for (int k = 0; k < shape[2]; k++) { + int cell_id = ravel_idx({i, j, k}, shape); + mres(i, j, k, 0, p) = cell_id; + } + } + } + + int dim = 1; + add_ghosts_reflecting_hi(ldims, mres, ib, p, dim, 0, 1); + + int bx = -ib[0]; + int by = -ib[1]; + int bz = -ib[2]; + for (int x = -bx; x < ldims[0] + bx; x++) { + for (int y = ldims[1] - by; y < ldims[1]; y++) { + for (int z = -bz; z < ldims[2] + bz; z++) { + int cell_id = ravel_idx({bx + x, by + y, bz + z}, shape); + + int reflected_y = 2 * ldims[1] - y - 1; + int reflected_cell_id = + ravel_idx({bx + x, by + reflected_y, bz + z}, shape); + EXPECT_EQ(mres(bx + x, by + y, bz + z, 0, p), + cell_id + reflected_cell_id) + << "xyz " << x << " " << y << " " << z; + } + } + } +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) +{ + // TODO + EXPECT_TRUE(false); +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) +{ + Grid_t* grid_ptr = setupGrid(); + auto& grid = *grid_ptr; + + Int3 ldims = grid.ldims; + Int3 shape = ldims + 2 * grid.ibn; + Int3 ib = -grid.ibn; + + auto mres = psc::mflds::zeros(grid, 1, ib); + + EXPECT_EQ(grid.n_patches(), 1); + int p = 0; + + for (int i = 0; i < shape[0]; i++) { + for (int j = 0; j < shape[1]; j++) { + for (int k = 0; k < shape[2]; k++) { + int cell_id = ravel_idx({i, j, k}, shape); + mres(i, j, k, 0, p) = cell_id; + } + } + } + + int dim = 2; + add_ghosts_reflecting_hi(ldims, mres, ib, p, dim, 0, 1); + + int bx = -ib[0]; + int by = -ib[1]; + int bz = -ib[2]; + for (int x = -bx; x < ldims[0] + bx; x++) { + for (int y = -by; y < ldims[1] + by; y++) { + for (int z = ldims[2] - bz; z < ldims[2]; z++) { + int cell_id = ravel_idx({bx + x, by + y, bz + z}, shape); + + int reflected_z = 2 * ldims[2] - z - 1; + int reflected_cell_id = + ravel_idx({bx + x, by + y, bz + reflected_z}, shape); + EXPECT_EQ(mres(bx + x, by + y, bz + z, 0, p), + cell_id + reflected_cell_id) + << "xyz " << x << " " << y << " " << z; + } + } + } +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingLowZ) +{ + // TODO + EXPECT_TRUE(false); +} + +// ====================================================================== +// main + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + MPI_Finalize(); + return rc; +}