From 0602ddd8a36a75eed0b5d202ec4be1e4ca7b8c0c Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 5 Dec 2019 19:50:31 +0000 Subject: [PATCH] refactor seed and seed aray variable names --- include/openmc/angle_energy.h | 2 +- include/openmc/distribution.h | 34 +++++++-------- include/openmc/distribution_angle.h | 4 +- include/openmc/distribution_energy.h | 26 +++++------ include/openmc/distribution_multi.h | 16 +++---- include/openmc/distribution_spatial.h | 18 ++++---- include/openmc/math_functions.h | 22 +++++----- include/openmc/mgxs.h | 8 ++-- include/openmc/particle.h | 4 +- include/openmc/photon.h | 10 ++--- include/openmc/physics.h | 6 +-- include/openmc/random_lcg.h | 18 ++++---- include/openmc/reaction_product.h | 4 +- include/openmc/scattdata.h | 14 +++--- include/openmc/secondary_correlated.h | 4 +- include/openmc/secondary_kalbach.h | 4 +- include/openmc/secondary_nbody.h | 4 +- include/openmc/secondary_thermal.h | 20 ++++----- include/openmc/secondary_uncorrelated.h | 4 +- include/openmc/source.h | 8 ++-- include/openmc/surface.h | 2 +- include/openmc/thermal.h | 8 ++-- src/distribution.cpp | 32 +++++++------- src/distribution_angle.cpp | 6 +-- src/distribution_energy.cpp | 22 +++++----- src/distribution_multi.cpp | 16 +++---- src/distribution_spatial.cpp | 18 ++++---- src/math_functions.cpp | 34 +++++++-------- src/mgxs.cpp | 10 ++--- src/particle.cpp | 4 +- src/particle_restart.cpp | 2 +- src/photon.cpp | 58 ++++++++++++------------- src/physics.cpp | 46 ++++++++++---------- src/plot.cpp | 8 ++-- src/random_lcg.cpp | 18 ++++---- src/reaction_product.cpp | 8 ++-- src/scattdata.cpp | 26 +++++------ src/secondary_correlated.cpp | 12 ++--- src/secondary_kalbach.cpp | 14 +++--- src/secondary_nbody.cpp | 26 +++++------ src/secondary_thermal.cpp | 30 ++++++------- src/secondary_uncorrelated.cpp | 8 ++-- src/simulation.cpp | 2 +- src/source.cpp | 14 +++--- src/surface.cpp | 6 +-- src/thermal.cpp | 12 ++--- 46 files changed, 336 insertions(+), 336 deletions(-) diff --git a/include/openmc/angle_energy.h b/include/openmc/angle_energy.h index c88acd569..e2f449c6a 100644 --- a/include/openmc/angle_energy.h +++ b/include/openmc/angle_energy.h @@ -15,7 +15,7 @@ namespace openmc { class AngleEnergy { public: virtual void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const = 0; + uint64_t* seed) const = 0; virtual ~AngleEnergy() = default; }; diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index d62ba1624..f9d718cdb 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -21,7 +21,7 @@ namespace openmc { class Distribution { public: virtual ~Distribution() = default; - virtual double sample(uint64_t* prn_seed) const = 0; + virtual double sample(uint64_t* seed) const = 0; }; //============================================================================== @@ -34,9 +34,9 @@ public: Discrete(const double* x, const double* p, int n); //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; // Properties const std::vector& x() const { return x_; } @@ -59,9 +59,9 @@ public: Uniform(double a, double b) : a_{a}, b_{b} {}; //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; private: double a_; //!< Lower bound of distribution double b_; //!< Upper bound of distribution @@ -77,9 +77,9 @@ public: Maxwell(double theta) : theta_{theta} { }; //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; private: double theta_; //!< Factor in exponential [eV] }; @@ -94,9 +94,9 @@ public: Watt(double a, double b) : a_{a}, b_{b} { }; //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; private: double a_; //!< Factor in exponential [eV] double b_; //!< Factor in square root [1/eV] @@ -112,9 +112,9 @@ public: Normal(double mean_value, double std_dev) : mean_value_{mean_value}, std_dev_{std_dev} { }; //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; private: double mean_value_; //!< middle of distribution [eV] double std_dev_; //!< standard deviation [eV] @@ -131,9 +131,9 @@ public: Muir(double e0, double m_rat, double kt) : e0_{e0}, m_rat_{m_rat}, kt_{kt} { }; //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; private: // example DT fusion m_rat = 5 (D = 2 + T = 3) // ion temp = 20000 eV @@ -154,9 +154,9 @@ public: const double* c=nullptr); //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; // x property std::vector& x() { return x_; } @@ -185,9 +185,9 @@ public: Equiprobable(const double* x, int n) : x_{x, x+n} { }; //! Sample a value from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled value - double sample(uint64_t* prn_seed) const; + double sample(uint64_t* seed) const; private: std::vector x_; //! Possible outcomes }; diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h index 23d7aa008..d62236c29 100644 --- a/include/openmc/distribution_angle.h +++ b/include/openmc/distribution_angle.h @@ -23,9 +23,9 @@ public: //! Sample an angle given an incident particle energy //! \param[in] E Particle energy in [eV] - //! \param[inout] prn_seed pseudorandom number seed pointer + //! \param[inout] seed pseudorandom number seed pointer //! \return Cosine of the angle in the range [-1,1] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; //! Determine whether angle distribution is empty //! \return Whether distribution is empty diff --git a/include/openmc/distribution_energy.h b/include/openmc/distribution_energy.h index dd59a2359..7b739ba94 100644 --- a/include/openmc/distribution_energy.h +++ b/include/openmc/distribution_energy.h @@ -22,7 +22,7 @@ namespace openmc { class EnergyDistribution { public: - virtual double sample(double E, uint64_t* prn_seed) const = 0; + virtual double sample(double E, uint64_t* seed) const = 0; virtual ~EnergyDistribution() = default; }; @@ -36,9 +36,9 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer //! \return Sampled energy in [eV] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; private: int primary_flag_; //!< Indicator of whether the photon is a primary or //!< non-primary photon. @@ -56,9 +56,9 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer //! \return Sampled energy in [eV] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; private: double threshold_; //!< Energy threshold in lab, (A + 1)/A * |Q| double mass_ratio_; //!< (A/(A+1))^2 @@ -76,9 +76,9 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer //! \return Sampled energy in [eV] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; private: //! Outgoing energy for a single incoming energy struct CTTable { @@ -106,9 +106,9 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer //! \return Sampled energy in [eV] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; private: Tabulated1D theta_; //!< Incoming energy dependent parameter double u_; //!< Restriction energy @@ -125,9 +125,9 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer //! \return Sampled energy in [eV] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; private: Tabulated1D theta_; //!< Incoming energy dependent parameter double u_; //!< Restriction energy @@ -144,9 +144,9 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer //! \return Sampled energy in [eV] - double sample(double E, uint64_t* prn_seed) const; + double sample(double E, uint64_t* seed) const; private: Tabulated1D a_; //!< Energy-dependent 'a' parameter Tabulated1D b_; //!< Energy-dependent 'b' parameter diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h index 08f1cd0b5..493ad85e0 100644 --- a/include/openmc/distribution_multi.h +++ b/include/openmc/distribution_multi.h @@ -23,9 +23,9 @@ public: virtual ~UnitSphereDistribution() = default; //! Sample a direction from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Direction sampled - virtual Direction sample(uint64_t* prn_seed) const = 0; + virtual Direction sample(uint64_t* seed) const = 0; Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction }; @@ -40,9 +40,9 @@ public: explicit PolarAzimuthal(pugi::xml_node node); //! Sample a direction from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Direction sampled - Direction sample(uint64_t* prn_seed) const; + Direction sample(uint64_t* seed) const; private: UPtrDist mu_; //!< Distribution of polar angle UPtrDist phi_; //!< Distribution of azimuthal angle @@ -57,9 +57,9 @@ public: Isotropic() { }; //! Sample a direction from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled direction - Direction sample(uint64_t* prn_seed) const; + Direction sample(uint64_t* seed) const; }; //============================================================================== @@ -72,9 +72,9 @@ public: explicit Monodirectional(pugi::xml_node node) : UnitSphereDistribution{node} { }; //! Sample a direction from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled direction - Direction sample(uint64_t* prn_seed) const; + Direction sample(uint64_t* seed) const; }; using UPtrAngle = std::unique_ptr; diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 4895dd6ae..53857fac3 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -17,7 +17,7 @@ public: virtual ~SpatialDistribution() = default; //! Sample a position from the distribution - virtual Position sample(uint64_t* prn_seed) const = 0; + virtual Position sample(uint64_t* seed) const = 0; }; //============================================================================== @@ -29,9 +29,9 @@ public: explicit CartesianIndependent(pugi::xml_node node); //! Sample a position from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled position - Position sample(uint64_t* prn_seed) const; + Position sample(uint64_t* seed) const; private: UPtrDist x_; //!< Distribution of x coordinates UPtrDist y_; //!< Distribution of y coordinates @@ -47,9 +47,9 @@ public: explicit SphericalIndependent(pugi::xml_node node); //! Sample a position from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled position - Position sample(uint64_t* prn_seed) const; + Position sample(uint64_t* seed) const; private: UPtrDist r_; //!< Distribution of r coordinates UPtrDist theta_; //!< Distribution of theta coordinates @@ -66,9 +66,9 @@ public: explicit SpatialBox(pugi::xml_node node, bool fission=false); //! Sample a position from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled position - Position sample(uint64_t* prn_seed) const; + Position sample(uint64_t* seed) const; // Properties bool only_fissionable() const { return only_fissionable_; } @@ -89,9 +89,9 @@ public: explicit SpatialPoint(pugi::xml_node node); //! Sample a position from the distribution - //! \param prn_seed Pseudorandom number seed pointer + //! \param seed Pseudorandom number seed pointer //! \return Sampled position - Position sample(uint64_t* prn_seed) const; + Position sample(uint64_t* seed) const; private: Position r_; //!< Single position at which sites are generated }; diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h index 1440d13b9..e2d776dc2 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -129,14 +129,14 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]); //! \param mu The cosine of angle in lab or CM //! \param phi The azimuthal angle; will randomly chosen angle if a nullptr //! is passed -//! \param prn_seed A pointer to the pseudorandom seed +//! \param seed A pointer to the pseudorandom seed //============================================================================== extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi, - uint64_t* prn_seed); + uint64_t* seed); Direction rotate_angle(Direction u, double mu, const double* phi, - uint64_t* prn_seed); + uint64_t* seed); //============================================================================== //! Samples an energy from the Maxwell fission distribution based on a direct @@ -147,11 +147,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi, //! rule C64 in the Monte Carlo Sampler LA-9721-MS. //! //! \param T The tabulated function of the incoming energy -//! \param prn_seed A pointer to the pseudorandom seed +//! \param seed A pointer to the pseudorandom seed //! \return The sampled outgoing energy //============================================================================== -extern "C" double maxwell_spectrum(double T, uint64_t* prn_seed); +extern "C" double maxwell_spectrum(double T, uint64_t* seed); //============================================================================== //! Samples an energy from a Watt energy-dependent fission distribution. @@ -163,11 +163,11 @@ extern "C" double maxwell_spectrum(double T, uint64_t* prn_seed); //! //! \param a Watt parameter a //! \param b Watt parameter b -//! \param prn_seed A pointer to the pseudorandom seed +//! \param seed A pointer to the pseudorandom seed //! \return The sampled outgoing energy //============================================================================== -extern "C" double watt_spectrum(double a, double b, uint64_t* prn_seed); +extern "C" double watt_spectrum(double a, double b, uint64_t* seed); //============================================================================== //! Samples an energy from the Gaussian energy-dependent fission distribution. @@ -180,11 +180,11 @@ extern "C" double watt_spectrum(double a, double b, uint64_t* prn_seed); //! //! @param mean mean of the Gaussian distribution //! @param std_dev standard deviation of the Gaussian distribution -//! @param prn_seed A pointer to the pseudorandom seed +//! @param seed A pointer to the pseudorandom seed //! @result The sampled outgoing energy //============================================================================== -extern "C" double normal_variate(double mean, double std_dev, uint64_t* prn_seed); +extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed); //============================================================================== //! Samples an energy from the Muir (Gaussian) energy-dependent distribution. @@ -196,12 +196,12 @@ extern "C" double normal_variate(double mean, double std_dev, uint64_t* prn_seed //! @param e0 peak neutron energy [eV] //! @param m_rat ratio of the fusion reactants to AMU //! @param kt the ion temperature of the reactants [eV] -//! @param prn_seed A pointer to the pseudorandom seed +//! @param seed A pointer to the pseudorandom seed //! @result The sampled outgoing energy //============================================================================== extern "C" double muir_spectrum(double e0, double m_rat, double kt, - uint64_t* prn_seed); + uint64_t* seed); //============================================================================== //! Doppler broadens the windowed multipole curvefit. diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 131d248e1..b7dc374fd 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -154,9 +154,9 @@ class Mgxs { //! @param gin Incoming energy group. //! @param dg Sampled delayed group index. //! @param gout Sampled outgoing energy group. - //! @param prn_seed Pseudorandom seed pointer + //! @param seed Pseudorandom seed pointer void - sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed); + sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed); //! \brief Samples the outgoing energy and angle from a scatter event. //! @@ -164,9 +164,9 @@ class Mgxs { //! @param gout Sampled outgoing energy group. //! @param mu Sampled cosine of the change-in-angle. //! @param wgt Weight of the particle to be adjusted. - //! @param prn_seed Pseudorandom seed pointer. + //! @param seed Pseudorandom seed pointer. void - sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed); + sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* seed); //! \brief Calculates cross section quantities needed for tracking. //! diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 00476e3ef..3888a847d 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -292,8 +292,8 @@ public: bool write_track_ {false}; // Current PRNG state - uint64_t prn_seeds_[N_STREAMS]; // current seeds - int stream_; // current RNG stream + uint64_t seeds_[N_STREAMS]; // current seeds + int stream_; // current RNG stream }; } // namespace openmc diff --git a/include/openmc/photon.h b/include/openmc/photon.h index c7bbeddd2..b652a4d8b 100644 --- a/include/openmc/photon.h +++ b/include/openmc/photon.h @@ -45,12 +45,12 @@ public: void calculate_xs(Particle& p) const; void compton_scatter(double alpha, bool doppler, double* alpha_out, - double* mu, int* i_shell, uint64_t* prn_seed) const; + double* mu, int* i_shell, uint64_t* seed) const; - double rayleigh_scatter(double alpha, uint64_t* prn_seed) const; + double rayleigh_scatter(double alpha, uint64_t* seed) const; void pair_production(double alpha, double* E_electron, double* E_positron, - double* mu_electron, double* mu_positron, uint64_t* prn_seed) const; + double* mu_electron, double* mu_positron, uint64_t* seed) const; void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; @@ -97,14 +97,14 @@ public: private: void compton_doppler(double alpha, double mu, double* E_out, int* i_shell, - uint64_t* prn_seed) const; + uint64_t* seed) const; }; //============================================================================== // Non-member functions //============================================================================== -std::pair klein_nishina(double alpha, uint64_t* prn_seed); +std::pair klein_nishina(double alpha, uint64_t* seed); void free_memory_photon(); diff --git a/include/openmc/physics.h b/include/openmc/physics.h index fce164cd6..fb43775a6 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -72,17 +72,17 @@ void sab_scatter(int i_nuclide, int i_sab, Particle* p); //! dependence of cross sections in treating resonance elastic scattering such //! as the DBRC and a new, accelerated scheme are also implemented here. Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT, uint64_t* prn_seed); + Direction v_neut, double xs_eff, double kT, uint64_t* seed); //! samples a target velocity based on the free gas scattering formulation, used //! by most Monte Carlo codes, in which cross section is assumed to be constant //! in energy. Excellent documentation for this method can be found in //! FRA-TM-123. Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT, - uint64_t* prn_seed); + uint64_t* seed); void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, - Particle::Bank* site, uint64_t* prn_seed); + Particle::Bank* site, uint64_t* seed); //! handles all reactions with a single secondary neutron (other than fission), //! i.e. level scattering, (n,np), (n,na), etc. diff --git a/include/openmc/random_lcg.h b/include/openmc/random_lcg.h index 5521f1df0..f5d5aa15b 100644 --- a/include/openmc/random_lcg.h +++ b/include/openmc/random_lcg.h @@ -21,11 +21,11 @@ constexpr int64_t DEFAULT_SEED {1}; //============================================================================== //! Generate a pseudo-random number using a linear congruential generator. -//! @param prn_seed Pseudorandom number seed pointer +//! @param seed Pseudorandom number seed pointer //! @return A random number between 0 and 1 //============================================================================== -double prn(uint64_t* prn_seed); +double prn(uint64_t* seed); //============================================================================== //! Generate a random number which is 'n' times ahead from the current seed. @@ -34,11 +34,11 @@ double prn(uint64_t* prn_seed); //! `prn()` 'n' times, though without the side effect of altering the RNG //! state. //! @param n The number of RNG seeds to skip ahead by -//! @param prn_seed Pseudorandom number seed +//! @param seed Pseudorandom number seed //! @return A random number between 0 and 1 //============================================================================== -double future_prn(int64_t n, uint64_t prn_seed); +double future_prn(int64_t n, uint64_t seed); //============================================================================== //! Set a RNG seed to a unique value based on a unique particle ID by striding @@ -55,21 +55,21 @@ uint64_t init_seed(int64_t id, int offset); //! Set the RNG seeds to unique values based on the ID of the particle. This //! function initializes the seeds for all RNG streams of the particle via //! striding. -//! @param prn_seeds Pseudorandom number seed array +//! @param seeds Pseudorandom number seed array //! @param id The particle ID //============================================================================== -void init_particle_seeds(int64_t id, uint64_t* prn_seeds); +void init_particle_seeds(int64_t id, uint64_t* seeds); //============================================================================== //! Advance the random number seed 'n' times from the current seed. This //! differs from the future_prn() function in that this function does alter //! the RNG state. -//! @param prn_seed Pseudorandom number seed pointer +//! @param seed Pseudorandom number seed pointer //! @param n The number of RNG seeds to skip ahead by //============================================================================== -void advance_prn_seed(int64_t n, uint64_t* prn_seed); +void advance_prn_seed(int64_t n, uint64_t* seed); //============================================================================== //! Advance a random number seed 'n' times. @@ -81,7 +81,7 @@ void advance_prn_seed(int64_t n, uint64_t* prn_seed); //! @param seed The starting to seed to advance from //============================================================================== -uint64_t future_seed(uint64_t n, uint64_t prn_seed); +uint64_t future_seed(uint64_t n, uint64_t seed); //============================================================================== // API FUNCTIONS diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index c61eec6b0..c842b6f96 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -42,8 +42,8 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom seed pointer - void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seed) const; + //! \param[inout] seed Pseudorandom seed pointer + void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const; Particle::Type particle_; //!< Particle type EmissionMode emission_mode_; //!< Emission mode diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index 8f3a79b77..f17e689bf 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -65,9 +65,9 @@ class ScattData { //! @param gout Sampled outgoing energy group. //! @param mu Sampled cosine of the change-in-angle. //! @param wgt Weight of the particle to be adjusted. - //! @param prn_seed Pseudorandom number seed pointer + //! @param seed Pseudorandom number seed pointer virtual void - sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed) = 0; + sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed) = 0; //! \brief Initializes the ScattData object from a given scatter and //! multiplicity matrix. @@ -110,9 +110,9 @@ class ScattData { //! @param gin Incoming energy group. //! @param gout Sampled outgoing energy group. //! @param i_gout Sampled outgoing energy group index. - //! @param prn_seed Pseudorandom number seed pointer + //! @param seed Pseudorandom number seed pointer void - sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seed); + sample_energy(int gin, int& gout, int& i_gout, uint64_t* seed); //! \brief Provides a cross section value given certain parameters //! @@ -163,7 +163,7 @@ class ScattDataLegendre: public ScattData { calc_f(int gin, int gout, double mu); void - sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed); + sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed); size_t get_order() {return dist[0][0].size() - 1;}; @@ -199,7 +199,7 @@ class ScattDataHistogram: public ScattData { calc_f(int gin, int gout, double mu); void - sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed); + sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed); size_t get_order() {return dist[0][0].size();}; @@ -240,7 +240,7 @@ class ScattDataTabular: public ScattData { calc_f(int gin, int gout, double mu); void - sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed); + sample(int gin, int& gout, double& mu, double& wgt, uint64_t* seed); size_t get_order() {return dist[0][0].size();}; diff --git a/include/openmc/secondary_correlated.h b/include/openmc/secondary_correlated.h index 6b7c36796..3fdfd8fad 100644 --- a/include/openmc/secondary_correlated.h +++ b/include/openmc/secondary_correlated.h @@ -38,9 +38,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; // energy property std::vector& energy() { return energy_; } diff --git a/include/openmc/secondary_kalbach.h b/include/openmc/secondary_kalbach.h index b9f190635..0864367c0 100644 --- a/include/openmc/secondary_kalbach.h +++ b/include/openmc/secondary_kalbach.h @@ -29,9 +29,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: //! Outgoing energy/angle at a single incoming energy struct KMTable { diff --git a/include/openmc/secondary_nbody.h b/include/openmc/secondary_nbody.h index b810d17d0..a0f6787ad 100644 --- a/include/openmc/secondary_nbody.h +++ b/include/openmc/secondary_nbody.h @@ -24,9 +24,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: int n_bodies_; //!< Number of particles distributed double mass_ratio_; //!< Total mass of particles [neutron mass] diff --git a/include/openmc/secondary_thermal.h b/include/openmc/secondary_thermal.h index 74bb3265f..6e1768f46 100644 --- a/include/openmc/secondary_thermal.h +++ b/include/openmc/secondary_thermal.h @@ -31,9 +31,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section }; @@ -53,9 +53,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: double debye_waller_; }; @@ -76,9 +76,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: const std::vector& energy_; //!< Energies at which cosines are tabulated xt::xtensor mu_out_; //!< Cosines for each incident energy @@ -100,9 +100,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: const std::vector& energy_; //!< Incident energies xt::xtensor energy_out_; //!< Outgoing energies for each incident energy @@ -125,9 +125,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom number seed pointer + //! \param[inout] seed Pseudorandom number seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; private: //! Secondary energy/angle distribution struct DistEnergySab { diff --git a/include/openmc/secondary_uncorrelated.h b/include/openmc/secondary_uncorrelated.h index 66ef834f1..79b8b5031 100644 --- a/include/openmc/secondary_uncorrelated.h +++ b/include/openmc/secondary_uncorrelated.h @@ -29,9 +29,9 @@ public: //! \param[in] E_in Incoming energy in [eV] //! \param[out] E_out Outgoing energy in [eV] //! \param[out] mu Outgoing cosine with respect to current direction - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const override; + uint64_t* seed) const override; // Accessors AngleDistribution& angle() { return angle_; } diff --git a/include/openmc/source.h b/include/openmc/source.h index 8a7c30a09..a177995ea 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -38,9 +38,9 @@ public: explicit SourceDistribution(pugi::xml_node node); //! Sample from the external source distribution - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer //! \return Sampled site - Particle::Bank sample(uint64_t* prn_seed) const; + Particle::Bank sample(uint64_t* seed) const; // Properties double strength() const { return strength_; } @@ -61,9 +61,9 @@ extern "C" void initialize_source(); //! Sample a site from all external source distributions in proportion to their //! source strength -//! \param[inout] prn_seed Pseudorandom seed pointer +//! \param[inout] seed Pseudorandom seed pointer //! \return Sampled source site -Particle::Bank sample_external_source(uint64_t* prn_seed); +Particle::Bank sample_external_source(uint64_t* seed); //! Fill source bank at end of generation for fixed source simulations void fill_source_bank_fixedsource(); diff --git a/include/openmc/surface.h b/include/openmc/surface.h index 35fa84bd1..2d5000d28 100644 --- a/include/openmc/surface.h +++ b/include/openmc/surface.h @@ -117,7 +117,7 @@ public: virtual Direction reflect(Position r, Direction u) const; virtual Direction diffuse_reflect(Position r, Direction u, - uint64_t* prn_seed) const; + uint64_t* seed) const; //! Evaluate the equation describing the surface. //! diff --git a/include/openmc/thermal.h b/include/openmc/thermal.h index 37ad8f4cb..16be335b6 100644 --- a/include/openmc/thermal.h +++ b/include/openmc/thermal.h @@ -64,9 +64,9 @@ public: //! \param[in] E_in Incident neutron energy in [eV] //! \param[out] E_out Outgoing neutron energy in [eV] //! \param[out] mu Outgoing scattering angle cosine - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void sample(const NuclideMicroXS& micro_xs, double E_in, - double* E_out, double* mu, uint64_t* prn_seed); + double* E_out, double* mu, uint64_t* seed); private: struct Reaction { // Default constructor @@ -101,9 +101,9 @@ public: //! \param[out] i_temp corresponding temperature index //! \param[out] elastic Thermal elastic scattering cross section //! \param[out] inelastic Thermal inelastic scattering cross section - //! \param[inout] prn_seed Pseudorandom seed pointer + //! \param[inout] seed Pseudorandom seed pointer void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic, - double* inelastic, uint64_t* prn_seed) const; + double* inelastic, uint64_t* seed) const; //! Determine whether table applies to a particular nuclide //! diff --git a/src/distribution.cpp b/src/distribution.cpp index e48ce8baa..92cabef0d 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -35,11 +35,11 @@ Discrete::Discrete(const double* x, const double* p, int n) normalize(); } -double Discrete::sample(uint64_t* prn_seed) const +double Discrete::sample(uint64_t* seed) const { int n = x_.size(); if (n > 1) { - double xi = prn(prn_seed); + double xi = prn(seed); double c = 0.0; for (int i = 0; i < n; ++i) { c += p_[i]; @@ -74,9 +74,9 @@ Uniform::Uniform(pugi::xml_node node) b_ = params.at(1); } -double Uniform::sample(uint64_t* prn_seed) const +double Uniform::sample(uint64_t* seed) const { - return a_ + prn(prn_seed)*(b_ - a_); + return a_ + prn(seed)*(b_ - a_); } //============================================================================== @@ -88,9 +88,9 @@ Maxwell::Maxwell(pugi::xml_node node) theta_ = std::stod(get_node_value(node, "parameters")); } -double Maxwell::sample(uint64_t* prn_seed) const +double Maxwell::sample(uint64_t* seed) const { - return maxwell_spectrum(theta_, prn_seed); + return maxwell_spectrum(theta_, seed); } //============================================================================== @@ -108,9 +108,9 @@ Watt::Watt(pugi::xml_node node) b_ = params.at(1); } -double Watt::sample(uint64_t* prn_seed) const +double Watt::sample(uint64_t* seed) const { - return watt_spectrum(a_, b_, prn_seed); + return watt_spectrum(a_, b_, seed); } //============================================================================== @@ -127,9 +127,9 @@ Normal::Normal(pugi::xml_node node) std_dev_ = params.at(1); } -double Normal::sample(uint64_t* prn_seed) const +double Normal::sample(uint64_t* seed) const { - return normal_variate(mean_value_, std_dev_, prn_seed); + return normal_variate(mean_value_, std_dev_, seed); } //============================================================================== @@ -147,9 +147,9 @@ Muir::Muir(pugi::xml_node node) kt_ = params.at(2); } -double Muir::sample(uint64_t* prn_seed) const +double Muir::sample(uint64_t* seed) const { - return muir_spectrum(e0_, m_rat_, kt_, prn_seed); + return muir_spectrum(e0_, m_rat_, kt_, seed); } //============================================================================== @@ -220,10 +220,10 @@ void Tabular::init(const double* x, const double* p, std::size_t n, const double } } -double Tabular::sample(uint64_t* prn_seed) const +double Tabular::sample(uint64_t* seed) const { // Sample value of CDF - double c = prn(prn_seed); + double c = prn(seed); // Find first CDF bin which is above the sampled value double c_i = c_[0]; @@ -263,11 +263,11 @@ double Tabular::sample(uint64_t* prn_seed) const // Equiprobable implementation //============================================================================== -double Equiprobable::sample(uint64_t* prn_seed) const +double Equiprobable::sample(uint64_t* seed) const { std::size_t n = x_.size(); - double r = prn(prn_seed); + double r = prn(seed); int i = std::floor((n - 1)*r); double xl = x_[i]; diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index 076f82470..970ea9040 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -62,7 +62,7 @@ AngleDistribution::AngleDistribution(hid_t group) } } -double AngleDistribution::sample(double E, uint64_t* prn_seed) const +double AngleDistribution::sample(double E, uint64_t* seed) const { // Determine number of incoming energies auto n = energy_.size(); @@ -83,10 +83,10 @@ double AngleDistribution::sample(double E, uint64_t* prn_seed) const } // Sample between the ith and (i+1)th bin - if (r > prn(prn_seed)) ++i; + if (r > prn(seed)) ++i; // Sample i-th distribution - double mu = distribution_[i]->sample(prn_seed); + double mu = distribution_[i]->sample(seed); // Make sure mu is in range [-1,1] and return if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu); diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 3e1995145..82754faf8 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -25,7 +25,7 @@ DiscretePhoton::DiscretePhoton(hid_t group) read_attribute(group, "atomic_weight_ratio", A_); } -double DiscretePhoton::sample(double E, uint64_t* prn_seed) const +double DiscretePhoton::sample(double E, uint64_t* seed) const { if (primary_flag_ == 2) { return energy_ + A_/(A_+ 1)*E; @@ -44,7 +44,7 @@ LevelInelastic::LevelInelastic(hid_t group) read_attribute(group, "mass_ratio", mass_ratio_); } -double LevelInelastic::sample(double E, uint64_t* prn_seed ) const +double LevelInelastic::sample(double E, uint64_t* seed) const { return mass_ratio_*(E - threshold_); } @@ -146,7 +146,7 @@ ContinuousTabular::ContinuousTabular(hid_t group) } // incoming energies } -double ContinuousTabular::sample(double E, uint64_t* prn_seed) const +double ContinuousTabular::sample(double E, uint64_t* seed) const { // Read number of interpolation regions and incoming energies bool histogram_interp; @@ -177,7 +177,7 @@ double ContinuousTabular::sample(double E, uint64_t* prn_seed) const if (histogram_interp) { l = i; } else { - l = r > prn(prn_seed) ? i + 1 : i; + l = r > prn(seed) ? i + 1 : i; } // Interpolation for energy E1 and EK @@ -197,7 +197,7 @@ double ContinuousTabular::sample(double E, uint64_t* prn_seed) const // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); n_discrete = distribution_[l].n_discrete; - double r1 = prn(prn_seed); + double r1 = prn(seed); double c_k = distribution_[l].c[0]; int k = 0; int end = n_energy_out - 2; @@ -275,14 +275,14 @@ MaxwellEnergy::MaxwellEnergy(hid_t group) close_dataset(dset); } -double MaxwellEnergy::sample(double E, uint64_t* prn_seed) const +double MaxwellEnergy::sample(double E, uint64_t* seed) const { // Get temperature corresponding to incoming energy double theta = theta_(E); while (true) { // Sample maxwell fission spectrum - double E_out = maxwell_spectrum(theta, prn_seed); + double E_out = maxwell_spectrum(theta, seed); // Accept energy based on restriction energy if (E_out <= E - u_) return E_out; @@ -301,7 +301,7 @@ Evaporation::Evaporation(hid_t group) close_dataset(dset); } -double Evaporation::sample(double E, uint64_t* prn_seed) const +double Evaporation::sample(double E, uint64_t* seed) const { // Get temperature corresponding to incoming energy double theta = theta_(E); @@ -313,7 +313,7 @@ double Evaporation::sample(double E, uint64_t* prn_seed) const // density function double x; while (true) { - x = -std::log((1.0 - v*prn(prn_seed))*(1.0 - v*prn(prn_seed))); + x = -std::log((1.0 - v*prn(seed))*(1.0 - v*prn(seed))); if (x <= y) break; } @@ -338,7 +338,7 @@ WattEnergy::WattEnergy(hid_t group) close_dataset(dset); } -double WattEnergy::sample(double E, uint64_t* prn_seed) const +double WattEnergy::sample(double E, uint64_t* seed) const { // Determine Watt parameters at incident energy double a = a_(E); @@ -346,7 +346,7 @@ double WattEnergy::sample(double E, uint64_t* prn_seed) const while (true) { // Sample energy-dependent Watt fission spectrum - double E_out = watt_spectrum(a, b, prn_seed); + double E_out = watt_spectrum(a, b, seed); // Accept energy based on restriction energy if (E_out <= E - u_) return E_out; diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp index 711234201..ea24e3ea3 100644 --- a/src/distribution_multi.cpp +++ b/src/distribution_multi.cpp @@ -53,31 +53,31 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node) } } -Direction PolarAzimuthal::sample(uint64_t* prn_seed) const +Direction PolarAzimuthal::sample(uint64_t* seed) const { // Sample cosine of polar angle - double mu = mu_->sample(prn_seed); + double mu = mu_->sample(seed); if (mu == 1.0) return u_ref_; // Sample azimuthal angle - double phi = phi_->sample(prn_seed); + double phi = phi_->sample(seed); // If the reference direction is along the z-axis, rotate the aziumthal angle // to match spherical coordinate conventions. // TODO: apply this change directly to rotate_angle if (u_ref_.x == 0 && u_ref_.y == 0) phi += 0.5*PI; - return rotate_angle(u_ref_, mu, &phi, prn_seed); + return rotate_angle(u_ref_, mu, &phi, seed); } //============================================================================== // Isotropic implementation //============================================================================== -Direction Isotropic::sample(uint64_t* prn_seed) const +Direction Isotropic::sample(uint64_t* seed) const { - double phi = 2.0*PI*prn(prn_seed); - double mu = 2.0*prn(prn_seed) - 1.0; + double phi = 2.0*PI*prn(seed); + double mu = 2.0*prn(seed) - 1.0; return {mu, std::sqrt(1.0 - mu*mu) * std::cos(phi), std::sqrt(1.0 - mu*mu) * std::sin(phi)}; } @@ -86,7 +86,7 @@ Direction Isotropic::sample(uint64_t* prn_seed) const // Monodirectional implementation //============================================================================== -Direction Monodirectional::sample(uint64_t* prn_seed) const +Direction Monodirectional::sample(uint64_t* seed) const { return u_ref_; } diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index ce2b7f34a..a5c45e4e9 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -46,9 +46,9 @@ CartesianIndependent::CartesianIndependent(pugi::xml_node node) } } -Position CartesianIndependent::sample(uint64_t* prn_seed) const +Position CartesianIndependent::sample(uint64_t* seed) const { - return {x_->sample(prn_seed), y_->sample(prn_seed), z_->sample(prn_seed)}; + return {x_->sample(seed), y_->sample(seed), z_->sample(seed)}; } //============================================================================== @@ -107,11 +107,11 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node) } -Position SphericalIndependent::sample(uint64_t* prn_seed) const +Position SphericalIndependent::sample(uint64_t* seed) const { - double r = r_->sample(prn_seed); - double theta = theta_->sample(prn_seed); - double phi = phi_->sample(prn_seed); + double r = r_->sample(seed); + double theta = theta_->sample(seed); + double phi = phi_->sample(seed); double x = r*sin(theta)*cos(phi) + origin_.x; double y = r*sin(theta)*sin(phi) + origin_.y; double z = r*cos(theta) + origin_.z; @@ -135,9 +135,9 @@ SpatialBox::SpatialBox(pugi::xml_node node, bool fission) upper_right_ = Position{params[3], params[4], params[5]}; } -Position SpatialBox::sample(uint64_t* prn_seed) const +Position SpatialBox::sample(uint64_t* seed) const { - Position xi {prn(prn_seed), prn(prn_seed), prn(prn_seed)}; + Position xi {prn(seed), prn(seed), prn(seed)}; return lower_left_ + xi*(upper_right_ - lower_left_); } @@ -157,7 +157,7 @@ SpatialPoint::SpatialPoint(pugi::xml_node node) r_ = Position{params.data()}; } -Position SpatialPoint::sample(uint64_t* prn_seed) const +Position SpatialPoint::sample(uint64_t* seed) const { return r_; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 21afd696c..860ff2e8c 100644 --- a/src/math_functions.cpp +++ b/src/math_functions.cpp @@ -630,22 +630,22 @@ void calc_zn_rad(int n, double rho, double zn_rad[]) { } -void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t* prn_seed) { - Direction u = rotate_angle({uvw}, mu, phi, prn_seed); +void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t* seed) { + Direction u = rotate_angle({uvw}, mu, phi, seed); uvw[0] = u.x; uvw[1] = u.y; uvw[2] = u.z; } -Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* prn_seed) +Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* seed) { // Sample azimuthal angle in [0,2pi) if none provided double phi_; if (phi != nullptr) { phi_ = (*phi); } else { - phi_ = 2.0*PI*prn(prn_seed); + phi_ = 2.0*PI*prn(seed); } // Precompute factors to save flops @@ -675,11 +675,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* prn_ } -double maxwell_spectrum(double T, uint64_t* prn_seed) { +double maxwell_spectrum(double T, uint64_t* seed) { // Set the random numbers - double r1 = prn(prn_seed); - double r2 = prn(prn_seed); - double r3 = prn(prn_seed); + double r1 = prn(seed); + double r2 = prn(seed); + double r3 = prn(seed); // determine cosine of pi/2*r double c = std::cos(PI / 2. * r3); @@ -691,33 +691,33 @@ double maxwell_spectrum(double T, uint64_t* prn_seed) { } -double normal_variate(double mean, double standard_deviation, uint64_t* prn_seed) { +double normal_variate(double mean, double standard_deviation, uint64_t* seed) { // perhaps there should be a limit to the number of resamples while ( true ) { - double v1 = 2 * prn(prn_seed) - 1.; - double v2 = 2 * prn(prn_seed) - 1.; + double v1 = 2 * prn(seed) - 1.; + double v2 = 2 * prn(seed) - 1.; double r = std::pow(v1, 2) + std::pow(v2, 2); double r2 = std::pow(r, 2); if (r2 < 1) { double z = std::sqrt(-2.0 * std::log(r2)/r2); - z *= (prn(prn_seed) <= 0.5) ? v1 : v2; + z *= (prn(seed) <= 0.5) ? v1 : v2; return mean + standard_deviation*z; } } } -double muir_spectrum(double e0, double m_rat, double kt, uint64_t* prn_seed) { +double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) { // note sigma here is a factor of 2 shy of equation // 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS double sigma = std::sqrt(2.*e0*kt/m_rat); - return normal_variate(e0, sigma, prn_seed); + return normal_variate(e0, sigma, seed); } -double watt_spectrum(double a, double b, uint64_t* prn_seed) { - double w = maxwell_spectrum(a, prn_seed); - double E_out = w + 0.25 * a * a * b + (2. * prn(prn_seed) - 1.) * std::sqrt(a * a * b * w); +double watt_spectrum(double a, double b, uint64_t* seed) { + double w = maxwell_spectrum(a, seed); + double E_out = w + 0.25 * a * a * b + (2. * prn(seed) - 1.) * std::sqrt(a * a * b * w); return E_out; } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index e2195d512..bba58f668 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -529,7 +529,7 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu, //============================================================================== void -Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed) +Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* seed) { // This method assumes that the temperature and angle indices are set #ifdef _OPENMP @@ -544,8 +544,8 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed) double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin); // sample random numbers - double xi_pd = prn(prn_seed) * nu_fission; - double xi_gout = prn(prn_seed); + double xi_pd = prn(seed) * nu_fission; + double xi_gout = prn(seed); // Select whether the neutron is prompt or delayed if (xi_pd <= prob_prompt) { @@ -585,7 +585,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed) //============================================================================== void -Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed) +Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* seed) { // This method assumes that the temperature and angle indices are set // Sample the data @@ -594,7 +594,7 @@ Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_ #else int tid = 0; #endif - xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, prn_seed); + xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, seed); } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index 6fb41f295..f7bce7367 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -678,7 +678,7 @@ Particle::write_restart() const } // #pragma omp critical } -uint64_t* Particle::current_seed() {return prn_seeds_ + stream_;} -const uint64_t* Particle::current_seed() const {return prn_seeds_ + stream_;} +uint64_t* Particle::current_seed() {return seeds_ + stream_;} +const uint64_t* Particle::current_seed() const {return seeds_ + stream_;} } // namespace openmc diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index b7d94dd80..053193207 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -98,7 +98,7 @@ void run_particle_restart() throw std::runtime_error{"Unexpected run mode: " + std::to_string(previous_run_mode)}; } - init_particle_seeds(particle_seed, p.prn_seeds_); + init_particle_seeds(particle_seed, p.seeds_); // Transport neutron p.transport(); diff --git a/src/photon.cpp b/src/photon.cpp index 066e182fc..ed096c396 100644 --- a/src/photon.cpp +++ b/src/photon.cpp @@ -290,12 +290,12 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element) } void PhotonInteraction::compton_scatter(double alpha, bool doppler, - double* alpha_out, double* mu, int* i_shell, uint64_t* prn_seed) const + double* alpha_out, double* mu, int* i_shell, uint64_t* seed) const { double form_factor_xmax = 0.0; while (true) { // Sample Klein-Nishina distribution for trial energy and angle - std::tie(*alpha_out, *mu) = klein_nishina(alpha, prn_seed); + std::tie(*alpha_out, *mu) = klein_nishina(alpha, seed); // Note that the parameter used here does not correspond exactly to the // momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the @@ -309,10 +309,10 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler, } // Perform rejection on form factor - if (prn(prn_seed) < form_factor_x / form_factor_xmax) { + if (prn(seed) < form_factor_x / form_factor_xmax) { if (doppler) { double E_out; - this->compton_doppler(alpha, *mu, &E_out, i_shell, prn_seed); + this->compton_doppler(alpha, *mu, &E_out, i_shell, seed); *alpha_out = E_out/MASS_ELECTRON_EV; } else { *i_shell = -1; @@ -323,14 +323,14 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler, } void PhotonInteraction::compton_doppler(double alpha, double mu, - double* E_out, int* i_shell, uint64_t* prn_seed) const + double* E_out, int* i_shell, uint64_t* seed) const { auto n = data::compton_profile_pz.size(); int shell; // index for shell while (true) { // Sample electron shell - double rn = prn(prn_seed); + double rn = prn(seed); double c = 0.0; for (shell = 0; shell < electron_pdf_.size(); ++shell) { c += electron_pdf_(shell); @@ -377,7 +377,7 @@ void PhotonInteraction::compton_doppler(double alpha, double mu, } // Sample value on bounded cdf - c = prn(prn_seed)*c_max; + c = prn(seed)*c_max; // Determine pz corresponding to sampled cdf value auto cdf_shell = xt::view(profile_cdf_, shell, xt::all()); @@ -418,7 +418,7 @@ void PhotonInteraction::compton_doppler(double alpha, double mu, if (E_out1 > 0.0) { if (E_out2 > 0.0) { // If both are positive, pick one at random - *E_out = prn(prn_seed) < 0.5 ? E_out1 : E_out2; + *E_out = prn(seed) < 0.5 ? E_out1 : E_out2; } else { *E_out = E_out1; } @@ -496,7 +496,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const xs.last_E = p.E_; } -double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seed) const +double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* seed) const { double mu; while (true) { @@ -507,7 +507,7 @@ double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seed) con double F_max = coherent_int_form_factor_(x2_max); // Sample cumulative distribution - double F = prn(prn_seed)*F_max; + double F = prn(seed)*F_max; // Determine x^2 corresponding to F const auto& x {coherent_int_form_factor_.x()}; @@ -519,14 +519,14 @@ double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seed) con // Calculate mu mu = 1.0 - 2.0*x2/x2_max; - if (prn(prn_seed) < 0.5*(1.0 + mu*mu)) break; + if (prn(seed) < 0.5*(1.0 + mu*mu)) break; } return mu; } void PhotonInteraction::pair_production(double alpha, double* E_electron, double* E_positron, double* mu_electron, double* mu_positron, - uint64_t* prn_seed) const + uint64_t* seed) const { constexpr double r[] { 122.81, 73.167, 69.228, 67.301, 64.696, 61.228, @@ -593,12 +593,12 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron, double u2 = phi2_max; double e; while (true) { - double rn = prn(prn_seed); + double rn = prn(seed); // Sample the index i in (1, 2) using the point probabilities // p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2) int i; - if (prn(prn_seed) < u1/(u1 + u2)) { + if (prn(seed) < u1/(u1 + u2)) { i = 1; // Sample e from pi_1 using the inverse transform method @@ -619,10 +619,10 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron, t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b))); if (i == 1) { double phi1 = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4; - if (prn(prn_seed) <= phi1/phi1_max) break; + if (prn(seed) <= phi1/phi1_max) break; } else { double phi2 = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4; - if (prn(prn_seed) <= phi2/phi2_max) break; + if (prn(seed) <= phi2/phi2_max) break; } } @@ -635,13 +635,13 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron, // p(mu) = C/(1 - beta*mu)^2 using the inverse transform method. double beta = std::sqrt(*E_electron*(*E_electron + 2.0*MASS_ELECTRON_EV)) / (*E_electron + MASS_ELECTRON_EV) ; - double rn = 2.0*prn(prn_seed) - 1.0; + double rn = 2.0*prn(seed) - 1.0; *mu_electron = (rn + beta)/(rn*beta + 1.0); // Sample the scattering angle of the positron beta = std::sqrt(*E_positron*(*E_positron + 2.0*MASS_ELECTRON_EV)) / (*E_positron + MASS_ELECTRON_EV); - rn = 2.0*prn(prn_seed) - 1.0; + rn = 2.0*prn(seed) - 1.0; *mu_positron = (rn + beta)/(rn*beta + 1.0); } @@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-member functions //============================================================================== -std::pair klein_nishina(double alpha, uint64_t* prn_seed) +std::pair klein_nishina(double alpha, uint64_t* seed) { double alpha_out, mu; double beta = 1.0 + 2.0*alpha; @@ -720,19 +720,19 @@ std::pair klein_nishina(double alpha, uint64_t* prn_seed) double t = beta/(beta + 8.0); double x; while (true) { - if (prn(prn_seed) < t) { + if (prn(seed) < t) { // Left branch of flow chart - double r = 2.0*prn(prn_seed); + double r = 2.0*prn(seed); x = 1.0 + alpha*r; - if (prn(prn_seed) < 4.0/x*(1.0 - 1.0/x)) { + if (prn(seed) < 4.0/x*(1.0 - 1.0/x)) { mu = 1 - r; break; } } else { // Right branch of flow chart - x = beta/(1.0 + 2.0*alpha*prn(prn_seed)); + x = beta/(1.0 + 2.0*alpha*prn(seed)); mu = 1.0 + (1.0 - x)/alpha; - if (prn(prn_seed) < 0.5*(mu*mu + 1.0/x)) break; + if (prn(seed) < 0.5*(mu*mu + 1.0/x)) break; } } alpha_out = alpha/x; @@ -740,24 +740,24 @@ std::pair klein_nishina(double alpha, uint64_t* prn_seed) } else { // Koblinger's direct method double gamma = 1.0 - std::pow(beta, -2); - double s = prn(prn_seed)*(4.0/alpha + 0.5*gamma + + double s = prn(seed)*(4.0/alpha + 0.5*gamma + (1.0 - (1.0 + beta)/(alpha*alpha))*std::log(beta)); if (s <= 2.0/alpha) { // For first term, x = 1 + 2ar // Therefore, a' = a/(1 + 2ar) - alpha_out = alpha/(1.0 + 2.0*alpha*prn(prn_seed)); + alpha_out = alpha/(1.0 + 2.0*alpha*prn(seed)); } else if (s <= 4.0/alpha) { // For third term, x = beta/(1 + 2ar) // Therefore, a' = a(1 + 2ar)/beta - alpha_out = alpha*(1.0 + 2.0*alpha*prn(prn_seed))/beta; + alpha_out = alpha*(1.0 + 2.0*alpha*prn(seed))/beta; } else if (s <= 4.0/alpha + 0.5*gamma) { // For fourth term, x = 1/sqrt(1 - gamma*r) // Therefore, a' = a*sqrt(1 - gamma*r) - alpha_out = alpha*std::sqrt(1.0 - gamma*prn(prn_seed)); + alpha_out = alpha*std::sqrt(1.0 - gamma*prn(seed)); } else { // For third term, x = beta^r // Therefore, a' = a/beta^r - alpha_out = alpha/std::pow(beta, prn(prn_seed)); + alpha_out = alpha/std::pow(beta, prn(seed)); } // Calculate cosine of scattering angle based on basic relation diff --git a/src/physics.cpp b/src/physics.cpp index 7b86259fe..4b80d1039 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -750,7 +750,7 @@ void sab_scatter(int i_nuclide, int i_sab, Particle* p) } Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT, uint64_t* prn_seed) + Direction v_neut, double xs_eff, double kT, uint64_t* seed) { // check if nuclide is a resonant scatterer ResScatMethod sampling_method; @@ -782,7 +782,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, case ResScatMethod::cxs: // sample target velocity with the constant cross section (cxs) approx. - return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed); + return sample_cxs_target_velocity(nuc->awr_, E, u, kT, seed); case ResScatMethod::dbrc: case ResScatMethod::rvs: { @@ -816,7 +816,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, if (i_E_up == i_E_low) { // Handle degenerate case -- if the upper/lower bounds occur for the same // index, then using cxs is probably a good approximation - return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed); + return sample_cxs_target_velocity(nuc->awr_, E, u, kT, seed); } if (sampling_method == ResScatMethod::dbrc) { @@ -840,7 +840,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, Direction v_target; while (true) { // sample target velocity with the constant cross section (cxs) approx. - v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed); + v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, seed); Direction v_rel = v_neut - v_target; E_rel = v_rel.dot(v_rel); if (E_rel < E_up) break; @@ -849,7 +849,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, // perform Doppler broadening rejection correction (dbrc) double xs_0K = nuc->elastic_xs_0K(E_rel); double R = xs_0K / xs_max; - if (prn(prn_seed) < R) return v_target; + if (prn(seed) < R) return v_target; } } else if (sampling_method == ResScatMethod::rvs) { @@ -869,10 +869,10 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, while (true) { // directly sample Maxwellian - double E_t = -kT * std::log(prn(prn_seed)); + double E_t = -kT * std::log(prn(seed)); // sample a relative energy using the xs cdf - double cdf_rel = cdf_low + prn(prn_seed)*(cdf_up - cdf_low); + double cdf_rel = cdf_low + prn(seed)*(cdf_up - cdf_low); int i_E_rel = lower_bound_index(&nuc->xs_cdf_[i_E_low-1], &nuc->xs_cdf_[i_E_up+1], cdf_rel); double E_rel = nuc->energy_0K_[i_E_low + i_E_rel]; @@ -890,7 +890,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, if (std::abs(mu) < 1.0) { // set and accept target velocity E_t /= nuc->awr_; - return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, prn_seed); + return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, seed); } } } @@ -901,7 +901,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, } Direction -sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t* prn_seed) +sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t* seed) { double beta_vn = std::sqrt(awr * E / kT); double alpha = 1.0/(1.0 + std::sqrt(PI)*beta_vn/2.0); @@ -910,10 +910,10 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_ double mu; while (true) { // Sample two random numbers - double r1 = prn(prn_seed); - double r2 = prn(prn_seed); + double r1 = prn(seed); + double r2 = prn(seed); - if (prn(prn_seed) < alpha) { + if (prn(seed) < alpha) { // With probability alpha, we sample the distribution p(y) = // y*e^(-y). This can be done with sampling scheme C45 frmo the Monte // Carlo sampler @@ -925,7 +925,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_ // e^(-y^2). This can be done with sampling scheme C61 from the Monte // Carlo sampler - double c = std::cos(PI/2.0 * prn(prn_seed)); + double c = std::cos(PI/2.0 * prn(seed)); beta_vt_sq = -std::log(r1) - std::log(r2)*c*c; } @@ -933,14 +933,14 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_ double beta_vt = std::sqrt(beta_vt_sq); // Sample cosine of angle between neutron and target velocity - mu = 2.0*prn(prn_seed) - 1.0; + mu = 2.0*prn(seed) - 1.0; // Determine rejection probability double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) / (beta_vn + beta_vt); // Perform rejection sampling on vt and mu - if (prn(prn_seed) < accept_prob) break; + if (prn(seed) < accept_prob) break; } // Determine speed of target nucleus @@ -948,19 +948,19 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_ // Determine velocity vector of target nucleus based on neutron's velocity // and the sampled angle between them - return vt * rotate_angle(u, mu, nullptr, prn_seed); + return vt * rotate_angle(u, mu, nullptr, seed); } -void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t* prn_seed) +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t* seed) { // Sample cosine of angle -- fission neutrons are always emitted // isotropically. Sometimes in ACE data, fission reactions actually have // an angular distribution listed, but for those that do, it's simply just // a uniform distribution in mu - double mu = 2.0 * prn(prn_seed) - 1.0; + double mu = 2.0 * prn(seed) - 1.0; // Sample azimuthal angle uniformly in [0,2*pi) - double phi = 2.0*PI*prn(prn_seed); + double phi = 2.0*PI*prn(seed); site->u.x = mu; site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi); site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi); @@ -971,12 +971,12 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed); double beta = nu_d / nu_t; - if (prn(prn_seed) < beta) { + if (prn(seed) < beta) { // ==================================================================== // DELAYED NEUTRON SAMPLED // sampled delayed precursor group - double xi = prn(prn_seed)*nu_d; + double xi = prn(seed)*nu_d; double prob = 0.0; int group; for (group = 1; group < nuc->n_precursor_; ++group) { @@ -1000,7 +1000,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part while (true) { // sample from energy/angle distribution -- note that mu has already been // sampled above and doesn't need to be resampled - rx->products_[group].sample(E_in, site->E, mu, prn_seed); + rx->products_[group].sample(E_in, site->E, mu, seed); // resample if energy is greater than maximum neutron energy constexpr int neutron = static_cast(Particle::Type::neutron); @@ -1025,7 +1025,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part // sample from prompt neutron energy distribution int n_sample = 0; while (true) { - rx->products_[0].sample(E_in, site->E, mu, prn_seed); + rx->products_[0].sample(E_in, site->E, mu, seed); // resample if energy is greater than maximum neutron energy constexpr int neutron = static_cast(Particle::Type::neutron); diff --git a/src/plot.cpp b/src/plot.cpp index dbc1ccb9e..31b771a47 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -80,7 +80,7 @@ namespace model { std::vector plots; std::unordered_map plot_map; -uint64_t plotter_prn_seed = 1; +uint64_t plotter_seed = 1; } // namespace model @@ -960,9 +960,9 @@ voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace) } RGBColor random_color(void) { - return {int(prn(&model::plotter_prn_seed)*255), - int(prn(&model::plotter_prn_seed)*255), - int(prn(&model::plotter_prn_seed)*255)}; + return {int(prn(&model::plotter_seed)*255), + int(prn(&model::plotter_seed)*255), + int(prn(&model::plotter_seed)*255)}; } extern "C" int openmc_id_map(const void* plot, int32_t* data_out) diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp index dbae0f8c8..8af093fb5 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -22,24 +22,24 @@ constexpr double prn_norm {1.0 / prn_mod}; // 2^-63 // PRN //============================================================================== -double prn(uint64_t* prn_seed) +double prn(uint64_t* seed) { // This algorithm uses bit-masking to find the next integer(8) value to be // used to calculate the random number. - *prn_seed = (prn_mult * (*prn_seed) + prn_add) & prn_mask; + *seed = (prn_mult * (*seed) + prn_add) & prn_mask; // Once the integer is calculated, we just need to divide by 2**m, // represented here as multiplying by a pre-calculated factor - return (*prn_seed) * prn_norm; + return (*seed) * prn_norm; } //============================================================================== // FUTURE_PRN //============================================================================== -double future_prn(int64_t n, uint64_t prn_seed) +double future_prn(int64_t n, uint64_t seed) { - return future_seed(static_cast(n), prn_seed) * prn_norm; + return future_seed(static_cast(n), seed) * prn_norm; } //============================================================================== @@ -66,16 +66,16 @@ void init_particle_seeds(int64_t id, uint64_t* prn_seeds) // ADVANCE_PRN_SEED //============================================================================== -void advance_prn_seed(int64_t n, uint64_t* prn_seed) +void advance_prn_seed(int64_t n, uint64_t* seed) { - *prn_seed = future_seed(static_cast(n), *prn_seed); + *seed = future_seed(static_cast(n), *seed); } //============================================================================== // FUTURE_SEED //============================================================================== -uint64_t future_seed(uint64_t n, uint64_t prn_seed) +uint64_t future_seed(uint64_t n, uint64_t seed) { // Make sure nskip is less than 2^M. n &= prn_mask; @@ -106,7 +106,7 @@ uint64_t future_seed(uint64_t n, uint64_t prn_seed) } // With G and C, we can now find the new seed. - return (g_new * prn_seed + c_new) & prn_mask; + return (g_new * seed + c_new) & prn_mask; } //============================================================================== diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index d92d14517..50870a3d1 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -77,25 +77,25 @@ ReactionProduct::ReactionProduct(hid_t group) } void ReactionProduct::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { auto n = applicability_.size(); if (n > 1) { double prob = 0.0; - double c = prn(prn_seed); + double c = prn(seed); for (int i = 0; i < n; ++i) { // Determine probability that i-th energy distribution is sampled prob += applicability_[i](E_in); // If i-th distribution is sampled, sample energy from the distribution if (c <= prob) { - distribution_[i]->sample(E_in, E_out, mu, prn_seed); + distribution_[i]->sample(E_in, E_out, mu, seed); break; } } } else { // If only one distribution is present, go ahead and sample it - distribution_[0]->sample(E_in, E_out, mu, prn_seed); + distribution_[0]->sample(E_in, E_out, mu, seed); } } diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 5d06375e5..c61d0249c 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -167,10 +167,10 @@ ScattData::base_combine(size_t max_order, //============================================================================== void -ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seed) +ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t* seed) { // Sample the outgoing group - double xi = prn(prn_seed); + double xi = prn(seed); double prob = 0.; i_gout = 0; for (gout = gmin[gin]; gout < gmax[gin]; ++gout) { @@ -348,21 +348,21 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu) void ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt, - uint64_t* prn_seed) + uint64_t* seed) { // Sample the outgoing energy using the base-class method int i_gout; - sample_energy(gin, gout, i_gout, prn_seed); + sample_energy(gin, gout, i_gout, seed); // Now we can sample mu using the scattering kernel using rejection // sampling from a rectangular bounding box double M = max_val[gin][i_gout]; int samples; for (samples = 0; samples < MAX_SAMPLE; ++samples) { - mu = 2. * prn(prn_seed) - 1.; + mu = 2. * prn(seed) - 1.; double f = calc_f(gin, gout, mu); if (f > 0.) { - double u = prn(prn_seed) * M; + double u = prn(seed) * M; if (u <= f) break; } } @@ -537,14 +537,14 @@ ScattDataHistogram::calc_f(int gin, int gout, double mu) void ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt, - uint64_t* prn_seed) + uint64_t* seed) { // Sample the outgoing energy using the base-class method int i_gout; - sample_energy(gin, gout, i_gout, prn_seed); + sample_energy(gin, gout, i_gout, seed); // Determine the outgoing cosine bin - double xi = prn(prn_seed); + double xi = prn(seed); int imu; if (xi < dist[gin][i_gout][0]) { @@ -556,7 +556,7 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt, } // Randomly select mu within the imu bin - mu = prn(prn_seed) * dmu + this->mu[imu]; + mu = prn(seed) * dmu + this->mu[imu]; if (mu < -1.) { mu = -1.; @@ -741,15 +741,15 @@ ScattDataTabular::calc_f(int gin, int gout, double mu) void ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt, - uint64_t* prn_seed) + uint64_t* seed) { // Sample the outgoing energy using the base-class method int i_gout; - sample_energy(gin, gout, i_gout, prn_seed); + sample_energy(gin, gout, i_gout, seed); // Determine the outgoing cosine bin int NP = this->mu.shape()[0]; - double xi = prn(prn_seed); + double xi = prn(seed); double c_k = dist[gin][i_gout][0]; int k; diff --git a/src/secondary_correlated.cpp b/src/secondary_correlated.cpp index 81e1ec1f9..8e0f23d52 100644 --- a/src/secondary_correlated.cpp +++ b/src/secondary_correlated.cpp @@ -153,14 +153,14 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group) } void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< // Before the secondary distribution refactor, an isotropic polar cosine was // always sampled but then overwritten with the polar cosine sampled from the // correlated distribution. To preserve the random number stream, we keep // this dummy sampling here but can remove it later (will change answers) - mu = 2.0*prn(prn_seed) - 1.0; + mu = 2.0*prn(seed) - 1.0; // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< // Find energy bin and calculate interpolation factor -- if the energy is @@ -180,7 +180,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, } // Sample between the ith and [i+1]th bin - int l = r > prn(prn_seed) ? i + 1 : i; + int l = r > prn(seed) ? i + 1 : i; // Interpolation for energy E1 and EK int n_energy_out = distribution_[i].e_out.size(); @@ -199,7 +199,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); n_discrete = distribution_[l].n_discrete; - double r1 = prn(prn_seed); + double r1 = prn(seed); double c_k = distribution_[l].c[0]; int k = 0; int end = n_energy_out - 2; @@ -260,9 +260,9 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, // Find correlated angular distribution for closest outgoing energy bin if (r1 - c_k < c_k1 - r1) { - mu = distribution_[l].angle[k]->sample(prn_seed); + mu = distribution_[l].angle[k]->sample(seed); } else { - mu = distribution_[l].angle[k + 1]->sample(prn_seed); + mu = distribution_[l].angle[k + 1]->sample(seed); } } diff --git a/src/secondary_kalbach.cpp b/src/secondary_kalbach.cpp index b0f9854ac..6a31a4376 100644 --- a/src/secondary_kalbach.cpp +++ b/src/secondary_kalbach.cpp @@ -113,14 +113,14 @@ KalbachMann::KalbachMann(hid_t group) } // incoming energies } -void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* prn_seed) const +void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* seed) const { // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< // Before the secondary distribution refactor, an isotropic polar cosine was // always sampled but then overwritten with the polar cosine sampled from the // correlated distribution. To preserve the random number stream, we keep // this dummy sampling here but can remove it later (will change answers) - mu = 2.0*prn(prn_seed) - 1.0; + mu = 2.0*prn(seed) - 1.0; // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< // Find energy bin and calculate interpolation factor -- if the energy is @@ -140,7 +140,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* prn_s } // Sample between the ith and [i+1]th bin - int l = r > prn(prn_seed) ? i + 1 : i; + int l = r > prn(seed) ? i + 1 : i; // Interpolation for energy E1 and EK int n_energy_out = distribution_[i].e_out.size(); @@ -159,7 +159,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* prn_s // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); n_discrete = distribution_[l].n_discrete; - double r1 = prn(prn_seed); + double r1 = prn(seed); double c_k = distribution_[l].c[0]; int k = 0; int end = n_energy_out - 2; @@ -229,11 +229,11 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* prn_s } // Sampled correlated angle from Kalbach-Mann parameters - if (prn(prn_seed) > km_r) { - double T = (2.0*prn(prn_seed) - 1.0) * std::sinh(km_a); + if (prn(seed) > km_r) { + double T = (2.0*prn(seed) - 1.0) * std::sinh(km_a); mu = std::log(T + std::sqrt(T*T + 1.0))/km_a; } else { - double r1 = prn(prn_seed); + double r1 = prn(seed); mu = std::log(r1*std::exp(km_a) + (1.0 - r1)*std::exp(-km_a))/km_a; } } diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp index 46e56145b..8f6ba80e2 100644 --- a/src/secondary_nbody.cpp +++ b/src/secondary_nbody.cpp @@ -22,38 +22,38 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group) } void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // By definition, the distribution of the angle is isotropic for an N-body // phase space distribution - mu = 2.0*prn(prn_seed) - 1.0; + mu = 2.0*prn(seed) - 1.0; // Determine E_max parameter double Ap = mass_ratio_; double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_); // x is essentially a Maxwellian distribution - double x = maxwell_spectrum(1.0, prn_seed); + double x = maxwell_spectrum(1.0, seed); double y; double r1, r2, r3, r4, r5, r6; switch (n_bodies_) { case 3: - y = maxwell_spectrum(1.0, prn_seed); + y = maxwell_spectrum(1.0, seed); break; case 4: - r1 = prn(prn_seed); - r2 = prn(prn_seed); - r3 = prn(prn_seed); + r1 = prn(seed); + r2 = prn(seed); + r3 = prn(seed); y = -std::log(r1*r2*r3); break; case 5: - r1 = prn(prn_seed); - r2 = prn(prn_seed); - r3 = prn(prn_seed); - r4 = prn(prn_seed); - r5 = prn(prn_seed); - r6 = prn(prn_seed); + r1 = prn(seed); + r2 = prn(seed); + r3 = prn(seed); + r4 = prn(seed); + r5 = prn(seed); + r6 = prn(seed); y = -std::log(r1*r2*r3*r4) - std::log(r5) * std::pow(std::cos(PI/2.0*r6), 2); break; default: diff --git a/src/secondary_thermal.cpp b/src/secondary_thermal.cpp index d7c5ee63a..f61922a82 100644 --- a/src/secondary_thermal.cpp +++ b/src/secondary_thermal.cpp @@ -33,7 +33,7 @@ CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) void CoherentElasticAE::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // Get index and interpolation factor for elastic grid int i; @@ -43,7 +43,7 @@ CoherentElasticAE::sample(double E_in, double& E_out, double& mu, // Sample a Bragg edge between 1 and i const auto& factors = xs_.factors(); - double prob = prn(prn_seed) * factors[i+1]; + double prob = prn(seed) * factors[i+1]; int k = 0; if (prob >= factors.front()) { k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob); @@ -67,11 +67,11 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group) void IncoherentElasticAE::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4 double c = 2 * E_in * debye_waller_; - mu = std::log(1.0 + prn(prn_seed)*(std::exp(2.0*c) - 1))/c - 1.0; + mu = std::log(1.0 + prn(seed)*(std::exp(2.0*c) - 1))/c - 1.0; // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4) E_out = E_in; @@ -90,7 +90,7 @@ IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group, void IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // Get index and interpolation factor for elastic grid int i; @@ -102,7 +102,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu, // Sample outgoing cosine bin int n_mu = mu_out_.shape()[1]; - int k = prn(prn_seed) * n_mu; + int k = prn(seed) * n_mu; // Rather than use the sampled discrete mu directly, it is smeared over // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the @@ -125,7 +125,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu, mu_out_(i, k+1) + f*(mu_out_(i+1, k+1) - mu_out_(i, k+1)); // Smear cosine - mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5); + mu += std::min(mu - mu_left, mu_right - mu)*(prn(seed) - 0.5); // Energy doesn't change in elastic scattering E_out = E_in; @@ -146,7 +146,7 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group, void IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // Get index and interpolation factor for inelastic grid int i; @@ -164,10 +164,10 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu, int n = energy_out_.shape()[1]; if (!skewed_) { // All bins equally likely - j = prn(prn_seed) * n; + j = prn(seed) * n; } else { // Distribution skewed away from edge points - double r = prn(prn_seed) * (n - 3); + double r = prn(seed) * (n - 3); if (r > 1.0) { // equally likely N-4 middle bins j = r + 1; @@ -195,7 +195,7 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu, // Sample outgoing cosine bin int m = mu_out_.shape()[2]; - int k = prn(prn_seed) * m; + int k = prn(seed) * m; // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1] double mu_ijk = mu_out_(i, j, k); @@ -250,7 +250,7 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group) void IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // Get index and interpolation factor for inelastic grid int i; @@ -263,7 +263,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu, // Determine outgoing energy bin // (First reset n_energy_out to the right value) auto n = distribution_[l].n_e_out; - double r1 = prn(prn_seed); + double r1 = prn(seed); double c_j = distribution_[l].e_out_cdf[0]; double c_j1; std::size_t j; @@ -303,7 +303,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu, // Sample outgoing cosine bin int n_mu = distribution_[l].mu.shape()[1]; - std::size_t k = prn(prn_seed) * n_mu; + std::size_t k = prn(seed) * n_mu; // Rather than use the sampled discrete mu directly, it is smeared over // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the @@ -328,7 +328,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu, mu_right = mu_l(j, k+1) + f*(mu_l(j+1, k+1) - mu_l(j, k+1)); // Smear cosine - mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5); + mu += std::min(mu - mu_left, mu_right - mu)*(prn(seed) - 0.5); } } // namespace openmc diff --git a/src/secondary_uncorrelated.cpp b/src/secondary_uncorrelated.cpp index 6f65e0ef2..1421fcf83 100644 --- a/src/secondary_uncorrelated.cpp +++ b/src/secondary_uncorrelated.cpp @@ -53,7 +53,7 @@ UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group) void UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, - uint64_t* prn_seed) const + uint64_t* seed) const { // Sample cosine of scattering angle if (fission_) { @@ -62,14 +62,14 @@ UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu, mu = 1.0; // <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<< } else if (!angle_.empty()) { - mu = angle_.sample(E_in, prn_seed); + mu = angle_.sample(E_in, seed); } else { // no angle distribution given => assume isotropic for all energies - mu = 2.0*prn(prn_seed) - 1.0; + mu = 2.0*prn(seed) - 1.0; } // Sample outgoing energy - E_out = energy_->sample(E_in, prn_seed); + E_out = energy_->sample(E_in, seed); } } // namespace openmc diff --git a/src/simulation.cpp b/src/simulation.cpp index 405102827..2f3230c6e 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -476,7 +476,7 @@ void initialize_history(Particle* p, int64_t index_source) // set random number seed int64_t particle_seed = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + p->id_; - init_particle_seeds(particle_seed, p->prn_seeds_); + init_particle_seeds(particle_seed, p->seeds_); // set particle trace simulation::trace = false; diff --git a/src/source.cpp b/src/source.cpp index 01d333649..628631fdb 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -142,7 +142,7 @@ SourceDistribution::SourceDistribution(pugi::xml_node node) } -Particle::Bank SourceDistribution::sample(uint64_t* prn_seed) const +Particle::Bank SourceDistribution::sample(uint64_t* seed) const { Particle::Bank site; @@ -158,7 +158,7 @@ Particle::Bank SourceDistribution::sample(uint64_t* prn_seed) const site.particle = particle_; // Sample spatial distribution - site.r = space_->sample(prn_seed); + site.r = space_->sample(seed); double xyz[] {site.r.x, site.r.y, site.r.z}; // Now search to see if location exists in geometry @@ -200,7 +200,7 @@ Particle::Bank SourceDistribution::sample(uint64_t* prn_seed) const ++n_accept; // Sample angle - site.u = angle_->sample(prn_seed); + site.u = angle_->sample(seed); // Check for monoenergetic source above maximum particle energy auto p = static_cast(particle_); @@ -218,7 +218,7 @@ Particle::Bank SourceDistribution::sample(uint64_t* prn_seed) const while (true) { // Sample energy spectrum - site.E = energy_->sample(prn_seed); + site.E = energy_->sample(seed); // Resample if energy falls outside minimum or maximum particle energy if (site.E < data::energy_max[p] && site.E > data::energy_min[p]) break; @@ -287,7 +287,7 @@ void initialize_source() } } -Particle::Bank sample_external_source(uint64_t* prn_seed) +Particle::Bank sample_external_source(uint64_t* seed) { // Determine total source strength double total_strength = 0.0; @@ -297,7 +297,7 @@ Particle::Bank sample_external_source(uint64_t* prn_seed) // Sample from among multiple source distributions int i = 0; if (model::external_sources.size() > 1) { - double xi = prn(prn_seed)*total_strength; + double xi = prn(seed)*total_strength; double c = 0.0; for (; i < model::external_sources.size(); ++i) { c += model::external_sources[i].strength(); @@ -306,7 +306,7 @@ Particle::Bank sample_external_source(uint64_t* prn_seed) } // Sample source site from i-th source distribution - Particle::Bank site {model::external_sources[i].sample(prn_seed)}; + Particle::Bank site {model::external_sources[i].sample(seed)}; // If running in MG, convert site.E to group if (!settings::run_CE) { diff --git a/src/surface.cpp b/src/surface.cpp index feea34d8d..abdaacbed 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -197,7 +197,7 @@ Surface::reflect(Position r, Direction u) const } Direction -Surface::diffuse_reflect(Position r, Direction u, uint64_t* prn_seed) const +Surface::diffuse_reflect(Position r, Direction u, uint64_t* seed) const { // Diffuse reflect direction according to the normal. // cosine distribution @@ -208,10 +208,10 @@ Surface::diffuse_reflect(Position r, Direction u, uint64_t* prn_seed) const // sample from inverse function, u=sqrt(rand) since p(u)=2u, so F(u)=u^2 const double mu = (projection>=0.0) ? - -std::sqrt(prn(prn_seed)) : std::sqrt(prn(prn_seed)); + -std::sqrt(prn(seed)) : std::sqrt(prn(seed)); // sample azimuthal distribution uniformly - u = rotate_angle(n, mu, nullptr, prn_seed); + u = rotate_angle(n, mu, nullptr, seed); // normalize the direction return u/u.norm(); diff --git a/src/thermal.cpp b/src/thermal.cpp index 93933e2c1..5b0ee0d20 100644 --- a/src/thermal.cpp +++ b/src/thermal.cpp @@ -151,7 +151,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector& tem void ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic, double* inelastic, - uint64_t* prn_seed) const + uint64_t* seed) const { // Determine temperature for S(a,b) table double kT = sqrtkT*sqrtkT; @@ -173,7 +173,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp, // Randomly sample between temperature i and i+1 double f = (kT - kTs_[i]) / (kTs_[i+1] - kTs_[i]); - if (f > prn(prn_seed)) ++i; + if (f > prn(seed)) ++i; } // Set temperature index @@ -266,13 +266,13 @@ ThermalData::calculate_xs(double E, double* elastic, double* inelastic) const void ThermalData::sample(const NuclideMicroXS& micro_xs, double E, - double* E_out, double* mu, uint64_t* prn_seed) + double* E_out, double* mu, uint64_t* seed) { // Determine whether inelastic or elastic scattering will occur - if (prn(prn_seed) < micro_xs.thermal_elastic / micro_xs.thermal) { - elastic_.distribution->sample(E, *E_out, *mu, prn_seed); + if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) { + elastic_.distribution->sample(E, *E_out, *mu, seed); } else { - inelastic_.distribution->sample(E, *E_out, *mu, prn_seed); + inelastic_.distribution->sample(E, *E_out, *mu, seed); } // Because of floating-point roundoff, it may be possible for mu to be