Fix #2522: Gaussian beam positions do not change with warpx.random_seed (#2523)

* Draw Gaussian beam position with amrex random engine

* Update benchmarks

* Update tolerance in space-charge tests

* Update benchmark for space charge initialization test

* Update benchmarks

* Update benchmark

* Clean-up code

* Update benchmarks

1 files changed, 5 insertions, 10 deletions

diff --git a/Source/Particles/PhysicalParticleContainer.cpp b/Source/Particles/PhysicalParticleContainer.cpp
index 9e795dfd4..c331b7a61 100644
--- a/Source/Particles/PhysicalParticleContainer.cpp
+++ b/Source/Particles/PhysicalParticleContainer.cpp
@@ -408,11 +408,6 @@ PhysicalParticleContainer::AddGaussianBeam (
     const Real q_tot, long npart,
     const int do_symmetrize) {
 
-    std::mt19937_64 mt(0451);
-    std::normal_distribution<double> distx(x_m, x_rms);
-    std::normal_distribution<double> disty(y_m, y_rms);
-    std::normal_distribution<double> distz(z_m, z_rms);
-
     // Declare temporary vectors on the CPU
     Gpu::HostVector<ParticleReal> particle_x;
     Gpu::HostVector<ParticleReal> particle_y;
@@ -432,14 +427,14 @@ PhysicalParticleContainer::AddGaussianBeam (
         for (long i = 0; i < npart; ++i) {
 #if (defined WARPX_DIM_3D) || (defined WARPX_DIM_RZ)
             const Real weight = q_tot/(npart*charge);
-            const Real x = distx(mt);
-            const Real y = disty(mt);
-            const Real z = distz(mt);
+            const Real x = amrex::RandomNormal(x_m, x_rms);
+            const Real y = amrex::RandomNormal(y_m, y_rms);
+            const Real z = amrex::RandomNormal(z_m, z_rms);
 #elif (defined WARPX_DIM_XZ)
             const Real weight = q_tot/(npart*charge*y_rms);
-            const Real x = distx(mt);
+            const Real x = amrex::RandomNormal(x_m, x_rms);
             constexpr Real y = 0._prt;
-            const Real z = distz(mt);
+            const Real z = amrex::RandomNormal(z_m, z_rms);
 #endif
             if (plasma_injector->insideBounds(x, y, z)  &&
                 std::abs( x - x_m ) < x_cut * x_rms     &&