diff --git a/src/include/boundary_injector.hxx b/src/include/boundary_injector.hxx deleted file mode 100644 index 69cf4fea4f..0000000000 --- a/src/include/boundary_injector.hxx +++ /dev/null @@ -1,168 +0,0 @@ -#pragma once - -#include - -#include "grid.hxx" -#include "rng.hxx" -#include "particle.h" -#include -#include "pushp.hxx" -#include "dim.hxx" -#include "setup_particles.hxx" -#include "kg/VecRange.hxx" -#include "../libpsc/psc_push_particles/inc_push.cxx" -#include "injector_base.hxx" - -/// @brief A particle generator for use with @ref BoundaryInjector. Samples -/// particles from a (possibly shifted) Maxwellian distribution. -class ParticleGeneratorMaxwellian -{ -public: - using Real = psc::particle::Inject::Real; - using Real3 = psc::particle::Inject::Real3; - - // FIXME would be nice to just pass 1 thing for kind-related info - ParticleGeneratorMaxwellian(int kind_idx, Grid_t::Kind kind, Real3 mean_u, - Real3 temperature, bool correct_gamma = false) - : kind_idx{kind_idx}, correct_gamma{correct_gamma} - { - for (int d = 0; d < 3; d++) { - Real stdev_u = sqrt(temperature[d] / kind.m); - vdfs[d] = VelocityDistributionFunction{mean_u[d], stdev_u}; - } - } - - psc::particle::Inject get(Real3 min_pos, Real3 pos_range) - { - Real3 x; - for (int d = 0; d < 3; d++) { - x[d] = min_pos[d] + uniform_dist.get() * pos_range[d]; - } - - Real3 u{vdfs[0].get(), vdfs[1].get(), vdfs[2].get()}; - if (correct_gamma) { - u = vel_to_4vel(u); - } - Real w = 1.0; - psc::particle::Tag tag = 0; - - return {x, u, w, kind_idx, tag}; - } - -private: - using VelocityDistributionFunction = rng::Normal; - Vec3 vdfs; - int kind_idx; - bool correct_gamma; - rng::Uniform uniform_dist{0.0, 1.0}; -}; - -/// @brief Injects particles on a given boundary, sampling from a given particle -/// generator. For precise control over multiple particle species, use one -/// BoundaryInjector per species. -/// @tparam PARTICLE_GENERATOR a type that defines `get(min_pos, pos_range)` and -/// returns an injectable particle within that range of positions (usually a -/// grid cell); see @ref ParticleGeneratorMaxwellian -/// @tparam PUSH_PARTICLES type that provides the types `Mparticles`, -/// `MfieldsState`, `Current`, `real_t`, etc. -template -class BoundaryInjector - : public InjectorBase -{ - static const int INJECT_DIM_IDX_ = 1; - -public: - using ParticleGenerator = PARTICLE_GENERATOR; - using PushParticles = PUSH_PARTICLES; - - using Mparticles = typename PushParticles::Mparticles; - using MfieldsState = typename PushParticles::MfieldsState; - using Current = typename PushParticles::Current; - using real_t = typename PushParticles::real_t; - using Real3 = Vec3; - - BoundaryInjector(ParticleGenerator particle_generator, Grid_t& grid) - : particle_generator_{particle_generator}, - advance_{grid.dt}, - prts_per_unit_density_{grid.norm.prts_per_unit_density} - {} - - /// Injects particles at the lower y-bound as if there were a population of - /// particles just beyond the edge. The imaginary particle population has unit - /// density, and individual particles from that population are sampled using - /// the given ParticleGenerator. - /// - /// Some of these limitations may be removed in the future. - void inject(Mparticles& mprts, MfieldsState& mflds) override - { - static_assert(INJECT_DIM_IDX_ == 1, - "only injection at lower bound of y is supported"); - - const Grid_t& grid = mprts.grid(); - auto injectors_by_patch = mprts.injector(); - - Real3 dxi = grid.domain.dx_inv; - Current current(grid); - - for (int p = 0; p < grid.n_patches(); p++) { - if (!grid.atBoundaryLo(p, INJECT_DIM_IDX_)) { - continue; - } - - Int3 ilo = {0, 0, 0}; - Int3 ihi = grid.ldims; - - ilo[INJECT_DIM_IDX_] = -1; - ihi[INJECT_DIM_IDX_] = 0; - - auto&& injector = injectors_by_patch[p]; - auto flds = mflds[p]; - typename Current::fields_t J(flds); - - for (Int3 initial_idx : VecRange(ilo, ihi)) { - Real3 cell_corner = Double3(initial_idx) * grid.domain.dx; - int n_prts_to_try_inject = - get_n_in_cell(1.0, prts_per_unit_density_, true); - - for (int prt_count = 0; prt_count < n_prts_to_try_inject; prt_count++) { - psc::particle::Inject prt = - particle_generator_.get(cell_corner, grid.domain.dx); - - Real3 v = advance_.calc_v(prt.u); - Real3 initial_x = prt.x; - advance_.push_x(prt.x, v); - - if (prt.x[INJECT_DIM_IDX_] < 0.0) { - // don't inject a particle that fails to enter the patch - continue; - } - - // GOTCHA: currently, injectors expect particle positions to be - // global, but current deposition expects patch-local - psc::particle::Inject prt_with_global_x = prt; - prt_with_global_x.x += grid.patches[p].xb; - injector(prt_with_global_x); - - // Update currents - // Taken from push_particles_1vb.hxx PushParticlesVb::push_mprts() - - Real3 initial_normalized_pos = initial_x * dxi; - Real3 final_normalized_pos = prt.x * dxi; - Int3 final_idx = final_normalized_pos.fint(); - - // CURRENT DENSITY BETWEEN (n+.5)*dt and (n+1.5)*dt - real_t qni_wni = grid.kinds[prt.kind].q * prt.w; - current.calc_j(J, initial_normalized_pos, final_normalized_pos, - final_idx, initial_idx, qni_wni, v); - } - } - } - } - -private: - ParticleGenerator particle_generator_; - // can't move along x or z, or else might leave patch - AdvanceParticle advance_; - real_t prts_per_unit_density_; -}; diff --git a/src/include/fields3d.hxx b/src/include/fields3d.hxx index 69e3b7d15b..c6cfb1ef4c 100644 --- a/src/include/fields3d.hxx +++ b/src/include/fields3d.hxx @@ -417,6 +417,7 @@ struct MfieldsStateFromMfields : MfieldsStateBase { using fields_view_t = typename Mfields::fields_view_t; using real_t = typename Mfields::real_t; + using Real3 = Vec3; using Real = typename Mfields::Real; using Storage = typename Mfields::Storage; using space = gt::space::host; diff --git a/src/include/fields3d.inl b/src/include/fields3d.inl index 42e32514b8..5e73dc3db8 100644 --- a/src/include/fields3d.inl +++ b/src/include/fields3d.inl @@ -39,8 +39,6 @@ public: // FIXME, should just check for consistency? (# ghosts might differ, too) // reader.get("ib", mflds.ib, launch); // reader.get("im", mflds.im, launch); - reader.beginStep(kg::io::StepMode::Read); - auto n_comps = mflds.n_comps(); auto shape = makeDims(n_comps, mflds.gdims()); assert(reader.variableShape() == shape); @@ -55,7 +53,6 @@ public: {}); //{ib, im}); } reader.performGets(); - reader.endStep(); for (int p = 0; p < mflds.n_patches(); p++) { auto h_flds = make_Fields3d(h_mflds[p]); diff --git a/src/include/injector_simple.hxx b/src/include/injector_simple.hxx index 5d0b3dce19..b31b1c7de9 100644 --- a/src/include/injector_simple.hxx +++ b/src/include/injector_simple.hxx @@ -36,6 +36,18 @@ struct InjectorSimple mprts_.push_back(p_, prt); } + void inject_local(const psc::particle::Inject& new_prt) + { + auto prt = + Particle{Real3(new_prt.x), + Real3(new_prt.u), + real_t(new_prt.w * mprts_.grid().kinds[new_prt.kind].q), + new_prt.kind, + mprts_.uid_gen(), + new_prt.tag}; + mprts_.push_back(p_, prt); + } + void reweight(const psc::particle::Inject& new_prt) { auto& grid = mprts_.grid(); diff --git a/src/include/output_fields.hxx b/src/include/output_fields.hxx index 1a39679d47..09c66a60a6 100644 --- a/src/include/output_fields.hxx +++ b/src/include/output_fields.hxx @@ -92,13 +92,17 @@ struct OutputTfieldItemParams : BaseOutputFieldItemParams // Returns whether to accumulate on this timestep. bool do_accum(int timestep) { + if (!enabled()) { + return false; + } + // next_out could be this timestep int n_intervals_elapsed = (timestep - 1) / out_interval; int next_out = out_interval * (n_intervals_elapsed + 1); bool in_averaging_range = next_out - timestep < average_length; bool on_averaging_step = (next_out - timestep) % sample_interval == 0; - return enabled() && in_averaging_range && on_averaging_step; + return in_averaging_range && on_averaging_step; } }; diff --git a/src/include/psc.hxx b/src/include/psc.hxx index 364f7d9cf6..ac1dda9fec 100644 --- a/src/include/psc.hxx +++ b/src/include/psc.hxx @@ -7,9 +7,11 @@ #include #include +#include "../libpsc/psc_bnd_fields/field_bc_base.hxx" +#include "../libpsc/psc_bnd_fields/conducting_wall.hxx" +#include "../libpsc/psc_particle_injectors/injector_base.hxx" #include "gauss_corrector_base.hxx" #include "diagnostic_base.hxx" -#include "injector_base.hxx" #include "external_current_base.hxx" #include #include @@ -117,6 +119,7 @@ struct Psc using BndFields = typename PscConfig::BndFields; using BndParticles = typename PscConfig::BndParticles; using Dim = typename PscConfig::Dim; + using FieldBcBaseT = FieldBcBase; using GaussCorrectorBaseT = GaussCorrectorBase; using DiagnosticBaseT = DiagnosticBase; using InjectorBaseT = InjectorBase; @@ -146,6 +149,20 @@ struct Psc } } +#ifndef USE_CUDA + for (int d = 0; d < 3; d++) { + using psc::bnd::LoHi; + using psc::bnd::field::ConductingWall; + + if (grid.bc.fld_lo[d] == BND_FLD_CONDUCTING_WALL) { + add_field_bc(new ConductingWall{d, LoHi::Lo}); + } + if (grid.bc.fld_hi[d] == BND_FLD_CONDUCTING_WALL) { + add_field_bc(new ConductingWall{d, LoHi::Hi}); + } + } +#endif + int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank); log_.open("mem-" + std::to_string(rank) + ".log"); @@ -162,6 +179,13 @@ struct Psc // TODO: improve ownership model: we should own these objects (i.e., use // unique_ptr), but don't want to burden the user with C++ boilerplate. + void add_field_bc(FieldBcBaseT* field_bc) + { + if (field_bc) { + field_bcs_.push_back(field_bc); + } + } + void add_gauss_corrector(GaussCorrectorBaseT* corrector) { if (corrector) { @@ -219,11 +243,17 @@ struct Psc void pre_first_step() { + for (auto field_bc : field_bcs_) { + field_bc->apply_h_bcs(mflds_); + } bndf.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); bnd_.fill_ghosts(mflds_, JXI, JXI + 3); + for (auto field_bc : field_bcs_) { + field_bc->apply_e_bcs(mflds_); + } bndf.fill_ghosts_E(mflds_); bnd_.fill_ghosts(mflds_, EX, EX + 3); @@ -286,7 +316,7 @@ struct Psc psc_stats_val[st_nr_particles] = mprts_.size(); - if (grid().timestep() % p_.stats_every == 0) { + if (p_.stats_every > 0 && grid().timestep() % p_.stats_every == 0) { print_status(); } @@ -414,6 +444,9 @@ struct Psc mpi_printf(comm, "***** Bnd fields J...\n"); prof_start(pr_bndf); + for (auto field_bc : field_bcs_) { + field_bc->apply_j_bcs(mflds_); + } bndf.add_ghosts_J(mflds_); bnd_.add_ghosts(mflds_, JXI, JXI + 3); bnd_.fill_ghosts(mflds_, JXI, JXI + 3); @@ -428,6 +461,9 @@ struct Psc mpi_printf(comm, "***** Bnd fields B (1 of 2)...\n"); prof_restart(pr_bndf); + for (auto field_bc : field_bcs_) { + field_bc->apply_h_bcs(mflds_); + } bndf.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); prof_stop(pr_bndf); @@ -441,6 +477,9 @@ struct Psc mpi_printf(comm, "***** Bnd fields E...\n"); prof_restart(pr_bndf); + for (auto field_bc : field_bcs_) { + field_bc->apply_e_bcs(mflds_); + } bndf.fill_ghosts_E(mflds_); bnd_.fill_ghosts(mflds_, EX, EX + 3); prof_stop(pr_bndf); @@ -463,6 +502,9 @@ struct Psc mpi_printf(comm, "***** Bnd fields B (2 of 2)...\n"); prof_restart(pr_bndf); + for (auto field_bc : field_bcs_) { + field_bc->apply_h_bcs(mflds_); + } bndf.fill_ghosts_H(mflds_); bnd_.fill_ghosts(mflds_, HX, HX + 3); prof_stop(pr_bndf); @@ -555,6 +597,7 @@ protected: Balance& balance_; Collision& collision_; Checks& checks_; + std::vector field_bcs_; std::vector gauss_correctors_; std::vector diagnostics_; std::vector injectors_; diff --git a/src/include/setup_particles.hxx b/src/include/setup_particles.hxx index f335c289a5..0dd10b22d4 100644 --- a/src/include/setup_particles.hxx +++ b/src/include/setup_particles.hxx @@ -29,19 +29,6 @@ struct psc_particle_np psc::particle::Tag tag; }; -/** - * @brief Calculates gamma * v for the given velocity v. - * @tparam Real real type - * @param v the actual velocity - * @return the spatial components of the corresponding 4-velocity - */ -template -Vec3 vel_to_4vel(Vec3 v) -{ - Real gamma = 1.0 / sqrt(1.0 - v.mag2()); - return v * gamma; -} - struct InitNptFunc { // Initialize particles according to a Maxwellian. @@ -121,6 +108,51 @@ int get_n_in_cell(real_t density, real_t prts_per_unit_density, return std::max(1, int(density * prts_per_unit_density + .5)); } +/** + * @brief Boosts velocities using cached intermediate values for a particular + * Lorentz frame. + */ +struct VelocityBooster +{ + VelocityBooster(Double3 frame_v) + : frame_gamma(1.0 / std::sqrt(1.0 - frame_v.mag2())), + frame_u(frame_v * frame_gamma), + frame_dir(frame_v / frame_v.mag()) + { + // FIXME kind of hacky + if (frame_v.mag2() == 0.0) { + frame_dir = {1, 0, 0}; + } + } + + /** + * @param prt_v a particle's "unprimed" proper velocity + * @return its "primed" proper velocity + */ + Double3 boost(Double3 prt_u) + { + double prt_gamma = std::sqrt(1.0 + prt_u.mag2()); + return prt_u + (frame_gamma - 1.0) * prt_u.dot(frame_dir) * frame_dir - + frame_u * prt_gamma; + } + + /** + * @param prt_v a particle's "unprimed" non-proper velocity + * @return its "primed" proper velocity + */ + Double3 boost_and_make_proper(Double3 prt_v) + { + double prt_gamma = 1.0 / std::sqrt(1.0 - prt_v.mag2()); + Double3 prt_u = prt_v * prt_gamma; + return prt_u + (frame_gamma - 1.0) * prt_u.dot(frame_dir) * frame_dir - + frame_u * prt_gamma; + } + + double frame_gamma; + Double3 frame_u; + Double3 frame_dir; +}; + // ====================================================================== // SetupParticles @@ -213,13 +245,27 @@ struct SetupParticles return [=]() { static rng::Normal dist; + if (initial_momentum_gamma_correction) { + // FIXME cache this (static doesn't work) + VelocityBooster booster{-npt.p}; + + Double3 prt_v; + for (int d = 0; d < 3; d++) { + // sample velocity in plasma frame + prt_v[d] = dist.get(0.0, std::sqrt(npt.T[d] / m)); + } + + // boost to lab frame + // FIXME should really sample from Maxwell-Juttner + // this hack interprests v as u to handle rare case when v>1 + // v<<1 => v~= u anyways + return booster.boost(prt_v); + } + Double3 p; for (int i = 0; i < 3; i++) p[i] = dist.get(npt.p[i], beta * std::sqrt(npt.T[i] / m)); - if (initial_momentum_gamma_correction) { - p = vel_to_4vel(p); - } return p; }; } @@ -287,8 +333,8 @@ struct SetupParticles "have the exact same initial position distribution. This results in " "a charge density of 0 if there are two species with opposite " "charges, but the resulting charge density is nonzero in general. In " - "the latter, case, take special care to ensure Gauss' law isn't " - "violated."); + "the latter case, take special care to ensure Gauss' law isn't " + "violated.\n"); } int seed = rng::detail::get_process_seed(); diff --git a/src/kg/include/kg/Vec3.h b/src/kg/include/kg/Vec3.h index b74c7db0c1..31b028cb50 100644 --- a/src/kg/include/kg/Vec3.h +++ b/src/kg/include/kg/Vec3.h @@ -295,6 +295,19 @@ struct Vec : gt::sarray return res; } + /** + * @brief Returns a copy of this vec but with one component value replaced. + * @param d component index + * @param val new component value + * @return the copy + */ + KG_INLINE Vec with_component(int d, T val) const + { + Vec res = *this; + res[d] = val; + return res; + } + // conversion to pointer KG_INLINE operator const T*() const { return this->data(); } diff --git a/src/libpsc/axis.hxx b/src/libpsc/axis.hxx new file mode 100644 index 0000000000..88200b6431 --- /dev/null +++ b/src/libpsc/axis.hxx @@ -0,0 +1,48 @@ +#pragma once + +namespace psc +{ + +struct Axis +{ + static Axis X; + static Axis Y; + static Axis Z; + + Axis(int axis) : axis{axis} {} + + operator int() const { return axis; } + + Axis& operator++() + { + this->axis += 1; + return *this; + } + + Axis operator++(int) + { + Axis temp = *this; + ++(*this); + return temp; + } + + Axis next() + { + int next = (axis + 1) % 3; + return next; + } + + Axis prev() + { + int next = (axis + 2) % 3; + return next; + } + + int axis; +}; + +Axis Axis::X = Axis(0); +Axis Axis::Y = Axis(1); +Axis Axis::Z = Axis(2); + +} // namespace psc \ No newline at end of file diff --git a/src/libpsc/psc_bnd/psc_bnd_util.hxx b/src/libpsc/psc_bnd/psc_bnd_util.hxx new file mode 100644 index 0000000000..67a4862c95 --- /dev/null +++ b/src/libpsc/psc_bnd/psc_bnd_util.hxx @@ -0,0 +1,15 @@ +#pragma once + +namespace psc +{ +namespace bnd +{ + +enum LoHi +{ + Lo, + Hi, +}; + +} // namespace bnd +} // namespace psc \ No newline at end of file diff --git a/src/libpsc/psc_bnd_fields/conducting_wall.hxx b/src/libpsc/psc_bnd_fields/conducting_wall.hxx new file mode 100644 index 0000000000..46444a89e4 --- /dev/null +++ b/src/libpsc/psc_bnd_fields/conducting_wall.hxx @@ -0,0 +1,254 @@ +#pragma once + +#include "psc.h" +#include "../axis.hxx" +#include "kg/Vec3.h" +#include "kg/VecRange.hxx" +#include "field_bc_base.hxx" +#include "field_bc_util.hxx" +#include "../psc_bnd/psc_bnd_util.hxx" + +namespace psc +{ +namespace bnd +{ +namespace field +{ + +/** + * @brief A perfect electrical conductor. This implementation assumes particles + * are specularly reflected. + * @tparam Dim dimension type + * @tparam MfieldsState fields type + */ +template +struct ConductingWall : FieldBcBase +{ + using dim_t = Dim; + + ConductingWall(Axis d, LoHi lohi) : d{d}, lohi{lohi} {} + + /** + * @brief Move currents deposited within the wall onto the domain. The normal + * component is flipped in accordance with specular reflection. + * @param mflds fields + */ + void apply_j_bcs(MfieldsState& mflds) override + { + const Grid_t& grid = mflds.grid(); + + const Int3 dhat = Int3::unit(d); + + const int J0 = JXI + d; + const int J1 = JXI + d.next(); + const int J2 = JXI + d.prev(); + + for (int p = 0; p < mflds.n_patches(); p++) { + if (lohi == Lo && grid.atBoundaryLo(p, d)) { + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = 0; + stop[d] = start[d] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + // 1. transverse components: wall is at i3 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(J1, i3 + j3) += F(J1, i3 - j3); + F(J1, i3 - j3) = 0.0; + F(J2, i3 + j3) += F(J2, i3 - j3); + F(J2, i3 - j3) = 0.0; + } + + // 2. normal component: wall is at i3-dhat/2 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(J0, i3 + j3 - dhat) -= F(J0, i3 - j3); + F(J0, i3 - j3) = 0.0; + } + } + } + + if (lohi == Hi && grid.atBoundaryHi(p, d)) { + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = grid.ldims[d]; + stop[d] = start[d] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + // 1. transverse components: wall is at i3, and there's one less ghost + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d] - 1; j3[d]++) { + F(J1, i3 - j3) += F(J1, i3 + j3); + F(J1, i3 + j3) = 0.0; + F(J2, i3 - j3) += F(J2, i3 + j3); + F(J2, i3 + j3) = 0.0; + } + + // 2. normal component: wall is at i3-dhat/2 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(J0, i3 - j3) -= F(J0, i3 + j3 - dhat); + F(J0, i3 + j3 - dhat) = 0.0; + } + } + } + } + } + + /** + * @brief Set transverse E at the wall's surface to 0. Nominally, E would be 0 + * at interior points too, but instead we mirror E so that reflecting + * particles feel the "right" electric force. + * + * Note that the normal force is reflected in the interior. To see why this is + * necessary, consider a particle with v=0 at the surface. If the normal E + * force on it is nonzero and towards the wall, the particle will + * spontaneously bounce away from the wall. It will then return to the + * wall—with more kinetic energy than before—and bounce again. Flipping normal + * E avoids this runaway effect. + * @param mflds fields + */ + void apply_e_bcs(MfieldsState& mflds) override + { + const Grid_t& grid = mflds.grid(); + + const Int3 dhat = Int3::unit(d); + + const int E0 = EX + d; + const int E1 = EX + d.next(); + const int E2 = EX + d.prev(); + + for (int p = 0; p < mflds.n_patches(); p++) { + if (lohi == Lo && grid.atBoundaryLo(p, d)) { + detail::set_lower_ghosts_to_nan(mflds, p, d, EX, true); + + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = 0; + stop[d] = start[d] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + // 1. transverse components: wall is at i3 + F(E1, i3) = 0.0; + F(E2, i3) = 0.0; + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(E1, i3 - j3) = F(E1, i3 + j3); + F(E2, i3 - j3) = F(E2, i3 + j3); + } + + // 2. normal component: wall is at i3-dhat/2 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(E0, i3 - j3) = -F(E0, i3 + j3 - dhat); + } + } + } + + if (lohi == Hi && grid.atBoundaryHi(p, d)) { + detail::set_upper_ghosts_to_nan(mflds, p, d, EX, true); + + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = grid.ldims[d]; + stop[d] = start[d] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + // 1. transverse components: wall is at i3, and there's one less ghost + F(E1, i3) = 0.0; + F(E2, i3) = 0.0; + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d] - 1; j3[d]++) { + F(E1, i3 + j3) = F(E1, i3 - j3); + F(E2, i3 + j3) = F(E2, i3 - j3); + } + + // 2. normal component: wall is at i3-dhat/2 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(E0, i3 + j3) = -F(E0, i3 - j3 + dhat); + } + } + } + } + } + + /** + * @brief Normal H at the wall's surface is a no-op, but set interior H such + * that reflecting particles feel the "right" magnetic force. + * + * Transverse H is flipped to ensure that a particle at the wall's surface + * with nominally-nonzero normal velocity—and thus, actually zero average + * normal velocity, since half of the cloud is reflected—experiences no + * magnetic force. + * @param mflds fields + */ + void apply_h_bcs(MfieldsState& mflds) override + { + const Grid_t& grid = mflds.grid(); + + const Int3 dhat = Int3::unit(d); + + const int H0 = HX + d; + const int H1 = HX + d.next(); + const int H2 = HX + d.prev(); + + for (int p = 0; p < mflds.n_patches(); p++) { + if (lohi == Lo && grid.atBoundaryLo(p, d)) { + detail::set_lower_ghosts_to_nan(mflds, p, d, HX, false); + + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = 0; + stop[d] = start[d] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + // 1. transverse components: wall is at i3-dhat/2 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(H1, i3 - j3) = -F(H1, i3 + j3 - dhat); + F(H2, i3 - j3) = -F(H2, i3 + j3 - dhat); + } + + // 2. normal component: wall is at i3 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(H0, i3 - j3) = F(H0, i3 + j3); + } + } + } + + if (lohi == Hi && grid.atBoundaryHi(p, d)) { + detail::set_upper_ghosts_to_nan(mflds, p, d, HX, false); + + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = grid.ldims[d]; + stop[d] = start[d] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + // 1. transverse components: wall is at i3-dhat/2 + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d]; j3[d]++) { + F(H1, i3 + j3) = -F(H1, i3 - j3 + dhat); + F(H2, i3 + j3) = -F(H2, i3 - j3 + dhat); + } + + // 2. normal component: wall is at i3, and there's one less ghost + for (Int3 j3 = dhat; j3[d] <= mflds.ibn()[d] - 1; j3[d]++) { + F(H0, i3 + j3) = F(H0, i3 - j3); + } + } + } + } + } + + Axis d; + LoHi lohi; +}; + +} // namespace field +} // namespace bnd +} // namespace psc diff --git a/src/libpsc/psc_bnd_fields/field_bc_base.hxx b/src/libpsc/psc_bnd_fields/field_bc_base.hxx new file mode 100644 index 0000000000..8f5d2e7e9a --- /dev/null +++ b/src/libpsc/psc_bnd_fields/field_bc_base.hxx @@ -0,0 +1,11 @@ +#pragma once + +template +struct FieldBcBase +{ + virtual ~FieldBcBase() {} + + virtual void apply_j_bcs(MfieldsState& mflds) = 0; + virtual void apply_e_bcs(MfieldsState& mflds) = 0; + virtual void apply_h_bcs(MfieldsState& mflds) = 0; +}; diff --git a/src/libpsc/psc_bnd_fields/field_bc_util.hxx b/src/libpsc/psc_bnd_fields/field_bc_util.hxx new file mode 100644 index 0000000000..c1701e260e --- /dev/null +++ b/src/libpsc/psc_bnd_fields/field_bc_util.hxx @@ -0,0 +1,136 @@ +#pragma once + +#include + +#include "kg/VecRange.hxx" + +namespace psc +{ +namespace bnd +{ +namespace field +{ +namespace detail +{ + +/** + * @brief Set E or B lower ghosts to the given constants (each component has + * its own constant). + * @param mflds mflds + * @param p patch index + * @param d which dimension to set the ghosts of + * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set + * @param val the constants + * @param include_edge whether or not values located on exact domain edges + * should be considered "ghosts" + */ +template +void set_lower_ghosts(MfieldsState& mflds, int p, int d, int mb, + typename MfieldsState::Real3 val, bool include_edge) +{ + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + stop[d] = 0; + + // TODO use gtensor views instead of VecRange + + for (int m = mb; m < mb + 3; m++) { + for (Int3 i3 : VecRange(start, stop)) { + F(m, i3) = val[m - mb]; + } + } + + if (!include_edge) { + return; + } + + Int3 edge_start = mflds.ib(); + Int3 edge_stop = mflds.ib() + mflds.im(); + edge_start[d] = 0; + edge_stop[d] = 1; + + for (int m = mb; m < mb + 3; m++) { + bool edge_ec = mb == EX && m - mb != d; + bool edge_fc = mb == HX && m - mb == d; + + if (edge_ec || edge_fc) { + for (Int3 i3 : VecRange(edge_start, edge_stop)) { + F(m, i3) = val[m - mb]; + } + } + } +} + +/** + * @brief Set E or B upper ghosts to the given constants (each component has + * its own constant). + * @param mflds mflds + * @param p patch index + * @param d which dimension to set the ghosts of + * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set + * @param val the constants + * @param include_edge whether or not values located on exact domain edges + * should be considered "ghosts" + */ +template +void set_upper_ghosts(MfieldsState& mflds, int p, int d, int mb, + typename MfieldsState::Real3 val, bool include_edge) +{ + auto F = make_Fields3d(mflds[p]); + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d] = mflds.grid().ldims[d] + 1; + + // TODO use gtensor views instead of VecRange + + for (int m = mb; m < mb + 3; m++) { + for (Int3 i3 : VecRange(start, stop)) { + F(m, i3) = val[m - mb]; + } + } + + Int3 edge_start = mflds.ib(); + Int3 edge_stop = mflds.ib() + mflds.im(); + edge_start[d] = mflds.grid().ldims[d]; + edge_stop[d] = mflds.grid().ldims[d] + 1; + + for (int m = mb; m < mb + 3; m++) { + bool not_edge_ec = mb == EX && m - mb == d; + bool not_edge_fc = mb == HX && m - mb != d; + + if (not_edge_ec || not_edge_fc || include_edge) + for (Int3 i3 : VecRange(edge_start, edge_stop)) { + F(m, i3) = val[m - mb]; + } + } +} + +template +void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, + bool include_edge) +{ +#ifndef DEBUG + return; +#endif + using real_t = typename MfieldsState::real_t; + real_t nan = std::numeric_limits::quiet_NaN(); + set_lower_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); +} + +template +void set_upper_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, + bool include_edge) +{ +#ifndef DEBUG + return; +#endif + using real_t = typename MfieldsState::real_t; + real_t nan = std::numeric_limits::quiet_NaN(); + set_upper_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); +} + +} // namespace detail +} // namespace field +} // namespace bnd +} // namespace psc \ No newline at end of file diff --git a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx index 7697893c75..bf226e9e41 100644 --- a/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx +++ b/src/libpsc/psc_bnd_fields/psc_bnd_fields_impl.hxx @@ -1,9 +1,10 @@ +#pragma once #include "psc.h" #include "kg/VecRange.hxx" #include "fields.hxx" #include "bnd_fields.hxx" -#include "radiating_bnd.hxx" +#include "field_bc_util.hxx" #include @@ -24,617 +25,17 @@ struct BndFields_ : BndFieldsBase // ---------------------------------------------------------------------- // fill_ghosts_E - void fill_ghosts_E(MfieldsState& mflds) - { - const auto& grid = mflds.grid(); - - for (int p = 0; p < mflds.n_patches(); p++) { - // lo - for (int d = 0; d < 3; d++) { - if (grid.atBoundaryLo(p, d)) { - switch (grid.bc.fld_lo[d]) { - case BND_FLD_PERIODIC: { - break; - } - case BND_FLD_CONDUCTING_WALL: { - conducting_wall_E_lo(mflds, p, d); - break; - } - case BND_FLD_OPEN: { - set_lower_ghosts(mflds, p, d, EX, background_e, false); - break; - } - default: { - assert(0); - } - } - } - } - - // hi - for (int d = 0; d < 3; d++) { - if (grid.atBoundaryHi(p, d)) { - switch (grid.bc.fld_hi[d]) { - case BND_FLD_PERIODIC: { - break; - } - case BND_FLD_CONDUCTING_WALL: { - conducting_wall_E_hi(mflds, p, d); - break; - } - case BND_FLD_OPEN: { - set_upper_ghosts(mflds, p, d, EX, background_e, false); - break; - } - default: { - assert(0); - } - } - } - } - } - } + void fill_ghosts_E(MfieldsState& mflds) {} // ---------------------------------------------------------------------- // fill_ghosts_H - void fill_ghosts_H(MfieldsState& mflds) - { - const auto& grid = mflds.grid(); - - for (int p = 0; p < mflds.n_patches(); p++) { - // lo - for (int d = 0; d < 3; d++) { - if (grid.bc.fld_lo[d] == BND_FLD_OPEN && radiation) { - radiation->update_cache_lower(grid.time(), d); - } - - if (grid.atBoundaryLo(p, d)) { - switch (grid.bc.fld_lo[d]) { - case BND_FLD_PERIODIC: { - break; - } - case BND_FLD_CONDUCTING_WALL: { - conducting_wall_H_lo(mflds, p, d); - break; - } - case BND_FLD_OPEN: { - radiative_H_lo(mflds, p, d); - break; - } - default: { - assert(0); - } - } - } - } - // hi - for (int d = 0; d < 3; d++) { - if (grid.bc.fld_hi[d] == BND_FLD_OPEN && radiation) { - radiation->update_cache_upper(grid.time(), d); - } - - if (grid.atBoundaryHi(p, d)) { - switch (grid.bc.fld_hi[d]) { - case BND_FLD_PERIODIC: { - break; - } - case BND_FLD_CONDUCTING_WALL: { - conducting_wall_H_hi(mflds, p, d); - break; - } - case BND_FLD_OPEN: { - radiative_H_hi(mflds, p, d); - break; - } - default: { - assert(0); - } - } - } - } - } - } + void fill_ghosts_H(MfieldsState& mflds) {} // ---------------------------------------------------------------------- // add_ghosts_J - void add_ghosts_J(MfieldsState& mflds) - { - const auto& grid = mflds.grid(); - - for (int p = 0; p < mflds.n_patches(); p++) { - // lo - for (int d = 0; d < 3; d++) { - if (grid.atBoundaryLo(p, d)) { - switch (grid.bc.fld_lo[d]) { - case BND_FLD_PERIODIC: { - break; - } - case BND_FLD_CONDUCTING_WALL: { - conducting_wall_J_lo(mflds, p, d); - break; - } - case BND_FLD_OPEN: { - break; - } - default: { - assert(0); - } - } - } - } - // hi - for (int d = 0; d < 3; d++) { - if (grid.atBoundaryHi(p, d)) { - switch (grid.bc.fld_hi[d]) { - case BND_FLD_PERIODIC: { - break; - } - case BND_FLD_CONDUCTING_WALL: { - conducting_wall_J_hi(mflds, p, d); - break; - } - case BND_FLD_OPEN: { - break; - } - default: { - assert(0); - } - } - } - } - } - } - - static void set_lower_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, - bool include_edge) - { -#ifndef DEBUG - return; -#endif - real_t nan = std::numeric_limits::quiet_NaN(); - set_lower_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); - } - - static void set_upper_ghosts_to_nan(MfieldsState& mflds, int p, int d, int mb, - bool include_edge) - { -#ifndef DEBUG - return; -#endif - real_t nan = std::numeric_limits::quiet_NaN(); - set_upper_ghosts(mflds, p, d, mb, {nan, nan, nan}, include_edge); - } - - /** - * @brief Set E or B lower ghosts to the given constants (each component has - * its own constant). - * @param mflds mflds - * @param p patch index - * @param d which dimension to set the ghosts of - * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set - * @param val the constants - * @param include_edge whether or not values located on exact domain edges - * should be considered "ghosts" - */ - static void set_lower_ghosts(MfieldsState& mflds, int p, int d, int mb, - Real3 val, bool include_edge) - { - auto F = make_Fields3d(mflds[p]); - Int3 start = mflds.ib(); - Int3 stop = mflds.ib() + mflds.im(); - stop[d] = 0; - - // TODO use gtensor views instead of VecRange - - for (int m = mb; m < mb + 3; m++) { - for (Int3 i3 : VecRange(start, stop)) { - F(m, i3) = val[m - mb]; - } - } - - if (!include_edge) { - return; - } - - Int3 edge_start = mflds.ib(); - Int3 edge_stop = mflds.ib() + mflds.im(); - edge_start[d] = 0; - edge_stop[d] = 1; - - for (int m = mb; m < mb + 3; m++) { - bool edge_ec = mb == EX && m - mb != d; - bool edge_fc = mb == HX && m - mb == d; - - if (edge_ec || edge_fc) { - for (Int3 i3 : VecRange(edge_start, edge_stop)) { - F(m, i3) = val[m - mb]; - } - } - } - } - - /** - * @brief Set E or B upper ghosts to the given constants (each component has - * its own constant). - * @param mflds mflds - * @param p patch index - * @param d which dimension to set the ghosts of - * @param mb `EX` or `HX`; note that `mb+1` and `mb+2` are also set - * @param val the constants - * @param include_edge whether or not values located on exact domain edges - * should be considered "ghosts" - */ - static void set_upper_ghosts(MfieldsState& mflds, int p, int d, int mb, - Real3 val, bool include_edge) - { - auto F = make_Fields3d(mflds[p]); - Int3 start = mflds.ib(); - Int3 stop = mflds.ib() + mflds.im(); - start[d] = mflds.grid().ldims[d] + 1; - - // TODO use gtensor views instead of VecRange - - for (int m = mb; m < mb + 3; m++) { - for (Int3 i3 : VecRange(start, stop)) { - F(m, i3) = val[m - mb]; - } - } - - Int3 edge_start = mflds.ib(); - Int3 edge_stop = mflds.ib() + mflds.im(); - edge_start[d] = mflds.grid().ldims[d]; - edge_stop[d] = mflds.grid().ldims[d] + 1; - - for (int m = mb; m < mb + 3; m++) { - bool not_edge_ec = mb == EX && m - mb == d; - bool not_edge_fc = mb == HX && m - mb != d; - - if (not_edge_ec || not_edge_fc || include_edge) - for (Int3 i3 : VecRange(edge_start, edge_stop)) { - F(m, i3) = val[m - mb]; - } - } - } - - void conducting_wall_E_lo(MfieldsState& mflds, int p, int d) - { - set_lower_ghosts_to_nan(mflds, p, d, EX, true); - - auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); - - if (d == 1) { - for (int iz = -2; iz < ldims[2] + 2; iz++) { - // FIXME, needs to be for other dir, too, and it's ugly - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(EX, ix, 0, iz) = 0.; - F(EX, ix, -1, iz) = F(EX, ix, 1, iz); - - F(EY, ix, -1, iz) = -F(EY, ix, 0, iz); - - F(EZ, ix, 0, iz) = 0.; - F(EZ, ix, -1, iz) = F(EZ, ix, 1, iz); - } - } - } else if (d == 2) { - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(EX, ix, iy, 0) = 0.; - F(EX, ix, iy, -1) = F(EX, ix, iy, 1); - - F(EY, ix, iy, 0) = 0.; - F(EY, ix, iy, -1) = F(EY, ix, iy, 1); - - F(EZ, ix, iy, -1) = -F(EZ, ix, iy, 0); - } - } - } else { - assert(0); - } - } - - void conducting_wall_E_hi(MfieldsState& mflds, int p, int d) - { - set_upper_ghosts_to_nan(mflds, p, d, EX, true); - - auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); - - if (d == 1) { - int my _mrc_unused = ldims[1]; - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(EX, ix, my, iz) = 0.; - F(EX, ix, my + 1, iz) = F(EX, ix, my - 1, iz); - - F(EY, ix, my, iz) = -F(EY, ix, my - 1, iz); - - F(EZ, ix, my, iz) = 0.; - F(EZ, ix, my + 1, iz) = F(EZ, ix, my - 1, iz); - } - } - } else if (d == 2) { - int mz = ldims[2]; - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(EX, ix, iy, mz) = 0.; - F(EX, ix, iy, mz + 1) = F(EX, ix, iy, mz - 1); - - F(EY, ix, iy, mz) = 0.; - F(EY, ix, iy, mz + 1) = F(EY, ix, iy, mz - 1); - - F(EZ, ix, iy, mz) = -F(EZ, ix, iy, mz - 1); - } - } - } else { - assert(0); - } - } - - void conducting_wall_H_lo(MfieldsState& mflds, int p, int d) - { - set_lower_ghosts_to_nan(mflds, p, d, HX, false); - - auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); - - if (d == 1) { - for (int iz = -1; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HX, ix, -1, iz) = -F(HX, ix, 0, iz); - - F(HY, ix, -1, iz) = F(HY, ix, 1, iz); - - F(HZ, ix, -1, iz) = -F(HZ, ix, 0, iz); - } - } - } else if (d == 2) { - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HX, ix, iy, -1) = -F(HX, ix, iy, 0); - - F(HY, ix, iy, -1) = -F(HY, ix, iy, 0); - - F(HZ, ix, iy, -1) = F(HZ, ix, iy, 1); - } - } - } else { - assert(0); - } - } - - void conducting_wall_H_hi(MfieldsState& mflds, int p, int d) - { - set_upper_ghosts_to_nan(mflds, p, d, HX, false); - - auto F = make_Fields3d(mflds[p]); - - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); - - if (d == 1) { - int my _mrc_unused = ldims[1]; - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HX, ix, my, iz) = -F(HX, ix, my - 1, iz); - - F(HY, ix, my + 1, iz) = F(HY, ix, my - 1, iz); - - F(HZ, ix, my, iz) = -F(HZ, ix, my - 1, iz); - } - } - } else if (d == 2) { - int mz = ldims[2]; - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(HX, ix, iy, mz) = -F(HX, ix, iy, mz - 1); - - F(HY, ix, iy, mz) = -F(HY, ix, iy, mz - 1); - - F(HZ, ix, iy, mz + 1) = F(HZ, ix, iy, mz - 1); - } - } - } else { - assert(0); - } - } - - void conducting_wall_J_lo(MfieldsState& mflds, int p, int d) - { - auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); - - if (d == 1) { - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(JXI, ix, 1, iz) += F(JXI, ix, -1, iz); - F(JXI, ix, -1, iz) = 0.; - - F(JYI, ix, 0, iz) -= F(JYI, ix, -1, iz); - F(JYI, ix, -1, iz) = 0.; - - F(JZI, ix, 1, iz) += F(JZI, ix, -1, iz); - F(JZI, ix, -1, iz) = 0.; - } - } - } else if (d == 2) { - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(JXI, ix, iy, 1) += F(JXI, ix, iy, -1); - F(JXI, ix, iy, -1) = 0.; - - F(JYI, ix, iy, 1) += F(JYI, ix, iy, -1); - F(JYI, ix, iy, -1) = 0.; - - F(JZI, ix, iy, 0) -= F(JZI, ix, iy, -1); - F(JZI, ix, iy, -1) = 0.; - } - } - } else { - assert(0); - } - } - - void conducting_wall_J_hi(MfieldsState& mflds, int p, int d) - { - auto F = make_Fields3d(mflds[p]); - const int* ldims = mflds.grid().ldims; - Int3 ib = mflds.ib(), im = mflds.im(); - - if (d == 1) { - int my _mrc_unused = ldims[1]; - for (int iz = -2; iz < ldims[2] + 2; iz++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(JXI, ix, my - 1, iz) += F(JXI, ix, my + 1, iz); - F(JXI, ix, my + 1, iz) = 0.; - - F(JYI, ix, my - 1, iz) -= F(JYI, ix, my, iz); - F(JYI, ix, my, iz) = 0.; - - F(JZI, ix, my - 1, iz) += F(JZI, ix, my + 1, iz); - F(JZI, ix, my + 1, iz) = 0.; - } - } - } else if (d == 2) { - int mz = ldims[2]; - for (int iy = -2; iy < ldims[1] + 2; iy++) { - for (int ix = std::max(-2, ib[0]); - ix < std::min(ldims[0] + 2, ib[0] + im[0]); ix++) { - F(JXI, ix, iy, mz - 1) += F(JXI, ix, iy, mz + 1); - F(JXI, ix, iy, mz + 1) = 0.; - - F(JYI, ix, iy, mz - 1) += F(JYI, ix, iy, mz + 1); - F(JYI, ix, iy, mz + 1) = 0.; - - F(JZI, ix, iy, mz - 1) -= F(JZI, ix, iy, mz); - F(JZI, ix, iy, mz) = 0.; - } - } - } else { - assert(0); - } - } - - void radiative_H_lo(MfieldsState& mflds, int p, int d) - { - set_lower_ghosts_to_nan(mflds, p, d, HX, false); - - auto F = make_Fields3d(mflds[p]); - const Grid_t& grid = mflds.grid(); - real_t dt = grid.dt; - Real3 dtdx = dt * Real3(grid.domain.dx_inv); - - int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3; - int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; - int E1 = EX + d1, E2 = EX + d2; - int J1 = JXI + d1, J2 = JXI + d2; - - Int3 start = mflds.ib(); - Int3 stop = mflds.ib() + mflds.im(); - start[d0] = -1; - stop[d0] = 0; - - for (Int3 i3 : VecRange(start, stop)) { - Int3 edge_idx = i3 + Int3::unit(d0); - - real_t s = 0.0; - real_t p = 0.0; - if (radiation) { - Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * - Real3(grid.domain.dx); - Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * - Real3(grid.domain.dx); - s = radiation->pulse_s_lower(grid.time(), d0, p, x3_s); - p = radiation->pulse_p_lower(grid.time(), d0, p, x3_p); - } - - F(H2, i3) = - (4.f * s - 2.f * (F(E1, edge_idx) - background_e[d1]) - - dtdx[d2] * (F(H0, edge_idx) - F(H0, edge_idx - Int3::unit(d2))) - - (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) + - dt * F(J1, edge_idx)) / - (1.f + dtdx[d0]) + - background_h[d2]; - F(H1, i3) = - (-4.f * p + 2.f * (F(E2, edge_idx) - background_e[d2]) - - dtdx[d1] * (F(H0, edge_idx) - F(H0, edge_idx - Int3::unit(d1))) - - (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) + - dt * F(J2, edge_idx)) / - (1.f + dtdx[d0]) + - background_h[d1]; - } - } - - void radiative_H_hi(MfieldsState& mflds, int p, int d) - { - set_upper_ghosts_to_nan(mflds, p, d, HX, false); - - auto F = make_Fields3d(mflds[p]); - const Grid_t& grid = mflds.grid(); - Int3 ldims = grid.ldims; - real_t dt = grid.dt; - Real3 dtdx = dt * Real3(grid.domain.dx_inv); - - int d0 = (d + 0) % 3, d1 = (d + 1) % 3, d2 = (d + 2) % 3; - int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; - int E1 = EX + d1, E2 = EX + d2; - int J1 = JXI + d1, J2 = JXI + d2; - - Int3 start = mflds.ib(); - Int3 stop = mflds.ib() + mflds.im(); - start[d0] = grid.ldims[d0]; - stop[d0] = grid.ldims[d0] + 1; - - for (Int3 i3 : VecRange(start, stop)) { - Int3 edge_idx = i3 - Int3::unit(d0); - - real_t s = 0.0; - real_t p = 0.0; - if (radiation) { - Real3 x3_s = (Real3(edge_idx) + Real3::unit(d1) * real_t(0.5)) * - Real3(grid.domain.dx); - Real3 x3_p = (Real3(edge_idx) + Real3::unit(d2) * real_t(0.5)) * - Real3(grid.domain.dx); - s = radiation->pulse_s_upper(grid.time(), d0, p, x3_s); - p = radiation->pulse_p_upper(grid.time(), d0, p, x3_p); - } - - F(H2, i3) = (-4.f * s + 2.f * (F(E1, i3) - background_e[d1]) + - dtdx[d2] * (F(H0, i3) - F(H0, i3 - Int3::unit(d2))) - - (1.f - dtdx[d0]) * (F(H2, edge_idx) - background_h[d2]) - - dt * F(J1, i3)) / - (1.f + dtdx[d0]) + - background_h[d2]; - F(H1, i3) = (4.f * p - 2.f * (F(E2, i3) - background_e[d2]) + - dtdx[d1] * (F(H0, i3) - F(H0, i3 - Int3::unit(d1))) - - (1.f - dtdx[d0]) * (F(H1, edge_idx) - background_h[d1]) - - dt * F(J2, i3)) / - (1.f + dtdx[d0]) + - background_h[d1]; - } - } - - Vec3 background_e = {0.0, 0.0, 0.0}; - Vec3 background_h = {0.0, 0.0, 0.0}; - - RadiatingBoundary* radiation = nullptr; + void add_ghosts_J(MfieldsState& mflds) {} }; // ====================================================================== diff --git a/src/libpsc/psc_bnd_fields/radiating.hxx b/src/libpsc/psc_bnd_fields/radiating.hxx new file mode 100644 index 0000000000..f61eabc656 --- /dev/null +++ b/src/libpsc/psc_bnd_fields/radiating.hxx @@ -0,0 +1,261 @@ +#pragma once + +#include "psc.h" +#include "../axis.hxx" +#include "kg/Vec3.h" +#include "kg/VecRange.hxx" +#include "field_bc_base.hxx" +#include "field_bc_util.hxx" +#include "../psc_bnd/psc_bnd_util.hxx" + +namespace psc +{ +namespace bnd +{ +namespace field +{ + +/** + * @brief "Radiating" open boundary. Prescribe an arbitrary inflowing + * electromagnetic pulse. See Ruhl 2006 for details. + * @tparam Dim dimension type + * @tparam MfieldsState fields type + * @tparam P pulse type (i.e. a type that implements `PulseBase`) + */ +template +struct Radiating : FieldBcBase +{ + using dim_t = Dim; + using Pulse = P; + using real_t = typename MfieldsState::real_t; + using Real3 = typename MfieldsState::Real3; + + Radiating(Pulse pulse, Axis d, LoHi lohi) : pulse{pulse}, d{d}, lohi{lohi} {} + + /** + * @brief Don't do anything to the current. + * @param mflds fields + */ + void apply_j_bcs(MfieldsState& mflds) override {} + + /** + * @brief Set the normal E to 0. They could be self-consistently evolved + * instead, but ghost corners aren't handled yet. Transverse components are + * deep enough to not affect 1st-order particles. + * @param mflds fields + */ + void apply_e_bcs(MfieldsState& mflds) override + { + const Grid_t& grid = mflds.grid(); + pulse.tick(grid.time()); + + int d0 = d, d1 = d.next(), d2 = d.prev(); + int E0 = EX + d0, E1 = EX + d1, E2 = EX + d2; + + Int3 d0hat = Int3::unit(d0); + Int3 d1hat = Int3::unit(d1); + Int3 d2hat = Int3::unit(d2); + + for (int p = 0; p < mflds.n_patches(); p++) { + if (lohi == Lo && grid.atBoundaryLo(p, d)) { + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d0] = -1; + stop[d0] = start[d0] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + Real3 x_e1 = + (Real3(i3) + Real3(d1hat) * real_t(0.5)) * Real3(grid.domain.dx); + Real3 x_e2 = + (Real3(i3) + Real3(d2hat) * real_t(0.5)) * Real3(grid.domain.dx); + + F(E0, i3 - d0hat) = 0.0; + F(E1, i3) = pulse.sample_exterior_field(E1, grid.time(), p, x_e1); + F(E2, i3) = pulse.sample_exterior_field(E2, grid.time(), p, x_e2); + } + } + if (lohi == Hi && grid.atBoundaryHi(p, d)) { + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d0] = grid.ldims[d0] + 1; + stop[d0] = start[d0] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + Real3 x_e1 = + (Real3(i3) + Real3(d1hat) * real_t(0.5)) * Real3(grid.domain.dx); + Real3 x_e2 = + (Real3(i3) + Real3(d2hat) * real_t(0.5)) * Real3(grid.domain.dx); + + F(E0, i3) = 0.0; + F(E1, i3) = pulse.sample_exterior_field(E1, grid.time(), p, x_e1); + F(E2, i3) = pulse.sample_exterior_field(E2, grid.time(), p, x_e2); + } + } + } + } + + /** + * @brief Set the first layer of transverse H ghosts such that the inflowing S + * and P waves are prescribed at the boundary. The definitions of S and P + * differ from Ruhl 2006 by a factor of 2, and Ruhl's definitions aren't + * invariant under an x->y-z->x rotation. + * @param mflds fields + */ + void apply_h_bcs(MfieldsState& mflds) override + { + const Grid_t& grid = mflds.grid(); + real_t dt = grid.dt; + Real3 dtdx = dt * Real3(grid.domain.dx_inv); + + pulse.tick(grid.time()); + + int d0 = d, d1 = d.next(), d2 = d.prev(); + int H0 = HX + d0, H1 = HX + d1, H2 = HX + d2; + int E1 = EX + d1, E2 = EX + d2; + int J1 = JXI + d1, J2 = JXI + d2; + + Int3 d0hat = Int3::unit(d0); + Int3 d1hat = Int3::unit(d1); + Int3 d2hat = Int3::unit(d2); + + for (int p = 0; p < mflds.n_patches(); p++) { + if (lohi == Lo && grid.atBoundaryLo(p, d)) { + psc::bnd::field::detail::set_lower_ghosts_to_nan(mflds, p, d, HX, + false); + + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d0] = 0; + stop[d0] = start[d0] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + Real3 x3_s = + (Real3(i3) + Real3(d1hat) * real_t(0.5)) * Real3(grid.domain.dx); + Real3 x3_p = + (Real3(i3) + Real3(d2hat) * real_t(0.5)) * Real3(grid.domain.dx); + + real_t s = pulse.sample_exterior_field(E1, grid.time(), p, x3_s) + + pulse.sample_exterior_field(H2, grid.time(), p, x3_s); + real_t p = pulse.sample_exterior_field(E2, grid.time(), p, x3_p) - + pulse.sample_exterior_field(H1, grid.time(), p, x3_p); + + F(H2, i3 - d0hat) = (2.f * s - 2.f * F(E1, i3) - + dtdx[d2] * (F(H0, i3) - F(H0, i3 - d2hat)) - + (1.f - dtdx[d0]) * F(H2, i3) + dt * F(J1, i3)) / + (1.f + dtdx[d0]); + F(H1, i3 - d0hat) = (-2.f * p + 2.f * F(E2, i3) - + dtdx[d1] * (F(H0, i3) - F(H0, i3 - d1hat)) - + (1.f - dtdx[d0]) * F(H1, i3) - dt * F(J2, i3)) / + (1.f + dtdx[d0]); + } + } + + if (lohi == Hi && grid.atBoundaryHi(p, d)) { + psc::bnd::field::detail::set_upper_ghosts_to_nan(mflds, p, d, HX, + false); + + auto F = make_Fields3d(mflds[p]); + + Int3 start = mflds.ib(); + Int3 stop = mflds.ib() + mflds.im(); + start[d0] = grid.ldims[d0]; + stop[d0] = start[d0] + 1; + + for (Int3 i3 : VecRange(start, stop)) { + Real3 x3_s = + (Real3(i3) + Real3(d1hat) * real_t(0.5)) * Real3(grid.domain.dx); + Real3 x3_p = + (Real3(i3) + Real3(d2hat) * real_t(0.5)) * Real3(grid.domain.dx); + + real_t s = pulse.sample_exterior_field(E1, grid.time(), p, x3_s) - + pulse.sample_exterior_field(H2, grid.time(), p, x3_s); + real_t p = pulse.sample_exterior_field(E2, grid.time(), p, x3_p) + + pulse.sample_exterior_field(H1, grid.time(), p, x3_p); + + F(H2, i3) = (-2.f * s + 2.f * F(E1, i3) + + dtdx[d2] * (F(H0, i3) - F(H0, i3 - d2hat)) - + (1.f - dtdx[d0]) * F(H2, i3 - d0hat) - dt * F(J1, i3)) / + (1.f + dtdx[d0]); + F(H1, i3) = (2.f * p - 2.f * F(E2, i3) + + dtdx[d1] * (F(H0, i3) - F(H0, i3 - d1hat)) - + (1.f - dtdx[d0]) * F(H1, i3 - d0hat) + dt * F(J2, i3)) / + (1.f + dtdx[d0]); + } + } + } + } + + Axis d; + LoHi lohi; + +private: + Pulse pulse; +}; + +/** + * @brief The archetypal Pulse type used by the `Radiating` boundary condition. + * This type isn't used polymorphically, so extending it isn't strictly + * required. + * @tparam real_t real type + */ +template +struct PulseBase +{ + using Real3 = Vec3; + + /** + * @brief Sample a component of an out-of-domain field at a location and time. + * The field value will be used as part of a boundary condition calculation. + * @param m field component (e.g. `EX`) + * @param t time + * @param p patch index + * @param x3 cell-normalized location within the patch + * @return the field value + */ + virtual real_t sample_exterior_field(int m, double t, int p, Real3 x3) = 0; + + /** + * @brief Perform any operations that occur once per time step. + * @param t time + */ + virtual void tick(double t) {} +}; + +/** + * @brief A constant pulse. Use this for open boundaries that have constant + * external fields. + * @tparam real_t + */ +template +struct ConstantPulse : PulseBase +{ + using Real3 = Vec3; + + ConstantPulse(Real3 e, Real3 h) : e{e}, h{h} {} + + real_t sample_exterior_field(int m, double t, int p, Real3 x3) override + { + switch (m) { + case EX: return e[0]; + case EY: return e[1]; + case EZ: return e[2]; + case HX: return h[0]; + case HY: return h[1]; + case HZ: return h[2]; + default: return 0.0; + } + } + + Real3 e; + Real3 h; +}; + +} // namespace field +} // namespace bnd +} // namespace psc diff --git a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx b/src/libpsc/psc_bnd_fields/radiating_bnd.hxx deleted file mode 100644 index c6a40de9f0..0000000000 --- a/src/libpsc/psc_bnd_fields/radiating_bnd.hxx +++ /dev/null @@ -1,17 +0,0 @@ -#include "kg/Vec3.h" - -template -struct RadiatingBoundary -{ - using Real3 = Vec3; - - virtual real_t pulse_s_lower(double t, int d, int p, Real3 x3) = 0; - virtual real_t pulse_p_lower(double t, int d, int p, Real3 x3) = 0; - - virtual real_t pulse_s_upper(double t, int d, int p, Real3 x3) = 0; - virtual real_t pulse_p_upper(double t, int d, int p, Real3 x3) = 0; - - // FIXME these are a hack for a specific subclass to work - virtual void update_cache_lower(double t, int d) {} - virtual void update_cache_upper(double t, int d) {} -}; \ No newline at end of file diff --git a/src/libpsc/psc_particle_injectors/boundary_injector.hxx b/src/libpsc/psc_particle_injectors/boundary_injector.hxx new file mode 100644 index 0000000000..84f4ee32ed --- /dev/null +++ b/src/libpsc/psc_particle_injectors/boundary_injector.hxx @@ -0,0 +1,198 @@ +#pragma once + +#include + +#include "grid.hxx" +#include "rng.hxx" +#include "particle.h" +#include +#include "pushp.hxx" +#include "dim.hxx" +#include "setup_particles.hxx" +#include "kg/VecRange.hxx" +#include "../psc_push_particles/inc_push.cxx" +#include "../psc_bnd/psc_bnd_util.hxx" +#include "injector_base.hxx" + +using psc::bnd::LoHi; + +/// @brief A particle generator for use with @ref BoundaryInjector. Samples +/// particles from a (possibly shifted) Maxwellian distribution. +class ParticleGeneratorMaxwellian +{ +public: + using Real = psc::particle::Inject::Real; + using Real3 = psc::particle::Inject::Real3; + + // FIXME would be nice to just pass 1 thing for kind-related info + ParticleGeneratorMaxwellian(int kind_idx, Grid_t::Kind kind, Real3 mean_v, + Real3 temperature) + : kind_idx{kind_idx}, prt_booster{-mean_v} + { + for (int d = 0; d < 3; d++) { + vdfs[d] = VDF{0.0, sqrt(temperature[d] / kind.m)}; + } + } + + psc::particle::Inject get(Real3 min_pos, Real3 pos_range) + { + Real3 x; + for (int d = 0; d < 3; d++) { + x[d] = min_pos[d] + uniform_dist.get() * pos_range[d]; + } + + Real3 v{vdfs[0].get(), vdfs[1].get(), vdfs[2].get()}; + // FIXME should really sample from Maxwell-Juttner + // this hack interprests v as u to handle rare case when v>1 + // v<<1 => v~= u anyways + Real3 u = prt_booster.boost(v); + + Real w = 1.0; + psc::particle::Tag tag = 0; + + return {x, u, w, kind_idx, tag}; + } + +private: + using VDF = rng::Normal; + Vec3 vdfs; + VelocityBooster prt_booster; + int kind_idx; + rng::Uniform uniform_dist{0.0, 1.0}; +}; + +/// @brief Injects particles on a given boundary, sampling from a given particle +/// generator. For precise control over multiple particle species, use one +/// BoundaryInjector per species. +/// @tparam LOHI whether to inject at the lower or upper boundary +/// @tparam PARTICLE_GENERATOR a type that defines `get(min_pos, pos_range)` and +/// returns an injectable particle within that range of positions (usually a +/// grid cell); see @ref ParticleGeneratorMaxwellian +/// @tparam PUSH_PARTICLES type that provides the types `Mparticles`, +/// `MfieldsState`, `Current`, `real_t`, etc. +template +class BoundaryInjector + : public InjectorBase +{ + static const int INJECT_DIM_IDX_ = 1; + +public: + using ParticleGenerator = PARTICLE_GENERATOR; + using PushParticles = PUSH_PARTICLES; + + using Mparticles = typename PushParticles::Mparticles; + using MfieldsState = typename PushParticles::MfieldsState; + using Current = typename PushParticles::Current; + using real_t = typename PushParticles::real_t; + using InterpolateEM_t = typename PushParticles::InterpolateEM_t; + using Real3 = Vec3; + + static const bool lo = LOHI == LoHi::Lo; + + BoundaryInjector(ParticleGenerator particle_generator, real_t density = 1.0) + : particle_generator_{particle_generator}, density{density} + {} + + /// Injects particles at specified y-bounds as if there were a population of + /// particles just beyond the edge. The imaginary particle population is + /// sampled using the given ParticleGenerator. + /// + /// The dimensional limitations may be removed in the future. + void inject(Mparticles& mprts, MfieldsState& mflds) override + { + static_assert(INJECT_DIM_IDX_ == 1, + "only injection at lower bound of y is supported"); + + const Grid_t& grid = mprts.grid(); + auto injectors_by_patch = mprts.injector(); + + Real3 dxi = grid.domain.dx_inv; + Current current(grid); + + bool preaccelerate = true; + real_t npp = 0.5; // number of plasma periods + + real_t plasma_freq_sq = 0.0; + for (Grid_t::Kind kind : grid.kinds) { + plasma_freq_sq += density / kind.m; + } + real_t plasma_period = 2.f * M_PI / sqrt(plasma_freq_sq); + real_t t_accel = npp * plasma_period; + + for (int p = 0; p < grid.n_patches(); p++) { + if (lo ? grid.atBoundaryLo(p, INJECT_DIM_IDX_) + : grid.atBoundaryHi(p, INJECT_DIM_IDX_)) { + auto&& injector = injectors_by_patch[p]; + auto flds = mflds[p]; + typename Current::fields_t J(flds); + typename InterpolateEM_t::fields_t EM(flds.storage(), flds.ib()); + InterpolateEM_t ip; + AdvanceParticle advance{grid.dt}; + + Int3 start = Int3{0, 0, 0}.with_component( + INJECT_DIM_IDX_, lo ? -1 : grid.ldims[INJECT_DIM_IDX_]); + Int3 stop = grid.ldims.with_component(INJECT_DIM_IDX_, + start[INJECT_DIM_IDX_] + 1); + + for (Int3 initial_idx : VecRange(start, stop)) { + Real3 cell_corner = Real3(initial_idx) * grid.domain.dx; + int n_prts_to_try_inject = + get_n_in_cell(density, grid.norm.prts_per_unit_density, true); + + for (int prt_count = 0; prt_count < n_prts_to_try_inject; + prt_count++) { + // sample position uniformly from first ghost layer, and velocity + // from vdf at x=infty + psc::particle::Inject prt = + particle_generator_.get(cell_corner, grid.domain.dx); + + real_t m = grid.kinds[prt.kind].m; + real_t q = grid.kinds[prt.kind].q; + + Real3 initial_normalized_pos = prt.x * dxi; + + if (preaccelerate) { + real_t E_interp; + + ip.set_coeffs(initial_normalized_pos.with_component( + INJECT_DIM_IDX_, start[INJECT_DIM_IDX_] + (lo ? 1 : 0))); + switch (INJECT_DIM_IDX_) { + case 0: E_interp = ip.ex(EM); break; + case 1: E_interp = ip.ey(EM); break; + case 2: E_interp = ip.ez(EM); break; + default: assert(false); + } + + prt.u[INJECT_DIM_IDX_] += t_accel * q * E_interp / m; + } + + // push normal x + Real3 v = advance.calc_v(prt.u); + advance.push_x(prt.x, v); + + Real3 final_normalized_pos = prt.x * dxi; + Int3 final_idx = final_normalized_pos.fint(); + + if (lo ? final_idx[INJECT_DIM_IDX_] <= start[INJECT_DIM_IDX_] + : final_idx[INJECT_DIM_IDX_] >= start[INJECT_DIM_IDX_]) { + // don't inject a particle that fails to enter the patch + continue; + } + + injector.inject_local(prt); + + current.calc_j(J, initial_normalized_pos, final_normalized_pos, + final_idx, initial_idx, q * prt.w, v); + } + } + } + } + } + +public: + real_t density; + +private: + ParticleGenerator particle_generator_; +}; diff --git a/src/include/injector_base.hxx b/src/libpsc/psc_particle_injectors/injector_base.hxx similarity index 100% rename from src/include/injector_base.hxx rename to src/libpsc/psc_particle_injectors/injector_base.hxx diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx index 4f5c005b6b..59e2176653 100644 --- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx +++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx @@ -26,7 +26,7 @@ struct PushParticlesVb static void push_mprts(Mparticles& mprts, MfieldsState& mflds) { - const auto& grid = mprts.grid(); + const Grid_t& grid = mprts.grid(); Real3 dxi = Real3(grid.domain.dx).inv(); real_t dq_kind[MAX_NR_KINDS]; auto& kinds = grid.kinds; diff --git a/src/libpsc/tests/test_boundary_injector.cxx b/src/libpsc/tests/test_boundary_injector.cxx index d8a679f953..3398df8362 100644 --- a/src/libpsc/tests/test_boundary_injector.cxx +++ b/src/libpsc/tests/test_boundary_injector.cxx @@ -2,29 +2,29 @@ #include "test_common.hxx" -#include "boundary_injector.hxx" - #include "psc.hxx" #include "output_fields.hxx" #include "../psc_config.hxx" +#include "../psc_particle_injectors/boundary_injector.hxx" TEST(BoundaryInjectorTest, ParticleGeneratorMaxwellianTest) { int kind_idx = 15; Grid_t::Kind kind{1.0, 1836.0, "ion"}; ParticleGeneratorMaxwellian::Real w = 1.0; - ParticleGeneratorMaxwellian::Real3 mean_u{0.0, 5.0, 15.0}; + ParticleGeneratorMaxwellian::Real3 mean_v{0.0, .05, .15}; ParticleGeneratorMaxwellian::Real3 temperature{0.0, 0.0, 0.0}; ParticleGeneratorMaxwellian::Real3 pos{1.0, 2.0, 5.0}; - ParticleGeneratorMaxwellian gen{kind_idx, kind, mean_u, temperature}; + ParticleGeneratorMaxwellian gen{kind_idx, kind, mean_v, temperature}; auto prt = gen.get(pos, {0.0, 0.0, 0.0}); ASSERT_EQ(prt.kind, kind_idx); ASSERT_EQ(prt.w, w); ASSERT_EQ(prt.tag, 0); - ASSERT_EQ(prt.u, mean_u); // zero temperature => exact velocity + ASSERT_EQ(prt.u / sqrt(1 + prt.u.mag2()), + mean_v); // zero temperature => exact velocity ASSERT_EQ(prt.x, pos); } @@ -130,9 +130,9 @@ TEST(BoundaryInjectorTest, Integration1Particle) auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, balance, collision, checks); - psc.add_injector( - new BoundaryInjector( - ParticleGenerator(1, 1), grid)); + psc.add_injector(new BoundaryInjector( + ParticleGenerator(1, 1))); // ---------------------------------------------------------------------- // set up initial conditions @@ -188,8 +188,8 @@ TEST(BoundaryInjectorTest, IntegrationManyParticles) balance, collision, checks); psc.add_injector( - new BoundaryInjector( - ParticleGenerator(-1, 1), grid)); + new BoundaryInjector( + ParticleGenerator(-1, 1))); // ---------------------------------------------------------------------- // set up initial conditions @@ -242,11 +242,11 @@ TEST(BoundaryInjectorTest, IntegrationManySpecies) Collision collision{grid, 0, 0.1}; auto inject_electrons = - BoundaryInjector{ - ParticleGenerator(-1, 0), grid}; + BoundaryInjector{ + ParticleGenerator(-1, 0)}; auto inject_ions = - BoundaryInjector{ - ParticleGenerator(-1, 1), grid}; + BoundaryInjector{ + ParticleGenerator(-1, 1)}; auto psc = makePscIntegrator(psc_params, grid, mflds, mprts, balance, collision, checks); diff --git a/src/libpsc/tests/test_mfields_io.cxx b/src/libpsc/tests/test_mfields_io.cxx index cd62490c63..e1c32c7d06 100644 --- a/src/libpsc/tests/test_mfields_io.cxx +++ b/src/libpsc/tests/test_mfields_io.cxx @@ -75,7 +75,9 @@ TYPED_TEST(MfieldsTest, WriteRead) auto mflds2 = Mfields{grid, NR_FIELDS, {}}; { auto reader = io.open("test.bp", kg::io::Mode::Read); + reader.beginStep(kg::io::StepMode::Read); reader.get("mflds", mflds2); + reader.endStep(); reader.close(); } @@ -106,7 +108,9 @@ TYPED_TEST(MfieldsTest, WriteWithGhostsRead) auto mflds2 = Mfields{grid, NR_FIELDS, {}}; { auto reader = io.open("test.bp", kg::io::Mode::Read); + reader.beginStep(kg::io::StepMode::Read); reader.get("mflds", mflds2); + reader.endStep(); reader.close(); } @@ -137,7 +141,9 @@ TYPED_TEST(MfieldsTest, WriteReadWithGhosts) auto mflds2 = Mfields{grid, NR_FIELDS, {2, 2, 2}}; { auto reader = io.open("test.bp", kg::io::Mode::Read); + reader.beginStep(kg::io::StepMode::Read); reader.get("mflds", mflds2); + reader.endStep(); reader.close(); } diff --git a/src/psc_shock.cxx b/src/psc_shock.cxx index 2ff164d5f9..0f39645f45 100644 --- a/src/psc_shock.cxx +++ b/src/psc_shock.cxx @@ -4,10 +4,12 @@ #include "output_fields.hxx" #include "psc_config.hxx" -#include "include/boundary_injector.hxx" #include "input_params.hxx" #include "kg/include/kg/VecRange.hxx" #include "libpsc/psc_output_particles/output_particles_adios2_impl.hxx" +#include "libpsc/psc_bnd_fields/radiating.hxx" +#include "libpsc/psc_particle_injectors/boundary_injector.hxx" +#include "libpsc/axis.hxx" // ====================================================================== // PSC configuration @@ -41,17 +43,27 @@ using Real3 = Vec3; PscParams psc_params; -double electron_temperature; -double ion_temperature; double electron_mass; double ion_mass; +double n_upstream; Double3 v_upstream; - -Real3 background_h_upstream; - -Real3 background_e; -Real3 background_h; +double te_upstream; +double ti_upstream; +Real3 h0_upstream; +Real3 e0; // e0 is constant + +// ---------------------------- +// standing shock params +double n_downstream; +Double3 v_downstream; +double te_downstream; +double ti_downstream; +Real3 h0_downstream; + +double transition_half_width; +double transition_steepness = 2.0; // at least sqrt(3)~1.74 +// ---------------------------- Int3 gdims; Double3 lengths; @@ -64,10 +76,13 @@ int out_interval; int marder_interval; std::string turb_method; +std::string shock_method; int nicell; int seed; +std::string checkpoint_filename; + // ====================================================================== // setupParameters @@ -79,57 +94,127 @@ void setupParameters(int argc, char** argv) std::string path_to_params(argv[1]); InputParams inputParams(path_to_params); - psc_params.stats_every = 1000; + shock_method = inputParams.getOrDefault("shock_method", "wall"); + + psc_params.stats_every = + inputParams.getOrDefault("stats_interval", 1000); psc_params.cfl = inputParams.getOrDefault("cfl", .75); - psc_params.write_checkpoint_every_step = 0; + checkpoint_filename = + inputParams.getOrDefault("checkpoint_filename", ""); + + electron_mass = inputParams.get("m_e"); + ion_mass = inputParams.get("m_i"); + + n_upstream = 1.0; + te_upstream = inputParams.get("T_e"); + ti_upstream = inputParams.get("T_i"); + + double v_shock = inputParams.get("v_shock"); + double theta_bn_deg = inputParams.get("θ_Bn_deg"); + double b0 = inputParams.get("B_0"); + + double theta_bn = theta_bn_deg * M_PI / 180.0; + double theta_xz = inputParams.get("θ_xz_deg") * M_PI / 180.0; + h0_upstream = Real3{sin(theta_bn) * cos(theta_xz), cos(theta_bn), + sin(theta_bn) * sin(theta_xz)} * + b0; + + if (shock_method == "wall" || shock_method == "none") { + v_upstream = {0.0, v_shock / 1.5, 0.0}; // factor is empirical + e0 = -v_upstream.cross(h0_upstream); + + // relativistic correction + double gamma = 1 / sqrt(1 - v_upstream.mag2()); + e0 *= gamma; + h0_upstream *= Real3{gamma, 1.0, gamma}; + } else if (shock_method == "relaxation") { + v_upstream = {0.0, v_shock, 0.0}; + e0 = -v_upstream.cross(h0_upstream); // upstream and downstream are the same + + // for perpendicular shock (2013 Balogh eq.3.36 and normalization in + // sec.3.3.1) + if (theta_bn_deg != 90.0) { + LOG_ERROR("θ_Bn must be 90° for relaxation method; got %f°\n", + theta_bn_deg); + } - electron_temperature = inputParams.get("electron_temperature"); - ion_temperature = inputParams.get("ion_temperature"); - electron_mass = inputParams.get("electron_mass"); - ion_mass = inputParams.get("ion_mass"); + // no relativistic correction, since these aren't relativistic RH conditions - inputParams.errIfPresentAndNotEqual("v_upstream_x", 0.0, ""); - v_upstream = {0.0, inputParams.get("v_upstream_y"), 0.0}; - inputParams.errIfPresentAndNotEqual("v_upstream_z", 0.0, ""); + double b_norm = sqrt(ion_mass * n_upstream * v_upstream.mag2()); + double t_norm = 0.5 * ion_mass * v_upstream.mag2(); + double beta1 = + 2.0 * n_upstream * (te_upstream + ti_upstream) / h0_upstream.mag2(); + Double3 B1 = h0_upstream / b_norm; + double T1 = (te_upstream + ti_upstream) / t_norm; - double b_angle_y_to_x_rad = inputParams.get("b_angle_y_to_x_rad"); - double b_mag = inputParams.get("b_mag"); - background_h_upstream = - b_mag * Real3{sin(b_angle_y_to_x_rad), cos(b_angle_y_to_x_rad), 0.0}; + double MA_sq = + v_upstream.mag2() * ion_mass * n_upstream / h0_upstream.mag2(); - double gamma = 1 / sqrt(1 - v_upstream.mag2()); - background_e = -gamma * v_upstream.cross(background_h_upstream); - // note: this only holds for vx=vz=0 - background_h = background_h_upstream * Real3{gamma, 0.0, gamma}; + double tmp = 1.0 + (1.0 + 2.5 * beta1) * B1.mag2(); + double r = 8.0 / (tmp + sqrt(sqr(tmp) + 2 * B1.mag2())); + r = inputParams.getOrDefault("r", r); + + double heating_factor = + 1.0 + + 4.0 / (5.0 * T1) * ((sqr(r) - 1.0) / (2.0 * sqr(r)) + (1.0 - r) / MA_sq); + + n_downstream = n_upstream * r; + v_downstream = v_upstream / r; + h0_downstream = h0_upstream * r; + + te_downstream = + inputParams.getOrDefault("T_e2", te_upstream * heating_factor); + ti_downstream = + inputParams.getOrDefault("T_i2", ti_upstream * heating_factor); + } gdims[0] = inputParams.get("nx"); gdims[1] = inputParams.get("ny"); gdims[2] = inputParams.get("nz"); psc_params.nmax = inputParams.get("nt"); - n_patches[0] = inputParams.get("n_patches_x"); - n_patches[1] = inputParams.get("n_patches_y"); - n_patches[2] = inputParams.get("n_patches_z"); + n_patches[0] = inputParams.get("npx"); + n_patches[1] = inputParams.get("npy"); + n_patches[2] = inputParams.get("npz"); - Double3 dx = {inputParams.get("dx"), inputParams.get("dy"), - inputParams.get("dz")}; - - lengths = Double3(gdims) * dx; - - if (inputParams.warnIfPresent("turb_dB^2", "set turb_dB instead")) { - turb_db2 = inputParams.get("turb_dB^2"); + if (inputParams.has("lx")) { + lengths[0] = inputParams.get("lx"); } else { - turb_db2 = sqr(inputParams.get("turb_dB")); + lengths[0] = inputParams.get("dx") * gdims[0]; + } + if (inputParams.has("ly")) { + lengths[1] = inputParams.get("ly"); + } else { + lengths[1] = inputParams.get("dy") * gdims[1]; + } + if (inputParams.has("lz")) { + lengths[2] = inputParams.get("lz"); + } else { + lengths[2] = inputParams.get("dz") * gdims[2]; + } + + transition_half_width = lengths[1] / 8.0; + + turb_db2 = sqr(inputParams.get("dB")); + turb_correlation_length = inputParams.get("L_c"); + + if (inputParams.has("checkpoint_interval")) { + psc_params.write_checkpoint_every_step = + inputParams.get("checkpoint_interval"); + inputParams.errIfPresent( + "n_checkpoints", + "n_checkpoints is mutually exclusive with checkpoint_interval"); + } else if (inputParams.has("n_checkpoints")) { + int n_checkpoints = inputParams.get("n_checkpoints"); + if (n_checkpoints > 0) { + psc_params.write_checkpoint_every_step = psc_params.nmax / n_checkpoints; + } } - turb_correlation_length = inputParams.get("turb_correlation_length"); int n_writes = inputParams.getOrDefault("n_writes", 100); out_interval = psc_params.nmax / n_writes; marder_interval = inputParams.getOrDefault("marder_interval", -1); - inputParams.errIfPresentAndEqual("mirror_domain", true, - "only 'false' is permitted"); - turb_method = inputParams.getOrDefault("turb_method", "alfven_dense"); @@ -141,6 +226,23 @@ void setupParameters(int argc, char** argv) dst << src.rdbuf(); } +template +T interpolate_across_shock(T upstream, T downstream, double y) +{ + if (transition_half_width == 0.0) { + return y > 0.0 ? downstream : upstream; + } else { + // smooth transition function + y /= transition_half_width; + double sigmoid_val = abs(y) >= 1.0 + ? (y > 0.0 ? 1.0 : -1.0) + : tanh(transition_steepness * y / (1.0 - y * y)); + double weight_downstream = (sigmoid_val + 1.0) / 2.0; + double weight_upstream = 1.0 - weight_downstream; + return upstream * weight_upstream + downstream * weight_downstream; + } +} + // ====================================================================== // setupGrid // @@ -152,14 +254,39 @@ void setupParameters(int argc, char** argv) Grid_t* setupGrid() { // FIXME add a check to catch mismatch between Dim and n grid points early - Double3 corner = {0.0, 0.0, 0.0}; + + Double3 corner; + int bnd_fld_lower; + int bnd_fld_upper; + int bnd_prt_lower; + int bnd_prt_upper; + + if (shock_method == "wall") { + corner = {0.0, 0.0, 0.0}; + bnd_fld_lower = BND_FLD_OPEN; + bnd_fld_upper = BND_FLD_CONDUCTING_WALL; + bnd_prt_lower = BND_PRT_OPEN; + bnd_prt_upper = BND_PRT_REFLECTING; + } else if (shock_method == "none") { + corner = {0.0, 0.0, 0.0}; + bnd_fld_lower = BND_FLD_PERIODIC; + bnd_fld_upper = BND_FLD_PERIODIC; + bnd_prt_lower = BND_PRT_PERIODIC; + bnd_prt_upper = BND_PRT_PERIODIC; + } else if (shock_method == "relaxation") { + corner = {0.0, -lengths[1] / 2.0, 0.0}; + bnd_fld_lower = BND_FLD_OPEN; + bnd_fld_upper = BND_FLD_OPEN; + bnd_prt_lower = BND_PRT_OPEN; + bnd_prt_upper = BND_PRT_OPEN; + } + auto domain = Grid_t::Domain{gdims, lengths, corner, n_patches}; - auto bc = - psc::grid::BC{{BND_FLD_PERIODIC, BND_FLD_OPEN, BND_FLD_PERIODIC}, - {BND_FLD_PERIODIC, BND_FLD_CONDUCTING_WALL, BND_FLD_PERIODIC}, - {BND_PRT_PERIODIC, BND_PRT_OPEN, BND_PRT_PERIODIC}, - {BND_PRT_PERIODIC, BND_PRT_REFLECTING, BND_PRT_PERIODIC}}; + auto bc = psc::grid::BC{{BND_FLD_PERIODIC, bnd_fld_lower, BND_FLD_PERIODIC}, + {BND_FLD_PERIODIC, bnd_fld_upper, BND_FLD_PERIODIC}, + {BND_PRT_PERIODIC, bnd_prt_lower, BND_PRT_PERIODIC}, + {BND_PRT_PERIODIC, bnd_prt_upper, BND_PRT_PERIODIC}}; auto kinds = Grid_t::Kinds(NR_KINDS); kinds[KIND_ELECTRON] = {-1.0, electron_mass, "e"}; @@ -188,21 +315,33 @@ void initializeParticles(Balance& balance, Grid_t*& grid_ptr, Mparticles& mprts) setup_particles.random_offsets = true; setup_particles.initial_momentum_gamma_correction = true; - 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, - v_upstream, - {temperature, temperature, temperature}, - np.tag}); - }; - - partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, - init_np); + if (shock_method == "wall" || shock_method == "none") { + auto init_np = [&](int kind, Double3 pos, int p, Int3 idx, + psc_particle_np& np) { + double t = np.kind == KIND_ION ? ti_upstream : te_upstream; + np.n = 1.0; + np.p = setup_particles.createMaxwellian( + {np.kind, np.n, v_upstream, {t, t, t}, np.tag}); + }; + + partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, + init_np); + } else if (shock_method == "relaxation") { + auto init_np = [&](int kind, Double3 pos, int p, Int3 idx, + psc_particle_np& np) { + np.n = interpolate_across_shock(n_upstream, n_downstream, pos[1]); + // interpolate v_thermal, not T itself + double t = sqr(interpolate_across_shock( + sqrt(np.kind == KIND_ION ? ti_upstream : te_upstream), + sqrt(np.kind == KIND_ION ? ti_downstream : te_downstream), pos[1])); + Double3 v = interpolate_across_shock(v_upstream, v_downstream, pos[1]); + np.p = + setup_particles.createMaxwellian({np.kind, np.n, v, {t, t, t}, np.tag}); + }; + + partitionAndSetupParticles(setup_particles, balance, grid_ptr, mprts, + init_np); + } } // ====================================================================== @@ -218,13 +357,22 @@ void add_background_fields(MfieldsState& mflds) int n_ghosts = mflds.ibn().max(); grid.Foreach_3d(n_ghosts, n_ghosts, [&](int jx, int jy, int jz) { - field_patch(HX, jx, jy, jz) += background_h[0]; - field_patch(HY, jx, jy, jz) += background_h[1]; - field_patch(HZ, jx, jy, jz) += background_h[2]; + Real3 h0; + if (shock_method == "wall" || shock_method == "none") { + h0 = h0_upstream; + } else if (shock_method == "relaxation") { + Double3 pos = centering::get_pos(patch, {jx, jy, jz}, centering::NC, 0); + h0 = interpolate_across_shock(h0_upstream, h0_downstream, pos[1]); + h0[1] = h0_upstream[1]; // parallel B isn't compressed + } + + field_patch(HX, jx, jy, jz) += h0[0]; + field_patch(HY, jx, jy, jz) += h0[1]; + field_patch(HZ, jx, jy, jz) += h0[2]; - field_patch(EX, jx, jy, jz) += background_e[0]; - field_patch(EY, jx, jy, jz) += background_e[1]; - field_patch(EZ, jx, jy, jz) += background_e[2]; + field_patch(EX, jx, jy, jz) += e0[0]; + field_patch(EY, jx, jy, jz) += e0[1]; + field_patch(EZ, jx, jy, jz) += e0[2]; }); } } @@ -275,24 +423,40 @@ void inject_b_from_potential(MfieldsState& mflds, PscConfig::Mfields& vector_potential) { const auto& grid = mflds.grid(); + Real3 dx = grid.domain.dx; for (int p = 0; p < mflds.n_patches(); ++p) { auto field_patch = make_Fields3d(mflds[p]); auto vector_potential_patch = make_Fields3d(vector_potential[p]); grid.Foreach_3d(2, 1, [&](int jx, int jy, int jz) { - field_patch(HX, jx, jy, jz) = vector_potential_patch(AZ, jx, jy + 1, jz) - - vector_potential_patch(AZ, jx, jy, jz) - - vector_potential_patch(AY, jx, jy, jz + 1) + - vector_potential_patch(AY, jx, jy, jz); - field_patch(HY, jx, jy, jz) = vector_potential_patch(AX, jx, jy, jz + 1) - - vector_potential_patch(AX, jx, jy, jz) - - vector_potential_patch(AZ, jx + 1, jy, jz) + - vector_potential_patch(AZ, jx, jy, jz); - field_patch(HZ, jx, jy, jz) = vector_potential_patch(AY, jx + 1, jy, jz) - - vector_potential_patch(AY, jx, jy, jz) - - vector_potential_patch(AX, jx, jy + 1, jz) + - vector_potential_patch(AX, jx, jy, jz); + field_patch(HX, jx, jy, jz) = + (Dim::is_invar(1) ? 0 + : vector_potential_patch(AZ, jx, jy + 1, jz) - + vector_potential_patch(AZ, jx, jy, jz)) / + dx[1] - + (Dim::is_invar(2) ? 0 + : vector_potential_patch(AY, jx, jy, jz + 1) - + vector_potential_patch(AY, jx, jy, jz)) / + dx[2]; + field_patch(HY, jx, jy, jz) = + (Dim::is_invar(2) ? 0 + : vector_potential_patch(AX, jx, jy, jz + 1) - + vector_potential_patch(AX, jx, jy, jz)) / + dx[2] - + (Dim::is_invar(0) ? 0 + : vector_potential_patch(AZ, jx + 1, jy, jz) - + vector_potential_patch(AZ, jx, jy, jz)) / + dx[0]; + field_patch(HZ, jx, jy, jz) = + (Dim::is_invar(0) ? 0 + : vector_potential_patch(AY, jx + 1, jy, jz) - + vector_potential_patch(AY, jx, jy, jz)) / + dx[0] - + (Dim::is_invar(1) ? 0 + : vector_potential_patch(AX, jx, jy + 1, jz) - + vector_potential_patch(AX, jx, jy, jz)) / + dx[1]; }); } } @@ -320,7 +484,7 @@ void inject_plane_alfven_wave(PscConfig::Mfields& vector_potential, double db, } Double3 xp_hat{cos_theta * cos_phi, cos_theta * sin_phi, -sin_theta}; - Double3 yp_hat{sin_phi, -cos_phi, 0}; + Double3 yp_hat{-sin_phi, cos_phi, 0}; Double3 a_vec = db * cos(polarization) / k2 * xp_hat.cross(k_vec); Double3 b_vec = db * sin(polarization) / k2 * yp_hat.cross(k_vec); @@ -495,7 +659,8 @@ void inject_turbulence_dense(MfieldsState& mflds) Int3 i3_min = (1 - gdims) / 2; Int3 i3_max = gdims / 2; - // inject in only half of k-space, since +k and -k modes are indistinguishable + // inject in only half of k-space, since +k and -k modes are + // indistinguishable for (int d = 0; d < 3; d++) { if (gdims[d] > 2) { i3_min[d] = 0; @@ -581,7 +746,7 @@ void inject_turbulence_dense(MfieldsState& mflds) set_mean_b2(mflds, turb_db2); } -void initializeFields(MfieldsState& mflds) +void initialize_turbulence(MfieldsState& mflds) { if (turb_db2 > 0.0) { if (turb_method == "alfven_dense") { @@ -593,11 +758,9 @@ void initializeFields(MfieldsState& mflds) LOG_ERROR("Unrecognized turbulence method: %s\n", turb_method.c_str()); } } - - add_background_fields(mflds); } -struct AdvectedPeriodicFields : RadiatingBoundary +struct AdvectedPeriodicFields : psc::bnd::field::PulseBase { static const int DIM_Y = 1; @@ -605,18 +768,12 @@ struct AdvectedPeriodicFields : RadiatingBoundary Real3 background_e, Real3 background_h) : v_advect(v_advect), grid(mflds.grid()) { + // mflds must NOT include background fields at this point // FIXME would be better to exclude J, but the interpolator uses EX, etc. auto&& e_b_fields = mflds.storage().view(_all, _all, _all, _all, _all); // FIXME this probably isn't the best way to copy a gtensor array cycled_fields = gt::zeros_like(e_b_fields); cycled_fields.view(_all, _all, _all, _all, _all) = e_b_fields; - - for (int d = 0; d < 3; d++) { - cycled_fields.view(_all, _all, _all, EX + d, _all) = - cycled_fields.view(_all, _all, _all, EX + d, _all) - background_e[d]; - cycled_fields.view(_all, _all, _all, HX + d, _all) = - cycled_fields.view(_all, _all, _all, HX + d, _all) - background_h[d]; - } } Real3 advect_x3(Real3 x3, double t) @@ -628,15 +785,16 @@ struct AdvectedPeriodicFields : RadiatingBoundary { real_t patch_size = grid.domain.length[DIM_Y] / grid.domain.np[DIM_Y]; int n_patches_to_the_left = 0; - while (x3_advected[DIM_Y] < grid.domain.corner[DIM_Y]) { + // note: input x3 is already patch-local; we are just shifting to a + // *different* patch + while (x3_advected[DIM_Y] < 0.0) { x3_advected[DIM_Y] += patch_size; n_patches_to_the_left += 1; } return n_patches_to_the_left; } - void calc_e_h(double t, int p, Real3 x3, int d_e, real_t& e, int d_h, - real_t& h) + real_t sample_exterior_field(int m, double t, int p, Real3 x3) override { Real3 x3_advected = advect_x3(x3, t); int n_patches_to_the_left = shift_to_patch_local(x3_advected); @@ -649,56 +807,19 @@ struct AdvectedPeriodicFields : RadiatingBoundary auto em = decltype(ip)::fields_t( cycled_fields.view(_all, _all, _all, _all, p), -grid.ibn); - switch (d_e) { - case 0: e = ip.ex(em); break; - case 1: e = ip.ey(em); break; - case 2: e = ip.ez(em); break; - } - switch (d_h) { - case 0: h = ip.hx(em); break; - case 1: h = ip.hy(em); break; - case 2: h = ip.hz(em); break; + switch (m) { + case EX: return ip.ex(em) + e0[0]; + case EY: return ip.ey(em) + e0[1]; + case EZ: return ip.ez(em) + e0[2]; + case HX: return ip.hx(em) + h0_upstream[0]; + case HY: return ip.hy(em) + h0_upstream[1]; + case HZ: return ip.hz(em) + h0_upstream[2]; + default: return 0.0; } } - real_t pulse_s_lower(double t, int d, int p, Real3 x3) override + void tick(double t) override { - int d1 = (d + 1) % 3; - int d2 = (d + 2) % 3; - - real_t e, h; - calc_e_h(t, p, x3, d1, e, d2, h); - - return (e + h) / 2.0; - } - - real_t pulse_p_lower(double t, int d, int p, Real3 x3) override - { - int d1 = (d + 1) % 3; - int d2 = (d + 2) % 3; - - real_t e, h; - calc_e_h(t, p, x3, d2, e, d1, h); - - return (e - h) / 2.0; - } - - real_t pulse_s_upper(double t, int d, int p, Real3 x3) override - { - return 0.0; - } - - real_t pulse_p_upper(double t, int d, int p, Real3 x3) override - { - return 0.0; - } - - void update_cache_lower(double t, int d) override - { - if (d != DIM_Y) { - return; - } - // TODO make this work with >1 patch per process? assert(grid.n_patches() == 1); @@ -710,15 +831,20 @@ struct AdvectedPeriodicFields : RadiatingBoundary return; } - LOG_INFO("cycling turbulence...\n"); + cycle_turbulence(n_patches_to_the_left - n_patch_cycles); + } + + void cycle_turbulence(int n_patches) + { + LOG_INFO("cycling turbulence... (t=%f)\n", grid.time()); // hack: guess the rank based on how mrc does it for simple domains // (can't use mrc, because it wouldn't apply periodicity) Int3 np = grid.domain.np; Int3 proc = grid.localPatchInfo(0).idx3; - Int3 dest_proc = (proc + Int3::unit(DIM_Y)) % np; + Int3 dest_proc = (proc + Int3::unit(DIM_Y) * n_patches) % np; int dest_rank = flatten_index(dest_proc.reverse(), np.reverse()); - Int3 source_proc = (proc - Int3::unit(DIM_Y) + np) % np; + Int3 source_proc = ((proc - Int3::unit(DIM_Y) * n_patches) % np + np) % np; int source_rank = flatten_index(source_proc.reverse(), np.reverse()); MPI_Status status; @@ -726,7 +852,7 @@ struct AdvectedPeriodicFields : RadiatingBoundary MpiDtypeTraits::value(), dest_rank, 0, source_rank, 0, grid.comm(), &status); - n_patch_cycles += 1; + n_patch_cycles += n_patches; } const Grid_t& grid; @@ -806,25 +932,41 @@ static void run(int argc, char** argv) int oute_interval = -100; DiagEnergies oute{grid.comm(), oute_interval}; - auto ion_injector = - BoundaryInjector( - ParticleGeneratorMaxwellian( - KIND_ION, grid.kinds[KIND_ION], v_upstream, - {ion_temperature, ion_temperature, ion_temperature}, true), - grid); - auto electron_injector = - BoundaryInjector( + auto ion_injector_lo = BoundaryInjector( + ParticleGeneratorMaxwellian(KIND_ION, grid.kinds[KIND_ION], v_upstream, + {ti_upstream, ti_upstream, ti_upstream}), + n_upstream); + auto electron_injector_lo = + BoundaryInjector( + ParticleGeneratorMaxwellian(KIND_ELECTRON, grid.kinds[KIND_ELECTRON], + v_upstream, + {te_upstream, te_upstream, te_upstream}), + n_upstream); + + auto ion_injector_hi = BoundaryInjector( + ParticleGeneratorMaxwellian(KIND_ION, grid.kinds[KIND_ION], v_downstream, + {ti_downstream, ti_downstream, ti_downstream}), + n_downstream); + auto electron_injector_hi = + BoundaryInjector( ParticleGeneratorMaxwellian( - KIND_ELECTRON, grid.kinds[KIND_ELECTRON], v_upstream, - {electron_temperature, electron_temperature, electron_temperature}, - true), - grid); + KIND_ELECTRON, grid.kinds[KIND_ELECTRON], v_downstream, + {te_downstream, te_downstream, te_downstream}), + n_downstream); // ---------------------------------------------------------------------- // set up initial conditions - initializeParticles(balance, grid_ptr, mprts); - initializeFields(mflds); + if (checkpoint_filename.empty()) { + initializeParticles(balance, grid_ptr, mprts); + initialize_turbulence(mflds); + } else { + read_checkpoint(checkpoint_filename, *grid_ptr, mprts, mflds); + } // ---------------------------------------------------------------------- // run the simulation @@ -834,18 +976,56 @@ static void run(int argc, char** argv) psc.add_gauss_corrector(&marder); - psc.bndf.background_e = background_e; - psc.bndf.background_h = background_h; - psc.bndf.radiation = new AdvectedPeriodicFields{mflds, v_upstream[1], - background_e, background_h}; - psc.add_diagnostic(&out_fields); psc.add_diagnostic(&out_moments); psc.add_diagnostic(&outp); psc.add_diagnostic(&oute); - psc.add_injector(&ion_injector); - psc.add_injector(&electron_injector); + using psc::Axis; + using psc::bnd::LoHi; + using ConstantPulse = psc::bnd::field::ConstantPulse; + + if (shock_method != "none") { + psc.add_injector(&ion_injector_lo); + psc.add_injector(&electron_injector_lo); + + if (turb_db2 > 0.0 && v_upstream[1] > 0.0) { + if (checkpoint_filename.empty()) { + // mflds is currently just the pure, initial turbulence + psc.add_field_bc(new psc::bnd::field::Radiating( + AdvectedPeriodicFields{mflds, v_upstream[1], e0, h0_upstream}, + Axis::Y, LoHi::Lo)); + } else { + // mflds is completely unrelated; need to re-initialize turbulence + MfieldsState mflds2{grid}; + initialize_turbulence(mflds2); + psc.add_field_bc(new psc::bnd::field::Radiating( + AdvectedPeriodicFields{mflds2, v_upstream[1], e0, h0_upstream}, + Axis::Y, LoHi::Lo)); + } + } else { + psc.add_field_bc( + new psc::bnd::field::Radiating( + ConstantPulse{e0, h0_upstream}, Axis::Y, LoHi::Lo)); + } + } + + if (shock_method == "relaxation") { + psc.add_injector(&ion_injector_hi); + psc.add_injector(&electron_injector_hi); + + psc.add_field_bc( + new psc::bnd::field::Radiating( + ConstantPulse{e0, h0_downstream}, Axis::Y, LoHi::Hi)); + } + + if (checkpoint_filename.empty()) { + // add background after initializing radiation inflow, which only wants + // the perturbations to B + add_background_fields(mflds); + } psc.integrate(); }