Allow negative mean velocity with NFluxPerCell (ECP-WarpX#4397)

* Allow negative drift velocity with Gaussian flux injection

* Simplify by using sign of u_m

* Add negative u_m test to flux_injection

* Update FluxInjection3D benchmark

* Save abs(u_m) as a temporary const

Examples/Tests/flux_injection/analysis_flux_injection_3d.py

@@ -57,11 +57,8 @@ def gaussian_dist(u, u_th):
     return 1./((2*np.pi)**.5*u_th) * np.exp(-u**2/(2*u_th**2) )
 
 def gaussian_flux_dist(u, u_th, u_m):
-    au_m = np.abs(u_m)
-    normalization_factor = u_th**2 * np.exp(-au_m**2/(2*u_th**2)) + (np.pi/2)**.5*au_m*u_th * (1 + erf(au_m/(2**.5*u_th)))
-    result = 1./normalization_factor * np.where( u>0, u * np.exp(-(u-au_m)**2/(2*u_th**2)), 0 )
-    if u_m < 0.:
-        result = result[::-1]
+    normalization_factor = u_th**2 * np.exp(-u_m**2/(2*u_th**2)) + (np.pi/2)**.5*u_m*u_th * (1 + erf(u_m/(2**.5*u_th)))
+    result = 1./normalization_factor * np.where( u>0, u * np.exp(-(u-u_m)**2/(2*u_th**2)), 0 )
     return result
 
 def compare_gaussian(u, w, u_th, label=''):
@@ -84,10 +81,10 @@ def compare_gaussian_flux(u, w, u_th, u_m, label=''):
 
 # Load data and perform check
 
-plt.figure(figsize=(5,7))
+plt.figure(figsize=(8,7))
 
-plt.subplot(211)
-plt.title('Electrons')
+plt.subplot(221)
+plt.title('Electrons u_m=0.07')
 
 ux = ad['electron','particle_momentum_x'].to_ndarray()/(m_e*c)
 uy = ad['electron','particle_momentum_y'].to_ndarray()/(m_e*c)
@@ -97,21 +94,46 @@ def compare_gaussian_flux(u, w, u_th, u_m, label=''):
 compare_gaussian(ux, w, u_th=0.1, label='u_x')
 compare_gaussian_flux(uy, w, u_th=0.1, u_m=0.07, label='u_y')
 compare_gaussian(uz, w, u_th=0.1, label='u_z')
-plt.legend(loc=0)
 
-plt.subplot(212)
-plt.title('Protons')
+plt.subplot(223)
+plt.title('Protons u_m=0.05')
 
 ux = ad['proton','particle_momentum_x'].to_ndarray()/(m_p*c)
 uy = ad['proton','particle_momentum_y'].to_ndarray()/(m_p*c)
 uz = ad['proton','particle_momentum_z'].to_ndarray()/(m_p*c)
 w = ad['proton', 'particle_weight'].to_ndarray()
 
-compare_gaussian_flux(ux, w, u_th=0.1, u_m=-0.05, label='u_x')
+compare_gaussian_flux(-ux, w, u_th=0.1, u_m=0.05, label='u_x')
 compare_gaussian(uy, w, u_th=0.1, label='u_y')
 compare_gaussian(uz, w, u_th=0.1, label='u_z')
-plt.legend(loc=0)
 
+plt.subplot(222)
+plt.title('Electrons u_m=-0.07')
+
+ux = ad['electron_negative','particle_momentum_x'].to_ndarray()/(m_e*c)
+uy = ad['electron_negative','particle_momentum_y'].to_ndarray()/(m_e*c)
+uz = ad['electron_negative','particle_momentum_z'].to_ndarray()/(m_e*c)
+w = ad['electron_negative', 'particle_weight'].to_ndarray()
+
+compare_gaussian(ux, w, u_th=0.1, label='u_x')
+compare_gaussian(uy, w, u_th=0.1, label='u_y')
+compare_gaussian_flux(uz, w, u_th=0.1, u_m=-0.07, label='u_z')
+plt.legend(loc=(1.02, 0.5))
+
+plt.subplot(224)
+plt.title('Protons u_m=-0.05')
+
+ux = ad['proton_negative','particle_momentum_x'].to_ndarray()/(m_p*c)
+uy = ad['proton_negative','particle_momentum_y'].to_ndarray()/(m_p*c)
+uz = ad['proton_negative','particle_momentum_z'].to_ndarray()/(m_p*c)
+w = ad['proton_negative', 'particle_weight'].to_ndarray()
+
+compare_gaussian(ux, w, u_th=0.1, label='u_x')
+compare_gaussian(uy, w, u_th=0.1, label='u_y')
+compare_gaussian_flux(-uz, w, u_th=0.1, u_m=-0.05, label='u_z')
+#plt.legend(loc=0)
+
+plt.tight_layout()
 plt.savefig('Distribution.png')
 
 # Verify checksum

Examples/Tests/flux_injection/inputs_3d

@@ -28,7 +28,7 @@ boundary.field_lo = periodic periodic periodic
 boundary.field_hi = periodic periodic periodic
 
 # particles
-particles.species_names = electron proton
+particles.species_names = electron proton electron_negative proton_negative
 algo.particle_shape = 3
 
 electron.charge = -q_e
@@ -61,6 +61,37 @@ proton.ux_m = 0.05
 proton.uy_th = 0.1
 proton.uz_th = 0.1
 
+# "negative" is negative u_m
+electron_negative.charge = -q_e
+electron_negative.mass = m_e
+electron_negative.injection_style = NFluxPerCell
+electron_negative.num_particles_per_cell = 100
+electron_negative.surface_flux_pos = -4.
+electron_negative.flux_normal_axis = z
+electron_negative.flux_direction = +1
+electron_negative.flux_profile = parse_flux_function
+electron_negative.flux_function(x,y,z,t) = "1."
+electron_negative.momentum_distribution_type = gaussianflux
+electron_negative.ux_th = 0.1
+electron_negative.uy_th = 0.1
+electron_negative.uz_th = 0.1
+electron_negative.uz_m = -0.07
+
+proton_negative.charge = +q_e
+proton_negative.mass = m_p
+proton_negative.injection_style = NFluxPerCell
+proton_negative.num_particles_per_cell = 100
+proton_negative.surface_flux_pos = 4.
+proton_negative.flux_normal_axis = z
+proton_negative.flux_direction = -1
+proton_negative.flux_profile = constant
+proton_negative.flux = 1.
+proton_negative.momentum_distribution_type = gaussianflux
+proton_negative.ux_th = 0.1
+proton_negative.uy_th = 0.1
+proton_negative.uz_th = 0.1
+proton_negative.uz_m = -0.05
+
 # Diagnostics
 diagnostics.diags_names = diag1
 diag1.intervals = 1000

Regression/Checksum/benchmarks_json/FluxInjection3D.json

@@ -1,21 +1,39 @@
 {
-  "electron": {
-    "particle_momentum_x": 1.1192116199394354e-18,
-    "particle_momentum_y": 2.238114590066897e-18,
-    "particle_momentum_z": 1.1156457989239732e-18,
-    "particle_position_x": 102495.14197173176,
-    "particle_position_y": 188132.22608016344,
-    "particle_position_z": 102423.13701045913,
+  "lev=0": {},
+  "electron_negative": {
+    "particle_momentum_x": 1.1222699783863554e-18,
+    "particle_momentum_y": 1.1202176725070554e-18,
+    "particle_momentum_z": 1.3925955132362978e-18,
+    "particle_position_x": 102352.09026544492,
+    "particle_position_y": 102418.88243172191,
+    "particle_position_z": 194298.7949373403,
     "particle_weight": 8.959999999999998e-07
   },
-  "lev=0": {},
   "proton": {
-    "particle_momentum_x": 3.835423016604918e-15,
-    "particle_momentum_y": 2.0468371931479925e-15,
-    "particle_momentum_z": 2.055186547721331e-15,
-    "particle_position_x": 189256.1546041931,
-    "particle_position_y": 102293.00576740496,
-    "particle_position_z": 102314.93877691089,
+    "particle_momentum_x": 3.8338884590187296e-15,
+    "particle_momentum_y": 2.0442156829943128e-15,
+    "particle_momentum_z": 2.045804260395492e-15,
+    "particle_position_x": 189238.69249885075,
+    "particle_position_y": 102242.91543133644,
+    "particle_position_z": 102297.92915049737,
+    "particle_weight": 8.959999999999998e-07
+  },
+  "electron": {
+    "particle_momentum_x": 1.1150196665556376e-18,
+    "particle_momentum_y": 2.2311586451156107e-18,
+    "particle_momentum_z": 1.1115069298757383e-18,
+    "particle_position_x": 102477.68142952351,
+    "particle_position_y": 188137.20906834095,
+    "particle_position_z": 102443.44417709563,
+    "particle_weight": 8.959999999999998e-07
+  },
+  "proton_negative": {
+    "particle_momentum_x": 2.051429985038282e-15,
+    "particle_momentum_y": 2.053711846305655e-15,
+    "particle_momentum_z": 2.7212135003240815e-15,
+    "particle_position_x": 102307.73649638034,
+    "particle_position_y": 102475.13878406698,
+    "particle_position_z": 193638.89296895845,
     "particle_weight": 8.959999999999998e-07
   }
 }

Source/Initialization/InjectorMomentum.H

@@ -104,37 +104,40 @@ namespace {
         // Momentum to be returned at the end of this function
         amrex::Real u = 0._rt;
 
+        const amrex::Real abs_u_m = std::abs(u_m);
+
         if (u_th == 0._rt) {
             u = u_m; // Trivial case ; avoids division by 0 in the rest of the code below
-        } else if (u_m < 0.6*u_th) {
-            // Mean velocity is lower than thermal velocity
-            // Use the distribution u*exp(-u**2*(1-u_m/u_th)/(2*u_th**2)) as an approximation
+        } else if (abs_u_m < 0.6*u_th) {
+            // Mean velocity magnitude is less than thermal velocity
+            // Use the distribution u*exp(-u**2*(1-abs(u_m)/u_th)/(2*u_th**2)) as an approximation
             // and then use the rejection method to correct it
-            // ( stop rejecting with probability exp(-u_m/(2*u_th**3)*(u-u_th)**2) )
+            // ( stop rejecting with probability exp(-abs(u_m)/(2*u_th**3)*(u-sign(u_m)*u_th)**2) )
             // Note that this is the method that is used in the common case u_m=0
-            const amrex::Real approx_u_th = u_th/std::sqrt( 1._rt - u_m/u_th );
-            const amrex::Real reject_prefactor = (u_m/u_th)/(2._rt*u_th*u_th); // To save computation
+            const amrex::Real umsign = std::copysign(1._rt, u_m);
+            const amrex::Real approx_u_th = u_th/std::sqrt( 1._rt - abs_u_m/u_th );
+            const amrex::Real reject_prefactor = (abs_u_m/u_th)/(2._rt*u_th*u_th); // To save computation
             bool reject = true;
             while (reject) {
-                // Generates u according to u*exp(-u**2/(2*approx_u_th**2),
+                // Generates u according to u*exp(-u**2/(2*approx_u_th**2)),
                 // using the method of the inverse cumulative function
                 amrex::Real xrand = 1._rt - amrex::Random(engine); // ensures urand > 0
                 u = approx_u_th * std::sqrt(2._rt*std::log(1._rt/xrand));
                 // Rejection method
                 xrand = amrex::Random(engine);
-                if (xrand < std::exp(-reject_prefactor*(u-u_th)*(u-u_th))) reject = false;
+                if (xrand < std::exp(-reject_prefactor*(u - umsign*u_th)*(u - umsign*u_th))) reject = false;
             }
         } else {
-            // Mean velocity is greater than thermal velocity
-            // Use the distribution exp(-(u-u_m-u_th**2/u_m)**2/(2*u_th**2)) as an approximation
+            // Mean velocity magnitude is greater than thermal velocity
+            // Use the distribution exp(-(u-u_m-u_th**2/abs(u_m))**2/(2*u_th**2)) as an approximation
             // and then use the rejection method to correct it
-            // ( stop rejecting with probability (u/u_m)*exp(1-(u/u_m)) ; note
+            // ( stop rejecting with probability (u/abs(u_m))*exp(1-(u/abs(u_m))) ; note
             // that this number is always between 0 and 1 )
             // Note that in the common case `u_m = 0`, this rejection method
             // is not used, and the above rejection method is used instead.
             bool reject = true;
-            const amrex::Real approx_u_m = u_m + u_th*u_th/u_m;
-            const amrex::Real inv_um = 1._rt/u_m; // To save computation
+            const amrex::Real approx_u_m = u_m + u_th*u_th/abs_u_m;
+            const amrex::Real inv_um = 1._rt/abs_u_m; // To save computation
             while (reject) {
                 // Approximate distribution: normal distribution, where we only retain positive u
                 u = -1._rt;
@@ -166,13 +169,6 @@ struct InjectorMomentumGaussianFlux
           m_flux_normal_axis(a_flux_normal_axis),
           m_flux_direction(a_flux_direction)
     {
-        // For now, do not allow negative `u_m` along the flux axis
-        bool raise_error = false;
-        if ((m_flux_normal_axis == 0) && (m_ux_m < 0)) raise_error = true;
-        if ((m_flux_normal_axis == 1) && (m_uy_m < 0)) raise_error = true;
-        if ((m_flux_normal_axis == 2) && (m_uz_m < 0)) raise_error = true;
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE( raise_error==false,
-            "When using the `gaussianflux` distribution, the central momentum along the flux axis must be positive or zero." );
     }
 
     AMREX_GPU_HOST_DEVICE