Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft: Use new rng in more places #317

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

Draft: Use new rng in more places #317

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2e6ad9a
psc_heating_impl: use rng::Normal
May 24, 2023
973aa05
cuda_heating: use rng::Normal
May 24, 2023
461ff81
binary_collision: use rng::Uniform
May 24, 2023
8270b02
psc_bnd_particles_open: use rng::Uniform
May 24, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
13 changes: 3 additions & 10 deletions src/include/binary_collision.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#pragma once

#include "cuda_compat.h"
#include "rng.hxx"

#include <cmath>

Expand All @@ -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()) / RAND_MAX;
} while (ran == real_t(0.f));

return ran;
}
real_t uniform() { return rng::Uniform<real_t>{}.get(); }
};

// ======================================================================
Expand Down
27 changes: 5 additions & 22 deletions src/libpsc/cuda/cuda_heating.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "heating_cuda_impl.hxx"
#include "balance.hxx"
#include "rng_state.hxx"
#include "rng.hxx"

#include <thrust/host_vector.h>
#include <thrust/device_vector.h>
Expand All @@ -20,23 +21,6 @@ extern std::size_t mem_heating;

using Float3 = Vec3<float>;

// ----------------------------------------------------------------------
// 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

Expand Down Expand Up @@ -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<float> 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();
}

// ----------------------------------------------------------------------
Expand Down
18 changes: 8 additions & 10 deletions src/libpsc/psc_bnd_particles/psc_bnd_particles_open.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@

#include "rng.hxx"
#include "fields.hxx"

using real_t = mparticles_t::real_t;
Expand All @@ -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);
Expand Down Expand Up @@ -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<double> uniform_distr;
do {
double sr = random1();
double sr = uniform_distr.get();

pxi[Z] = 0.;
for (int k = 0; k < nvdx - 1; k++) {
Expand All @@ -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] *
Expand All @@ -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;

Expand Down
22 changes: 5 additions & 17 deletions src/libpsc/psc_heating/psc_heating_impl.hxx
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@

#include "heating.hxx"
#include "rng.hxx"

#include <functional>
#include <stdlib.h>
Expand Down Expand Up @@ -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<real_t> 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)
Expand Down
Loading