diff options
Diffstat (limited to 'Source')
-rw-r--r-- | Source/Initialization/PlasmaInjector.H | 1 | ||||
-rw-r--r-- | Source/Initialization/PlasmaInjector.cpp | 1 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.H | 5 | ||||
-rw-r--r-- | Source/Particles/PhysicalParticleContainer.cpp | 97 |
4 files changed, 72 insertions, 32 deletions
diff --git a/Source/Initialization/PlasmaInjector.H b/Source/Initialization/PlasmaInjector.H index cd5802a55..f998e217e 100644 --- a/Source/Initialization/PlasmaInjector.H +++ b/Source/Initialization/PlasmaInjector.H @@ -293,6 +293,7 @@ public: amrex::Real z_rms; amrex::Real q_tot; long npart; + int do_symmetrize = 0; bool radially_weighted = true; diff --git a/Source/Initialization/PlasmaInjector.cpp b/Source/Initialization/PlasmaInjector.cpp index 52b5d8b61..f9642d1b6 100644 --- a/Source/Initialization/PlasmaInjector.cpp +++ b/Source/Initialization/PlasmaInjector.cpp @@ -286,6 +286,7 @@ PlasmaInjector::PlasmaInjector(int ispecies, const std::string& name) pp.get("z_rms", z_rms); pp.get("q_tot", q_tot); pp.get("npart", npart); + pp.query("do_symmetrize", do_symmetrize); gaussian_beam = true; parseMomentum(pp); } diff --git a/Source/Particles/PhysicalParticleContainer.H b/Source/Particles/PhysicalParticleContainer.H index df1e0296a..bd8cfb47e 100644 --- a/Source/Particles/PhysicalParticleContainer.H +++ b/Source/Particles/PhysicalParticleContainer.H @@ -94,7 +94,10 @@ public: void AddGaussianBeam(amrex::Real x_m, amrex::Real y_m, amrex::Real z_m, amrex::Real x_rms, amrex::Real y_rms, amrex::Real z_rms, - amrex::Real q_tot, long npart); + amrex::Real q_tot, long npart, int do_symmetrize); + + void CheckAndAddParticle(amrex::Real x, amrex::Real y, amrex::Real z, + std::array<amrex::Real, 3> u, amrex::Real weight); virtual void GetParticleSlice(const int direction, const amrex::Real z_old, const amrex::Real z_new, const amrex::Real t_boost, diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp index ab37bf6ff..1f517fccb 100644 --- a/Source/Particles/PhysicalParticleContainer.cpp +++ b/Source/Particles/PhysicalParticleContainer.cpp @@ -181,7 +181,8 @@ void PhysicalParticleContainer::MapParticletoBoostedFrame(Real& x, Real& y, Real void PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real x_rms, Real y_rms, Real z_rms, - Real q_tot, long npart) { + Real q_tot, long npart, + int do_symmetrize) { const Geometry& geom = m_gdb->Geom(0); RealBox containing_bx = geom.ProbDomain(); @@ -191,13 +192,15 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, std::normal_distribution<double> disty(y_m, y_rms); std::normal_distribution<double> distz(z_m, z_rms); - std::array<Real,PIdx::nattribs> attribs; - attribs.fill(0.0); - if (ParallelDescriptor::IOProcessor()) { - std::array<Real, 3> u; - Real weight; - for (long i = 0; i < npart; ++i) { + std::array<Real, 3> u; + Real weight; + // If do_symmetrize, create 4x fewer particles, and + // Replicate each particle 4 times (x,y) (-x,y) (x,-y) (-x,-y) + if (do_symmetrize){ + npart /= 4; + } + for (long i = 0; i < npart; ++i) { #if ( AMREX_SPACEDIM == 3 | WARPX_RZ) weight = q_tot/npart/charge; Real x = distx(mt); @@ -209,29 +212,27 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, Real y = 0.; Real z = distz(mt); #endif - if (plasma_injector->insideBounds(x, y, z)) { - plasma_injector->getMomentum(u, x, y, z); - if (WarpX::gamma_boost > 1.) { - MapParticletoBoostedFrame(x, y, z, u); - } - attribs[PIdx::ux] = u[0]; - attribs[PIdx::uy] = u[1]; - attribs[PIdx::uz] = u[2]; - attribs[PIdx::w ] = weight; - - if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) - { - auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); - particle_tile.push_back_real(particle_comps["xold"], x); - particle_tile.push_back_real(particle_comps["yold"], y); - particle_tile.push_back_real(particle_comps["zold"], z); - - particle_tile.push_back_real(particle_comps["uxold"], u[0]); - particle_tile.push_back_real(particle_comps["uyold"], u[1]); - particle_tile.push_back_real(particle_comps["uzold"], u[2]); - } - - AddOneParticle(0, 0, 0, x, y, z, attribs); + if (plasma_injector->insideBounds(x, y, z)) { + plasma_injector->getMomentum(u, x, y, z); + if (do_symmetrize){ + std::array<Real, 3> u_tmp; + Real x_tmp, y_tmp; + // Add four particles to the beam: + // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy) + for (int ix=0; ix<2; ix++){ + for (int iy=0; iy<2; iy++){ + u_tmp = u; + x_tmp = x*std::pow(-1,ix); + u_tmp[0] *= std::pow(-1,ix); + y_tmp = y*std::pow(-1,iy); + u_tmp[1] *= std::pow(-1,iy); + CheckAndAddParticle(x_tmp, y_tmp, z, + u_tmp, weight/4); + } + } + } else { + CheckAndAddParticle(x, y, z, u, weight); + } } } } @@ -239,6 +240,39 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m, } void +PhysicalParticleContainer::CheckAndAddParticle(Real x, Real y, Real z, + std::array<Real, 3> u, + Real weight) +{ + std::array<Real,PIdx::nattribs> attribs; + attribs.fill(0.0); + + // update attribs with input arguments + if (WarpX::gamma_boost > 1.) { + MapParticletoBoostedFrame(x, y, z, u); + } + attribs[PIdx::ux] = u[0]; + attribs[PIdx::uy] = u[1]; + attribs[PIdx::uz] = u[2]; + attribs[PIdx::w ] = weight; + + if (WarpX::do_boosted_frame_diagnostic && do_boosted_frame_diags) + { + // need to create old values + auto& particle_tile = DefineAndReturnParticleTile(0, 0, 0); + particle_tile.push_back_real(particle_comps["xold"], x); + particle_tile.push_back_real(particle_comps["yold"], y); + particle_tile.push_back_real(particle_comps["zold"], z); + + particle_tile.push_back_real(particle_comps["uxold"], u[0]); + particle_tile.push_back_real(particle_comps["uyold"], u[1]); + particle_tile.push_back_real(particle_comps["uzold"], u[2]); + } + // add particle + AddOneParticle(0, 0, 0, x, y, z, attribs); +} + +void PhysicalParticleContainer::AddParticles (int lev) { BL_PROFILE("PhysicalParticleContainer::AddParticles()"); @@ -263,7 +297,8 @@ PhysicalParticleContainer::AddParticles (int lev) plasma_injector->y_rms, plasma_injector->z_rms, plasma_injector->q_tot, - plasma_injector->npart); + plasma_injector->npart, + plasma_injector->do_symmetrize); return; |