diff --git a/src/include/binary_collision.hxx b/src/include/binary_collision.hxx index d8e995b221..f8dc5a738d 100644 --- a/src/include/binary_collision.hxx +++ b/src/include/binary_collision.hxx @@ -2,6 +2,7 @@ #pragma once #include "cuda_compat.h" +#include "rng.hxx" #include @@ -14,17 +15,9 @@ struct RngC // ---------------------------------------------------------------------- // uniform // - // returns random number in ]0:1] + // returns random number in [0:1[ - real_t uniform() - { - real_t ran; - do { - ran = real_t(random()) / real_t(RAND_MAX); - } while (ran == real_t(0.f)); - - return ran; - } + real_t uniform() { return rng::Uniform{}.get(); } }; // ====================================================================== diff --git a/src/libpsc/cuda/cuda_heating.cxx b/src/libpsc/cuda/cuda_heating.cxx index c64802c99c..f892578ff5 100644 --- a/src/libpsc/cuda/cuda_heating.cxx +++ b/src/libpsc/cuda/cuda_heating.cxx @@ -8,6 +8,7 @@ #include "heating_cuda_impl.hxx" #include "balance.hxx" #include "rng_state.hxx" +#include "rng.hxx" #include #include @@ -20,23 +21,6 @@ extern std::size_t mem_heating; using Float3 = Vec3; -// ---------------------------------------------------------------------- -// bm_normal2 - -static inline float2 bm_normal2(void) -{ - float u1, u2; - do { - u1 = random() * (1.f / RAND_MAX); - u2 = random() * (1.f / RAND_MAX); - } while (u1 <= 0.f); - - float2 rv; - rv.x = sqrtf(-2.f * logf(u1)) * cosf(2.f * M_PI * u2); - rv.y = sqrtf(-2.f * logf(u1)) * sinf(2.f * M_PI * u2); - return rv; -} - // ---------------------------------------------------------------------- // d_particle_kick @@ -176,14 +160,13 @@ struct cuda_heating_foil __host__ void particle_kick(DParticleCuda& prt, float H, float heating_dt) { - float2 r01 = bm_normal2(); - float2 r23 = bm_normal2(); + rng::Normal standard_normal_distr; float Dp = sqrtf(H * heating_dt); - prt.x[0] += Dp * r01.x; - prt.x[1] += Dp * r01.y; - prt.x[2] += Dp * r23.x; + prt.x[0] += Dp * standard_normal_distr.get(); + prt.x[1] += Dp * standard_normal_distr.get(); + prt.x[2] += Dp * standard_normal_distr.get(); } // ---------------------------------------------------------------------- diff --git a/src/libpsc/psc_bnd_particles/psc_bnd_particles_open.cxx b/src/libpsc/psc_bnd_particles/psc_bnd_particles_open.cxx index dc1dd68347..78976a820d 100644 --- a/src/libpsc/psc_bnd_particles/psc_bnd_particles_open.cxx +++ b/src/libpsc/psc_bnd_particles/psc_bnd_particles_open.cxx @@ -1,4 +1,5 @@ +#include "rng.hxx" #include "fields.hxx" using real_t = mparticles_t::real_t; @@ -8,8 +9,6 @@ static const int debug_every_step = 10; static inline bool at_lo_boundary(int p, int d); static inline bool at_hi_boundary(int p, int d); -static inline double random1() { return random() / (double)RAND_MAX; } - static void copy_to_mrc_fld(struct mrc_fld* m3, struct psc_mfields* mflds) { mfields_t mf(mflds); @@ -472,8 +471,9 @@ static inline double inject_particles(int p, struct psc_mparticles* mprts, do { nnm++; // FIXME, shouldn't have to loop here + rng::Uniform uniform_distr; do { - double sr = random1(); + double sr = uniform_distr.get(); pxi[Z] = 0.; for (int k = 0; k < nvdx - 1; k++) { @@ -485,26 +485,24 @@ static inline double inject_particles(int p, struct psc_mparticles* mprts, } } while (pxi[Z] == 0); - double sr = random1(); double yya = 0.; double yy0; int icount = 0; do { icount++; yy0 = yya; - yya = yy0 - (erf(yy0) - (2. * sr - 1.)) / + yya = yy0 - (erf(yy0) - (2. * uniform_distr.get() - 1.)) / (2. / sqrt(M_PI) * exp(-sqr(yy0))); } while (fabs(yya - yy0) > 1.e-15 && icount != 100); pxi[X] = v[X] + yya * sqrt(W[YY] / (W[XX] * W[YY] - sqr(W[XY]))) + (pxi[Z] - v[Z]) * vv[ZX] / vv[ZZ]; - sr = random1(); yya = 0.0; icount = 0; do { icount++; yy0 = yya; - yya = yy0 - (erf(yy0) - (2. * sr - 1.)) / + yya = yy0 - (erf(yy0) - (2. * uniform_distr.get() - 1.)) / (2. / sqrt(M_PI) * exp(-sqr(yy0))); } while (fabs(yya - yy0) > 1.e-15 && icount != 100); pxi[Y] = v[Y] + 1. / W[YY] * @@ -517,11 +515,11 @@ static inline double inject_particles(int p, struct psc_mparticles* mprts, } } while (sqr(pxi[X]) + sqr(pxi[Y]) + sqr(pxi[Z]) > 1.); - double xr = random1(); + double xr = uniform_distr.get(); xi[X] = pos[X] + xr * grid.dx[X]; - double yr = random1(); + double yr = uniform_distr.get(); xi[Y] = pos[Y] + yr * grid.dx[Y]; - double zr = random1(); + double zr = uniform_distr.get(); double dz = dir * zr * pxi[Z] * ppsc->dt; xi[Z] = pos[Z] + dz; diff --git a/src/libpsc/psc_heating/psc_heating_impl.hxx b/src/libpsc/psc_heating/psc_heating_impl.hxx index f82c1a1142..a9b590e7af 100644 --- a/src/libpsc/psc_heating/psc_heating_impl.hxx +++ b/src/libpsc/psc_heating/psc_heating_impl.hxx @@ -1,5 +1,6 @@ #include "heating.hxx" +#include "rng.hxx" #include #include @@ -28,28 +29,15 @@ struct Heating__ : HeatingBase void kick_particle(Particle& prt, real_t H) { - float ran1, ran2, ran3, ran4, ran5, ran6; - do { - ran1 = random() / ((float)RAND_MAX + 1); - ran2 = random() / ((float)RAND_MAX + 1); - ran3 = random() / ((float)RAND_MAX + 1); - ran4 = random() / ((float)RAND_MAX + 1); - ran5 = random() / ((float)RAND_MAX + 1); - ran6 = random() / ((float)RAND_MAX + 1); - } while (ran1 >= 1.f || ran2 >= 1.f || ran3 >= 1.f || ran4 >= 1.f || - ran5 >= 1.f || ran6 >= 1.f); - - real_t ranx = sqrtf(-2.f * logf(1.0 - ran1)) * cosf(2.f * M_PI * ran2); - real_t rany = sqrtf(-2.f * logf(1.0 - ran3)) * cosf(2.f * M_PI * ran4); - real_t ranz = sqrtf(-2.f * logf(1.0 - ran5)) * cosf(2.f * M_PI * ran6); + rng::Normal standard_normal_distr; real_t Dpxi = sqrtf(H * heating_dt_); real_t Dpyi = sqrtf(H * heating_dt_); real_t Dpzi = sqrtf(H * heating_dt_); - prt.u[0] += Dpxi * ranx; - prt.u[1] += Dpyi * rany; - prt.u[2] += Dpzi * ranz; + prt.u[0] += Dpxi * standard_normal_distr.get(); + prt.u[1] += Dpyi * standard_normal_distr.get(); + prt.u[2] += Dpzi * standard_normal_distr.get(); } void operator()(Mparticles& mprts)