partially done with refactoring

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<double>& 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<double>& 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<double> x_; //! Possible outcomes
 };

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

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

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<UnitSphereDistribution>;

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
 };

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.

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

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<double, 2> 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<double, double> klein_nishina(double alpha);
+std::pair<double, double> klein_nishina(double alpha, uint64_t * prn_seeds, int stream);
 
 void free_memory_photon();
 

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.

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

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

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();};