diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index ad70b0b53..0ccc42478 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -25,7 +25,6 @@ using PscConfig = PscConfig1vbecDouble; // ---------------------------------------------------------------------- -using BgkMfields = PscConfig::Mfields; using MfieldsState = PscConfig::MfieldsState; using Mparticles = PscConfig::Mparticles; using Balance = PscConfig::Balance; @@ -73,6 +72,8 @@ double len_z; int out_interval; +bool mirror_domain; + } // namespace // ====================================================================== @@ -131,6 +132,8 @@ void setupParameters(int argc, char** argv) int n_writes = parsedParams.getOrDefault("n_writes", 100); out_interval = psc_params.nmax / n_writes; + mirror_domain = parsedParams.getOrDefault("mirror_domain", false); + std::ifstream src(path_to_params, std::ios::binary); std::ofstream dst("params_record.txt", std::ios::binary); dst << src.rdbuf(); @@ -147,17 +150,19 @@ 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 - {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}}; + int y_mult = mirror_domain ? 2 : 1; + auto domain = Grid_t::Domain{ + {nx, y_mult * ny, nz}, // n grid points + {len_x, y_mult * len_y, len_z}, // physical lengths + {0, 0, 0}, // location of lower corner + {n_patches_x, y_mult * n_patches_y, n_patches_z}}; // n patches + + auto BND_FLD_Y = mirror_domain ? BND_FLD_PERIODIC : BND_FLD_CONDUCTING_WALL; + auto BND_PRT_Y = mirror_domain ? BND_PRT_PERIODIC : BND_PRT_REFLECTING; + auto bc = psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_Y, BND_FLD_PERIODIC}, + {BND_FLD_PERIODIC, BND_FLD_Y, BND_FLD_PERIODIC}, + {BND_PRT_PERIODIC, BND_PRT_Y, BND_PRT_PERIODIC}, + {BND_PRT_PERIODIC, BND_PRT_Y, BND_PRT_PERIODIC}}; auto kinds = Grid_t::Kinds(NR_KINDS); kinds[KIND_ELECTRON] = {-1.0, electron_mass, "e"}; @@ -198,10 +203,11 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) double temperature = np.kind == KIND_ION ? ion_temperature : electron_temperature; np.n = 1.0; + double vy_coef = mirror_domain && crd[1] > len_y ? -1 : 1; np.p = setup_particles.createMaxwellian( {np.kind, np.n, - {v_upstream_x, v_upstream_y, v_upstream_z}, + {v_upstream_x, vy_coef * v_upstream_y, v_upstream_z}, {temperature, temperature, temperature}, np.tag}); }; @@ -226,10 +232,11 @@ void fillGhosts(MF& mfld, int compBegin, int compEnd) void initializeFields(MfieldsState& mflds) { setupFields(mflds, [&](int component, double coords[3]) { + double e_coef = mirror_domain && coords[1] > len_y ? -1 : 1; switch (component) { - case EX: return e_x; - case EY: return e_y; - case EZ: return e_z; + case EX: return e_coef * e_x; + case EY: return e_coef * e_y; + case EZ: return e_coef * e_z; case HX: return b_x; case HY: return b_y; case HZ: return b_z; @@ -257,9 +264,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 @@ -279,8 +283,7 @@ static void run(int argc, char** argv) // -- Checks ChecksParams checks_params{}; checks_params.gauss.check_interval = out_interval; - // checks_params.gauss.dump_always = true; - checks_params.gauss.err_threshold = 1e-5; + checks_params.continuity.check_interval = out_interval; Checks checks{grid, MPI_COMM_WORLD, checks_params};