From de449ddc8818a8af84c3236e88fb6d52cf9251b1 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 31 Mar 2025 13:08:51 -0400 Subject: [PATCH 01/56] +test_reflective_bcs: copy from psc_shock for now --- src/libpsc/tests/CMakeLists.txt | 1 + src/libpsc/tests/test_reflective_bcs.cxx | 238 +++++++++++++++++++++++ 2 files changed, 239 insertions(+) create mode 100644 src/libpsc/tests/test_reflective_bcs.cxx 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..92637f63d --- /dev/null +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -0,0 +1,238 @@ +#include "gtest/gtest.h" + +#include +#include +#include + +#include "DiagnosticsDefault.h" +#include "OutputFieldsDefault.h" +#include "../psc_config.hxx" + +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 +{ +// General PSC parameters +PscParams psc_params; + +int nx = 1; +int ny = 4; +int nz = 2; +int nt = 1; + +int n_patches_x = 1; +int n_patches_y = 1; +int n_patches_z = 1; + +double dx = 10; +double dy = 10; +double dz = 10; + +double len_x = nx * dx; +double len_y = ny * dy; +double len_z = nz * dz; + +int n_writes = 1; +int out_interval = nt / n_writes; + +} // namespace + +// ====================================================================== +// setupParameters + +void setupParameters(int argc, char** argv) +{ + psc_params.nmax = nt; + psc_params.stats_every = out_interval; + psc_params.cfl = .75; + + psc_params.write_checkpoint_every_step = 0; +} + +// ====================================================================== +// setupGrid +// +// This helper function is responsible for setting up the "Grid", +// which is really more than just the domain and its decomposition, it +// also encompasses PC normalization parameters, information about the +// particle kinds, etc. + +Grid_t* setupGrid() +{ + // FIXME add a check to catch mismatch between Dim and n grid points early + auto domain = + Grid_t::Domain{{nx, ny, nz}, // n grid points + {len_x, len_y, len_z}, // physical lengths + {0, 0, 0}, // location of lower corner + {n_patches_x, n_patches_y, n_patches_z}}; // 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 = 100; + + 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}; +} + +// ====================================================================== +// initializeParticles + +void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) +{ + auto injector = mprts.injector(); + injector[0]({{0.0, 5.0, 5.0}, {0.0, 0.0, 0.0}, 1.0, 0}); +} + +// ====================================================================== +// run + +static void run(int argc, char** argv) +{ + mpi_printf(MPI_COMM_WORLD, "*** Setting up...\n"); + + // ---------------------------------------------------------------------- + // setup various parameters first + + setupParameters(argc, argv); + + // ---------------------------------------------------------------------- + // Set up grid, state fields, particles + + auto grid_ptr = setupGrid(); + auto& grid = *grid_ptr; + MfieldsState mflds{grid}; + Mparticles mprts{grid}; + BgkMfields phi{grid, 1, mflds.ibn()}; + BgkMfields gradPhi{grid, 3, mflds.ibn()}; + BgkMfields divGradPhi{grid, 1, mflds.ibn()}; + + // ---------------------------------------------------------------------- + // 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 = out_interval; + outf_params.moments.pfield.out_interval = out_interval; + OutputFields outf{grid, outf_params}; + + // -- output particles + OutputParticlesParams outp_params{}; + outp_params.every_step = out_interval; + outp_params.data_dir = "."; + outp_params.basename = "prt"; + 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 + + initializeParticles(balance, grid_ptr, mprts); + + // ---------------------------------------------------------------------- + // run the simulation + + auto psc = + makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, + collision, checks, marder, diagnostics); + + psc.integrate(); +} + +// ====================================================================== +// main + +// ====================================================================== +// Test1 + +TEST(ReflectiveBcsTest, Test1) +{ + // psc_init(argc, argv); + // FIXME restore whatever previous functionality there was with options + int argc = 1; + char** argv = nullptr; + psc_init(argc, argv); + + run(argc, argv); + + psc_finalize(); +} + +int main(int argc, char** argv) +{ + // MPI_Init(&argc, &argv); + ::testing::InitGoogleTest(&argc, argv); + int rc = RUN_ALL_TESTS(); + // MPI_Finalize(); + return rc; +} From 0af729de98de2ff0b63aa556e5b043a6096521be Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 31 Mar 2025 13:16:41 -0400 Subject: [PATCH 02/56] test_reflective_bcs: bounce particle off wall --- src/libpsc/tests/test_reflective_bcs.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 92637f63d..b73f306d1 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -33,7 +33,7 @@ PscParams psc_params; int nx = 1; int ny = 4; int nz = 2; -int nt = 1; +int nt = 10; int n_patches_x = 1; int n_patches_y = 1; @@ -47,7 +47,7 @@ double len_x = nx * dx; double len_y = ny * dy; double len_z = nz * dz; -int n_writes = 1; +int n_writes = nt; int out_interval = nt / n_writes; } // namespace @@ -118,7 +118,7 @@ Grid_t* setupGrid() void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) { auto injector = mprts.injector(); - injector[0]({{0.0, 5.0, 5.0}, {0.0, 0.0, 0.0}, 1.0, 0}); + injector[0]({{0.0, 36.0, 5.0}, {0.0, 3.0, 0.0}, 1.0, 0}); } // ====================================================================== From a1a0f70cdf63d5490a7d8ab0bdee4a2cd40808cd Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 10:17:43 -0400 Subject: [PATCH 03/56] test_reflective_bcs: add slight z velocity --- src/libpsc/tests/test_reflective_bcs.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index b73f306d1..62e30d90e 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -118,7 +118,7 @@ Grid_t* setupGrid() void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) { auto injector = mprts.injector(); - injector[0]({{0.0, 36.0, 5.0}, {0.0, 3.0, 0.0}, 1.0, 0}); + injector[0]({{0.0, 36.0, 5.0}, {0.0, 3.0, 1.0}, 1.0, 0}); } // ====================================================================== From 62ef3522b4e72ef6e43fb24abf9309adac9515c4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 11:02:26 -0400 Subject: [PATCH 04/56] test_reflective_bcs: bring back full shocks setup that fails it has to run on debug mode, because otherwise the asserts aren't compiled --- src/libpsc/tests/test_reflective_bcs.cxx | 80 +++++++++++++++++++++--- 1 file changed, 71 insertions(+), 9 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 62e30d90e..24cdb26df 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -8,8 +8,20 @@ #include "OutputFieldsDefault.h" #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; +#ifdef USE_CUDA +using PscConfig = PscConfig1vbecCuda; +#else using PscConfig = PscConfig1vbecDouble; +#endif // ---------------------------------------------------------------------- @@ -30,24 +42,33 @@ namespace // General PSC parameters PscParams psc_params; +// matching: heliospheric .1 AU 3, shock tiny 2, real c and M, parallel B + +double electron_temperature = 7.8278E-05; +double ion_temperature = 7.8278E-05; +double ion_mass = 1.8360E+03; +double v_flow = 2.6685E-03; +double b_x = 0; +double b_y = 2.4000E-02; + int nx = 1; -int ny = 4; +int ny = 128; int nz = 2; -int nt = 10; +int nt = 100; int n_patches_x = 1; int n_patches_y = 1; int n_patches_z = 1; -double dx = 10; -double dy = 10; -double dz = 10; +double dx = 4.2849E+01; +double dy = 3.3475E-01; +double dz = 3.3475E-01; double len_x = nx * dx; double len_y = ny * dy; double len_z = nz * dz; -int n_writes = nt; +int n_writes = 10; int out_interval = nt / n_writes; } // namespace @@ -89,7 +110,7 @@ Grid_t* setupGrid() auto kinds = Grid_t::Kinds(NR_KINDS); kinds[KIND_ELECTRON] = {-1.0, 1.0, "e"}; - kinds[KIND_ION] = {1.0, 1.0, "i"}; + kinds[KIND_ION] = {1.0, ion_mass, "i"}; // --- generic setup auto norm_params = Grid_t::NormalizationParams::dimensionless(); @@ -117,8 +138,48 @@ Grid_t* setupGrid() void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) { - auto injector = mprts.injector(); - injector[0]({{0.0, 36.0, 5.0}, {0.0, 3.0, 1.0}, 1.0, 0}); + SetupParticles setup_particles(*grid_ptr); + setup_particles.centerer = Centering::Centerer(Centering::CC); + + auto init_np = [&](int kind, Double3 crd, int p, Int3 idx, + psc_particle_np& np) { + double temperature = + np.kind == KIND_ION ? ion_temperature : electron_temperature; + np.n = 1.0; + np.p = + setup_particles.createMaxwellian({np.kind, + np.n, + {0, v_flow, 0}, + {temperature, temperature, temperature}, + np.tag}); + }; + + partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, + init_np); +} + +// ====================================================================== +// fillGhosts + +template +void fillGhosts(MF& mfld, int compBegin, int compEnd) +{ + Bnd_ bnd{}; + bnd.fill_ghosts(mfld, compBegin, compEnd); +} + +// ====================================================================== +// initializeFields + +void initializeFields(MfieldsState& mflds) +{ + setupFields(mflds, [&](int component, double coords[3]) { + switch (component) { + case HX: return b_x; + case HY: return b_y; + default: return 0.0; + } + }); } // ====================================================================== @@ -198,6 +259,7 @@ static void run(int argc, char** argv) // set up initial conditions initializeParticles(balance, grid_ptr, mprts); + initializeFields(mflds); // ---------------------------------------------------------------------- // run the simulation From f89c88bd64874d0e45b2a756b2e48e5246711938 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 16:00:33 -0400 Subject: [PATCH 05/56] test_reflective_bcs: assert err(psc_params, *grid_ptr, mflds, mprts, balance, collision, checks, marder, diagnostics); - psc.integrate(); + for (; grid.timestep_ < nt; 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); + } } // ====================================================================== From d508832f3f8e16485d9f6f396caa4c50dcf09862 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 16:08:36 -0400 Subject: [PATCH 06/56] test_reflective_bcs: cleanup --- src/libpsc/tests/test_reflective_bcs.cxx | 36 ++---------------------- 1 file changed, 3 insertions(+), 33 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 93a6a96a4..b8fa38ffe 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -17,11 +17,7 @@ // EDIT to change order / floating point type / cuda / 2d/3d using Dim = dim_yz; -#ifdef USE_CUDA -using PscConfig = PscConfig1vbecCuda; -#else using PscConfig = PscConfig1vbecDouble; -#endif // ---------------------------------------------------------------------- @@ -39,11 +35,8 @@ using OutputParticles = PscConfig::OutputParticles; namespace { -// General PSC parameters PscParams psc_params; -// matching: heliospheric .1 AU 3, shock tiny 2, real c and M, parallel B - double electron_temperature = 7.8278E-05; double ion_temperature = 7.8278E-05; double ion_mass = 1.8360E+03; @@ -87,15 +80,9 @@ void setupParameters(int argc, char** argv) // ====================================================================== // setupGrid -// -// This helper function is responsible for setting up the "Grid", -// which is really more than just the domain and its decomposition, it -// also encompasses PC normalization parameters, information about the -// particle kinds, etc. Grid_t* setupGrid() { - // FIXME add a check to catch mismatch between Dim and n grid points early auto domain = Grid_t::Domain{{nx, ny, nz}, // n grid points {len_x, len_y, len_z}, // physical lengths @@ -158,16 +145,6 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) init_np); } -// ====================================================================== -// fillGhosts - -template -void fillGhosts(MF& mfld, int compBegin, int compEnd) -{ - Bnd_ bnd{}; - bnd.fill_ghosts(mfld, compBegin, compEnd); -} - // ====================================================================== // initializeFields @@ -201,9 +178,6 @@ static void run(int argc, char** argv) auto& grid = *grid_ptr; MfieldsState mflds{grid}; Mparticles mprts{grid}; - BgkMfields phi{grid, 1, mflds.ibn()}; - BgkMfields gradPhi{grid, 3, mflds.ibn()}; - BgkMfields divGradPhi{grid, 1, mflds.ibn()}; // ---------------------------------------------------------------------- // Set up various objects needed to run this case @@ -275,16 +249,11 @@ static void run(int argc, char** argv) } } -// ====================================================================== -// main - // ====================================================================== // Test1 TEST(ReflectiveBcsTest, Test1) { - // psc_init(argc, argv); - // FIXME restore whatever previous functionality there was with options int argc = 1; char** argv = nullptr; psc_init(argc, argv); @@ -294,11 +263,12 @@ TEST(ReflectiveBcsTest, Test1) psc_finalize(); } +// ====================================================================== +// main + int main(int argc, char** argv) { - // MPI_Init(&argc, &argv); ::testing::InitGoogleTest(&argc, argv); int rc = RUN_ALL_TESTS(); - // MPI_Finalize(); return rc; } From b17f487bcbea687aa76e8da3e078bd1707f1e0c5 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 4 Apr 2025 14:23:27 -0400 Subject: [PATCH 07/56] test_reflective_bcs: simplify while still failing --- src/libpsc/tests/test_reflective_bcs.cxx | 40 +++++------------------- 1 file changed, 7 insertions(+), 33 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index b8fa38ffe..cc91d8e8c 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -37,15 +37,8 @@ namespace { PscParams psc_params; -double electron_temperature = 7.8278E-05; -double ion_temperature = 7.8278E-05; -double ion_mass = 1.8360E+03; -double v_flow = 2.6685E-03; -double b_x = 0; -double b_y = 2.4000E-02; - int nx = 1; -int ny = 128; +int ny = 8; // test success if this goes to 4 int nz = 2; int nt = 100; @@ -97,11 +90,11 @@ Grid_t* setupGrid() auto kinds = Grid_t::Kinds(NR_KINDS); kinds[KIND_ELECTRON] = {-1.0, 1.0, "e"}; - kinds[KIND_ION] = {1.0, ion_mass, "i"}; + kinds[KIND_ION] = {1.0, 1.0, "i"}; // --- generic setup auto norm_params = Grid_t::NormalizationParams::dimensionless(); - norm_params.nicell = 100; + norm_params.nicell = 1; double dt = psc_params.cfl * courant_length(domain); Grid_t::Normalization norm{norm_params}; @@ -130,35 +123,17 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) auto init_np = [&](int kind, Double3 crd, int p, Int3 idx, psc_particle_np& np) { - double temperature = - np.kind == KIND_ION ? ion_temperature : electron_temperature; np.n = 1.0; - np.p = - setup_particles.createMaxwellian({np.kind, - np.n, - {0, v_flow, 0}, - {temperature, temperature, temperature}, - np.tag}); + + // how can i make this fail without setting some t=0? + np.p = setup_particles.createMaxwellian( + {np.kind, np.n, {0, 0, 0}, {0, 0, 0.0000006}, np.tag}); }; partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, init_np); } -// ====================================================================== -// initializeFields - -void initializeFields(MfieldsState& mflds) -{ - setupFields(mflds, [&](int component, double coords[3]) { - switch (component) { - case HX: return b_x; - case HY: return b_y; - default: return 0.0; - } - }); -} - // ====================================================================== // run @@ -233,7 +208,6 @@ static void run(int argc, char** argv) // set up initial conditions initializeParticles(balance, grid_ptr, mprts); - initializeFields(mflds); // ---------------------------------------------------------------------- // run the simulation From 4050d24d750561f5eaed3e33114958f452cce962 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 11 Apr 2025 11:27:44 -0400 Subject: [PATCH 08/56] test_reflective_bcs: use just 1 particle --- src/libpsc/tests/test_reflective_bcs.cxx | 33 ++++++------------------ 1 file changed, 8 insertions(+), 25 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index cc91d8e8c..9fcbcb26a 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -46,9 +46,9 @@ int n_patches_x = 1; int n_patches_y = 1; int n_patches_z = 1; -double dx = 4.2849E+01; -double dy = 3.3475E-01; -double dz = 3.3475E-01; +double dx = 10; +double dy = 10; +double dz = 10; double len_x = nx * dx; double len_y = ny * dy; @@ -113,27 +113,6 @@ Grid_t* setupGrid() return new Grid_t{domain, bc, kinds, norm, dt, -1, ibn}; } -// ====================================================================== -// initializeParticles - -void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) -{ - SetupParticles setup_particles(*grid_ptr); - setup_particles.centerer = Centering::Centerer(Centering::CC); - - auto init_np = [&](int kind, Double3 crd, int p, Int3 idx, - psc_particle_np& np) { - np.n = 1.0; - - // how can i make this fail without setting some t=0? - np.p = setup_particles.createMaxwellian( - {np.kind, np.n, {0, 0, 0}, {0, 0, 0.0000006}, np.tag}); - }; - - partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, - init_np); -} - // ====================================================================== // run @@ -207,7 +186,11 @@ static void run(int argc, char** argv) // ---------------------------------------------------------------------- // set up initial conditions - initializeParticles(balance, grid_ptr, mprts); + { + auto injector = mprts.injector(); + auto inj = injector[0]; + inj({{0, 15, 5}, {0, -.001, 0}, 1, 0}); + } // ---------------------------------------------------------------------- // run the simulation From 0666140234229ea899a4775ccbb77086ba2e2c8c Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 11 Apr 2025 11:29:13 -0400 Subject: [PATCH 09/56] fields_item, +add_ghosts: extract add_ghosts_reflecting_hi --- src/include/add_ghosts_reflecting.hxx | 38 +++++++++++++++++++++++++++ src/include/fields_item.hxx | 35 +----------------------- 2 files changed, 39 insertions(+), 34 deletions(-) create mode 100644 src/include/add_ghosts_reflecting.hxx diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx new file mode 100644 index 000000000..907d94197 --- /dev/null +++ b/src/include/add_ghosts_reflecting.hxx @@ -0,0 +1,38 @@ +#pragma once + +#include "../kg/include/kg/Vec3.h" +#include "grid.hxx" + +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); + } +} \ No newline at end of file diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index cbcbab3ba..96eb36bc3 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 @@ -87,40 +88,6 @@ private: } } - 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) From 23ba92fb6dc8d3bded58664d489e6231d2f69080 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 15 Apr 2025 15:02:24 -0400 Subject: [PATCH 10/56] test_reflective_bcs: don't dump --- src/libpsc/tests/test_reflective_bcs.cxx | 13 ++++--------- 1 file changed, 4 insertions(+), 9 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 9fcbcb26a..5f93115b2 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -54,9 +54,6 @@ double len_x = nx * dx; double len_y = ny * dy; double len_z = nz * dz; -int n_writes = 10; -int out_interval = nt / n_writes; - } // namespace // ====================================================================== @@ -65,7 +62,7 @@ int out_interval = nt / n_writes; void setupParameters(int argc, char** argv) { psc_params.nmax = nt; - psc_params.stats_every = out_interval; + psc_params.stats_every = 1; psc_params.cfl = .75; psc_params.write_checkpoint_every_step = 0; @@ -167,15 +164,13 @@ static void run(int argc, char** argv) // -- output fields OutputFieldsParams outf_params{}; - outf_params.fields.pfield.out_interval = out_interval; - outf_params.moments.pfield.out_interval = out_interval; + 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 = out_interval; - outp_params.data_dir = "."; - outp_params.basename = "prt"; + outp_params.every_step = -1; OutputParticles outp{grid, outp_params}; int oute_interval = -100; From 780a39f2af51cd411798c86a0408fe20b0bc7f49 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 15 Apr 2025 16:10:25 -0400 Subject: [PATCH 11/56] checks_impl: -MP template from ChecksCommon --- src/libpsc/psc_checks/checks_impl.hxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) 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>; From 9b5e5898e30b473fcbfe1214a688552ebc399812 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 12:31:12 -0400 Subject: [PATCH 12/56] test_reflective_bcs: rename test1 --- src/libpsc/tests/test_reflective_bcs.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 5f93115b2..0a04f2546 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -202,9 +202,9 @@ static void run(int argc, char** argv) } // ====================================================================== -// Test1 +// Integration test (pun not intended) -TEST(ReflectiveBcsTest, Test1) +TEST(ReflectiveBcsTest, Integration) { int argc = 1; char** argv = nullptr; From 4ee752f8100a43c4b077d5590bc6419e5a5921e5 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 12:35:41 -0400 Subject: [PATCH 13/56] test_reflective_bcs: add unit tests that fail --- src/libpsc/tests/test_reflective_bcs.cxx | 27 ++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 0a04f2546..8beaeb9f7 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -215,6 +215,33 @@ TEST(ReflectiveBcsTest, Integration) psc_finalize(); } +// ====================================================================== +// Unit tests + +TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) +{ + // TODO + EXPECT_TRUE(false); +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) +{ + // TODO + EXPECT_TRUE(false); +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) +{ + // TODO + EXPECT_TRUE(false); +} + +TEST(ReflectiveBcsTest, AddGhostsReflectingLowZ) +{ + // TODO + EXPECT_TRUE(false); +} + // ====================================================================== // main From c7a6874249149690a7738318664a9725b03847b6 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 13:44:55 -0400 Subject: [PATCH 14/56] fields_item, add_ghosts_reflecting: mv lo --- src/include/add_ghosts_reflecting.hxx | 34 +++++++++++++++++++++++++++ src/include/fields_item.hxx | 34 --------------------------- 2 files changed, 34 insertions(+), 34 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 907d94197..47cc6de4d 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -3,6 +3,40 @@ #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) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 96eb36bc3..854ce4c7a 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -54,40 +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_boundary(const Grid_t& grid, FE& mres_gt, const Int3& ib, int p, int mb, int me) From d2f8123eec6f28e2e2e508b58f0c46563bf52e3f Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 13:55:47 -0400 Subject: [PATCH 15/56] test_reflective_bcs: some magic numbers --- src/libpsc/tests/test_reflective_bcs.cxx | 31 +++++------------------- 1 file changed, 6 insertions(+), 25 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 8beaeb9f7..2f6c09e90 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -36,24 +36,6 @@ using OutputParticles = PscConfig::OutputParticles; namespace { PscParams psc_params; - -int nx = 1; -int ny = 8; // test success if this goes to 4 -int nz = 2; -int nt = 100; - -int n_patches_x = 1; -int n_patches_y = 1; -int n_patches_z = 1; - -double dx = 10; -double dy = 10; -double dz = 10; - -double len_x = nx * dx; -double len_y = ny * dy; -double len_z = nz * dz; - } // namespace // ====================================================================== @@ -61,7 +43,7 @@ double len_z = nz * dz; void setupParameters(int argc, char** argv) { - psc_params.nmax = nt; + psc_params.nmax = 100; psc_params.stats_every = 1; psc_params.cfl = .75; @@ -73,11 +55,10 @@ void setupParameters(int argc, char** argv) Grid_t* setupGrid() { - auto domain = - Grid_t::Domain{{nx, ny, nz}, // n grid points - {len_x, len_y, len_z}, // physical lengths - {0, 0, 0}, // location of lower corner - {n_patches_x, n_patches_y, n_patches_z}}; // n patches + 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}, @@ -194,7 +175,7 @@ static void run(int argc, char** argv) makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, collision, checks, marder, diagnostics); - for (; grid.timestep_ < nt; grid.timestep_++) { + 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); From 577cb56beba93cced0697bf68a4e5965ced3efff Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 14:00:47 -0400 Subject: [PATCH 16/56] test_reflective_bcs: inline setup_params --- src/libpsc/tests/test_reflective_bcs.cxx | 16 +++------------- 1 file changed, 3 insertions(+), 13 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 2f6c09e90..ae128f306 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -38,18 +38,6 @@ namespace PscParams psc_params; } // namespace -// ====================================================================== -// setupParameters - -void setupParameters(int argc, char** argv) -{ - psc_params.nmax = 100; - psc_params.stats_every = 1; - psc_params.cfl = .75; - - psc_params.write_checkpoint_every_step = 0; -} - // ====================================================================== // setupGrid @@ -101,7 +89,9 @@ static void run(int argc, char** argv) // ---------------------------------------------------------------------- // setup various parameters first - setupParameters(argc, argv); + psc_params.nmax = 100; + psc_params.stats_every = 1; + psc_params.cfl = .75; // ---------------------------------------------------------------------- // Set up grid, state fields, particles From ad0c34a812ee570b64b9319059138438625c7c66 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 15:48:36 -0400 Subject: [PATCH 17/56] test_reflective_bcs: mv mpi init --- src/libpsc/tests/test_reflective_bcs.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index ae128f306..d5cedc24a 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -179,11 +179,8 @@ TEST(ReflectiveBcsTest, Integration) { int argc = 1; char** argv = nullptr; - psc_init(argc, argv); run(argc, argv); - - psc_finalize(); } // ====================================================================== @@ -218,7 +215,9 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowZ) int main(int argc, char** argv) { + MPI_Init(&argc, &argv); ::testing::InitGoogleTest(&argc, argv); int rc = RUN_ALL_TESTS(); + MPI_Finalize(); return rc; } From 5bd37a0800feaa83495f9e9f0b538ad354c25c5d Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:31:04 -0400 Subject: [PATCH 18/56] test_reflective_bcs: impl AddGhostsReflectingHighY --- src/libpsc/tests/test_reflective_bcs.cxx | 48 +++++++++++++++++++++++- 1 file changed, 46 insertions(+), 2 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index d5cedc24a..2c6809091 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -186,10 +186,54 @@ TEST(ReflectiveBcsTest, Integration) // ====================================================================== // Unit tests +int ravel_idx(Int3 idx, Int3 shape) +{ + return idx[0] + shape[0] * (idx[1] + shape[1] * (idx[2])); +} + TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) { - // TODO - EXPECT_TRUE(false); + 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) From 2d2c3aef6964c092d7860243c835fac5accf4a3e Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:36:11 -0400 Subject: [PATCH 19/56] add_ghosts_reflecting: fix sign errors in y idx --- src/include/add_ghosts_reflecting.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 47cc6de4d..1700f3c23 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -49,7 +49,7 @@ void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, { 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); + mres_gt(ix - ib[0], iy + 1 - ib[1], iz - ib[2], m, p); } } } From 001e20c98258afed12fff5a7e42a65c69af2375d Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:36:44 -0400 Subject: [PATCH 20/56] add_ghosts_reflecting: scan entire ghost region --- src/include/add_ghosts_reflecting.hxx | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 1700f3c23..34a1010de 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -43,13 +43,15 @@ void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, { 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; - { + // FIXME only need to scan 1 cell into the ghost region for 1st-order + // moments + for (int iz = ib[2]; iz < ldims[2] - ib[2]; iz++) { + for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { + for (int iy = ldims[1] + ib[1]; iy < ldims[1]; iy++) { for (int m = mb; m < me; m++) { + int iy_reflected = 2 * ldims[1] - iy - 1; 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); + mres_gt(ix - ib[0], iy_reflected - ib[1], iz - ib[2], m, p); } } } From 1b108a59fe1e54a73a802b69498989afc2c33d6f Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:41:40 -0400 Subject: [PATCH 21/56] test_reflective_bcs: impl AddGhostsReflectingHighZ --- src/libpsc/tests/test_reflective_bcs.cxx | 44 ++++++++++++++++++++++-- 1 file changed, 42 insertions(+), 2 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 2c6809091..5f839c6d1 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -6,6 +6,7 @@ #include "DiagnosticsDefault.h" #include "OutputFieldsDefault.h" +#include "add_ghosts_reflecting.hxx" #include "../psc_config.hxx" // ====================================================================== @@ -244,8 +245,47 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) { - // TODO - EXPECT_TRUE(false); + 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) From 4bfe629f867169d4d409f392f0599d193df0bfd0 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:44:55 -0400 Subject: [PATCH 22/56] add_ghosts_reflecting: also fix z --- src/include/add_ghosts_reflecting.hxx | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 34a1010de..265aa7090 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -41,10 +41,8 @@ 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; + // FIXME only need to scan 1 cell into the ghost region for 1st-order moments if (d == 1) { - // FIXME only need to scan 1 cell into the ghost region for 1st-order - // moments for (int iz = ib[2]; iz < ldims[2] - ib[2]; iz++) { for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { for (int iy = ldims[1] + ib[1]; iy < ldims[1]; iy++) { @@ -57,13 +55,13 @@ void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, } } } 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 iz = ldims[2] + ib[2]; iz < ldims[2]; iz++) { + for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { + for (int iy = ib[1]; iy < ldims[1] - ib[1]; iy++) { for (int m = mb; m < me; m++) { + int iz_reflected = 2 * ldims[2] - iz - 1; 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); + mres_gt(ix - ib[0], iy - ib[1], iz_reflected - ib[2], m, p); } } } From 05ac0109127299fc0cb93810631336046dd12e3e Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:52:25 -0400 Subject: [PATCH 23/56] add_ghosts_reflecting: reorder iter --- src/include/add_ghosts_reflecting.hxx | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 265aa7090..73a35484a 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -42,11 +42,12 @@ 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 + if (d == 1) { - for (int iz = ib[2]; iz < ldims[2] - ib[2]; iz++) { - for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { + for (int m = mb; m < me; m++) { + for (int iz = ib[2]; iz < ldims[2] - ib[2]; iz++) { for (int iy = ldims[1] + ib[1]; iy < ldims[1]; iy++) { - for (int m = mb; m < me; m++) { + for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { int iy_reflected = 2 * ldims[1] - iy - 1; mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += mres_gt(ix - ib[0], iy_reflected - ib[1], iz - ib[2], m, p); @@ -55,10 +56,10 @@ void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, } } } else if (d == 2) { - for (int iz = ldims[2] + ib[2]; iz < ldims[2]; iz++) { - for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { + for (int m = mb; m < me; m++) { + for (int iz = ldims[2] + ib[2]; iz < ldims[2]; iz++) { for (int iy = ib[1]; iy < ldims[1] - ib[1]; iy++) { - for (int m = mb; m < me; m++) { + for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { int iz_reflected = 2 * ldims[2] - iz - 1; mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += mres_gt(ix - ib[0], iy - ib[1], iz_reflected - ib[2], m, p); From e035122fb743ae4dc284a9877348d9c20b745e8d Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 17:13:33 -0400 Subject: [PATCH 24/56] add_ghosts_reflecting: generalize revert this if it's horribly inefficient --- src/include/add_ghosts_reflecting.hxx | 44 +++++++++++++-------------- 1 file changed, 21 insertions(+), 23 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 73a35484a..c318ba33b 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -43,31 +43,29 @@ void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, { // FIXME only need to scan 1 cell into the ghost region for 1st-order moments - if (d == 1) { - for (int m = mb; m < me; m++) { - for (int iz = ib[2]; iz < ldims[2] - ib[2]; iz++) { - for (int iy = ldims[1] + ib[1]; iy < ldims[1]; iy++) { - for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { - int iy_reflected = 2 * ldims[1] - iy - 1; - mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += - mres_gt(ix - ib[0], iy_reflected - ib[1], iz - ib[2], m, p); - } - } - } - } - } else if (d == 2) { - for (int m = mb; m < me; m++) { - for (int iz = ldims[2] + ib[2]; iz < ldims[2]; iz++) { - for (int iy = ib[1]; iy < ldims[1] - ib[1]; iy++) { - for (int ix = ib[0]; ix < ldims[0] - ib[0]; ix++) { - int iz_reflected = 2 * ldims[2] - iz - 1; - mres_gt(ix - ib[0], iy - ib[1], iz - ib[2], m, p) += - mres_gt(ix - ib[0], iy - ib[1], iz_reflected - ib[2], m, p); - } + 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; } } } - } else { - assert(0); } } \ No newline at end of file From ad03f7deabf20cdb39cebdb0b6d7325eb4bca1c9 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 18 Apr 2025 10:37:19 -0400 Subject: [PATCH 25/56] test_reflective_bcs: +init_mres --- src/libpsc/tests/test_reflective_bcs.cxx | 31 ++++++++++++------------ 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 5f839c6d1..899c0c721 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -192,6 +192,19 @@ 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, AddGhostsReflectingHighY) { Grid_t* grid_ptr = setupGrid(); @@ -206,14 +219,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) 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; - } - } - } + init_mres(mres, shape, p); int dim = 1; add_ghosts_reflecting_hi(ldims, mres, ib, p, dim, 0, 1); @@ -257,14 +263,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) 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; - } - } - } + init_mres(mres, shape, p); int dim = 2; add_ghosts_reflecting_hi(ldims, mres, ib, p, dim, 0, 1); From 6fb51052a6e1986d946aa879d0475eb9bc696b14 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 18 Apr 2025 10:45:22 -0400 Subject: [PATCH 26/56] test_reflective_bcs: impl lo tests --- src/libpsc/tests/test_reflective_bcs.cxx | 72 ++++++++++++++++++++++-- 1 file changed, 68 insertions(+), 4 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 899c0c721..2c6d3bf5a 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -245,8 +245,40 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) { - // TODO - EXPECT_TRUE(false); + 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(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, AddGhostsReflectingHighZ) @@ -289,8 +321,40 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) TEST(ReflectiveBcsTest, AddGhostsReflectingLowZ) { - // TODO - EXPECT_TRUE(false); + 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(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; + } + } + } } // ====================================================================== From 2532a677723c8647a416b955f4164d0e007e6506 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 18 Apr 2025 10:46:02 -0400 Subject: [PATCH 27/56] add_ghosts_reflecting: rewrite lo --- src/include/add_ghosts_reflecting.hxx | 47 +++++++++++++-------------- 1 file changed, 23 insertions(+), 24 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index c318ba33b..7b77b9f53 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -7,33 +7,32 @@ 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); - } + // 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] = 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; } } } - } else { - assert(0); } } From 3ea024c58acbc8581a8c7deadb92f8907c4185ae Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 18 Apr 2025 11:08:57 -0400 Subject: [PATCH 28/56] test_reflective_bcs: ensure particle reflects --- src/libpsc/tests/test_reflective_bcs.cxx | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 2c6d3bf5a..14083730e 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -150,12 +150,15 @@ static void run(int argc, char** argv) auto diagnostics = makeDiagnosticsDefault(outf, outp, oute); + EXPECT_EQ(grid.n_patches(), 1); + int p = 0; + // ---------------------------------------------------------------------- // set up initial conditions { auto injector = mprts.injector(); - auto inj = injector[0]; + auto inj = injector[p]; inj({{0, 15, 5}, {0, -.001, 0}, 1, 0}); } @@ -166,11 +169,19 @@ static void run(int argc, char** argv) makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, collision, checks, marder, diagnostics); + auto prts = mprts.accessor()[p]; + + EXPECT_EQ(prts.size(), 1); + EXPECT_LT(prts[0].u()[1], 0.0); + 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); } + + EXPECT_EQ(prts.size(), 1); + EXPECT_GT(prts[0].u()[1], 0.0); } // ====================================================================== From ce9e67ba08d0b124934a1fd0e260923499e23ddb Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 18 Apr 2025 11:09:40 -0400 Subject: [PATCH 29/56] test_reflective_bcs: inline run() --- src/libpsc/tests/test_reflective_bcs.cxx | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 14083730e..31b6da7d5 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -81,9 +81,9 @@ Grid_t* setupGrid() } // ====================================================================== -// run +// Integration test (pun not intended) -static void run(int argc, char** argv) +TEST(ReflectiveBcsTest, Integration) { mpi_printf(MPI_COMM_WORLD, "*** Setting up...\n"); @@ -184,17 +184,6 @@ static void run(int argc, char** argv) EXPECT_GT(prts[0].u()[1], 0.0); } -// ====================================================================== -// Integration test (pun not intended) - -TEST(ReflectiveBcsTest, Integration) -{ - int argc = 1; - char** argv = nullptr; - - run(argc, argv); -} - // ====================================================================== // Unit tests From 7fe83afee995b2edca60a8a801aa678f68cea631 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 18 Apr 2025 14:10:07 -0400 Subject: [PATCH 30/56] test_reflective_bcs: condense setup --- src/libpsc/tests/test_reflective_bcs.cxx | 67 +++++------------------- 1 file changed, 14 insertions(+), 53 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 31b6da7d5..f4caa26b4 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -85,77 +85,42 @@ Grid_t* setupGrid() TEST(ReflectiveBcsTest, Integration) { - mpi_printf(MPI_COMM_WORLD, "*** Setting up...\n"); - // ---------------------------------------------------------------------- - // setup various parameters first + // setup 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_params.gauss.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}; + 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); - EXPECT_EQ(grid.n_patches(), 1); - int p = 0; + 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]; @@ -165,10 +130,6 @@ TEST(ReflectiveBcsTest, Integration) // ---------------------------------------------------------------------- // run the simulation - auto psc = - makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, - collision, checks, marder, diagnostics); - auto prts = mprts.accessor()[p]; EXPECT_EQ(prts.size(), 1); From 1c9228ca0d8553df7f5067b4e141adce8680d502 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 29 Apr 2025 11:13:17 -0400 Subject: [PATCH 31/56] test_reflective_bcs: inj 2 prts --- src/libpsc/tests/test_reflective_bcs.cxx | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index f4caa26b4..8b4406556 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -124,7 +124,11 @@ TEST(ReflectiveBcsTest, Integration) { auto injector = mprts.injector(); auto inj = injector[p]; - inj({{0, 15, 5}, {0, -.001, 0}, 1, 0}); + double vy = .1; + 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}, 0, 0}); // electron + inj({{0, y_center, 5}, {0, vy, 0}, 1, 0}); // positron } // ---------------------------------------------------------------------- @@ -132,8 +136,9 @@ TEST(ReflectiveBcsTest, Integration) auto prts = mprts.accessor()[p]; - EXPECT_EQ(prts.size(), 1); + EXPECT_EQ(prts.size(), 2); EXPECT_LT(prts[0].u()[1], 0.0); + EXPECT_GT(prts[1].u()[1], 0.0); for (; grid.timestep_ < psc_params.nmax; grid.timestep_++) { psc.step(); @@ -141,8 +146,9 @@ TEST(ReflectiveBcsTest, Integration) ASSERT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); } - EXPECT_EQ(prts.size(), 1); + EXPECT_EQ(prts.size(), 2); EXPECT_GT(prts[0].u()[1], 0.0); + EXPECT_LT(prts[1].u()[1], 0.0); } // ====================================================================== From deaa8ddea79e7647fe3938d97c2f6305cc6b9a14 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 29 Apr 2025 11:23:29 -0400 Subject: [PATCH 32/56] test_reflective_bcs: switch reverse/expects --- src/libpsc/tests/test_reflective_bcs.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 8b4406556..f0966a38e 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -136,19 +136,19 @@ TEST(ReflectiveBcsTest, Integration) auto prts = mprts.accessor()[p]; - EXPECT_EQ(prts.size(), 2); - EXPECT_LT(prts[0].u()[1], 0.0); - EXPECT_GT(prts[1].u()[1], 0.0); + ASSERT_EQ(prts.size(), 2); + ASSERT_LT(prts[0].u()[1], 0.0); + ASSERT_GT(prts[1].u()[1], 0.0); 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); + EXPECT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); + EXPECT_LT(checks.gauss.last_max_err, checks.gauss.err_threshold); } - EXPECT_EQ(prts.size(), 2); - EXPECT_GT(prts[0].u()[1], 0.0); - EXPECT_LT(prts[1].u()[1], 0.0); + ASSERT_EQ(prts.size(), 2); + ASSERT_GT(prts[0].u()[1], 0.0); + ASSERT_LT(prts[1].u()[1], 0.0); } // ====================================================================== From 09c14d3cea8f35379192c08283a56ed22f3ec127 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 29 Apr 2025 11:27:55 -0400 Subject: [PATCH 33/56] test_reflective_bcs: more asserts --- src/libpsc/tests/test_reflective_bcs.cxx | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index f0966a38e..97d3a6295 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -140,6 +140,12 @@ TEST(ReflectiveBcsTest, Integration) 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; grid.timestep_++) { psc.step(); EXPECT_LT(checks.continuity.last_max_err, checks.continuity.err_threshold); From 86bfabe035a83a8b961b35bb31f1f6e924e94451 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 29 Apr 2025 11:29:36 -0400 Subject: [PATCH 34/56] test_reflective_bcs: fix weight and kind --- src/libpsc/tests/test_reflective_bcs.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 97d3a6295..153735ee9 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -127,8 +127,8 @@ TEST(ReflectiveBcsTest, Integration) double vy = .1; 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}, 0, 0}); // electron - inj({{0, y_center, 5}, {0, vy, 0}, 1, 0}); // positron + inj({{0, y_center, 5}, {0, -vy, 0}, 1, 0}); // electron + inj({{0, y_center, 5}, {0, vy, 0}, 1, 1}); // positron } // ---------------------------------------------------------------------- From 73cbb1ff138bf72fd68782cb1a04a416659c652a Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 29 Apr 2025 11:35:06 -0400 Subject: [PATCH 35/56] test_reflective_bcs: go faster --- src/libpsc/tests/test_reflective_bcs.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 153735ee9..9ca38cc7b 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -124,7 +124,7 @@ TEST(ReflectiveBcsTest, Integration) { auto injector = mprts.injector(); auto inj = injector[p]; - double vy = .1; + 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 From 36f7064fc9b8c407aec2784e59f292a508974285 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 1 May 2025 10:45:54 -0400 Subject: [PATCH 36/56] add_ghosts_reflecting; *: rename to _cc --- src/include/add_ghosts_reflecting.hxx | 8 ++++---- src/include/fields_item.hxx | 4 ++-- src/libpsc/tests/test_reflective_bcs.cxx | 16 ++++++++-------- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 7b77b9f53..db64f2019 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -4,8 +4,8 @@ #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) +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 @@ -37,8 +37,8 @@ void add_ghosts_reflecting_lo(const Int3& ldims, FE& mres_gt, const Int3& ib, } template -void add_ghosts_reflecting_hi(const Int3& ldims, FE& mres_gt, const Int3& ib, - int p, int d, int mb, int me) +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 diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 854ce4c7a..cf8783787 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -63,7 +63,7 @@ private: 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); + add_ghosts_reflecting_lo_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } } } @@ -72,7 +72,7 @@ private: 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); + add_ghosts_reflecting_hi_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } } } diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 9ca38cc7b..48d1c5f1a 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -178,7 +178,7 @@ void init_mres(FE& mres_gt, Int3 shape, int p) } } -TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) +TEST(ReflectiveBcsTest, AddGhostsReflectingHighCcY) { Grid_t* grid_ptr = setupGrid(); auto& grid = *grid_ptr; @@ -195,7 +195,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) init_mres(mres, shape, p); int dim = 1; - add_ghosts_reflecting_hi(ldims, mres, ib, p, dim, 0, 1); + add_ghosts_reflecting_hi_cc(ldims, mres, ib, p, dim, 0, 1); int bx = -ib[0]; int by = -ib[1]; @@ -216,7 +216,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighY) } } -TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) +TEST(ReflectiveBcsTest, AddGhostsReflectingLowCcY) { Grid_t* grid_ptr = setupGrid(); auto& grid = *grid_ptr; @@ -233,7 +233,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) init_mres(mres, shape, p); int dim = 1; - add_ghosts_reflecting_lo(ldims, mres, ib, p, dim, 0, 1); + add_ghosts_reflecting_lo_cc(ldims, mres, ib, p, dim, 0, 1); int bx = -ib[0]; int by = -ib[1]; @@ -254,7 +254,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowY) } } -TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) +TEST(ReflectiveBcsTest, AddGhostsReflectingHighCcZ) { Grid_t* grid_ptr = setupGrid(); auto& grid = *grid_ptr; @@ -271,7 +271,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) init_mres(mres, shape, p); int dim = 2; - add_ghosts_reflecting_hi(ldims, mres, ib, p, dim, 0, 1); + add_ghosts_reflecting_hi_cc(ldims, mres, ib, p, dim, 0, 1); int bx = -ib[0]; int by = -ib[1]; @@ -292,7 +292,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighZ) } } -TEST(ReflectiveBcsTest, AddGhostsReflectingLowZ) +TEST(ReflectiveBcsTest, AddGhostsReflectingLowCcZ) { Grid_t* grid_ptr = setupGrid(); auto& grid = *grid_ptr; @@ -309,7 +309,7 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowZ) init_mres(mres, shape, p); int dim = 2; - add_ghosts_reflecting_lo(ldims, mres, ib, p, dim, 0, 1); + add_ghosts_reflecting_lo_cc(ldims, mres, ib, p, dim, 0, 1); int bx = -ib[0]; int by = -ib[1]; From 1e814c50d0bdc0a59ae045f6d4f2b7dece8ba7e9 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 1 May 2025 10:53:39 -0400 Subject: [PATCH 37/56] add_ghosts_reflecting: +nc stubs --- src/include/add_ghosts_reflecting.hxx | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index db64f2019..4cfa0a1da 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -67,4 +67,14 @@ void add_ghosts_reflecting_hi_cc(const Int3& ldims, FE& mres_gt, const Int3& ib, } } } -} \ No newline at end of file +} + +template +void add_ghosts_reflecting_lo_nc(const Int3& ldims, FE& mres_gt, const Int3& ib, + int p, int d, int mb, int me) +{} + +template +void add_ghosts_reflecting_hi_nc(const Int3& ldims, FE& mres_gt, const Int3& ib, + int p, int d, int mb, int me) +{} From 23af57f67b867379c12f47fca614d9888c28fc51 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 1 May 2025 13:17:38 -0400 Subject: [PATCH 38/56] test_reflective_bcs: +HighNcY test --- src/libpsc/tests/test_reflective_bcs.cxx | 41 ++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 48d1c5f1a..10b1c0a01 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -330,6 +330,47 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowCcZ) } } +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; + } + } + } +} + // ====================================================================== // main From c5c06388b8c1b9c47cad8855b0ec92d4b4c368d8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 11:26:41 -0400 Subject: [PATCH 39/56] add_ghosts_reflecting: fix order --- src/include/add_ghosts_reflecting.hxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 4cfa0a1da..f9f565912 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -9,8 +9,8 @@ void add_ghosts_reflecting_lo_cc(const Int3& ldims, FE& mres_gt, const Int3& ib, { // FIXME only need to scan 1 cell into the ghost region for 1st-order moments - Int3 idx_end = ldims - ib; Int3 idx_begin = ib; + Int3 idx_end = ldims - ib; idx_begin[d] = 0; idx_end[d] = -ib[d]; @@ -42,8 +42,8 @@ void add_ghosts_reflecting_hi_cc(const Int3& ldims, FE& mres_gt, const Int3& ib, { // FIXME only need to scan 1 cell into the ghost region for 1st-order moments - Int3 idx_end = ldims - ib; Int3 idx_begin = ib; + Int3 idx_end = ldims - ib; idx_begin[d] = ldims[d] + ib[d]; idx_end[d] = ldims[d]; From f30cbc427a74af01005d6a23a7226ac88e22ec97 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 11:31:22 -0400 Subject: [PATCH 40/56] add_ghosts_reflecting: impl hi_nc --- src/include/add_ghosts_reflecting.hxx | 37 ++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index f9f565912..557aeda52 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -77,4 +77,39 @@ void add_ghosts_reflecting_lo_nc(const Int3& ldims, FE& mres_gt, const Int3& ib, 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; + } + } + } + } +} From cd039ed35a494bc2929b7ac412c1c0889977bea1 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 11:35:22 -0400 Subject: [PATCH 41/56] test_reflective_bcs: +LowNcY test --- src/libpsc/tests/test_reflective_bcs.cxx | 41 ++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index 10b1c0a01..fd6483357 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -371,6 +371,47 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingHighNcY) } } +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; + } + } + } +} + // ====================================================================== // main From 20f787416d6abd5f50c4296692a5ce6c3d06e1ab Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 11:41:38 -0400 Subject: [PATCH 42/56] add_ghosts_reflecting: impl lo_nc --- src/include/add_ghosts_reflecting.hxx | 37 ++++++++++++++++++++++++++- 1 file changed, 36 insertions(+), 1 deletion(-) diff --git a/src/include/add_ghosts_reflecting.hxx b/src/include/add_ghosts_reflecting.hxx index 557aeda52..1509dfcbb 100644 --- a/src/include/add_ghosts_reflecting.hxx +++ b/src/include/add_ghosts_reflecting.hxx @@ -72,7 +72,42 @@ void add_ghosts_reflecting_hi_cc(const Int3& ldims, FE& mres_gt, const Int3& ib, 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, From 113c81dab3a81d79a1e2f4aa63c20b10cf60082c Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 12:04:12 -0400 Subject: [PATCH 43/56] test_reflective_bcs: +NcZ tests --- src/libpsc/tests/test_reflective_bcs.cxx | 82 ++++++++++++++++++++++++ 1 file changed, 82 insertions(+) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index fd6483357..ca27cd17b 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -412,6 +412,88 @@ TEST(ReflectiveBcsTest, AddGhostsReflectingLowNcY) } } +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 From be28d6c7f94bf125e6f6c4ef8cc1b740a6665cb5 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 12:04:26 -0400 Subject: [PATCH 44/56] test_reflective_bcs: assert reflection happened --- src/libpsc/tests/test_reflective_bcs.cxx | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/src/libpsc/tests/test_reflective_bcs.cxx b/src/libpsc/tests/test_reflective_bcs.cxx index ca27cd17b..7b570f3d2 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -146,12 +146,25 @@ TEST(ReflectiveBcsTest, Integration) 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); From 644cc12ce7fedeff325ad02df701ae01df5ba9b4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 13:38:03 -0400 Subject: [PATCH 45/56] fields_item: de-nest ifs --- src/include/fields_item.hxx | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index cf8783787..df65750ce 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -60,20 +60,18 @@ private: { // 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_cc(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_cc(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_cc(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_cc(grid.ldims, mres_gt, ib, p, d, mb, me); } } } From 183778807e686a25e892c1804bd3f78d9797ea11 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 13:49:42 -0400 Subject: [PATCH 46/56] fields_item:+DomainBoundary --- src/include/fields_item.hxx | 46 ++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 21 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index df65750ce..35ee9b2dd 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -33,30 +33,12 @@ inline std::vector addKindSuffix( // ====================================================================== // ItemMomentBnd -template -class ItemMomentBnd +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_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++) { @@ -75,6 +57,28 @@ private: } } } +}; + +// ====================================================================== +// 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_; From 43eb8d4278df9bf376895f3e18abf5d6e7ecc148 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:17:01 -0400 Subject: [PATCH 47/56] fields_item: template by Centering --- src/include/fields_item.hxx | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 35ee9b2dd..e70698e97 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -3,6 +3,7 @@ #include "../libpsc/psc_bnd/psc_bnd_impl.hxx" #include "bnd.hxx" +#include "centering.hxx" #include "fields.hxx" #include "fields3d.hxx" #include "particles.hxx" @@ -33,8 +34,18 @@ inline std::vector addKindSuffix( // ====================================================================== // ItemMomentBnd +template class DomainBoundary { +public: + template + 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 static void add_ghosts_boundary(const Grid_t& grid, FE& mres_gt, @@ -73,8 +84,8 @@ public: 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)); + 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)); From c8fc874fc3b6ce659448c38937620fc36d6810f0 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:19:02 -0400 Subject: [PATCH 48/56] fields_item: impl NC specialization --- src/include/fields_item.hxx | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index e70698e97..effb09ab0 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -70,6 +70,33 @@ public: } }; +template <> +class DomainBoundary +{ +public: + template + 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++) { + // 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++) { + // 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 From d877429d01c616dfb071ae97ab1fda5a64e17fa8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:22:31 -0400 Subject: [PATCH 49/56] fields_item: template ItemMomentBnd by centering --- src/include/fields_item.hxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index effb09ab0..0402392fa 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -100,7 +100,7 @@ public: // ====================================================================== // ItemMomentBnd -template +template class ItemMomentBnd { public: @@ -111,8 +111,8 @@ public: 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)); + 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)); @@ -162,6 +162,6 @@ public: } private: - ItemMomentBnd bnd_; + ItemMomentBnd bnd_; std::vector comp_names_; }; From 4a3f307f12d493309f303e5a8ced0899599ab67d Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:33:13 -0400 Subject: [PATCH 50/56] deposit: +DepositCentering to add centering info --- src/include/psc/deposit.hxx | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/src/include/psc/deposit.hxx b/src/include/psc/deposit.hxx index cd386a099..a394d0c78 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"; } From edabbdc6fbda6e46fce291e51524c3a0082b2f32 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:38:13 -0400 Subject: [PATCH 51/56] deposit: using DepositNorm --- src/include/psc/deposit.hxx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/include/psc/deposit.hxx b/src/include/psc/deposit.hxx index a394d0c78..dabe9da5f 100644 --- a/src/include/psc/deposit.hxx +++ b/src/include/psc/deposit.hxx @@ -277,8 +277,9 @@ public: using real_t = R; using dim_t = D; using real3_t = gt::sarray; + using DepositNorm = DEPOSIT_NORM; - 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} {} @@ -290,7 +291,7 @@ public: real3_t x = xi * dxi_; real_t value = fnqs_ * val; - DepositNorm deposit; + DepositNorm deposit; deposit(flds, ib, x, value); } From 9fbd31ed30cbc5dd2e132725a13b61a78ce0dfb3 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:40:49 -0400 Subject: [PATCH 52/56] deposit: expose CENTERING in Deposit --- src/include/psc/deposit.hxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/include/psc/deposit.hxx b/src/include/psc/deposit.hxx index dabe9da5f..af0a89f3f 100644 --- a/src/include/psc/deposit.hxx +++ b/src/include/psc/deposit.hxx @@ -270,7 +270,7 @@ namespace code // deposition in code units // units template class DepositNorm> + template class DEPOSIT_NORM> class Deposit { public: @@ -278,6 +278,7 @@ public: 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(); } From 40b340b3eed8dfe6c42496c4d0520770b2d0a4bf Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:48:20 -0400 Subject: [PATCH 53/56] moment: expose CENTERING --- src/include/psc/moment.hxx | 14 ++++++++++++++ 1 file changed, 14 insertions(+) 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) From dd1d74d76641be053c4b762d632aea0e073f7a20 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 14:49:57 -0400 Subject: [PATCH 54/56] fields_item: get centering from moment_type --- src/include/fields_item.hxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/include/fields_item.hxx b/src/include/fields_item.hxx index 0402392fa..2a4cbe857 100644 --- a/src/include/fields_item.hxx +++ b/src/include/fields_item.hxx @@ -162,6 +162,6 @@ public: } private: - ItemMomentBnd bnd_; + ItemMomentBnd bnd_; std::vector comp_names_; }; From b4fc341e437acf0cb5531b7e44133b98664a8c9d Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 2 May 2025 16:42:17 -0400 Subject: [PATCH 55/56] +test_reflective_bcs_integration: try using psc_init --- src/libpsc/tests/CMakeLists.txt | 11 +- src/libpsc/tests/test_reflective_bcs.cxx | 90 --------- .../tests/test_reflective_bcs_integration.cxx | 183 ++++++++++++++++++ 3 files changed, 191 insertions(+), 93 deletions(-) create mode 100644 src/libpsc/tests/test_reflective_bcs_integration.cxx diff --git a/src/libpsc/tests/CMakeLists.txt b/src/libpsc/tests/CMakeLists.txt index 35f990dc5..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}) @@ -22,6 +24,7 @@ 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) @@ -33,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 index 7b570f3d2..2a321e066 100644 --- a/src/libpsc/tests/test_reflective_bcs.cxx +++ b/src/libpsc/tests/test_reflective_bcs.cxx @@ -80,96 +80,6 @@ Grid_t* setupGrid() 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 prts = mprts.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); -} - // ====================================================================== // Unit tests 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..4b3226b73 --- /dev/null +++ b/src/libpsc/tests/test_reflective_bcs_integration.cxx @@ -0,0 +1,183 @@ +#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 prts = mprts.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; +} From ad22f61d834270be3ef1a65719a8d744fe5cf2e0 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 5 May 2025 14:49:17 -0400 Subject: [PATCH 56/56] test_reflective_bcs_integration: persist accessor --- src/libpsc/tests/test_reflective_bcs_integration.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/libpsc/tests/test_reflective_bcs_integration.cxx b/src/libpsc/tests/test_reflective_bcs_integration.cxx index 4b3226b73..1dd70e6d2 100644 --- a/src/libpsc/tests/test_reflective_bcs_integration.cxx +++ b/src/libpsc/tests/test_reflective_bcs_integration.cxx @@ -134,7 +134,8 @@ TEST(ReflectiveBcsTest, Integration) // ---------------------------------------------------------------------- // run the simulation - auto prts = mprts.accessor()[p]; + auto accessor = mprts.accessor(); + auto prts = accessor[p]; ASSERT_EQ(prts.size(), 2); ASSERT_LT(prts[0].u()[1], 0.0);