Random numbers

Random is a static class for generating pseudo-random numbers, and exist so that all random number operations can be aggregated into a single class. Internally, Random use a Mersenne-Twister random number generation supplied by std::random.

To use the Random class, simply include <CD_Random.H>, e.g.,

#include <CD_Random.H>

See the Random API for further details.

Drawing random numbers

General approach

The general routine for drawing a random number is

/*!
  @brief For getting a random number from a user-supplied distribution. T must be a distribution for which we can call T(s_rng)
  @param[in] a_distribution Distribution. Must have object of type
*/
template <typename T>
inline static Real
get(T& a_distribution);

Here, the template parameter T is some distribution that follows the appropriate C++ template constraints of <random>. For example:

std::uniform_real_distribution<Real> dist(0.0, 100.0);

const Real randomNumber = Random::get(dist);

Pre-defined distributions

Pre-defined distributions exist for performing the following operations:

/*!
  @brief Get Poisson distributed number
  @param[in] a_mean Poisson mean value (inverse rate)
*/
template <typename T, typename = std::enable_if_t<std::is_integral<T>::value>>
inline static T
getPoisson(const Real a_mean);

/*!
  @brief Get Poisson distributed number. 
  @param[in] a_N Number of trials
  @param[in] a_p The usual success probability in binomial distributions
*/
template <typename T, typename = std::enable_if_t<std::is_integral<T>::value>>
inline static T
getBinomial(const T a_N, const Real a_p) noexcept;

/*!
  @brief Get a uniform real number on the interval [0,1]
*/
inline static Real
getUniformReal01();

/*!
  @brief Get a uniform real number on the interval [-1,1]
*/
inline static Real
getUniformReal11();

/*!
  @brief Get a number from a normal distribution centered on zero and variance 1
*/
inline static Real
getNormal01();

/*!
  @brief Get a random direction in space.
  @details Uses Marsaglia algorithm. 
*/
inline static RealVect
getDirection();

Seeding the RNG

By default, the random number generator is seeded using the MPI rank. It is necessary to seed MPI ranks using different seeds to avoid them producing the same number sequences. If using MPI+OpenMP, additional steps are taken to ensure that each OpenMP thread obtains a unique random number generator.

Driver (see Driver) will seed the random number generator, and user can override the seed by setting Random.seed = <number> in the input script. If the user sets <number> < 0 then a random seed will be produced based on the elapsed CPU clock time. If running with MPI, this seed is obtained by only one of the MPI ranks, and this seed is then broadcast to all the other ranks. The other ranks will then increment the seed by their own MPI rank number so that each MPI rank gets a unique seed.