From e2bdf69f0312596112dbed04ab3b36706bde767b Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 14 Feb 2025 12:05:43 -0500 Subject: [PATCH 01/22] +psc_shock currently just a maxwellian --- src/CMakeLists.txt | 12 ++- src/psc_shock.cxx | 258 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 265 insertions(+), 5 deletions(-) create mode 100644 src/psc_shock.cxx diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 402396b11..de0e7a3ec 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -8,9 +8,10 @@ macro(add_psc_executable name) target_link_libraries(${name} psc) endmacro(add_psc_executable) -if (NOT USE_CUDA) +if(NOT USE_CUDA) add_psc_executable(psc_bgk) endif() + add_psc_executable(psc_bubble_yz) add_psc_executable(psc_flatfoil_yz) add_psc_executable(psc_flatfoil_yz_small) @@ -19,15 +20,16 @@ add_psc_executable(psc_harris_xz) add_psc_executable(psc_harris_yz) add_psc_executable(psc_2d_shock) add_psc_executable(psc_radiation) +add_psc_executable(psc_shock) -if (NOT USE_CUDA) +if(NOT USE_CUDA) install( TARGETS psc_bgk RUNTIME DESTINATION bin - ) + ) endif() install( - TARGETS psc_bubble_yz psc_flatfoil_yz psc_whistler psc_harris_xz + TARGETS psc_bubble_yz psc_flatfoil_yz psc_whistler psc_harris_xz psc_shock RUNTIME DESTINATION bin - ) +) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx new file mode 100644 index 000000000..c9593e0f1 --- /dev/null +++ b/src/psc_shock.cxx @@ -0,0 +1,258 @@ +#include +#include +#include + +#include "DiagnosticsDefault.h" +#include "OutputFieldsDefault.h" +#include "psc_config.hxx" + +#include "psc_bgk_util/bgk_params.hxx" +#include "psc_bgk_util/table.hxx" +#include "psc_bgk_util/params_parser.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 + +// ---------------------------------------------------------------------- + +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; + +double v_thermal = 1.0; +double box_size = 1.0; +} // namespace + +// ====================================================================== +// setupParameters + +void setupParameters(int argc, char** argv) +{ + psc_params.nmax = 100; + psc_params.stats_every = 10; + 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{ + {1, 16, 16}, // n grid points + {box_size, box_size, box_size}, // physical lengths + {-box_size / 2.0, -box_size / 2.0, -box_size / 2.0}, // origin offset + {1, 1, 1}}; // n patches + + auto bc = + psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_PERIODIC}, + {BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_PERIODIC}, + {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_PERIODIC}, + {BND_PRT_PERIODIC, BND_PRT_PERIODIC, 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) +{ + SetupParticles setup_particles(*grid_ptr); + setup_particles.centerer = Centering::Centerer(Centering::NC); + + auto init_np = [&](int kind, Double3 crd, int p, Int3 idx, + psc_particle_np& np) { + double temperature = sqr(v_thermal); + + np.n = 1.0; + np.p = + setup_particles.createMaxwellian({np.kind, + np.n, + {0, 0, 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]) { return 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 = 10; + + // -- Collision + int collision_interval = 0; + double collision_nu = .1; + Collision collision{grid, collision_interval, collision_nu}; + + // -- Checks + ChecksParams checks_params{}; + checks_params.gauss_every_step = 50; + // checks_params.gauss_dump_always = true; + checks_params.gauss_threshold = 1e-5; + + 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 = 5; + 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; + 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); + initializeFields(mflds); + + // ---------------------------------------------------------------------- + // run the simulation + + auto psc = + makePscIntegrator(psc_params, *grid_ptr, mflds, mprts, balance, + collision, checks, marder, diagnostics); + + psc.integrate(); +} + +// ====================================================================== +// main + +int main(int argc, char** argv) +{ + // psc_init(argc, argv); + // FIXME restore whatever previous functionality there was with options + int temp = 1; + psc_init(temp, argv); + + run(argc, argv); + + psc_finalize(); + return 0; +} From 2acaaf3a0bd56a00c55a05464c1e18030db7d22c Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 18 Feb 2025 11:22:08 -0500 Subject: [PATCH 02/22] psc_shock: fix superluminal temp --- src/psc_shock.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index c9593e0f1..5ef2c5c1f 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -44,7 +44,7 @@ namespace // General PSC parameters PscParams psc_params; -double v_thermal = 1.0; +double v_thermal = 1e-3; double box_size = 1.0; } // namespace From 81c9380ecf103646a0fa70a7fafe799c43cc8083 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 18 Feb 2025 11:22:22 -0500 Subject: [PATCH 03/22] psc_shock: add bulk flow --- src/psc_shock.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 5ef2c5c1f..158bf1839 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -124,7 +124,7 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) np.p = setup_particles.createMaxwellian({np.kind, np.n, - {0, 0, 0}, + {0, 10 * v_thermal, 0}, {temperature, temperature, temperature}, np.tag}); }; From f49754f8735c12dec230bc8b92994d8c02e826e6 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 18 Feb 2025 11:28:48 -0500 Subject: [PATCH 04/22] psc_shock: make y bounds reflective --- src/psc_shock.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 158bf1839..a74a4b8af 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -78,10 +78,10 @@ Grid_t* setupGrid() {1, 1, 1}}; // n patches auto bc = - psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_PERIODIC}, - {BND_FLD_PERIODIC, BND_FLD_PERIODIC, BND_FLD_PERIODIC}, - {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_PERIODIC}, - {BND_PRT_PERIODIC, BND_PRT_PERIODIC, BND_PRT_PERIODIC}}; + 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"}; From bd0c0cdebfaee0223c4716739f45a6f7861cd48e Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 18 Feb 2025 11:28:56 -0500 Subject: [PATCH 05/22] psc_shock: run for longer --- src/psc_shock.cxx | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index a74a4b8af..a8b74d1e3 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -53,8 +53,8 @@ double box_size = 1.0; void setupParameters(int argc, char** argv) { - psc_params.nmax = 100; - psc_params.stats_every = 10; + psc_params.nmax = 10000; + psc_params.stats_every = 100; psc_params.cfl = .75; psc_params.write_checkpoint_every_step = 0; @@ -182,7 +182,7 @@ static void run(int argc, char** argv) Balance balance{.1}; // -- Sort - psc_params.sort_interval = 10; + psc_params.sort_interval = 100; // -- Collision int collision_interval = 0; @@ -201,7 +201,7 @@ static void run(int argc, char** argv) double marder_diffusion = 0.9; int marder_loop = 3; bool marder_dump = false; - psc_params.marder_interval = 5; + psc_params.marder_interval = 50; Marder marder(grid, marder_diffusion, marder_loop, marder_dump); // ---------------------------------------------------------------------- @@ -209,13 +209,13 @@ static void run(int argc, char** argv) // -- output fields OutputFieldsParams outf_params{}; - outf_params.fields.pfield.out_interval = 1; - outf_params.moments.pfield.out_interval = 1; + outf_params.fields.pfield.out_interval = 100; + outf_params.moments.pfield.out_interval = 100; OutputFields outf{grid, outf_params}; // -- output particles OutputParticlesParams outp_params{}; - outp_params.every_step = 1; + outp_params.every_step = 100; outp_params.data_dir = "."; outp_params.basename = "prt"; OutputParticles outp{grid, outp_params}; From daa5ca39aa30699015a9a9c405a57e0210580565 Mon Sep 17 00:00:00 2001 From: James McClung Date: Thu, 6 Mar 2025 14:10:47 -0500 Subject: [PATCH 06/22] psc_shock: add a background B field --- src/psc_shock.cxx | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index a8b74d1e3..92e805479 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -148,7 +148,12 @@ void fillGhosts(MF& mfld, int compBegin, int compEnd) void initializeFields(MfieldsState& mflds) { - setupFields(mflds, [&](int component, double coords[3]) { return 0.0; }); + setupFields(mflds, [&](int component, double coords[3]) { + if (component == HZ) { + return 0.01; + } + return 0.0; + }); } // ====================================================================== From 76989d4bb1d5234af85c9323be21c1718a7612b7 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 7 Mar 2025 14:22:42 -0500 Subject: [PATCH 07/22] psc_shock: init to "tiny" .1 AU helio case --- src/psc_shock.cxx | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 92e805479..e28f80e43 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -44,8 +44,12 @@ namespace // General PSC parameters PscParams psc_params; -double v_thermal = 1e-3; -double box_size = 1.0; +double v_thermal = 8.85E-03; +double v_flow = 2.67E-03; +double b_background = 6.9713E-02; +double len_x = 1.0; +double len_y = 500.0; +double len_z = 20.0; } // namespace // ====================================================================== @@ -53,7 +57,7 @@ double box_size = 1.0; void setupParameters(int argc, char** argv) { - psc_params.nmax = 10000; + psc_params.nmax = 34000; psc_params.stats_every = 100; psc_params.cfl = .75; @@ -71,11 +75,10 @@ void setupParameters(int argc, char** argv) Grid_t* setupGrid() { // FIXME add a check to catch mismatch between Dim and n grid points early - auto domain = Grid_t::Domain{ - {1, 16, 16}, // n grid points - {box_size, box_size, box_size}, // physical lengths - {-box_size / 2.0, -box_size / 2.0, -box_size / 2.0}, // origin offset - {1, 1, 1}}; // n patches + auto domain = Grid_t::Domain{{1, 1000, 40}, // n grid points + {len_x, len_y, len_z}, // physical lengths + {-len_x / 2.0, 0, 0}, // origin offset + {1, 1, 1}}; // n patches auto bc = psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL, BND_FLD_PERIODIC}, @@ -85,11 +88,12 @@ 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, 100.0, "i"}; // --- generic setup - auto norm_params = Grid_t::NormalizationParams::dimensionless(); + auto norm_params = Grid_t::NormalizationParams(); norm_params.nicell = 100; + norm_params.n0 = 5.00E+08; double dt = psc_params.cfl * courant_length(domain); Grid_t::Normalization norm{norm_params}; @@ -124,7 +128,7 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) np.p = setup_particles.createMaxwellian({np.kind, np.n, - {0, 10 * v_thermal, 0}, + {0, v_flow, 0}, {temperature, temperature, temperature}, np.tag}); }; @@ -149,8 +153,8 @@ void fillGhosts(MF& mfld, int compBegin, int compEnd) void initializeFields(MfieldsState& mflds) { setupFields(mflds, [&](int component, double coords[3]) { - if (component == HZ) { - return 0.01; + if (component == HX) { + return b_background; } return 0.0; }); From b2332a5cfb69c194121e2746fc372a0aecbeddae Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 11 Mar 2025 14:14:40 -0400 Subject: [PATCH 08/22] psc_shock: x corner to 0 --- src/psc_shock.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index e28f80e43..e62d29d53 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -77,8 +77,8 @@ Grid_t* setupGrid() // FIXME add a check to catch mismatch between Dim and n grid points early auto domain = Grid_t::Domain{{1, 1000, 40}, // n grid points {len_x, len_y, len_z}, // physical lengths - {-len_x / 2.0, 0, 0}, // origin offset - {1, 1, 1}}; // n patches + {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}, From 04148f2f44ac5fa5910efdeaf6f66af8486bce73 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 11 Mar 2025 14:27:31 -0400 Subject: [PATCH 09/22] psc_shock: output less frequently --- src/psc_shock.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index e62d29d53..b414ef59c 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -218,13 +218,13 @@ static void run(int argc, char** argv) // -- output fields OutputFieldsParams outf_params{}; - outf_params.fields.pfield.out_interval = 100; - outf_params.moments.pfield.out_interval = 100; + outf_params.fields.pfield.out_interval = 500; + outf_params.moments.pfield.out_interval = 500; OutputFields outf{grid, outf_params}; // -- output particles OutputParticlesParams outp_params{}; - outp_params.every_step = 100; + outp_params.every_step = 1000; outp_params.data_dir = "."; outp_params.basename = "prt"; OutputParticles outp{grid, outp_params}; From ebda6ee9cd335518de256163dc91f0775d4d3327 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 11 Mar 2025 14:33:29 -0400 Subject: [PATCH 10/22] psc_shock: set temperature directly v_thermal is different for ions and electrons --- src/psc_shock.cxx | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index b414ef59c..1c7cad714 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -44,7 +44,7 @@ namespace // General PSC parameters PscParams psc_params; -double v_thermal = 8.85E-03; +double temperature = 7.8278E-05; double v_flow = 2.67E-03; double b_background = 6.9713E-02; double len_x = 1.0; @@ -122,8 +122,6 @@ 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 = sqr(v_thermal); - np.n = 1.0; np.p = setup_particles.createMaxwellian({np.kind, From f05ceeaf2edce3fbb9bf23b404b52a217fb13b58 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 11 Mar 2025 15:05:25 -0400 Subject: [PATCH 11/22] psc_shock: go back to dimensionless I'm uneasy about this being correct --- src/psc_shock.cxx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 1c7cad714..902988d11 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -91,9 +91,8 @@ Grid_t* setupGrid() kinds[KIND_ION] = {1.0, 100.0, "i"}; // --- generic setup - auto norm_params = Grid_t::NormalizationParams(); + auto norm_params = Grid_t::NormalizationParams::dimensionless(); norm_params.nicell = 100; - norm_params.n0 = 5.00E+08; double dt = psc_params.cfl * courant_length(domain); Grid_t::Normalization norm{norm_params}; From 63004d7c864acf7ab50f32e7cdbfa8e8c74b9eca Mon Sep 17 00:00:00 2001 From: James McClung Date: Wed, 12 Mar 2025 12:22:32 -0400 Subject: [PATCH 12/22] psc_shock: fix particle centering --- src/psc_shock.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 902988d11..ba8d2e445 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -117,7 +117,7 @@ Grid_t* setupGrid() void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) { SetupParticles setup_particles(*grid_ptr); - setup_particles.centerer = Centering::Centerer(Centering::NC); + setup_particles.centerer = Centering::Centerer(Centering::CC); auto init_np = [&](int kind, Double3 crd, int p, Int3 idx, psc_particle_np& np) { From c56174b18467e105bd878cf262e761e8e90dd18c Mon Sep 17 00:00:00 2001 From: James McClung Date: Mon, 17 Mar 2025 15:34:36 -0400 Subject: [PATCH 13/22] psc_shock: configure to Egedal 2012 ish --- src/psc_shock.cxx | 46 ++++++++++++++++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 14 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index ba8d2e445..f2f6021ed 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -44,12 +44,28 @@ namespace // General PSC parameters PscParams psc_params; -double temperature = 7.8278E-05; -double v_flow = 2.67E-03; -double b_background = 6.9713E-02; -double len_x = 1.0; -double len_y = 500.0; -double len_z = 20.0; +double electron_temperature = 2.0833E-02; +double ion_temperature = 1.0417E-01; +double ion_mass = 400.0; +double v_flow = 6.2460E-02; +double b_background = 5.0002E-01; + +int nx = 1; +int ny = 40960; // Egedal 2012 uses 40960 +int nz = 2; // Egedal 2012 uses 3840 +int nt = 3700100 / 30; // need 3700100 to match Egedal 2012 + +double dx = 1.0; +double dy = 1.5625E-01; +double dz = dy; + +double len_x = nx * dx; +double len_y = ny * dy; +double len_z = nz * dz; + +int n_writes = 100; +int out_interval = nt / n_writes; + } // namespace // ====================================================================== @@ -57,8 +73,8 @@ double len_z = 20.0; void setupParameters(int argc, char** argv) { - psc_params.nmax = 34000; - psc_params.stats_every = 100; + psc_params.nmax = nt; + psc_params.stats_every = out_interval; psc_params.cfl = .75; psc_params.write_checkpoint_every_step = 0; @@ -75,7 +91,7 @@ void setupParameters(int argc, char** argv) Grid_t* setupGrid() { // FIXME add a check to catch mismatch between Dim and n grid points early - auto domain = Grid_t::Domain{{1, 1000, 40}, // n grid points + 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 {1, 1, 1}}; // n patches @@ -88,7 +104,7 @@ Grid_t* setupGrid() auto kinds = Grid_t::Kinds(NR_KINDS); kinds[KIND_ELECTRON] = {-1.0, 1.0, "e"}; - kinds[KIND_ION] = {1.0, 100.0, "i"}; + kinds[KIND_ION] = {1.0, ion_mass, "i"}; // --- generic setup auto norm_params = Grid_t::NormalizationParams::dimensionless(); @@ -121,6 +137,8 @@ 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, @@ -197,7 +215,7 @@ static void run(int argc, char** argv) // -- Checks ChecksParams checks_params{}; - checks_params.gauss_every_step = 50; + checks_params.gauss_every_step = out_interval; // checks_params.gauss_dump_always = true; checks_params.gauss_threshold = 1e-5; @@ -215,13 +233,13 @@ static void run(int argc, char** argv) // -- output fields OutputFieldsParams outf_params{}; - outf_params.fields.pfield.out_interval = 500; - outf_params.moments.pfield.out_interval = 500; + 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 = 1000; + outp_params.every_step = out_interval; outp_params.data_dir = "."; outp_params.basename = "prt"; OutputParticles outp{grid, outp_params}; From 7869a7cab027a5d1ddd81957af1a4d0ca3e05f9b Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 21 Mar 2025 15:26:06 -0400 Subject: [PATCH 14/22] psc_shock: configure to parallel heliospheric case --- src/psc_shock.cxx | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index f2f6021ed..c233999f0 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -44,20 +44,23 @@ namespace // General PSC parameters PscParams psc_params; -double electron_temperature = 2.0833E-02; -double ion_temperature = 1.0417E-01; -double ion_mass = 400.0; -double v_flow = 6.2460E-02; -double b_background = 5.0002E-01; +// 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 = 40960; // Egedal 2012 uses 40960 -int nz = 2; // Egedal 2012 uses 3840 -int nt = 3700100 / 30; // need 3700100 to match Egedal 2012 +int ny = 12800; +int nz = 2; +int nt = 5505294; -double dx = 1.0; -double dy = 1.5625E-01; -double dz = dy; +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; @@ -168,10 +171,11 @@ void fillGhosts(MF& mfld, int compBegin, int compEnd) void initializeFields(MfieldsState& mflds) { setupFields(mflds, [&](int component, double coords[3]) { - if (component == HX) { - return b_background; + switch (component) { + case HX: return b_x; + case HY: return b_y; + default: return 0.0; } - return 0.0; }); } From c35ba20ebdea8e01aab4e914a732aa90df8b8d3c Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 21 Mar 2025 15:39:59 -0400 Subject: [PATCH 15/22] psc_shock: +n_patches --- src/psc_shock.cxx | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index c233999f0..4e42a7fd1 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -58,6 +58,10 @@ int ny = 12800; int nz = 2; int nt = 5505294; +int n_patches_x = 1; +int n_patches_y = 128; +int n_patches_z = 1; + double dx = 4.2849E+01; double dy = 3.3475E-01; double dz = 3.3475E-01; @@ -94,10 +98,11 @@ void setupParameters(int argc, char** argv) 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 - {1, 1, 1}}; // n patches + 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}, From 06a6ba01d804bdc4291ab7d7d3fc2403f5acb2d8 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 21 Mar 2025 15:46:02 -0400 Subject: [PATCH 16/22] psc_shock: 128 cells per patch --- src/psc_shock.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 4e42a7fd1..462abf2aa 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -54,7 +54,7 @@ double b_x = 0; double b_y = 2.4000E-02; int nx = 1; -int ny = 12800; +int ny = 128 * 128; int nz = 2; int nt = 5505294; From e645b4ed24e3f3b2ef8470289020685091788e86 Mon Sep 17 00:00:00 2001 From: James McClung Date: Fri, 4 Apr 2025 12:26:28 -0400 Subject: [PATCH 17/22] psc_shock: add electric field --- src/psc_shock.cxx | 34 +++++++++++++++++++++++++--------- 1 file changed, 25 insertions(+), 9 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 462abf2aa..c4ef8be8c 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -49,9 +49,21 @@ 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; + +double v_upstream_x = 0.0; +double v_upstream_y = 2.6685E-03; +double v_upstream_z = 0.0; + +double b_angle_y_to_x_deg = 15; +double b_angle_y_to_x = b_angle_y_to_x_deg * M_PI / 180; +double b_mag = 2.4000E-02; +double b_x = b_mag * sin(b_angle_y_to_x); +double b_y = b_mag * cos(b_angle_y_to_x); +double b_z = 0.0; + +double e_x = -(v_upstream_y * b_z - v_upstream_z * b_y); +double e_y = -(v_upstream_z * b_x - v_upstream_x * b_z); +double e_z = -(v_upstream_x * b_y - v_upstream_y * b_x); int nx = 1; int ny = 128 * 128; @@ -148,12 +160,12 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) 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}); + np.p = setup_particles.createMaxwellian( + {np.kind, + np.n, + {v_upstream_x, v_upstream_y, v_upstream_z}, + {temperature, temperature, temperature}, + np.tag}); }; partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, @@ -177,8 +189,12 @@ void initializeFields(MfieldsState& mflds) { setupFields(mflds, [&](int component, double coords[3]) { switch (component) { + case EX: return e_x; + case EY: return e_y; + case EZ: return e_z; case HX: return b_x; case HY: return b_y; + case HZ: return b_z; default: return 0.0; } }); From 90a866bb822ebd09bc8ac4bb0fae6bba900cb83d Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 8 Apr 2025 10:36:02 -0400 Subject: [PATCH 18/22] psc_shock: disable marder --- src/psc_shock.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index c4ef8be8c..1d1ccf98d 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -250,7 +250,7 @@ static void run(int argc, char** argv) double marder_diffusion = 0.9; int marder_loop = 3; bool marder_dump = false; - psc_params.marder_interval = 50; + psc_params.marder_interval = -1; Marder marder(grid, marder_diffusion, marder_loop, marder_dump); // ---------------------------------------------------------------------- From 7c995270a6c40b3c92d87b479aac05d6da315810 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 8 Apr 2025 10:26:18 -0400 Subject: [PATCH 19/22] psc_shock: reconfigure to #76 --- src/psc_shock.cxx | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 1d1ccf98d..413f7285f 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -44,11 +44,9 @@ 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 ion_mass = 1.0000E+02; double v_upstream_x = 0.0; double v_upstream_y = 2.6685E-03; @@ -74,7 +72,7 @@ int n_patches_x = 1; int n_patches_y = 128; int n_patches_z = 1; -double dx = 4.2849E+01; +double dx = 3.3475E-01; double dy = 3.3475E-01; double dz = 3.3475E-01; From b18b93f8ce2a1560cfc86351d3723d072ec99a56 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 8 Apr 2025 14:13:12 -0400 Subject: [PATCH 20/22] psc_shock: switch to ParsedParams --- src/psc_shock.cxx | 114 +++++++++++++++++++++++++++++++--------------- 1 file changed, 77 insertions(+), 37 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 413f7285f..20616cc39 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -6,8 +6,6 @@ #include "OutputFieldsDefault.h" #include "psc_config.hxx" -#include "psc_bgk_util/bgk_params.hxx" -#include "psc_bgk_util/table.hxx" #include "psc_bgk_util/params_parser.hxx" // ====================================================================== @@ -44,44 +42,36 @@ namespace // General PSC parameters PscParams psc_params; -double electron_temperature = 7.8278E-05; -double ion_temperature = 7.8278E-05; -double ion_mass = 1.0000E+02; +double electron_temperature; +double ion_temperature; +double electron_mass; +double ion_mass; -double v_upstream_x = 0.0; -double v_upstream_y = 2.6685E-03; -double v_upstream_z = 0.0; +double v_upstream_x; +double v_upstream_y; +double v_upstream_z; -double b_angle_y_to_x_deg = 15; -double b_angle_y_to_x = b_angle_y_to_x_deg * M_PI / 180; -double b_mag = 2.4000E-02; -double b_x = b_mag * sin(b_angle_y_to_x); -double b_y = b_mag * cos(b_angle_y_to_x); -double b_z = 0.0; +double b_x; +double b_y; +double b_z; -double e_x = -(v_upstream_y * b_z - v_upstream_z * b_y); -double e_y = -(v_upstream_z * b_x - v_upstream_x * b_z); -double e_z = -(v_upstream_x * b_y - v_upstream_y * b_x); +double e_x; +double e_y; +double e_z; -int nx = 1; -int ny = 128 * 128; -int nz = 2; -int nt = 5505294; +int nx; +int ny; +int nz; -int n_patches_x = 1; -int n_patches_y = 128; -int n_patches_z = 1; +int n_patches_x; +int n_patches_y; +int n_patches_z; -double dx = 3.3475E-01; -double dy = 3.3475E-01; -double dz = 3.3475E-01; +double len_x; +double len_y; +double len_z; -double len_x = nx * dx; -double len_y = ny * dy; -double len_z = nz * dz; - -int n_writes = 100; -int out_interval = nt / n_writes; +int out_interval; } // namespace @@ -90,11 +80,61 @@ int out_interval = nt / n_writes; void setupParameters(int argc, char** argv) { - psc_params.nmax = nt; - psc_params.stats_every = out_interval; - psc_params.cfl = .75; + if (argc != 2) { + std::cout << "Usage: " << argv[0] << " path/to/params\nExiting." + << std::endl; + exit(1); + } + std::string path_to_params(argv[1]); + ParsedParams parsedParams(path_to_params); + psc_params.stats_every = 1000; + psc_params.cfl = parsedParams.getOrDefault("cfl", .75); psc_params.write_checkpoint_every_step = 0; + + electron_temperature = parsedParams.get("electron_temperature"); + ion_temperature = parsedParams.get("ion_temperature"); + electron_mass = parsedParams.get("electron_mass"); + ion_mass = parsedParams.get("ion_mass"); + + v_upstream_x = parsedParams.get("v_upstream_x"); + v_upstream_y = parsedParams.get("v_upstream_y"); + v_upstream_z = parsedParams.get("v_upstream_z"); + + double b_angle_y_to_x_deg = parsedParams.get("b_angle_y_to_x_deg"); + double b_angle_y_to_x = b_angle_y_to_x_deg * M_PI / 180; + double b_mag = parsedParams.get("b_mag"); + b_x = b_mag * sin(b_angle_y_to_x); + b_y = b_mag * cos(b_angle_y_to_x); + b_z = 0.0; + + e_x = -(v_upstream_y * b_z - v_upstream_z * b_y); + e_y = -(v_upstream_z * b_x - v_upstream_x * b_z); + e_z = -(v_upstream_x * b_y - v_upstream_y * b_x); + + nx = parsedParams.get("nx"); + ny = parsedParams.get("ny"); + nz = parsedParams.get("nz"); + psc_params.nmax = parsedParams.get("nt"); + + n_patches_x = parsedParams.get("n_patches_x"); + n_patches_y = parsedParams.get("n_patches_y"); + n_patches_z = parsedParams.get("n_patches_z"); + + double dx = parsedParams.get("dx"); + double dy = parsedParams.get("dy"); + double dz = parsedParams.get("dz"); + + len_x = nx * dx; + len_y = ny * dy; + len_z = nz * dz; + + int n_writes = parsedParams.getOrDefault("n_writes", 100); + out_interval = psc_params.nmax / n_writes; + + std::ifstream src(path_to_params, std::ios::binary); + std::ofstream dst("params_record.txt", std::ios::binary); + dst << src.rdbuf(); } // ====================================================================== @@ -121,7 +161,7 @@ Grid_t* setupGrid() {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_ELECTRON] = {-1.0, electron_mass, "e"}; kinds[KIND_ION] = {1.0, ion_mass, "i"}; // --- generic setup From be5f852f2dc485e04221c4b506dded0ab281e832 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 8 Apr 2025 14:20:01 -0400 Subject: [PATCH 21/22] psc_shock: use initial gamma correction --- src/psc_shock.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 20616cc39..46370ed51 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -192,6 +192,7 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) { SetupParticles setup_particles(*grid_ptr); setup_particles.centerer = Centering::Centerer(Centering::CC); + setup_particles.initial_momentum_gamma_correction = true; auto init_np = [&](int kind, Double3 crd, int p, Int3 idx, psc_particle_np& np) { From 568547b6f8fb6e625e0ead9ee8421f2cf3c9b1a7 Mon Sep 17 00:00:00 2001 From: James McClung Date: Tue, 8 Apr 2025 14:21:32 -0400 Subject: [PATCH 22/22] psc_shock: use radians for B angle --- src/psc_shock.cxx | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 46370ed51..d90053dad 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -101,11 +101,10 @@ void setupParameters(int argc, char** argv) v_upstream_y = parsedParams.get("v_upstream_y"); v_upstream_z = parsedParams.get("v_upstream_z"); - double b_angle_y_to_x_deg = parsedParams.get("b_angle_y_to_x_deg"); - double b_angle_y_to_x = b_angle_y_to_x_deg * M_PI / 180; + double b_angle_y_to_x_rad = parsedParams.get("b_angle_y_to_x_rad"); double b_mag = parsedParams.get("b_mag"); - b_x = b_mag * sin(b_angle_y_to_x); - b_y = b_mag * cos(b_angle_y_to_x); + b_x = b_mag * sin(b_angle_y_to_x_rad); + b_y = b_mag * cos(b_angle_y_to_x_rad); b_z = 0.0; e_x = -(v_upstream_y * b_z - v_upstream_z * b_y);