From c254885439444526e053d2058df04f3844c1a251 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 31 Mar 2025 13:08:51 -0400 Subject: [PATCH 01/24] +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 44fbb09314be9bed90396dcb97cc806361549119 Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 31 Mar 2025 13:16:41 -0400 Subject: [PATCH 02/24] 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 4865b06f8a8c1646115399e29c2e06795e5c6ccf Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 10:17:43 -0400 Subject: [PATCH 03/24] 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 37901a3a3789d1869931b8f23f8c77f6482336c7 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 11:02:26 -0400 Subject: [PATCH 04/24] 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 f105198337e88be6922fa8048d67f6227ca349c9 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 16:00:33 -0400 Subject: [PATCH 05/24] 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 dab1b6b9ebe9e177eaf69176d711fef1e1f08621 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 3 Apr 2025 16:08:36 -0400 Subject: [PATCH 06/24] 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 951e0b88f98c31798408ada7106f410045ffcfe0 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 4 Apr 2025 14:23:27 -0400 Subject: [PATCH 07/24] 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 a1b6ed98f7eb417479be468b3530d722a30e273b Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 11 Apr 2025 11:27:44 -0400 Subject: [PATCH 08/24] 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 523bf7ce5dcadbb017dbb15f0d8cb1a32d4e7208 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 11 Apr 2025 11:29:13 -0400 Subject: [PATCH 09/24] 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 3b66b725eb6cc7705a7d7d84770744c616166639 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 15 Apr 2025 15:02:24 -0400 Subject: [PATCH 10/24] 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 4b5a7a95232d4fe5a98a28d821fd7539773b6a27 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 15 Apr 2025 16:10:25 -0400 Subject: [PATCH 11/24] 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 11cad3bba6ce5e59274e8deffa020de1caebdc4f Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 12:31:12 -0400 Subject: [PATCH 12/24] 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 3d549bc6f0f8c3d86542d35d2050d7c10b6232be Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 12:35:41 -0400 Subject: [PATCH 13/24] 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 06c02f8bbb9bb322663e77342e090a61e270c5a4 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 13:44:55 -0400 Subject: [PATCH 14/24] 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 f23e9665f9fa065ebd53630aeea311b2c5568497 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 13:55:47 -0400 Subject: [PATCH 15/24] 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 21f48bb5d2fecb14242d829c1cbdfceab1a5af18 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 14:00:47 -0400 Subject: [PATCH 16/24] 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 604003c09a694f70e93d7218e111f5c098b9284c Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 15:48:36 -0400 Subject: [PATCH 17/24] 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 963943000643ef1172320e9cc7e3933a78b8a35a Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:31:04 -0400 Subject: [PATCH 18/24] 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 55b5444a5c52a2dd4f509bbe3f447fc2d6444186 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:36:11 -0400 Subject: [PATCH 19/24] 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 4b7f8c4134836032f0ebdf5aaa03016920667177 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:36:44 -0400 Subject: [PATCH 20/24] 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 0737f05ac26216f5e8af052374ecb786f546d10c Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:41:40 -0400 Subject: [PATCH 21/24] 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 546496e114b919d4c1c247a17dc36eb245e7ffe6 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:44:55 -0400 Subject: [PATCH 22/24] 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 94cd195b9cc7f3702840e0991fbffbd136ff9697 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 16:52:25 -0400 Subject: [PATCH 23/24] 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 4131f1ffab03007db1c4b7ee07905f2ccd04371a Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 17 Apr 2025 17:13:33 -0400 Subject: [PATCH 24/24] 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