option do symmetrize gaussian particle beam

1 files changed, 60 insertions, 31 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 212084e64..38ea55d6e 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,24 @@ 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){
+                    // Add four particles to the beam:
+                    // (x,ux,y,uy) (-x,-ux,y,uy) (x,ux,-y,-uy) (-x,-ux,-y,-uy)
+                    std::array<Real, 3> u_tmp = u;
+                    for (int ix=0; ix<2; ix++){
+                        for (int iy=0; iy<2; iy++){
+                            x        *= std::pow(-1,ix);
+                            u_tmp[0] *= std::pow(-1,ix);
+                            y        *= std::pow(-1,iy);
+                            u_tmp[1] *= std::pow(-1,iy);
+                            CheckAndAddParticle(x, y, z, u_tmp, weight);
+                        }
+                    }
+                } else {
+                    CheckAndAddParticle(x, y, z, u, weight);
+                }
             }
         }
     }
@@ -239,6 +237,36 @@ 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);
+
+    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);
+}
+
+void
 PhysicalParticleContainer::AddParticles (int lev)
 {
     BL_PROFILE("PhysicalParticleContainer::AddParticles()");
@@ -263,7 +291,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;