From 5502f38580071e4bcf106eafeba5c8251e5fd747 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 21 Nov 2019 22:19:15 +0000 Subject: [PATCH] partially done with refactoring --- include/openmc/distribution.h | 16 ++-- include/openmc/distribution_angle.h | 2 +- include/openmc/distribution_energy.h | 14 +-- include/openmc/distribution_multi.h | 8 +- include/openmc/distribution_spatial.h | 10 +- include/openmc/math_functions.h | 26 ++++-- include/openmc/mgxs.h | 9 +- include/openmc/photon.h | 11 ++- include/openmc/physics.h | 10 +- include/openmc/plot.h | 2 +- include/openmc/reaction_product.h | 5 +- include/openmc/scattdata.h | 15 ++- src/bremsstrahlung.cpp | 6 +- src/distribution.cpp | 32 +++---- src/distribution_angle.cpp | 7 +- src/distribution_energy.cpp | 18 ++-- src/distribution_multi.cpp | 14 +-- src/distribution_spatial.cpp | 18 ++-- src/eigenvalue.cpp | 10 +- src/math_functions.cpp | 34 +++---- src/mgxs.cpp | 11 ++- src/nuclide.cpp | 8 +- src/particle.cpp | 6 +- src/photon.cpp | 69 +++++++------- src/physics.cpp | 128 +++++++++++++------------- src/physics_common.cpp | 2 +- src/physics_mg.cpp | 15 +-- src/plot.cpp | 12 ++- src/random_lcg.cpp | 1 - src/reaction_product.cpp | 9 +- src/scattdata.cpp | 29 +++--- 31 files changed, 301 insertions(+), 256 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index fecddc346..28d9b9e81 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -35,7 +35,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; // Properties const std::vector& x() const { return x_; } @@ -59,7 +59,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; private: double a_; //!< Lower bound of distribution double b_; //!< Upper bound of distribution @@ -76,7 +76,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; private: double theta_; //!< Factor in exponential [eV] }; @@ -92,7 +92,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; private: double a_; //!< Factor in exponential [eV] double b_; //!< Factor in square root [1/eV] @@ -109,7 +109,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; private: double mean_value_; //!< middle of distribution [eV] double std_dev_; //!< standard deviation [eV] @@ -127,7 +127,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; private: // example DT fusion m_rat = 5 (D = 2 + T = 3) // ion temp = 20000 eV @@ -149,7 +149,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; // x property std::vector& x() { return x_; } @@ -179,7 +179,7 @@ public: //! Sample a value from the distribution //! \return Sampled value - double sample() const; + double sample(uint64_t * prn_seeds, int stream) const; private: std::vector x_; //! Possible outcomes }; diff --git a/include/openmc/distribution_angle.h b/include/openmc/distribution_angle.h index 4344eb60d..7d8af1fb1 100644 --- a/include/openmc/distribution_angle.h +++ b/include/openmc/distribution_angle.h @@ -24,7 +24,7 @@ public: //! Sample an angle given an incident particle energy //! \param[in] E Particle energy in [eV] //! \return Cosine of the angle in the range [-1,1] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) 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 f13bd5de4..81ea03995 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) const = 0; + virtual double sample(double E, uint64_t * prn_seeds, int stream) const = 0; virtual ~EnergyDistribution() = default; }; @@ -37,7 +37,7 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] //! \return Sampled energy in [eV] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) const; private: int primary_flag_; //!< Indicator of whether the photon is a primary or //!< non-primary photon. @@ -56,7 +56,7 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] //! \return Sampled energy in [eV] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) const; private: double threshold_; //!< Energy threshold in lab, (A + 1)/A * |Q| double mass_ratio_; //!< (A/(A+1))^2 @@ -75,7 +75,7 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] //! \return Sampled energy in [eV] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) const; private: //! Outgoing energy for a single incoming energy struct CTTable { @@ -104,7 +104,7 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] //! \return Sampled energy in [eV] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) const; private: Tabulated1D theta_; //!< Incoming energy dependent parameter double u_; //!< Restriction energy @@ -122,7 +122,7 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] //! \return Sampled energy in [eV] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) const; private: Tabulated1D theta_; //!< Incoming energy dependent parameter double u_; //!< Restriction energy @@ -140,7 +140,7 @@ public: //! Sample energy distribution //! \param[in] E Incident particle energy in [eV] //! \return Sampled energy in [eV] - double sample(double E) const; + double sample(double E, uint64_t * prn_seeds, int stream) 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 487f52831..65e7ae005 100644 --- a/include/openmc/distribution_multi.h +++ b/include/openmc/distribution_multi.h @@ -24,7 +24,7 @@ public: //! Sample a direction from the distribution //! \return Direction sampled - virtual Direction sample() const = 0; + virtual Direction sample(uint64_t * prn_seeds, int stream) const = 0; Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction }; @@ -40,7 +40,7 @@ public: //! Sample a direction from the distribution //! \return Direction sampled - Direction sample() const; + Direction sample(uint64_t * prn_seeds, int stream) const; private: UPtrDist mu_; //!< Distribution of polar angle UPtrDist phi_; //!< Distribution of azimuthal angle @@ -56,7 +56,7 @@ public: //! Sample a direction from the distribution //! \return Sampled direction - Direction sample() const; + Direction sample(uint64_t * prn_seeds, int stream) const; }; //============================================================================== @@ -70,7 +70,7 @@ public: //! Sample a direction from the distribution //! \return Sampled direction - Direction sample() const; + Direction sample(uint64_t * prn_seeds, int stream) const; }; using UPtrAngle = std::unique_ptr; diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 6cf6f52d6..8642fd453 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() const = 0; + virtual Position sample(uint64_t * prn_seeds, int stream) const = 0; }; //============================================================================== @@ -30,7 +30,7 @@ public: //! Sample a position from the distribution //! \return Sampled position - Position sample() const; + Position sample(uint64_t * prn_seeds, int stream) const; private: UPtrDist x_; //!< Distribution of x coordinates UPtrDist y_; //!< Distribution of y coordinates @@ -47,7 +47,7 @@ public: //! Sample a position from the distribution //! \return Sampled position - Position sample() const; + Position sample(uint64_t * prn_seeds, int stream) const; private: UPtrDist r_; //!< Distribution of r coordinates UPtrDist theta_; //!< Distribution of theta coordinates @@ -65,7 +65,7 @@ public: //! Sample a position from the distribution //! \return Sampled position - Position sample() const; + Position sample(uint64_t * prn_seeds, int stream) const; // Properties bool only_fissionable() const { return only_fissionable_; } @@ -87,7 +87,7 @@ public: //! Sample a position from the distribution //! \return Sampled position - Position sample() const; + Position sample(uint64_t * prn_seeds, int stream) 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 39b703ab5..4246e354d 100644 --- a/include/openmc/math_functions.h +++ b/include/openmc/math_functions.h @@ -129,11 +129,15 @@ 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_seeds A pointer to the array of pseudorandom seeds +//! \param stream The pseudorandom stream index with which to sample from //============================================================================== -extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi); +extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi, + uint64_t * prn_seeds, int stream); -Direction rotate_angle(Direction u, double mu, const double* phi); +Direction rotate_angle(Direction u, double mu, const double* phi, + uint64_t * prn_streams, int stream); //============================================================================== //! Samples an energy from the Maxwell fission distribution based on a direct @@ -144,10 +148,12 @@ 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_seeds A pointer to the array of pseudorandom seeds +//! \param stream The pseudorandom stream index with which to sample from //! \return The sampled outgoing energy //============================================================================== -extern "C" double maxwell_spectrum(double T); +extern "C" double maxwell_spectrum(double T, uint64_t * prn_seeds, int stream); //============================================================================== //! Samples an energy from a Watt energy-dependent fission distribution. @@ -159,10 +165,12 @@ extern "C" double maxwell_spectrum(double T); //! //! \param a Watt parameter a //! \param b Watt parameter b +//! \param prn_seeds A pointer to the array of pseudorandom seeds +//! \param stream The pseudorandom stream index with which to sample from //! \return The sampled outgoing energy //============================================================================== -extern "C" double watt_spectrum(double a, double b); +extern "C" double watt_spectrum(double a, double b, uint64_t * prn_seeds, int stream); //============================================================================== //! Samples an energy from the Gaussian energy-dependent fission distribution. @@ -175,10 +183,13 @@ extern "C" double watt_spectrum(double a, double b); //! //! @param mean mean of the Gaussian distribution //! @param std_dev standard deviation of the Gaussian distribution +//! @param prn_seeds A pointer to the array of pseudorandom seeds +//! @param stream The pseudorandom stream index with which to sample from //! @result The sampled outgoing energy //============================================================================== -extern "C" double normal_variate(double mean, double std_dev); +extern "C" double normal_variate(double mean, double std_dev, uint64_t * prn_seeds, + int stream); //============================================================================== //! Samples an energy from the Muir (Gaussian) energy-dependent distribution. @@ -190,10 +201,13 @@ extern "C" double normal_variate(double mean, double std_dev); //! @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_seeds A pointer to the array of pseudorandom seeds +//! @param stream The pseudorandom stream index with which to sample from //! @result The sampled outgoing energy //============================================================================== -extern "C" double muir_spectrum(double e0, double m_rat, double kt); +extern "C" double muir_spectrum(double e0, double m_rat, double kt, uint64_t * prn_seeds, + int stream); //============================================================================== //! Doppler broadens the windowed multipole curvefit. diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index 3dfd6ba4c..cfb696a62 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -154,8 +154,10 @@ class Mgxs { //! @param gin Incoming energy group. //! @param dg Sampled delayed group index. //! @param gout Sampled outgoing energy group. + //! @param prn_seeds Pseudorandom seed array. + //! @param stream Pseudorandom stream index. void - sample_fission_energy(int gin, int& dg, int& gout); + sample_fission_energy(int gin, int& dg, int& gout, uint64_t * prn_seeds, int stream); //! \brief Samples the outgoing energy and angle from a scatter event. //! @@ -163,8 +165,11 @@ 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_seeds Array of pseudorandom stream seeds + //! @param stream Pseudorandom stream index void - sample_scatter(int gin, int& gout, double& mu, double& wgt); + sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, + int stream); //! \brief Calculates cross section quantities needed for tracking. //! diff --git a/include/openmc/photon.h b/include/openmc/photon.h index 0fbdbc0c9..f5b79455c 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) const; + double* mu, int* i_shell, uint64_t * prn_seeds, int stream) const; - double rayleigh_scatter(double alpha) const; + double rayleigh_scatter(double alpha, uint64_t * prn_seeds, int stream) const; void pair_production(double alpha, double* E_electron, double* E_positron, - double* mu_electron, double* mu_positron) const; + double* mu_electron, double* mu_positron, uint64_t * prn_seeds, int stream) const; void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const; @@ -96,14 +96,15 @@ public: xt::xtensor dcs_; private: - void compton_doppler(double alpha, double mu, double* E_out, int* i_shell) const; + void compton_doppler(double alpha, double mu, double* E_out, int* i_shell, + uint64_t * prn_seeds, int stream) const; }; //============================================================================== // Non-member functions //============================================================================== -std::pair klein_nishina(double alpha); +std::pair klein_nishina(double alpha, uint64_t * prn_seeds, int stream); void free_memory_photon(); diff --git a/include/openmc/physics.h b/include/openmc/physics.h index 107380eb8..8dc2650ae 100644 --- a/include/openmc/physics.h +++ b/include/openmc/physics.h @@ -59,7 +59,7 @@ void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_p void absorption(Particle* p, int i_nuclide); -void scatter(Particle*, int i_nuclide); +void scatter(Particle* p, int i_nuclide); //! Treats the elastic scattering of a neutron with a target. void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, @@ -72,15 +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); + Direction v_neut, double xs_eff, double kT, uint64_t * prn_seeds, int stream); //! 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); +Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT, + uint64_t * prn_seeds, int stream); -void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site); +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, + uint64_t * prn_seeds, int stream); //! 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/plot.h b/include/openmc/plot.h index da83f42aa..088699883 100644 --- a/include/openmc/plot.h +++ b/include/openmc/plot.h @@ -285,7 +285,7 @@ void create_voxel(Plot pl); //! Create a randomly generated RGB color //! \return RGBColor with random value -RGBColor random_color(); +RGBColor random_color(uint64_t * prn_seeds, int stream); } // namespace openmc diff --git a/include/openmc/reaction_product.h b/include/openmc/reaction_product.h index 9377b5b6a..23a4799f3 100644 --- a/include/openmc/reaction_product.h +++ b/include/openmc/reaction_product.h @@ -42,7 +42,10 @@ 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 - void sample(double E_in, double& E_out, double& mu) const; + //! \param[inout] prn_seeds Pseudorandom array of stream seeds + //! \param[in] stream Psuedorandom stream index + void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, + int stream) const; Particle::Type particle_; //!< Particle type EmissionMode emission_mode_; //!< Emission mode diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index 1a3a252a5..458a4c24d 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -65,8 +65,11 @@ 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_seeds Array of pseudorandom stream seeds + //! @param stream Pseudorandom stream index virtual void - sample(int gin, int& gout, double& mu, double& wgt) = 0; + sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, + int stream) = 0; //! \brief Initializes the ScattData object from a given scatter and //! multiplicity matrix. @@ -109,8 +112,10 @@ class ScattData { //! @param gin Incoming energy group. //! @param gout Sampled outgoing energy group. //! @param i_gout Sampled outgoing energy group index. + //! @param prn_seeds Array of pseudorandom stream seeds + //! @param stream Pseudorandom stream index void - sample_energy(int gin, int& gout, int& i_gout); + sample_energy(int gin, int& gout, int& i_gout, uint64_t * prn_seeds, int stream); //! \brief Provides a cross section value given certain parameters //! @@ -161,7 +166,7 @@ class ScattDataLegendre: public ScattData { calc_f(int gin, int gout, double mu); void - sample(int gin, int& gout, double& mu, double& wgt); + sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, int stream); size_t get_order() {return dist[0][0].size() - 1;}; @@ -197,7 +202,7 @@ class ScattDataHistogram: public ScattData { calc_f(int gin, int gout, double mu); void - sample(int gin, int& gout, double& mu, double& wgt); + sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, int stream); size_t get_order() {return dist[0][0].size();}; @@ -238,7 +243,7 @@ class ScattDataTabular: public ScattData { calc_f(int gin, int gout, double mu); void - sample(int gin, int& gout, double& mu, double& wgt); + sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, int stream); size_t get_order() {return dist[0][0].size();}; diff --git a/src/bremsstrahlung.cpp b/src/bremsstrahlung.cpp index e6b76d19d..40cd28f7a 100644 --- a/src/bremsstrahlung.cpp +++ b/src/bremsstrahlung.cpp @@ -65,7 +65,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) double y = std::exp(y_l + (y_r - y_l)*f); // Sample number of secondary bremsstrahlung photons - int n = y + prn(); + int n = y + prn(p.prn_seeds, p.stream); *E_lost = 0.0; if (n == 0) return; @@ -73,7 +73,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) // Sample index of the tabulated PDF in the energy grid, j or j+1 double c_max; int i_e; - if (prn() <= f || j == 0) { + if (prn(p.prn_seeds, p.stream) <= f || j == 0) { i_e = j + 1; // Interpolate the maximum value of the CDF at the incoming particle @@ -94,7 +94,7 @@ void thick_target_bremsstrahlung(Particle& p, double* E_lost) for (int i = 0; i < n; ++i) { // Generate a random number r and determine the index i for which // cdf(i) <= r*cdf,max <= cdf(i+1) - double c = prn()*c_max; + double c = prn(p.prn_seeds, p.stream)*c_max; int i_w = lower_bound_index(&mat->cdf(i_e, 0), &mat->cdf(i_e, 0) + i_e, c); // Sample the photon energy diff --git a/src/distribution.cpp b/src/distribution.cpp index 8b77d86b8..0b56f9857 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() const +double Discrete::sample(uint64_t * prn_seeds, int stream) const { int n = x_.size(); if (n > 1) { - double xi = prn(); + double xi = prn(prn_seeds, stream); 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() const +double Uniform::sample(uint64_t * prn_seeds, int stream) const { - return a_ + prn()*(b_ - a_); + return a_ + prn(prn_seeds, stream)*(b_ - a_); } //============================================================================== @@ -88,9 +88,9 @@ Maxwell::Maxwell(pugi::xml_node node) theta_ = std::stod(get_node_value(node, "parameters")); } -double Maxwell::sample() const +double Maxwell::sample(uint64_t * prn_seeds, int stream) const { - return maxwell_spectrum(theta_); + return maxwell_spectrum(theta_, prn_seeds, stream); } //============================================================================== @@ -108,9 +108,9 @@ Watt::Watt(pugi::xml_node node) b_ = params.at(1); } -double Watt::sample() const +double Watt::sample(uint64_t * prn_seeds, int stream) const { - return watt_spectrum(a_, b_); + return watt_spectrum(a_, b_, prn_seeds, stream); } //============================================================================== @@ -127,9 +127,9 @@ Normal::Normal(pugi::xml_node node) std_dev_ = params.at(1); } -double Normal::sample() const +double Normal::sample(uint64_t * prn_seeds, int stream) const { - return normal_variate(mean_value_, std_dev_); + return normal_variate(mean_value_, std_dev_, prn_seeds, stream); } //============================================================================== @@ -147,9 +147,9 @@ Muir::Muir(pugi::xml_node node) kt_ = params.at(2); } -double Muir::sample() const +double Muir::sample(uint64_t * prn_seeds, int stream) const { - return muir_spectrum(e0_, m_rat_, kt_); + return muir_spectrum(e0_, m_rat_, kt_, prn_seeds, stream); } //============================================================================== @@ -220,10 +220,10 @@ void Tabular::init(const double* x, const double* p, std::size_t n, const double } } -double Tabular::sample() const +double Tabular::sample(uint64_t * prn_seeds, int stream) const { // Sample value of CDF - double c = prn(); + double c = prn(prn_seeds, stream); // Find first CDF bin which is above the sampled value double c_i = c_[0]; @@ -263,11 +263,11 @@ double Tabular::sample() const // Equiprobable implementation //============================================================================== -double Equiprobable::sample() const +double Equiprobable::sample(uint64_t * prn_seeds, int stream) const { std::size_t n = x_.size(); - double r = prn(); + double r = prn(prn_seeds, stream); int i = std::floor((n - 1)*r); double xl = x_[i]; diff --git a/src/distribution_angle.cpp b/src/distribution_angle.cpp index 615f84f68..41a316c11 100644 --- a/src/distribution_angle.cpp +++ b/src/distribution_angle.cpp @@ -62,7 +62,8 @@ AngleDistribution::AngleDistribution(hid_t group) } } -double AngleDistribution::sample(double E) const +//double AngleDistribution::sample(double E) const +double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) const { // Determine number of incoming energies auto n = energy_.size(); @@ -83,10 +84,10 @@ double AngleDistribution::sample(double E) const } // Sample between the ith and (i+1)th bin - if (r > prn()) ++i; + if (r > prn(prn_seeds, stream)) ++i; // Sample i-th distribution - double mu = distribution_[i]->sample(); + double mu = distribution_[i]->sample(prn_seeds, stream); // 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 1b9be7346..1db293b0a 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -146,7 +146,7 @@ ContinuousTabular::ContinuousTabular(hid_t group) } // incoming energies } -double ContinuousTabular::sample(double E) const +double ContinuousTabular::sample(double E, uint64_t * prn_seeds, int stream) const { // Read number of interpolation regions and incoming energies bool histogram_interp; @@ -177,7 +177,7 @@ double ContinuousTabular::sample(double E) const if (histogram_interp) { l = i; } else { - l = r > prn() ? i + 1 : i; + l = r > prn(prn_seeds, stream) ? i + 1 : i; } // Interpolation for energy E1 and EK @@ -197,7 +197,7 @@ double ContinuousTabular::sample(double E) const // Determine outgoing energy bin n_energy_out = distribution_[l].e_out.size(); n_discrete = distribution_[l].n_discrete; - double r1 = prn(); + double r1 = prn(prn_seeds, stream); 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) const +double MaxwellEnergy::sample(double E, uint64_t * prn_seeds, int stream) const { // Get temperature corresponding to incoming energy double theta = theta_(E); while (true) { // Sample maxwell fission spectrum - double E_out = maxwell_spectrum(theta); + double E_out = maxwell_spectrum(theta, prn_seeds, stream); // 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) const +double Evaporation::sample(double E, uint64_t * prn_seeds, int stream) const { // Get temperature corresponding to incoming energy double theta = theta_(E); @@ -313,7 +313,7 @@ double Evaporation::sample(double E) const // density function double x; while (true) { - x = -std::log((1.0 - v*prn())*(1.0 - v*prn())); + x = -std::log((1.0 - v*prn(prn_seeds, stream))*(1.0 - v*prn(prn_seeds, stream))); if (x <= y) break; } @@ -338,7 +338,7 @@ WattEnergy::WattEnergy(hid_t group) close_dataset(dset); } -double WattEnergy::sample(double E) const +double WattEnergy::sample(double E, uint64_t * prn_seeds, int stream) const { // Determine Watt parameters at incident energy double a = a_(E); @@ -346,7 +346,7 @@ double WattEnergy::sample(double E) const while (true) { // Sample energy-dependent Watt fission spectrum - double E_out = watt_spectrum(a, b); + double E_out = watt_spectrum(a, b, prn_seeds, stream); // 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 c388254b6..4468febd3 100644 --- a/src/distribution_multi.cpp +++ b/src/distribution_multi.cpp @@ -53,14 +53,14 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node) } } -Direction PolarAzimuthal::sample() const +Direction PolarAzimuthal::sample(uint64_t * prn_seeds, int stream) const { // Sample cosine of polar angle - double mu = mu_->sample(); + double mu = mu_->sample(prn_seeds, stream); if (mu == 1.0) return u_ref_; // Sample azimuthal angle - double phi = phi_->sample(); + double phi = phi_->sample(prn_seeds, stream); // If the reference direction is along the z-axis, rotate the aziumthal angle // to match spherical coordinate conventions. @@ -74,10 +74,10 @@ Direction PolarAzimuthal::sample() const // Isotropic implementation //============================================================================== -Direction Isotropic::sample() const +Direction Isotropic::sample(uint64_t * prn_seeds, int stream) const { - double phi = 2.0*PI*prn(); - double mu = 2.0*prn() - 1.0; + double phi = 2.0*PI*prn(prn_seeds, stream); + double mu = 2.0*prn(prn_seeds, stream) - 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() const // Monodirectional implementation //============================================================================== -Direction Monodirectional::sample() const +Direction Monodirectional::sample(uint64_t * prn_seeds, int stream) const { return u_ref_; } diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index e92601016..9dd264107 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -46,9 +46,9 @@ CartesianIndependent::CartesianIndependent(pugi::xml_node node) } } -Position CartesianIndependent::sample() const +Position CartesianIndependent::sample(uint64_t * prn_seeds, int stream) const { - return {x_->sample(), y_->sample(), z_->sample()}; + return {x_->sample(prn_seeds, stream), y_->sample(prn_seeds, stream), z_->sample(prn_seeds, stream)}; } //============================================================================== @@ -107,11 +107,11 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node) } -Position SphericalIndependent::sample() const +Position SphericalIndependent::sample(uint64_t * prn_seeds, int stream) const { - double r = r_->sample(); - double theta = theta_->sample(); - double phi = phi_->sample(); + double r = r_->sample(prn_seeds, stream); + double theta = theta_->sample(prn_seeds, stream); + double phi = phi_->sample(prn_seeds, stream); 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() const +Position SpatialBox::sample(uint64_t * prn_seeds, int stream) const { - Position xi {prn(), prn(), prn()}; + Position xi {prn(prn_seeds, stream), prn(prn_seeds, stream), prn(prn_seeds, stream)}; 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() const +Position SpatialPoint::sample(uint64_t * prn_seeds, int stream) const { return r_; } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index 894157c5b..9bf197378 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -114,12 +114,16 @@ void synchronize_bank() fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank)); } + // Create pseudorandom number streams + uint64_t prn_seeds[N_STREAMS]; // seeds + int stream = STREAM_TRACKING; // current RNG stream + // Make sure all processors start at the same point for random sampling. Then // skip ahead in the sequence using the starting index in the 'global' // fission bank for each processor. - set_particle_seed(simulation::total_gen + overall_generation()); - advance_prn_seed(start); + set_particle_seed(simulation::total_gen + overall_generation(), prn_seeds); + advance_prn_seed(start, prn_seeds, stream); // TODO: What is starting stream here? Tracking? It doesn't appear to be enforced anywhere? // Determine how many fission sites we need to sample from the source bank // and the probability for selecting a site. @@ -154,7 +158,7 @@ void synchronize_bank() } // Randomly sample sites needed - if (prn() < p_sample) { + if (prn(prn_seeds, stream) < p_sample) { temp_sites[index_temp] = site; ++index_temp; } diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 7525a7ddd..6e639432b 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) { - Direction u = rotate_angle({uvw}, mu, phi); +void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t * prn_seeds, int stream) { + Direction u = rotate_angle({uvw}, mu, phi, prn_seeds, stream); uvw[0] = u.x; uvw[1] = u.y; uvw[2] = u.z; } -Direction rotate_angle(Direction u, double mu, const double* phi) +Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t * prn_seeds, int stream) { // Sample azimuthal angle in [0,2pi) if none provided double phi_; if (phi != nullptr) { phi_ = (*phi); } else { - phi_ = 2.0*PI*prn(); + phi_ = 2.0*PI*prn(prn_seeds, stream); } // Precompute factors to save flops @@ -675,11 +675,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi) } -double maxwell_spectrum(double T) { +double maxwell_spectrum(double T, uint64_t * prn_seeds, int stream) { // Set the random numbers - double r1 = prn(); - double r2 = prn(); - double r3 = prn(); + double r1 = prn(prn_seeds, stream); + double r2 = prn(prn_seeds, stream); + double r3 = prn(prn_seeds, stream); // determine cosine of pi/2*r double c = std::cos(PI / 2. * r3); @@ -691,33 +691,33 @@ double maxwell_spectrum(double T) { } -double normal_variate(double mean, double standard_deviation) { +double normal_variate(double mean, double standard_deviation, uint64_t * prn_seeds, int stream) { // perhaps there should be a limit to the number of resamples while ( true ) { - double v1 = 2 * prn() - 1.; - double v2 = 2 * prn() - 1.; + double v1 = 2 * prn(prn_seeds, stream) - 1.; + double v2 = 2 * prn(prn_seeds, stream) - 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() <= 0.5) ? v1 : v2; + z *= (prn(prn_seeds, stream) <= 0.5) ? v1 : v2; return mean + standard_deviation*z; } } } -double muir_spectrum(double e0, double m_rat, double kt) { +double muir_spectrum(double e0, double m_rat, double kt, uint64_t * prn_seeds, int stream) { // 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); + return normal_variate(e0, sigma, prn_seeds, stream); } -double watt_spectrum(double a, double b) { - double w = maxwell_spectrum(a); - double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w); +double watt_spectrum(double a, double b, uint64_t * prn_seeds, int stream) { + double w = maxwell_spectrum(a, prn_seeds, stream); + double E_out = w + 0.25 * a * a * b + (2. * prn(prn_seeds, stream) - 1.) * std::sqrt(a * a * b * w); return E_out; } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 0dd8c73cc..9dceb5127 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) +Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t * prn_seeds, int stream) { // 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) double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin); // sample random numbers - double xi_pd = prn() * nu_fission; - double xi_gout = prn(); + double xi_pd = prn(prn_seeds, stream) * nu_fission; + double xi_gout = prn(prn_seeds, stream); // Select whether the neutron is prompt or delayed if (xi_pd <= prob_prompt) { @@ -585,7 +585,8 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout) //============================================================================== void -Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt) +Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, + int stream) { // This method assumes that the temperature and angle indices are set // Sample the data @@ -594,7 +595,7 @@ Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt) #else int tid = 0; #endif - xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt); + xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, prn_seeds, stream); } //============================================================================== diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 3b092ddf4..1066b900c 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -581,7 +581,7 @@ void Nuclide::calculate_xs(int i_sab, int i_log_union, double sab_frac, Particle // Randomly sample between temperature i and i+1 f = (kT - kTs_[i_temp]) / (kTs_[i_temp + 1] - kTs_[i_temp]); - if (f > prn()) ++i_temp; + if (f > prn(p.prn_seeds, p.stream)) ++i_temp; break; } @@ -756,11 +756,11 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const // This guarantees the randomness and, at the same time, makes sure we // reuse random numbers for the same nuclide at different temperatures, // therefore preserving correlation of temperature in probability tables. - prn_set_stream(STREAM_URR_PTABLE); + p.stream = STREAM_URR_PTABLE; //TODO: to maintain the same random number stream as the Fortran code this //replaces, the seed is set with i_nuclide_ + 1 instead of i_nuclide_ - double r = future_prn(static_cast(i_nuclide_ + 1)); - prn_set_stream(STREAM_TRACKING); + double r = future_prn(static_cast(i_nuclide_ + 1), p.prn_seeds, p.stream); + p.stream = STREAM_TRACKING; int i_low = 0; while (urr.prob_(i_energy, URR_CUM_PROB, i_low) <= r) {++i_low;}; diff --git a/src/particle.cpp b/src/particle.cpp index b430d3a06..ab058a6c3 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -157,9 +157,9 @@ Particle::transport() while (true) { // Set the random number stream if (type_ == Particle::Type::neutron) { - prn_set_stream(STREAM_TRACKING); + stream = STREAM_TRACKING; } else { - prn_set_stream(STREAM_PHOTON); + stream = STREAM_PHOTON; } // Store pre-collision particle properties @@ -228,7 +228,7 @@ Particle::transport() } else if (macro_xs_.total == 0.0) { d_collision = INFINITY; } else { - d_collision = -std::log(prn()) / macro_xs_.total; + d_collision = -std::log(prn(prn_seeds, stream)) / macro_xs_.total; } // Select smaller of the two distances diff --git a/src/photon.cpp b/src/photon.cpp index ddde4bd6c..cbb4133b2 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) const + double* alpha_out, double* mu, int* i_shell, uint64_t * prn_seeds, int stream) 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); + std::tie(*alpha_out, *mu) = klein_nishina(alpha, prn_seeds, stream); // 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() < form_factor_x / form_factor_xmax) { + if (prn(prn_seeds, stream) < form_factor_x / form_factor_xmax) { if (doppler) { double E_out; - this->compton_doppler(alpha, *mu, &E_out, i_shell); + this->compton_doppler(alpha, *mu, &E_out, i_shell, prn_seeds, stream); *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) const + double* E_out, int* i_shell, uint64_t * prn_seeds, int stream) const { auto n = data::compton_profile_pz.size(); int shell; // index for shell while (true) { // Sample electron shell - double rn = prn(); + double rn = prn(prn_seeds, stream); 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()*c_max; + c = prn(prn_seeds, stream)*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() < 0.5 ? E_out1 : E_out2; + *E_out = prn(prn_seeds, stream) < 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) const +double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t * prn_seeds, int stream) const { double mu; while (true) { @@ -507,7 +507,7 @@ double PhotonInteraction::rayleigh_scatter(double alpha) const double F_max = coherent_int_form_factor_(x2_max); // Sample cumulative distribution - double F = prn()*F_max; + double F = prn(prn_seeds, stream)*F_max; // Determine x^2 corresponding to F const auto& x {coherent_int_form_factor_.x()}; @@ -519,13 +519,14 @@ double PhotonInteraction::rayleigh_scatter(double alpha) const // Calculate mu mu = 1.0 - 2.0*x2/x2_max; - if (prn() < 0.5*(1.0 + mu*mu)) break; + if (prn(prn_seeds, stream) < 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) const + double* E_positron, double* mu_electron, double* mu_positron, uint64_t * prn_seeds, + int stream) const { constexpr double r[] { 122.81, 73.167, 69.228, 67.301, 64.696, 61.228, @@ -592,12 +593,12 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron, double u2 = phi2_max; double e; while (true) { - double rn = prn(); + double rn = prn(prn_seeds, stream); // 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() < u1/(u1 + u2)) { + if (prn(prn_seeds, stream) < u1/(u1 + u2)) { i = 1; // Sample e from pi_1 using the inverse transform method @@ -618,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() <= phi1/phi1_max) break; + if (prn(prn_seeds, stream) <= phi1/phi1_max) break; } else { double phi2 = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4; - if (prn() <= phi2/phi2_max) break; + if (prn(prn_seeds, stream) <= phi2/phi2_max) break; } } @@ -634,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() - 1.0; + double rn = 2.0*prn(prn_seeds, stream) - 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() - 1.0; + rn = 2.0*prn(prn_seeds, stream) - 1.0; *mu_positron = (rn + beta)/(rn*beta + 1.0); } @@ -648,8 +649,8 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl { // If no transitions, assume fluorescent photon from captured free electron if (shell.n_transitions == 0) { - double mu = 2.0*prn() - 1.0; - double phi = 2.0*PI*prn(); + double mu = 2.0*prn(p.prn_seeds, p.stream) - 1.0; + double phi = 2.0*PI*prn(p.prn_seeds, p.stream); Direction u; u.x = mu; u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); @@ -660,7 +661,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl } // Sample transition - double rn = prn(); + double rn = prn(p.prn_seeds, p.stream); double c = 0.0; int i_transition; for (i_transition = 0; i_transition < shell.n_transitions; ++i_transition) { @@ -673,8 +674,8 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl int secondary = shell.transition_subshells(i_transition, 1); // Sample angle isotropically - double mu = 2.0*prn() - 1.0; - double phi = 2.0*PI*prn(); + double mu = 2.0*prn(p.prn_seeds, p.stream) - 1.0; + double phi = 2.0*PI*prn(p.prn_seeds, p.stream); Direction u; u.x = mu; u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); @@ -710,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl // Non-member functions //============================================================================== -std::pair klein_nishina(double alpha) +std::pair klein_nishina(double alpha, uint64_t * prn_seeds, int stream) { double alpha_out, mu; double beta = 1.0 + 2.0*alpha; @@ -719,19 +720,19 @@ std::pair klein_nishina(double alpha) double t = beta/(beta + 8.0); double x; while (true) { - if (prn() < t) { + if (prn(prn_seeds, stream) < t) { // Left branch of flow chart - double r = 2.0*prn(); + double r = 2.0*prn(prn_seeds, stream); x = 1.0 + alpha*r; - if (prn() < 4.0/x*(1.0 - 1.0/x)) { + if (prn(prn_seeds, stream) < 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()); + x = beta/(1.0 + 2.0*alpha*prn(prn_seeds, stream)); mu = 1.0 + (1.0 - x)/alpha; - if (prn() < 0.5*(mu*mu + 1.0/x)) break; + if (prn(prn_seeds, stream) < 0.5*(mu*mu + 1.0/x)) break; } } alpha_out = alpha/x; @@ -739,24 +740,24 @@ std::pair klein_nishina(double alpha) } else { // Koblinger's direct method double gamma = 1.0 - std::pow(beta, -2); - double s = prn()*(4.0/alpha + 0.5*gamma + + double s = prn(prn_seeds, stream)*(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()); + alpha_out = alpha/(1.0 + 2.0*alpha*prn(prn_seeds, stream)); } 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())/beta; + alpha_out = alpha*(1.0 + 2.0*alpha*prn(prn_seeds, stream))/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()); + alpha_out = alpha*std::sqrt(1.0 - gamma*prn(prn_seeds, stream)); } else { // For third term, x = beta^r // Therefore, a' = a/beta^r - alpha_out = alpha/std::pow(beta, prn()); + alpha_out = alpha/std::pow(beta, prn(prn_seeds, stream)); } // Calculate cosine of scattering angle based on basic relation diff --git a/src/physics.cpp b/src/physics.cpp index 0eb8c494e..165975ea7 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -112,9 +112,9 @@ void sample_neutron_reaction(Particle* p) // Create secondary photons if (settings::photon_transport) { - prn_set_stream(STREAM_PHOTON); + p->stream = STREAM_PHOTON; sample_secondary_photons(p, i_nuclide); - prn_set_stream(STREAM_TRACKING); + p->stream = STREAM_TRACKING; } // If survival biasing is being used, the following subroutine adjusts the @@ -133,9 +133,9 @@ void sample_neutron_reaction(Particle* p) // Advance URR seed stream 'N' times after energy changes if (p->E_ != p->E_last_) { - prn_set_stream(STREAM_URR_PTABLE); - advance_prn_seed(data::nuclides.size()); - prn_set_stream(STREAM_TRACKING); + p->stream = STREAM_URR_PTABLE; + advance_prn_seed(data::nuclides.size(), p->prn_seeds, p->stream); + p->stream = STREAM_TRACKING; } // Play russian roulette if survival biasing is turned on @@ -159,7 +159,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, // Sample the number of neutrons produced int nu = static_cast(nu_t); - if (prn() <= (nu_t - nu)) ++nu; + if (prn(p->prn_seeds, p->stream) <= (nu_t - nu)) ++nu; // Begin banking the source neutrons // First, if our bank is full then don't continue @@ -181,7 +181,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx, site.wgt = 1. / weight; // Sample delayed group and angle/energy for fission reaction - sample_fission_neutron(i_nuclide, rx, p->E_, &site); + sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->prn_seeds, p->stream); // Set the delayed group on the particle as well p->delayed_group_ = site.delayed_group; @@ -223,13 +223,13 @@ void sample_photon_reaction(Particle* p) // For tallying purposes, this routine might be called directly. In that // case, we need to sample a reaction via the cutoff variable double prob = 0.0; - double cutoff = prn() * micro.total; + double cutoff = prn(p->prn_seeds, p->stream) * micro.total; // Coherent (Rayleigh) scattering prob += micro.coherent; if (prob > cutoff) { - double mu = element.rayleigh_scatter(alpha); - p->u() = rotate_angle(p->u(), mu, nullptr); + double mu = element.rayleigh_scatter(alpha, p->prn_seeds, p->stream); + p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds, p->stream); p->event_ = EVENT_SCATTER; p->event_mt_ = COHERENT; return; @@ -240,7 +240,7 @@ void sample_photon_reaction(Particle* p) if (prob > cutoff) { double alpha_out, mu; int i_shell; - element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell); + element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p->prn_seeds, p->stream); // Determine binding energy of shell. The binding energy is 0.0 if // doppler broadening is not used. @@ -252,13 +252,13 @@ void sample_photon_reaction(Particle* p) } // Create Compton electron - double phi = 2.0*PI*prn(); + double phi = 2.0*PI*prn(p->prn_seeds, p->stream); double E_electron = (alpha - alpha_out)*MASS_ELECTRON_EV - e_b; int electron = static_cast(Particle::Type::electron); if (E_electron >= settings::energy_cutoff[electron]) { double mu_electron = (alpha - alpha_out*mu) / std::sqrt(alpha*alpha + alpha_out*alpha_out - 2.0*alpha*alpha_out*mu); - Direction u = rotate_angle(p->u(), mu_electron, &phi); + Direction u = rotate_angle(p->u(), mu_electron, &phi, p->prn_seeds, p->stream); p->create_secondary(u, E_electron, Particle::Type::electron); } @@ -272,7 +272,7 @@ void sample_photon_reaction(Particle* p) phi += PI; p->E_ = alpha_out*MASS_ELECTRON_EV; - p->u() = rotate_angle(p->u(), mu, &phi); + p->u() = rotate_angle(p->u(), mu, &phi, p->prn_seeds, p->stream); p->event_ = EVENT_SCATTER; p->event_mt_ = INCOHERENT; return; @@ -304,8 +304,8 @@ void sample_photon_reaction(Particle* p) // model in Serpent 2" by Toni Kaltiaisenaho double mu; while (true) { - double r = prn(); - if (4.0*(1.0 - r)*r >= prn()) { + double r = prn(p->prn_seeds, p->stream); + if (4.0*(1.0 - r)*r >= prn(p->prn_seeds, p->stream)) { double rel_vel = std::sqrt(E_electron * (E_electron + 2.0*MASS_ELECTRON_EV)) / (E_electron + MASS_ELECTRON_EV); mu = (2.0*r + rel_vel - 1.0) / (2.0*rel_vel*r - rel_vel + 1.0); @@ -313,7 +313,7 @@ void sample_photon_reaction(Particle* p) } } - double phi = 2.0*PI*prn(); + double phi = 2.0*PI*prn(p->prn_seeds, p->stream); Direction u; u.x = mu; u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); @@ -341,14 +341,14 @@ void sample_photon_reaction(Particle* p) double E_electron, E_positron; double mu_electron, mu_positron; element.pair_production(alpha, &E_electron, &E_positron, - &mu_electron, &mu_positron); + &mu_electron, &mu_positron, p->prn_seeds, p->stream); // Create secondary electron - Direction u = rotate_angle(p->u(), mu_electron, nullptr); + Direction u = rotate_angle(p->u(), mu_electron, nullptr, p->prn_seeds, p->stream); p->create_secondary(u, E_electron, Particle::Type::electron); // Create secondary positron - u = rotate_angle(p->u(), mu_positron, nullptr); + u = rotate_angle(p->u(), mu_positron, nullptr, p->prn_seeds, p->stream); p->create_secondary(u, E_positron, Particle::Type::positron); p->event_ = EVENT_ABSORB; @@ -382,8 +382,8 @@ void sample_positron_reaction(Particle* p) } // Sample angle isotropically - double mu = 2.0*prn() - 1.0; - double phi = 2.0*PI*prn(); + double mu = 2.0*prn(p->prn_seeds, p->stream) - 1.0; + double phi = 2.0*PI*prn(p->prn_seeds, p->stream); Direction u; u.x = mu; u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi); @@ -401,7 +401,7 @@ void sample_positron_reaction(Particle* p) int sample_nuclide(const Particle* p) { // Sample cumulative distribution function - double cutoff = prn() * p->macro_xs_.total; + double cutoff = prn(p->prn_seeds, p->stream) * p->macro_xs_.total; // Get pointers to nuclide/density arrays const auto& mat {model::materials[p->material_]}; @@ -426,7 +426,7 @@ int sample_nuclide(const Particle* p) int sample_element(Particle* p) { // Sample cumulative distribution function - double cutoff = prn() * p->macro_xs_.total; + double cutoff = prn(p->prn_seeds, p->stream) * p->macro_xs_.total; // Get pointers to elements, densities const auto& mat {model::materials[p->material_]}; @@ -479,7 +479,7 @@ Reaction* sample_fission(int i_nuclide, const Particle* p) int i_temp = p->neutron_xs_[i_nuclide].index_temp; int i_grid = p->neutron_xs_[i_nuclide].index_grid; double f = p->neutron_xs_[i_nuclide].interp_factor; - double cutoff = prn() * p->neutron_xs_[i_nuclide].fission; + double cutoff = prn(p->prn_seeds, p->stream) * p->neutron_xs_[i_nuclide].fission; double prob = 0.0; // Loop through each partial fission reaction type @@ -506,7 +506,7 @@ void sample_photon_product(int i_nuclide, const Particle* p, int* i_rx, int* i_p int i_temp = p->neutron_xs_[i_nuclide].index_temp; int i_grid = p->neutron_xs_[i_nuclide].index_grid; double f = p->neutron_xs_[i_nuclide].interp_factor; - double cutoff = prn() * p->neutron_xs_[i_nuclide].photon_prod; + double cutoff = prn(p->prn_seeds, p->stream) * p->neutron_xs_[i_nuclide].photon_prod; double prob = 0.0; // Loop through each reaction type @@ -554,7 +554,7 @@ void absorption(Particle* p, int i_nuclide) } else { // See if disappearance reaction happens if (p->neutron_xs_[i_nuclide].absorption > - prn() * p->neutron_xs_[i_nuclide].total) { + prn(p->prn_seeds, p->stream) * p->neutron_xs_[i_nuclide].total) { // Score absorption estimate of keff if (settings::run_mode == RUN_MODE_EIGENVALUE) { global_tally_absorption += p->wgt_ * p->neutron_xs_[ @@ -582,7 +582,7 @@ void scatter(Particle* p, int i_nuclide) // For tallying purposes, this routine might be called directly. In that // case, we need to sample a reaction via the cutoff variable - double cutoff = prn() * (micro.total - micro.absorption); + double cutoff = prn(p->prn_seeds, p->stream) * (micro.total - micro.absorption); bool sampled = false; // Calculate elastic cross section if it wasn't precalculated @@ -656,8 +656,8 @@ void scatter(Particle* p, int i_nuclide) int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide]; if (mat->p0_[i_nuc_mat]) { // Sample isotropic-in-lab outgoing direction - double mu = 2.0*prn() - 1.0; - double phi = 2.0*PI*prn(); + double mu = 2.0*prn(p->prn_seeds, p->stream) - 1.0; + double phi = 2.0*PI*prn(p->prn_seeds, p->stream); // Change direction of particle p->u().x = mu; @@ -684,7 +684,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Direction v_t {}; if (!p->neutron_xs_[i_nuclide].use_ptable) { v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n, - p->neutron_xs_[i_nuclide].elastic, kT); + p->neutron_xs_[i_nuclide].elastic, kT, p->prn_seeds, p->stream); } // Velocity of center-of-mass @@ -702,9 +702,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, auto& d = rx.products_[0].distribution_[0]; auto d_ = dynamic_cast(d.get()); if (d_) { - mu_cm = d_->angle().sample(p->E_); + mu_cm = d_->angle().sample(p->E_, p->prn_seeds, p->stream); } else { - mu_cm = 2.0*prn() - 1.0; + mu_cm = 2.0*prn(p->prn_seeds, p->stream) - 1.0; } // Determine direction cosines in CM @@ -713,7 +713,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, // Rotate neutron velocity vector to new angle -- note that the speed of the // neutron in CM does not change in elastic scattering. However, the speed // will change when we convert back to LAB - v_n = vel * rotate_angle(u_cm, mu_cm, nullptr); + v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->prn_seeds, p->stream); // Transform back to LAB frame v_n += v_cm; @@ -742,15 +742,15 @@ void sab_scatter(int i_nuclide, int i_sab, Particle* p) // Sample energy and angle double E_out; - data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p->E_, &E_out, &p->mu_); + data::thermal_scatt[i_sab]->data_[i_temp].sample(micro, p->E_, &E_out, &p->mu_, p->prn_seeds, p->stream); // Set energy to outgoing, change direction of particle p->E_ = E_out; - p->u() = rotate_angle(p->u(), p->mu_, nullptr); + p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->prn_seeds, p->stream); } Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u, - Direction v_neut, double xs_eff, double kT) + Direction v_neut, double xs_eff, double kT, uint64_t * prn_seeds, int stream) { // 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); + return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream); 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); + return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream); } 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); + v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream); 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() < R) return v_target; + if (prn(prn_seeds, stream) < 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()); + double E_t = -kT * std::log(prn(prn_seeds, stream)); // sample a relative energy using the xs cdf - double cdf_rel = cdf_low + prn()*(cdf_up - cdf_low); + double cdf_rel = cdf_low + prn(prn_seeds, stream)*(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); + return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, prn_seeds, stream); } } } @@ -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) +sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t * prn_seeds, int stream) { 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) double mu; while (true) { // Sample two random numbers - double r1 = prn(); - double r2 = prn(); + double r1 = prn(prn_seeds, stream); + double r2 = prn(prn_seeds, stream); - if (prn() < alpha) { + if (prn(prn_seeds, stream) < 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) // e^(-y^2). This can be done with sampling scheme C61 from the Monte // Carlo sampler - double c = std::cos(PI/2.0 * prn()); + double c = std::cos(PI/2.0 * prn(prn_seeds, stream)); 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) double beta_vt = std::sqrt(beta_vt_sq); // Sample cosine of angle between neutron and target velocity - mu = 2.0*prn() - 1.0; + mu = 2.0*prn(prn_seeds, stream) - 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() < accept_prob) break; + if (prn(prn_seeds, stream) < accept_prob) break; } // Determine speed of target nucleus @@ -948,19 +948,19 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT) // Determine velocity vector of target nucleus based on neutron's velocity // and the sampled angle between them - return vt * rotate_angle(u, mu, nullptr); + return vt * rotate_angle(u, mu, nullptr, prn_seeds, stream); } -void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site) +void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t * prn_seeds, int stream) { // 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() - 1.0; + double mu = 2.0 * prn(prn_seeds, stream) - 1.0; // Sample azimuthal angle uniformly in [0,2*pi) - double phi = 2.0*PI*prn(); + double phi = 2.0*PI*prn(prn_seeds, stream); 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() < beta) { + if (prn(prn_seeds, stream) < beta) { // ==================================================================== // DELAYED NEUTRON SAMPLED // sampled delayed precursor group - double xi = prn()*nu_d; + double xi = prn(prn_seeds, stream)*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); + rx->products_[group].sample(E_in, site->E, mu, prn_seeds, stream); // 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); + rx->products_[0].sample(E_in, site->E, mu, prn_seeds, stream); // resample if energy is greater than maximum neutron energy constexpr int neutron = static_cast(Particle::Type::neutron); @@ -1050,7 +1050,7 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) // sample outgoing energy and scattering cosine double E; double mu; - rx->products_[0].sample(E_in, E, mu); + rx->products_[0].sample(E_in, E, mu, p->prn_seeds, p->stream); // if scattering system is in center-of-mass, transfer cosine of scattering // angle and outgoing energy from CM to LAB @@ -1076,7 +1076,7 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p) p->mu_ = mu; // change direction of particle - p->u() = rotate_angle(p->u(), mu, nullptr); + p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds, p->stream); // evaluate yield double yield = (*rx->products_[0].yield_)(E_in); @@ -1097,7 +1097,7 @@ void sample_secondary_photons(Particle* p, int i_nuclide) double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod / p->neutron_xs_[i_nuclide].total; int y = static_cast(y_t); - if (prn() <= y_t - y) ++y; + if (prn(p->prn_seeds, p->stream) <= y_t - y) ++y; // Sample each secondary photon for (int i = 0; i < y; ++i) { @@ -1110,10 +1110,10 @@ void sample_secondary_photons(Particle* p, int i_nuclide) auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx]; double E; double mu; - rx->products_[i_product].sample(p->E_, E, mu); + rx->products_[i_product].sample(p->E_, E, mu, p->prn_seeds, p->stream); // Sample the new direction - Direction u = rotate_angle(p->u(), mu, nullptr); + Direction u = rotate_angle(p->u(), mu, nullptr, p->prn_seeds, p->stream); // Create the secondary photon p->create_secondary(u, E, Particle::Type::photon); diff --git a/src/physics_common.cpp b/src/physics_common.cpp index 42ddc2ecb..cfb4c01fe 100644 --- a/src/physics_common.cpp +++ b/src/physics_common.cpp @@ -12,7 +12,7 @@ namespace openmc { void russian_roulette(Particle* p) { if (p->wgt_ < settings::weight_cutoff) { - if (prn() < p->wgt_ / settings::weight_survive) { + if (prn(p->prn_seeds, p->stream) < p->wgt_ / settings::weight_survive) { p->wgt_ = settings::weight_survive; p->wgt_last_ = p->wgt_; } else { diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 8808adbb2..1d22a886b 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -78,10 +78,10 @@ void scatter(Particle* p) { data::mg.macro_xs_[p->material_].sample_scatter(p->g_last_, p->g_, p->mu_, - p->wgt_); + p->wgt_, p->prn_seeds, p->stream); // Rotate the angle - p->u() = rotate_angle(p->u(), p->mu_, nullptr); + p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->prn_seeds, p->stream); // Update energy value for downstream compatability (in tallying) p->E_ = data::mg.energy_bin_avg_[p->g_]; @@ -103,7 +103,7 @@ create_fission_sites(Particle* p, std::vector& bank) // Sample the number of neutrons produced int nu = static_cast(nu_t); - if (prn() <= (nu_t - int(nu_t))) { + if (prn(p->prn_seeds, p->stream) <= (nu_t - int(nu_t))) { nu++; } @@ -128,10 +128,10 @@ create_fission_sites(Particle* p, std::vector& bank) // Sample the cosine of the angle, assuming fission neutrons are emitted // isotropically - double mu = 2.*prn() - 1.; + double mu = 2.*prn(p->prn_seeds, p->stream) - 1.; // Sample the azimuthal angle uniformly in [0, 2.pi) - double phi = 2. * PI * prn(); + double phi = 2. * PI * prn(p->prn_seeds, p->stream ); site.u.x = mu; site.u.y = std::sqrt(1. - mu * mu) * std::cos(phi); site.u.z = std::sqrt(1. - mu * mu) * std::sin(phi); @@ -139,7 +139,8 @@ create_fission_sites(Particle* p, std::vector& bank) // Sample secondary energy distribution for the fission reaction int dg; int gout; - data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout); + data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout, p->prn_seeds, + p->stream); // Store the energy and delayed groups on the fission bank site.E = gout; // We add 1 to the delayed_group bc in MG, -1 is prompt, but in the rest @@ -179,7 +180,7 @@ absorption(Particle* p) global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission / p->macro_xs_.absorption; } else { - if (p->macro_xs_.absorption > prn() * p->macro_xs_.total) { + if (p->macro_xs_.absorption > prn(p->prn_seeds, p->stream) * p->macro_xs_.total) { #pragma omp atomic global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission / p->macro_xs_.absorption; diff --git a/src/plot.cpp b/src/plot.cpp index ba63da16b..cd8aa2c27 100644 --- a/src/plot.cpp +++ b/src/plot.cpp @@ -396,11 +396,15 @@ Plot::set_default_colors(pugi::xml_node plot_node) fatal_error(err_msg); } + // Initialize a random number stream + uint64_t prn_seed = 1; + int stream = 0; + for (auto& c : colors_) { - c = random_color(); + c = random_color(&prn_seed, stream); // make sure we don't interfere with some default colors while (c == RED || c == WHITE) { - c = random_color(); + c = random_color(&prn_seed, stream); } } } @@ -958,8 +962,8 @@ voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace) H5Sclose(memspace); } -RGBColor random_color() { - return {int(prn()*255), int(prn()*255), int(prn()*255)}; +RGBColor random_color(uint64_t * prn_seeds, int stream) { + return {int(prn(prn_seeds, stream)*255), int(prn(prn_seeds, stream)*255), int(prn(prn_seeds, stream)*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 fdd1547c6..62e8e17a3 100644 --- a/src/random_lcg.cpp +++ b/src/random_lcg.cpp @@ -141,7 +141,6 @@ prn_set_stream(int i) extern "C" int64_t openmc_get_seed() {return master_seed;} -// TODO: The idea of this function no longer works -- we'll need to update it somehow extern "C" void openmc_set_seed(int64_t new_seed) { diff --git a/src/reaction_product.cpp b/src/reaction_product.cpp index 34680a773..e047be8de 100644 --- a/src/reaction_product.cpp +++ b/src/reaction_product.cpp @@ -76,25 +76,26 @@ ReactionProduct::ReactionProduct(hid_t group) } } -void ReactionProduct::sample(double E_in, double& E_out, double& mu) const +void ReactionProduct::sample(double E_in, double& E_out, double& mu, + uint64_t * prn_seeds, int stream) const { auto n = applicability_.size(); if (n > 1) { double prob = 0.0; - double c = prn(); + double c = prn(prn_seeds, stream); 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); + distribution_[i]->sample(E_in, E_out, mu, prn_seeds, stream); break; } } } else { // If only one distribution is present, go ahead and sample it - distribution_[0]->sample(E_in, E_out, mu); + distribution_[0]->sample(E_in, E_out, mu, prn_seeds, stream); } } diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 5fc19c41e..a5eeae621 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) +ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t * prn_seeds, int stream) { // Sample the outgoing group - double xi = prn(); + double xi = prn(prn_seeds, stream); double prob = 0.; i_gout = 0; for (gout = gmin[gin]; gout < gmax[gin]; ++gout) { @@ -347,21 +347,22 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu) //============================================================================== void -ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt) +ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt, + uint64_t * prn_seeds, int stream) { // Sample the outgoing energy using the base-class method int i_gout; - sample_energy(gin, gout, i_gout); + sample_energy(gin, gout, i_gout, prn_seeds, stream); // 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() - 1.; + mu = 2. * prn(prn_seeds, stream) - 1.; double f = calc_f(gin, gout, mu); if (f > 0.) { - double u = prn() * M; + double u = prn(prn_seeds, stream) * M; if (u <= f) break; } } @@ -535,14 +536,15 @@ ScattDataHistogram::calc_f(int gin, int gout, double mu) //============================================================================== void -ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) +ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt, + uint64_t * prn_seeds, int stream) { // Sample the outgoing energy using the base-class method int i_gout; - sample_energy(gin, gout, i_gout); + sample_energy(gin, gout, i_gout, prn_seeds, stream); // Determine the outgoing cosine bin - double xi = prn(); + double xi = prn(prn_seeds, stream); int imu; if (xi < dist[gin][i_gout][0]) { @@ -554,7 +556,7 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt) } // Randomly select mu within the imu bin - mu = prn() * dmu + this->mu[imu]; + mu = prn(prn_seeds, stream) * dmu + this->mu[imu]; if (mu < -1.) { mu = -1.; @@ -738,15 +740,16 @@ ScattDataTabular::calc_f(int gin, int gout, double mu) //============================================================================== void -ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt) +ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt, + uint64_t * prn_seeds, int stream) { // Sample the outgoing energy using the base-class method int i_gout; - sample_energy(gin, gout, i_gout); + sample_energy(gin, gout, i_gout, prn_seeds, stream); // Determine the outgoing cosine bin int NP = this->mu.shape()[0]; - double xi = prn(); + double xi = prn(prn_seeds, stream); double c_k = dist[gin][i_gout][0]; int k;