Cut gaussian beam (#868)

* add option to cut gaussian

* cut gaussian distrubution at injection

* add test for cut gaussian beam

* increase tolerance

* doc: mention that cutting the gaussian results in a lower total charge

* more explicit

* beware the const-able

Author: MaxThevenet

Docs/source/running_cpp/parameters.rst

@@ -276,6 +276,8 @@ Particle initialization
       ``<species_name>.npart`` (number of particles in the beam),
       ``<species_name>.x/y/z_m`` (average position in `x/y/z`),
       ``<species_name>.x/y/z_rms`` (standard deviation in `x/y/z`),
+      ``<species_name>.x/y/z_rms`` (standard deviation in `x/y/z`),
+      ``<species_name>.x/y/z_cut`` (optional, particles with ``abs(x-x_m) > x_cut*x_rms`` are not injected, same for y and z. ``<species_name>.q_tot`` is the charge of the un-cut beam, so that cutting the distribution is likely to result in a lower total charge),
       and optional argument ``<species_name>.do_symmetrize`` (whether to
       symmetrize the beam in the x and y directions).
 

Examples/Tests/initial_distribution/analysis_distribution.py

@@ -19,7 +19,7 @@
 from read_raw_data import read_reduced_diags_histogram
 
 # print tolerance
-tolerance = 0.007
+tolerance = 0.01
 print('Tolerance:', tolerance)
 
 #===============================
@@ -96,17 +96,34 @@
 
 # parameters of theory
 x_rms = 0.25
+z_cut = 2.
 q_tot = -1.0e-20
 q_e   = -1.602176634e-19
 npart = q_tot/q_e
-db    = 0.04
+db    = bin_value[1]-bin_value[0]
 
 # compute the analytical solution
-f = npart*db * np.exp(-0.5*(bin_value/x_rms)**2)/(x_rms*np.sqrt(2.0*scc.pi))
-f_peak = np.amax(f)
+f_xy = npart*db * np.exp(-0.5*(bin_value/x_rms)**2)/(x_rms*np.sqrt(2.0*scc.pi)) * scs.erf(z_cut/np.sqrt(2.))
+f_z = npart*db * np.exp(-0.5*(bin_value/x_rms)**2)/(x_rms*np.sqrt(2.0*scc.pi))
+f_z[ np.absolute(bin_value) > z_cut * x_rms ] = 0.
+f_peak = np.amax(f_z)
 
 # compute error
-f4_error = np.sum(np.abs(f-h4x)+np.abs(f-h4y)+np.abs(f-h4z))/bin_value.size / f_peak
+f4_error = np.sum(np.abs(f_xy-h4x)+np.abs(f_xy-h4y)+np.abs(f_z-h4z))/bin_value.size / f_peak
+
+do_plot = False
+if do_plot:
+    import matplotlib.pyplot as plt
+    plt.figure()
+    plt.subplot(121)
+    plt.plot(bin_value, f_xy, '+-', label='ref')
+    plt.plot(bin_value, h4x, '+--', label='sim')
+    plt.legend()
+    plt.subplot(122)
+    plt.plot(bin_value, f_z, '+-', label='ref')
+    plt.plot(bin_value, h4z, '+--', label='sim')
+    plt.legend()
+    plt.savefig('toto.pdf', bbox_inches='tight')
 
 print('Gaussian position distribution difference:', f4_error)
 

Examples/Tests/initial_distribution/inputs

@@ -60,6 +60,7 @@ beam.injection_style            = "gaussian_beam"
 beam.x_rms                      = 0.25
 beam.y_rms                      = 0.25
 beam.z_rms                      = 0.25
+beam.z_cut                      = 2.
 beam.x_m                        = 0.0
 beam.y_m                        = 0.0
 beam.z_m                        = 0.0

Source/Initialization/PlasmaInjector.H

@@ -64,6 +64,9 @@ public:
     amrex::Real x_rms;
     amrex::Real y_rms;
     amrex::Real z_rms;
+    amrex::Real x_cut = std::numeric_limits<amrex::Real>::max();
+    amrex::Real y_cut = std::numeric_limits<amrex::Real>::max();
+    amrex::Real z_cut = std::numeric_limits<amrex::Real>::max();
     amrex::Real q_tot;
     long npart;
     int do_symmetrize = 0;

Source/Initialization/PlasmaInjector.cpp

@@ -176,6 +176,9 @@ PlasmaInjector::PlasmaInjector (int ispecies, const std::string& name)
         pp.get("x_rms", x_rms);
         pp.get("y_rms", y_rms);
         pp.get("z_rms", z_rms);
+        pp.query("x_cut", x_cut);
+        pp.query("y_cut", y_cut);
+        pp.query("z_cut", z_cut);
         pp.get("q_tot", q_tot);
         pp.get("npart", npart);
         pp.query("do_symmetrize", do_symmetrize);

Source/Particles/PhysicalParticleContainer.H

@@ -199,9 +199,10 @@ public:
 
     void MapParticletoBoostedFrame(amrex::Real& x, amrex::Real& y, amrex::Real& z, std::array<amrex::Real, 3>& u);
 
-    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, int do_symmetrize);
+    void AddGaussianBeam(const amrex::Real x_m, const amrex::Real y_m, const amrex::Real z_m,
+                         const amrex::Real x_rms, const amrex::Real y_rms, const amrex::Real z_rms,
+                         const amrex::Real x_cut, const amrex::Real y_cut, const amrex::Real z_cut,
+                         const amrex::Real q_tot, long npart, const int do_symmetrize);
 
     /** Load a particle beam from an external file
      *

Source/Particles/PhysicalParticleContainer.cpp

@@ -228,10 +228,12 @@ 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,
-                                           int do_symmetrize) {
+PhysicalParticleContainer::AddGaussianBeam (
+    const Real x_m, const Real y_m, const Real z_m,
+    const Real x_rms, const Real y_rms, const Real z_rms,
+    const Real x_cut, const Real y_cut, const Real z_cut,
+    const Real q_tot, long npart,
+    const int do_symmetrize) {
 
     std::mt19937_64 mt(0451);
     std::normal_distribution<double> distx(x_m, x_rms);
@@ -256,17 +258,20 @@ PhysicalParticleContainer::AddGaussianBeam(Real x_m, Real y_m, Real z_m,
         }
         for (long i = 0; i < npart; ++i) {
 #if (defined WARPX_DIM_3D) || (WARPX_DIM_RZ)
-            Real weight = q_tot/npart/charge;
-            Real x = distx(mt);
-            Real y = disty(mt);
-            Real z = distz(mt);
+            const Real weight = q_tot/(npart*charge);
+            const Real x = distx(mt);
+            const Real y = disty(mt);
+            const Real z = distz(mt);
 #elif (defined WARPX_DIM_XZ)
-            Real weight = q_tot/npart/charge/y_rms;
-            Real x = distx(mt);
-            Real y = 0.;
-            Real z = distz(mt);
+            const Real weight = q_tot/(npart*charge*y_rms);
+            const Real x = distx(mt);
+            constexpr Real y = 0.;
+            const Real z = distz(mt);
 #endif
-            if (plasma_injector->insideBounds(x, y, z)) {
+            if (plasma_injector->insideBounds(x, y, z)  &&
+                std::abs( x - x_m ) < x_cut * x_rms     &&
+                std::abs( y - y_m ) < y_cut * y_rms     &&
+                std::abs( z - z_m ) < z_cut * z_rms   ) {
                 XDim3 u = plasma_injector->getMomentum(x, y, z);
                 u.x *= PhysConst::c;
                 u.y *= PhysConst::c;
@@ -388,6 +393,9 @@ PhysicalParticleContainer::AddParticles (int lev)
                         plasma_injector->x_rms,
                         plasma_injector->y_rms,
                         plasma_injector->z_rms,
+                        plasma_injector->x_cut,
+                        plasma_injector->y_cut,
+                        plasma_injector->z_cut,
                         plasma_injector->q_tot,
                         plasma_injector->npart,
                         plasma_injector->do_symmetrize);