diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx new file mode 100644 index 000000000..1509dfcbb --- /dev/null +++ b/src/include/add_ghosts_reflecting.hxx @@ -0,0 +1,150 @@ +#pragma once + +#include "../kg/include/kg/Vec3.h" +#include "grid.hxx" + +template +void add_ghosts_reflecting_lo_cc(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_begin = ib; + Int3 idx_end = ldims - ib; + + idx_begin[d] = 0; + idx_end[d] = -ib[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] = -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; + } + } + } + } +} + +template +void add_ghosts_reflecting_hi_cc(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_begin = ib; + Int3 idx_end = ldims - 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; + } + } + } + } +} + +template +void add_ghosts_reflecting_lo_nc(const Int3& ldims, FE& mres_gt, const Int3& ib, + int p, int d, int mb, int me) +{ + // When lower boundary has 2 ghosts, upper boundary of nc grid only has 1 + // because nc has +1 point compared to cc, and takes it from upper ghost + // region. Thus, ignore a ghost from lower ghost region to balance it out. + + Int3 unused_ghost = {!!ib[0], !!ib[1], !!ib[2]}; + + Int3 idx_begin = ib + unused_ghost; + Int3 idx_end = ldims - ib; + + // Also note that the boundary nc point is the point about which reflections + // occur, and isn't itself reflected. + + idx_begin[d] = 1; + idx_end[d] = 1 - ib[d] - unused_ghost[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] = -idx[d]; + + 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; + } + } + } + } +} + +template +void add_ghosts_reflecting_hi_nc(const Int3& ldims, FE& mres_gt, const Int3& ib, + int p, int d, int mb, int me) +{ + // When lower boundary has 2 ghosts, upper boundary of nc grid only has 1 + // because nc has +1 point compared to cc, and takes it from upper ghost + // region. Thus, ignore a ghost from lower ghost region to balance it out. + + Int3 unused_ghost = {!!ib[0], !!ib[1], !!ib[2]}; + + Int3 idx_begin = ib + unused_ghost; + Int3 idx_end = ldims - ib; + + // Also note that the extra nc point is the point about which reflections + // occur, and isn't itself reflected. + + idx_begin[d] = ldims[d] + ib[d] + unused_ghost[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]; + + 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; + } + } + } + } +} diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index cbcbab3ba..2a4cbe857 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -3,10 +3,12 @@ #include "../libpsc/psc_bnd/psc_bnd_impl.hxx" #include "bnd.hxx" +#include "centering.hxx" #include "fields.hxx" #include "fields3d.hxx" #include "particles.hxx" #include "psc_fields_c.h" +#include "add_ghosts_reflecting.hxx" #include @@ -32,118 +34,89 @@ inline std::vector addKindSuffix( // ====================================================================== // ItemMomentBnd -template -class ItemMomentBnd +template +class DomainBoundary { public: - using storage_type = S; - - ItemMomentBnd(const Grid_t& grid) {} - - void add_ghosts(const Grid_t& grid, storage_type& mres_gt, const Int3& ib) - { - for (int p = 0; p < mres_gt.shape(4); p++) { - add_ghosts_boundary(grid, mres_gt, ib, p, 0, mres_gt.shape(3)); - } - - bnd_.add_ghosts(grid, mres_gt, ib, 0, mres_gt.shape(3)); - } - -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); - } - } + static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, + const Int3& ib, int p, int mb, int me); +}; +template <> +class DomainBoundary +{ +public: template - void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, - int p, int d, int mb, int me) + static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, + const Int3& ib, int p, 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); - } - } - } + // lo + for (int d = 0; d < 3; d++) { + // FIXME why reflect for open BCs? + if (grid.atBoundaryLo(p, d) && (grid.bc.prt_lo[d] == BND_PRT_REFLECTING || + grid.bc.prt_lo[d] == BND_PRT_OPEN)) { + add_ghosts_reflecting_lo_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } - } 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); - } - } - } + } + // hi + for (int d = 0; d < 3; d++) { + // FIXME why reflect for open BCs? + if (grid.atBoundaryHi(p, d) && (grid.bc.prt_hi[d] == BND_PRT_REFLECTING || + grid.bc.prt_hi[d] == BND_PRT_OPEN)) { + add_ghosts_reflecting_hi_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } - } else { - assert(0); } } +}; +template <> +class DomainBoundary +{ +public: template - void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, - int p, int mb, int me) + static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, + const Int3& ib, int p, int mb, int me) { // lo for (int d = 0; d < 3; d++) { - if (grid.atBoundaryLo(p, d)) { - if (grid.bc.prt_lo[d] == BND_PRT_REFLECTING || - grid.bc.prt_lo[d] == BND_PRT_OPEN) { - add_ghosts_reflecting_lo(grid.ldims, mres_gt, ib, p, d, mb, me); - } + // FIXME why reflect for open BCs? + if (grid.atBoundaryLo(p, d) && (grid.bc.prt_lo[d] == BND_PRT_REFLECTING || + grid.bc.prt_lo[d] == BND_PRT_OPEN)) { + add_ghosts_reflecting_lo_nc(grid.ldims, mres_gt, ib, p, d, mb, me); } } // hi for (int d = 0; d < 3; d++) { - if (grid.atBoundaryHi(p, d)) { - if (grid.bc.prt_hi[d] == BND_PRT_REFLECTING || - grid.bc.prt_hi[d] == BND_PRT_OPEN) { - add_ghosts_reflecting_hi(grid.ldims, mres_gt, ib, p, d, mb, me); - } + // FIXME why reflect for open BCs? + if (grid.atBoundaryHi(p, d) && (grid.bc.prt_hi[d] == BND_PRT_REFLECTING || + grid.bc.prt_hi[d] == BND_PRT_OPEN)) { + add_ghosts_reflecting_hi_nc(grid.ldims, mres_gt, ib, p, d, mb, me); } } } +}; + +// ====================================================================== +// ItemMomentBnd + +template +class ItemMomentBnd +{ +public: + using storage_type = S; + + ItemMomentBnd(const Grid_t& grid) {} + + void add_ghosts(const Grid_t& grid, storage_type& mres_gt, const Int3& ib) + { + for (int p = 0; p < mres_gt.shape(4); p++) { + DomainBoundary::add_ghosts_boundary(grid, mres_gt, ib, p, 0, + mres_gt.shape(3)); + } + + bnd_.add_ghosts(grid, mres_gt, ib, 0, mres_gt.shape(3)); + } private: Bnd bnd_; @@ -189,6 +162,6 @@ public: } private: - ItemMomentBnd bnd_; + ItemMomentBnd bnd_; std::vector comp_names_; }; diff --git a/src/include/psc/deposit.hxx b/src/include/psc/deposit.hxx index cd386a099..af0a89f3f 100644 --- a/src/include/psc/deposit.hxx +++ b/src/include/psc/deposit.hxx @@ -1,6 +1,8 @@ #pragma once +#include "centering.hxx" + #include #include #include @@ -156,12 +158,21 @@ public: } }; +// ---------------------------------------------------------------------------- + +template +class DepositCentering +{ +public: + static const centering::Centering CENTERING = C; +}; + // ---------------------------------------------------------------------------- // deposit directly, assuming grid spacing == 1, pass in patch-relative position // x template -class Deposit1stNc +class Deposit1stNc : public DepositCentering { public: static std::string suffix() { return "_1st_nc"; } @@ -184,7 +195,7 @@ public: }; template -class Deposit1stCc +class Deposit1stCc : public DepositCentering { public: static std::string suffix() { return "_1st_cc"; } @@ -211,7 +222,7 @@ public: // x template -class Deposit2ndNc +class Deposit2ndNc : public DepositCentering { public: static std::string suffix() { return "_2nd_nc"; } @@ -259,15 +270,17 @@ namespace code // deposition in code units // units template class DepositNorm> + template class DEPOSIT_NORM> class Deposit { public: using real_t = R; using dim_t = D; using real3_t = gt::sarray; + using DepositNorm = DEPOSIT_NORM; + static const centering::Centering CENTERING = DepositNorm::CENTERING; - static std::string suffix() { return DepositNorm::suffix(); } + static std::string suffix() { return DepositNorm::suffix(); } Deposit(const real3_t& dx, real_t fnqs) : dxi_{real_t(1.) / dx}, fnqs_{fnqs} {} @@ -279,7 +292,7 @@ public: real3_t x = xi * dxi_; real_t value = fnqs_ * val; - DepositNorm deposit; + DepositNorm deposit; deposit(flds, ib, x, value); } diff --git a/src/include/psc/moment.hxx b/src/include/psc/moment.hxx index 70f4a0564..39672f7ed 100644 --- a/src/include/psc/moment.hxx +++ b/src/include/psc/moment.hxx @@ -1,6 +1,8 @@ #pragma once +#include "centering.hxx" + #include "psc_bits.h" #include @@ -120,6 +122,8 @@ class moment_n { public: using dim_t = D; + static const centering::Centering CENTERING = + DepositCode::CENTERING; static std::string name() { return "n" + DepositCode::suffix(); } static std::vector comp_names(const Grid_t::Kinds& kinds) @@ -146,6 +150,8 @@ class moment_rho { public: using dim_t = D; + static const centering::Centering CENTERING = + DepositCode::CENTERING; static std::string name() { return "rho" + DepositCode::suffix(); } static std::vector comp_names(const Grid_t::Kinds& kinds) @@ -171,6 +177,8 @@ class moment_v { public: using dim_t = D; + static const centering::Centering CENTERING = + DepositCode::CENTERING; static std::string name() { return "v" + DepositCode::suffix(); } static std::vector comp_names(const Grid_t::Kinds& kinds) @@ -201,6 +209,8 @@ class moment_p { public: using dim_t = D; + static const centering::Centering CENTERING = + DepositCode::CENTERING; static std::string name() { return "p" + DepositCode::suffix(); } static std::vector comp_names(const Grid_t::Kinds& kinds) @@ -228,6 +238,8 @@ class moment_T { public: using dim_t = D; + static const centering::Centering CENTERING = + DepositCode::CENTERING; static std::string name() { return "T" + DepositCode::suffix(); } static std::vector comp_names(const Grid_t::Kinds& kinds) @@ -261,6 +273,8 @@ class moment_all { public: using dim_t = D; + static const centering::Centering CENTERING = + DepositCode::CENTERING; static std::string name() { return "all" + DepositCode::suffix(); } static std::vector comp_names(const Grid_t::Kinds& kinds) 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..3b6b89533 100644 --- a/src/libpsc/tests/CMakeLists.txt +++ b/src/libpsc/tests/CMakeLists.txt @@ -1,14 +1,16 @@ set(MPIEXEC_COMMAND ${MPIEXEC_EXECUTABLE} -# ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} + + # ${MPIEXEC_NUMPROC_FLAG} ${MPIEXEC_MAX_NUMPROCS} ) macro(add_psc_test name) add_executable(${name}) target_gtensor_sources(${name} PRIVATE ${name}.cxx) target_link_libraries(${name} psc GTest::gtest) - if (USE_GTEST_DISCOVER_TESTS) + + if(USE_GTEST_DISCOVER_TESTS) gtest_discover_tests(${name} EXEC_WRAPPER ${MPIEXEC_COMMAND}) else() gtest_add_tests(TARGET ${name} EXEC_WRAPPER ${MPIEXEC_COMMAND}) @@ -21,6 +23,8 @@ 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_reflective_bcs_integration) add_psc_test(test_mfields) add_psc_test(test_mfields_cuda) add_psc_test(test_bnd) @@ -32,9 +36,11 @@ add_psc_test(test_push_particles_2) add_psc_test(test_push_fields) add_psc_test(test_moments) add_psc_test(test_collision) -if (USE_CUDA AND NOT USE_VPIC) + +if(USE_CUDA AND NOT USE_VPIC) add_psc_test(test_collision_cuda) endif() + add_psc_test(test_balance) add_psc_test(TestUniqueIdGenerator) add_psc_test(test_mfields_io) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx new file mode 100644 index 000000000..2a321e066 --- /dev/null +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -0,0 +1,430 @@ +#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}; +} + +// ====================================================================== +// Unit tests + +int ravel_idx(Int3 idx, Int3 shape) +{ + return idx[0] + shape[0] * (idx[1] + shape[1] * (idx[2])); +} + +template +void init_mres(FE& mres_gt, Int3 shape, int p) +{ + 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_gt(i, j, k, 0, p) = cell_id; + } + } + } +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingHighCcY) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 1; + add_ghosts_reflecting_hi_cc(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, AddGhostsReflectingLowCcY) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 1; + add_ghosts_reflecting_lo_cc(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 = 0; y < by; 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 = -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, AddGhostsReflectingHighCcZ) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 2; + add_ghosts_reflecting_hi_cc(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, AddGhostsReflectingLowCcZ) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 2; + add_ghosts_reflecting_lo_cc(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 = 0; z < bz; z++) { + int cell_id = ravel_idx({bx + x, by + y, bz + z}, shape); + + int reflected_z = -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, AddGhostsReflectingHighNcY) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 1; + add_ghosts_reflecting_hi_nc(ldims, mres, ib, p, dim, 0, 1); + + int bx_cc = -ib[0] - !grid.isInvar(0); + int by_cc = -ib[1] - !grid.isInvar(1); + int bz_cc = -ib[2] - !grid.isInvar(2); + Int3 ldims_cc = + ldims + Int3{!grid.isInvar(0), !grid.isInvar(1), !grid.isInvar(2)}; + + for (int x = -bx_cc; x < ldims_cc[0] + bx_cc; x++) { + for (int y = ldims_cc[1] - by_cc - 1; y < ldims_cc[1] - 1; y++) { + for (int z = -bz_cc; z < ldims_cc[2] + bz_cc; z++) { + int cell_id = ravel_idx({x - ib[0], y - ib[1], z - ib[2]}, shape); + + int reflected_y = 2 * (ldims_cc[1] - 1) - y; + int reflected_cell_id = + ravel_idx({x - ib[0], reflected_y - ib[1], z - ib[2]}, shape); + EXPECT_EQ(mres(x - ib[0], y - ib[1], z - ib[2], 0, p), + cell_id + reflected_cell_id) + << "xyz " << x << " " << y << " " << z; + } + } + } +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingLowNcY) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 1; + add_ghosts_reflecting_lo_nc(ldims, mres, ib, p, dim, 0, 1); + + int bx_cc = -ib[0] - !grid.isInvar(0); + int by_cc = -ib[1] - !grid.isInvar(1); + int bz_cc = -ib[2] - !grid.isInvar(2); + Int3 ldims_cc = + ldims + Int3{!grid.isInvar(0), !grid.isInvar(1), !grid.isInvar(2)}; + + for (int x = -bx_cc; x < ldims_cc[0] + bx_cc; x++) { + for (int y = 1; y < 1 + by_cc; y++) { + for (int z = -bz_cc; z < ldims_cc[2] + bz_cc; z++) { + int cell_id = ravel_idx({x - ib[0], y - ib[1], z - ib[2]}, shape); + + int reflected_y = -y; + int reflected_cell_id = + ravel_idx({x - ib[0], reflected_y - ib[1], z - ib[2]}, shape); + EXPECT_EQ(mres(x - ib[0], y - ib[1], z - ib[2], 0, p), + cell_id + reflected_cell_id) + << "xyz " << x << " " << y << " " << z; + } + } + } +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingHighNcZ) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 2; + add_ghosts_reflecting_hi_nc(ldims, mres, ib, p, dim, 0, 1); + + int bx_cc = -ib[0] - !grid.isInvar(0); + int by_cc = -ib[1] - !grid.isInvar(1); + int bz_cc = -ib[2] - !grid.isInvar(2); + Int3 ldims_cc = + ldims + Int3{!grid.isInvar(0), !grid.isInvar(1), !grid.isInvar(2)}; + + for (int x = -bx_cc; x < ldims_cc[0] + bx_cc; x++) { + for (int y = -by_cc; y < ldims_cc[2] + by_cc; y++) { + for (int z = ldims_cc[2] - bz_cc - 1; z < ldims_cc[2] - 1; z++) { + int cell_id = ravel_idx({x - ib[0], y - ib[1], z - ib[2]}, shape); + + int reflected_z = 2 * (ldims_cc[2] - 1) - z; + int reflected_cell_id = + ravel_idx({x - ib[0], y - ib[1], reflected_z - ib[2]}, shape); + EXPECT_EQ(mres(x - ib[0], y - ib[1], z - ib[2], 0, p), + cell_id + reflected_cell_id) + << "xyz " << x << " " << y << " " << z; + } + } + } +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingLowNcZ) +{ + 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; + + init_mres(mres, shape, p); + + int dim = 2; + add_ghosts_reflecting_lo_nc(ldims, mres, ib, p, dim, 0, 1); + + int bx_cc = -ib[0] - !grid.isInvar(0); + int by_cc = -ib[1] - !grid.isInvar(1); + int bz_cc = -ib[2] - !grid.isInvar(2); + Int3 ldims_cc = + ldims + Int3{!grid.isInvar(0), !grid.isInvar(1), !grid.isInvar(2)}; + + for (int x = -bx_cc; x < ldims_cc[0] + bx_cc; x++) { + for (int y = -by_cc; y < ldims_cc[1] + by_cc; y++) { + for (int z = 1; z < 1 + bz_cc; z++) { + int cell_id = ravel_idx({x - ib[0], y - ib[1], z - ib[2]}, shape); + + int reflected_z = -z; + int reflected_cell_id = + ravel_idx({x - ib[0], y - ib[1], reflected_z - ib[2]}, shape); + EXPECT_EQ(mres(x - ib[0], y - ib[1], z - ib[2], 0, p), + cell_id + reflected_cell_id) + << "xyz " << x << " " << y << " " << z; + } + } + } +} + +// ====================================================================== +// main + +int main(int argc, char** argv) +{ + MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + MPI_Finalize(); + return rc; +} diff --git a/src/libpsc/tests/test_reflective_bcs_integration.cxx b/src/libpsc/tests/test_reflective_bcs_integration.cxx new file mode 100644 index 000000000..1dd70e6d2 --- /dev/null +++ b/src/libpsc/tests/test_reflective_bcs_integration.cxx @@ -0,0 +1,184 @@ +#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}; +} + +// ====================================================================== +// Integration test (pun not intended) + +TEST(ReflectiveBcsTest, Integration) +{ + // ---------------------------------------------------------------------- + // setup + + psc_params.nmax = 100; + 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_center = grid.domain.length[1] / 2.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 + } + + // ---------------------------------------------------------------------- + // 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); + + bool about_to_reflect = false; + bool reflected = false; + + for (; grid.timestep_ < psc_params.nmax; grid.timestep_++) { + psc.step(); + EXPECT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + EXPECT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); + + if (prts[0].u()[1] > 0.0) { + reflected = about_to_reflect; + break; + } + + about_to_reflect = + prts[0].x()[1] < grid.domain.dx[1] / 2.0 && prts[0].u()[1] < 0.0; + } + + 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); +} + +// ====================================================================== +// main + +int main(int argc, char** argv) +{ + psc_init(argc, argv); + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + psc_finalize(); + return rc; +}