CD_RandomImplem.H

Go to the documentation of this file.

/*
 * SPDX-FileCopyrightText: 2021-2026 SINTEF Energy Research
 *
 * SPDX-License-Identifier: GPL-3.0-or-later
 */
 
#ifndef CD_RANDOMIMPLEM_H
#define CD_RANDOMIMPLEM_H
 
// Std includes
#include <chrono>
#ifdef _OPENMP
#include <omp.h>
#endif
 
// Chombo includes
#include <SPMD.H>
#include <CH_Timer.H>
#include <ParmParse.H>
 
// Our includes
#include <CD_Random.H>
#include <CD_NamespaceHeader.H>
 
inline void
Random::seed()
{
  if (!s_seeded) {
    ParmParse pp("Random");
 
    if (pp.contains("seed")) {
      int seed;
      pp.get("seed", seed);
      if (seed < 0) {
        Random::setRandomSeed();
      }
      else {
        Random::setSeed(seed);
      }
    }
    else {
      Random::setSeed(0);
    }
 
    s_seeded = true;
  }
}
 
inline void
Random::setSeed(const int a_seed)
{
#ifdef CH_MPI
#ifdef _OPENMP
#pragma omp parallel
  {
    const int seed = a_seed + procID() * omp_get_num_threads() + omp_get_thread_num();
 
    s_rng = std::mt19937_64(seed);
  }
#else
  const int seed = a_seed + procID();
 
  s_rng = std::mt19937_64(seed);
#endif
#else
#ifdef _OPENMP
#pragma omp parallel
  {
    const int seed = a_seed + omp_get_thread_num();
 
    s_rng = std::mt19937_64(seed);
  }
#else
  const int seed = a_seed;
 
  s_rng = std::mt19937_64(seed);
#endif
#endif
 
  s_seeded = true;
}
 
inline void
Random::setRandomSeed()
{
  auto seed = static_cast<int>(std::chrono::system_clock::now().time_since_epoch().count());
 
  // Special hook for MPI -- master rank broadcasts the seed and everyone increments by their processor ID.
#ifdef CH_MPI
  MPI_Bcast(&seed, 1, MPI_INT, 0, Chombo_MPI::comm);
#endif
 
  Random::setSeed(seed);
}
 
template <typename T, typename>
inline T
Random::getPoisson(const Real a_mean)
{
  CH_assert(s_seeded);
 
  T ret = (T)0;
 
  if (a_mean <= 0.0) {
    return ret;
  }
 
  if (a_mean < 250.0) {
    std::poisson_distribution<T> poisson(a_mean);
 
    ret = poisson(s_rng);
  }
  else {
    // Normal approximation for large means. The continuous draw is rounded to the nearest integer;
    // truncating it would bias every sample by -1/2.
    std::normal_distribution<Real> normal(a_mean, sqrt(a_mean));
 
    ret = (T)llround(std::max(normal(s_rng), (Real)0.0));
  }
 
  return ret;
}
 
template <typename T, typename>
inline T
Random::getBinomial(const T a_N, const Real a_p) noexcept
{
 
  CH_assert(s_seeded);
 
  T ret = 0;
 
  const bool useNormalApprox = (a_N > (T)9.0 * ((1.0 - a_p) / a_p)) && (a_N > (T)9.0 * a_p / (1.0 - a_p));
 
  if (useNormalApprox) {
    const Real mean = a_N * a_p;
 
    std::normal_distribution<Real> normalDist(mean, mean * (1.0 - a_p));
 
    ret = (T)std::max(normalDist(s_rng), (Real)0.0);
  }
  else {
    std::binomial_distribution<T> binomDist(a_N, a_p);
 
    ret = binomDist(s_rng);
  }
 
  return ret;
}
 
inline Real
Random::getUniformReal01()
{
  CH_assert(s_seeded);
 
  return s_uniform01(s_rng);
}
 
inline Real
Random::getUniformReal11()
{
  CH_assert(s_seeded);
 
  return s_uniform11(s_rng);
}
 
inline Real
Random::getNormal01()
{
  CH_assert(s_seeded);
 
  return s_normal01(s_rng);
}
 
inline RealVect
Random::getDirection()
{
  CH_assert(s_seeded);
 
  constexpr Real safety = 1.E-12;
 
#if CH_SPACEDIM == 2
  Real x1 = 2.0;
  Real x2 = 2.0;
  Real r  = x1 * x1 + x2 * x2;
  while (r >= 1.0 || r < safety) {
    x1 = Random::getUniformReal11();
    x2 = Random::getUniformReal11();
    r  = x1 * x1 + x2 * x2;
  }
 
  return RealVect(x1, x2) / sqrt(r);
#elif CH_SPACEDIM == 3
  Real x1 = 2.0;
  Real x2 = 2.0;
  Real r  = x1 * x1 + x2 * x2;
  while (r >= 1.0 || r < safety) {
    x1 = Random::getUniformReal11();
    x2 = Random::getUniformReal11();
    r  = x1 * x1 + x2 * x2;
  }
 
  const Real x = 2 * x1 * sqrt(1 - r);
  const Real y = 2 * x2 * sqrt(1 - r);
  const Real z = 1 - 2 * r;
 
  return RealVect(x, y, z);
#endif
}
 
template <typename T>
inline Real
Random::get(T& a_distribution)
{
  CH_assert(s_seeded);
 
  return a_distribution(s_rng);
}
 
template <typename T>
inline size_t
Random::getDiscrete(T& a_distribution)
{
  CH_assert(s_seeded);
 
  return a_distribution(s_rng);
}
 
inline RealVect
Random::randomPosition(const RealVect& a_cellPos,
                       const RealVect& a_lo,
                       const RealVect& a_hi,
                       const RealVect& a_bndryCentroid,
                       const RealVect& a_bndryNormal,
                       const Real      a_dx,
                       const Real      a_kappa) noexcept
{
 
  RealVect pos;
 
  if (a_kappa < 1.0) { // Rejection sampling.
    pos = Random::randomPosition(a_lo, a_hi, a_bndryCentroid, a_bndryNormal);
  }
  else { // Regular cell. Get a position.
    pos = Random::randomPosition(a_lo, a_hi);
  }
 
  // Convert from unit cell coordinates to physical coordinates.
  pos = a_cellPos + pos * a_dx;
 
  return pos;
}
 
inline RealVect
Random::randomPosition(const RealVect& a_lo,
                       const RealVect& a_hi,
                       const RealVect& a_bndryCentroid,
                       const RealVect& a_bndryNormal) noexcept
{
  constexpr int maxIter = 100;
 
  RealVect pos   = Random::randomPosition(a_lo, a_hi);
  bool     valid = (pos - a_bndryCentroid).dotProduct(a_bndryNormal) >= 0.0;
 
  for (int iter = 0; !valid && iter < maxIter; ++iter) {
    pos   = Random::randomPosition(a_lo, a_hi);
    valid = (pos - a_bndryCentroid).dotProduct(a_bndryNormal) >= 0.0;
  }
 
  // Fall back to the boundary centroid for degenerate cut-cells where the
  // valid half-space and the bounding box do not overlap.
  if (!valid) {
    pos = a_bndryCentroid;
  }
 
  return pos;
}
 
inline RealVect
Random::randomPosition(const RealVect& a_lo, const RealVect& a_hi) noexcept
{
 
  RealVect pos = RealVect::Zero;
 
  for (int dir = 0; dir < SpaceDim; dir++) {
    pos[dir] = a_lo[dir] + Random::getUniformReal01() * (a_hi[dir] - a_lo[dir]);
  }
 
  return pos;
}
 
#include <CD_NamespaceFooter.H>
 
#endif