crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Merge pull request #9 from jtramm/particle_owned_seeds2

Browse source 
CI
This commit is contained in:
John Tramm 2019-12-05 11:57:07 -06:00 • committed by GitHub
commit c061500aeb
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
59 changed files with 535 additions and 756 deletions
Download patch file Download diff file Expand all files Collapse all files

4
examples/xml/pincell/settings.xml
View file

@ -3,8 +3,8 @@
<!-- Define how many particles to run and for how many batches -->
<run_mode>eigenvalue</run_mode>
<batches>10</batches>
<inactive>5</inactive>
<batches>100</batches>
<inactive>10</inactive>
<particles>1000</particles>
<!-- The starting source is a uniform distribution over the entire pin

4
include/openmc/angle_energy.h
View file

@ -14,8 +14,8 @@ namespace openmc {
class AngleEnergy {
public:
virtual void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const = 0;
virtual void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const = 0;
virtual ~AngleEnergy() = default;
};

42
include/openmc/distribution.h
View file

@ -21,7 +21,7 @@ namespace openmc {
class Distribution {
public:
virtual ~Distribution() = default;
virtual double sample(uint64_t * prn_seeds, int stream) const = 0;
virtual double sample(uint64_t* prn_seed) const = 0;
};
//==============================================================================
@ -34,10 +34,9 @@ public:
Discrete(const double* x, const double* p, int n);
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
// Properties
const std::vector<double>& x() const { return x_; }
@ -60,10 +59,9 @@ public:
Uniform(double a, double b) : a_{a}, b_{b} {};
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
private:
double a_; //!< Lower bound of distribution
double b_; //!< Upper bound of distribution
@ -79,10 +77,9 @@ public:
Maxwell(double theta) : theta_{theta} { };
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
private:
double theta_; //!< Factor in exponential [eV]
};
@ -97,10 +94,9 @@ public:
Watt(double a, double b) : a_{a}, b_{b} { };
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
private:
double a_; //!< Factor in exponential [eV]
double b_; //!< Factor in square root [1/eV]
@ -116,10 +112,9 @@ public:
Normal(double mean_value, double std_dev) : mean_value_{mean_value}, std_dev_{std_dev} { };
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
private:
double mean_value_; //!< middle of distribution [eV]
double std_dev_; //!< standard deviation [eV]
@ -136,10 +131,9 @@ public:
Muir(double e0, double m_rat, double kt) : e0_{e0}, m_rat_{m_rat}, kt_{kt} { };
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
private:
// example DT fusion m_rat = 5 (D = 2 + T = 3)
// ion temp = 20000 eV
@ -160,10 +154,9 @@ public:
const double* c=nullptr);
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
// x property
std::vector<double>& x() { return x_; }
@ -192,10 +185,9 @@ public:
Equiprobable(const double* x, int n) : x_{x, x+n} { };
//! Sample a value from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t * prn_seeds, int stream) const;
double sample(uint64_t* prn_seed) const;
private:
std::vector<double> x_; //! Possible outcomes
};

5
include/openmc/distribution_angle.h
View file

@ -23,10 +23,9 @@ public:
//! Sample an angle given an incident particle energy
//! \param[in] E Particle energy in [eV]
//! \param[inout] prn_seeds Array of pseudorandom number seeds
//! \param[in] stream Pseudorandom stream index
//! \param[inout] prn_seed pseudorandom number seed pointer
//! \return Cosine of the angle in the range [-1,1]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
//! Determine whether angle distribution is empty
//! \return Whether distribution is empty

32
include/openmc/distribution_energy.h
View file

@ -22,7 +22,7 @@ namespace openmc {
class EnergyDistribution {
public:
virtual double sample(double E, uint64_t * prn_seeds, int stream) const = 0;
virtual double sample(double E, uint64_t* prn_seed) const = 0;
virtual ~EnergyDistribution() = default;
};
@ -36,10 +36,9 @@ public:
//! Sample energy distribution
//! \param[in] E Incident particle energy in [eV]
//! \param[inout] prn_seeds Pseudorandom number genererator seed array
//! \param[in] stream Pseudorandom number genererator stream index
//! \param[inout] prn_seed Pseudorandom number seed pointer
//! \return Sampled energy in [eV]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
private:
int primary_flag_; //!< Indicator of whether the photon is a primary or
//!< non-primary photon.
@ -57,10 +56,9 @@ public:
//! Sample energy distribution
//! \param[in] E Incident particle energy in [eV]
//! \param[inout] prn_seeds Pseudorandom number genererator seed array
//! \param[in] stream Pseudorandom number genererator stream index
//! \param[inout] prn_seed Pseudorandom number seed pointer
//! \return Sampled energy in [eV]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
private:
double threshold_; //!< Energy threshold in lab, (A + 1)/A * |Q|
double mass_ratio_; //!< (A/(A+1))^2
@ -78,10 +76,9 @@ public:
//! Sample energy distribution
//! \param[in] E Incident particle energy in [eV]
//! \param[inout] prn_seeds Pseudorandom number genererator seed array
//! \param[in] stream Pseudorandom number genererator stream index
//! \param[inout] prn_seed Pseudorandom number seed pointer
//! \return Sampled energy in [eV]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
private:
//! Outgoing energy for a single incoming energy
struct CTTable {
@ -109,10 +106,9 @@ public:
//! Sample energy distribution
//! \param[in] E Incident particle energy in [eV]
//! \param[inout] prn_seeds Pseudorandom number genererator seed array
//! \param[in] stream Pseudorandom number genererator stream index
//! \param[inout] prn_seed Pseudorandom number seed pointer
//! \return Sampled energy in [eV]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
private:
Tabulated1D theta_; //!< Incoming energy dependent parameter
double u_; //!< Restriction energy
@ -129,10 +125,9 @@ public:
//! Sample energy distribution
//! \param[in] E Incident particle energy in [eV]
//! \param[inout] prn_seeds Pseudorandom number genererator seed array
//! \param[in] stream Pseudorandom number genererator stream index
//! \param[inout] prn_seed Pseudorandom number seed pointer
//! \return Sampled energy in [eV]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
private:
Tabulated1D theta_; //!< Incoming energy dependent parameter
double u_; //!< Restriction energy
@ -149,10 +144,9 @@ public:
//! Sample energy distribution
//! \param[in] E Incident particle energy in [eV]
//! \param[inout] prn_seeds Pseudorandom number genererator seed array
//! \param[in] stream Pseudorandom number genererator stream index
//! \param[inout] prn_seed Pseudorandom number seed pointer
//! \return Sampled energy in [eV]
double sample(double E, uint64_t * prn_seeds, int stream) const;
double sample(double E, uint64_t* prn_seed) const;
private:
Tabulated1D a_; //!< Energy-dependent 'a' parameter
Tabulated1D b_; //!< Energy-dependent 'b' parameter

20
include/openmc/distribution_multi.h
View file

@ -23,10 +23,9 @@ public:
virtual ~UnitSphereDistribution() = default;
//! Sample a direction from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Direction sampled
virtual Direction sample(uint64_t * prn_seeds, int stream) const = 0;
virtual Direction sample(uint64_t* prn_seed) const = 0;
Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
};
@ -41,10 +40,9 @@ public:
explicit PolarAzimuthal(pugi::xml_node node);
//! Sample a direction from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Direction sampled
Direction sample(uint64_t * prn_seeds, int stream) const;
Direction sample(uint64_t* prn_seed) const;
private:
UPtrDist mu_; //!< Distribution of polar angle
UPtrDist phi_; //!< Distribution of azimuthal angle
@ -59,10 +57,9 @@ public:
Isotropic() { };
//! Sample a direction from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled direction
Direction sample(uint64_t * prn_seeds, int stream) const;
Direction sample(uint64_t* prn_seed) const;
};
//==============================================================================
@ -75,10 +72,9 @@ public:
explicit Monodirectional(pugi::xml_node node) : UnitSphereDistribution{node} { };
//! Sample a direction from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled direction
Direction sample(uint64_t * prn_seeds, int stream) const;
Direction sample(uint64_t* prn_seed) const;
};
using UPtrAngle = std::unique_ptr<UnitSphereDistribution>;

22
include/openmc/distribution_spatial.h
View file

@ -17,7 +17,7 @@ public:
virtual ~SpatialDistribution() = default;
//! Sample a position from the distribution
virtual Position sample(uint64_t * prn_seeds, int stream) const = 0;
virtual Position sample(uint64_t* prn_seed) const = 0;
};
//==============================================================================
@ -29,10 +29,9 @@ public:
explicit CartesianIndependent(pugi::xml_node node);
//! Sample a position from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t * prn_seeds, int stream) const;
Position sample(uint64_t* prn_seed) const;
private:
UPtrDist x_; //!< Distribution of x coordinates
UPtrDist y_; //!< Distribution of y coordinates
@ -48,10 +47,9 @@ public:
explicit SphericalIndependent(pugi::xml_node node);
//! Sample a position from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t * prn_seeds, int stream) const;
Position sample(uint64_t* prn_seed) const;
private:
UPtrDist r_; //!< Distribution of r coordinates
UPtrDist theta_; //!< Distribution of theta coordinates
@ -68,10 +66,9 @@ public:
explicit SpatialBox(pugi::xml_node node, bool fission=false);
//! Sample a position from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t * prn_seeds, int stream) const;
Position sample(uint64_t* prn_seed) const;
// Properties
bool only_fissionable() const { return only_fissionable_; }
@ -92,10 +89,9 @@ public:
explicit SpatialPoint(pugi::xml_node node);
//! Sample a position from the distribution
//! \param prn_seeds Array of pseudorandom number seeds
//! \param stream Pseudorandom stream index
//! \param prn_seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t * prn_seeds, int stream) const;
Position sample(uint64_t* prn_seed) const;
private:
Position r_; //!< Single position at which sites are generated
};

2
include/openmc/geometry.h
View file

@ -28,14 +28,12 @@ extern std::vector<int64_t> overlap_check_count;
// Information about nearest boundary crossing
//==============================================================================
/*
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface_index {0}; //!< if boundary is surface, index in surfaces vector
int coord_level; //!< coordinate level after crossing boundary
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
*/
//==============================================================================
//! Check two distances by coincidence tolerance

28
include/openmc/math_functions.h
View file

@ -129,15 +129,14 @@ extern "C" void calc_zn_rad(int n, double rho, double zn_rad[]);
//! \param mu The cosine of angle in lab or CM
//! \param phi The azimuthal angle; will randomly chosen angle if a nullptr
//! is passed
//! \param prn_seeds A pointer to the array of pseudorandom seeds
//! \param stream The pseudorandom stream index with which to sample from
//! \param prn_seed A pointer to the pseudorandom seed
//==============================================================================
extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi,
uint64_t * prn_seeds, int stream);
uint64_t* prn_seed);
Direction rotate_angle(Direction u, double mu, const double* phi,
uint64_t * prn_streams, int stream);
uint64_t* prn_seed);
//==============================================================================
//! Samples an energy from the Maxwell fission distribution based on a direct
@ -148,12 +147,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi,
//! rule C64 in the Monte Carlo Sampler LA-9721-MS.
//!
//! \param T The tabulated function of the incoming energy
//! \param prn_seeds A pointer to the array of pseudorandom seeds
//! \param stream The pseudorandom stream index with which to sample from
//! \param prn_seed A pointer to the pseudorandom seed
//! \return The sampled outgoing energy
//==============================================================================
extern "C" double maxwell_spectrum(double T, uint64_t * prn_seeds, int stream);
extern "C" double maxwell_spectrum(double T, uint64_t* prn_seed);
//==============================================================================
//! Samples an energy from a Watt energy-dependent fission distribution.
@ -165,12 +163,11 @@ extern "C" double maxwell_spectrum(double T, uint64_t * prn_seeds, int stream);
//!
//! \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
//! \param prn_seed A pointer to the pseudorandom seed
//! \return The sampled outgoing energy
//==============================================================================
extern "C" double watt_spectrum(double a, double b, uint64_t * prn_seeds, int stream);
extern "C" double watt_spectrum(double a, double b, uint64_t* prn_seed);
//==============================================================================
//! Samples an energy from the Gaussian energy-dependent fission distribution.
@ -183,13 +180,11 @@ extern "C" double watt_spectrum(double a, double b, uint64_t * prn_seeds, int st
//!
//! @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
//! @param prn_seed A pointer to the pseudorandom seed
//! @result The sampled outgoing energy
//==============================================================================
extern "C" double normal_variate(double mean, double std_dev, uint64_t * prn_seeds,
int stream);
extern "C" double normal_variate(double mean, double std_dev, uint64_t* prn_seed);
//==============================================================================
//! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
@ -201,13 +196,12 @@ extern "C" double normal_variate(double mean, double std_dev, uint64_t * prn_see
//! @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
//! @param prn_seed A pointer to the pseudorandom seed
//! @result The sampled outgoing energy
//==============================================================================
extern "C" double muir_spectrum(double e0, double m_rat, double kt,
uint64_t * prn_seeds, int stream);
uint64_t* prn_seed);
//==============================================================================
//! Doppler broadens the windowed multipole curvefit.

11
include/openmc/mgxs.h
View file

@ -154,10 +154,9 @@ 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.
//! @param prn_seed Pseudorandom seed pointer
void
sample_fission_energy(int gin, int& dg, int& gout, uint64_t * prn_seeds, int stream);
sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed);
//! \brief Samples the outgoing energy and angle from a scatter event.
//!
@ -165,11 +164,9 @@ class Mgxs {
//! @param gout Sampled outgoing energy group.
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
//! @param prn_seeds Array of pseudorandom stream seeds
//! @param stream Pseudorandom stream index
//! @param prn_seed Pseudorandom seed pointer.
void
sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds,
int stream);
sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
//! \brief Calculates cross section quantities needed for tracking.
//!

24
include/openmc/particle.h
View file

@ -129,21 +129,8 @@ struct MacroXS {
double incoherent; //!< macroscopic incoherent xs
double photoelectric; //!< macroscopic photoelectric xs
double pair_production; //!< macroscopic pair production xs
int i_grid;
};
//==============================================================================
// Information about nearest boundary crossing
//==============================================================================
struct BoundaryInfo {
double distance {INFINITY}; //!< distance to nearest boundary
int surface_index {0}; //!< if boundary is surface, index in surfaces vector
int coord_level; //!< coordinate level after crossing boundary
std::array<int, 3> lattice_translation {}; //!< which way lattice indices will change
};
//============================================================================
//! State of a particle being transported through geometry
//============================================================================
@ -288,7 +275,6 @@ public:
// Indices for various arrays
int surface_ {0}; //!< index for surface particle is on
BoundaryInfo boundary_;
int cell_born_ {-1}; //!< index for cell particle was born in
int material_ {-1}; //!< index for current material
int material_last_ {-1}; //!< index for last material
@ -313,14 +299,4 @@ public:
} // namespace openmc
extern template class std::vector<openmc::Particle::Bank>;
namespace openmc{
extern std::vector<Particle> particle_bank;
#pragma omp threadprivate(particle_bank)
}
#endif // OPENMC_PARTICLE_H

10
include/openmc/photon.h
View file

@ -45,12 +45,12 @@ public:
void calculate_xs(Particle& p) const;
void compton_scatter(double alpha, bool doppler, double* alpha_out,
double* mu, int* i_shell, uint64_t * prn_seeds, int stream) const;
double* mu, int* i_shell, uint64_t* prn_seed) const;
double rayleigh_scatter(double alpha, uint64_t * prn_seeds, int stream) const;
double rayleigh_scatter(double alpha, uint64_t* prn_seed) const;
void pair_production(double alpha, double* E_electron, double* E_positron,
double* mu_electron, double* mu_positron, uint64_t * prn_seeds, int stream) const;
double* mu_electron, double* mu_positron, uint64_t* prn_seed) const;
void atomic_relaxation(const ElectronSubshell& shell, Particle& p) const;
@ -97,14 +97,14 @@ public:
private:
void compton_doppler(double alpha, double mu, double* E_out, int* i_shell,
uint64_t * prn_seeds, int stream) const;
uint64_t* prn_seed) const;
};
//==============================================================================
// Non-member functions
//==============================================================================
std::pair<double, double> klein_nishina(double alpha, uint64_t * prn_seeds, int stream);
std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seed);
void free_memory_photon();

6
include/openmc/physics.h
View file

@ -72,17 +72,17 @@ void sab_scatter(int i_nuclide, int i_sab, Particle* p);
//! dependence of cross sections in treating resonance elastic scattering such
//! as the DBRC and a new, accelerated scheme are also implemented here.
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
Direction v_neut, double xs_eff, double kT, uint64_t * prn_seeds, int stream);
Direction v_neut, double xs_eff, double kT, uint64_t* prn_seed);
//! samples a target velocity based on the free gas scattering formulation, used
//! by most Monte Carlo codes, in which cross section is assumed to be constant
//! in energy. Excellent documentation for this method can be found in
//! FRA-TM-123.
Direction sample_cxs_target_velocity(double awr, double E, Direction u, double kT,
uint64_t * prn_seeds, int stream);
uint64_t* prn_seed);
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in,
Particle::Bank* site, uint64_t * prn_seeds, int stream);
Particle::Bank* site, uint64_t* prn_seed);
//! handles all reactions with a single secondary neutron (other than fission),
//! i.e. level scattering, (n,np), (n,na), etc.

53
include/openmc/random_lcg.h
View file

@ -10,53 +10,64 @@ namespace openmc {
// Module constants.
//==============================================================================
constexpr int N_STREAMS = 6;
extern "C" const int STREAM_TRACKING;
extern "C" const int STREAM_TALLIES;
extern "C" const int STREAM_SOURCE;
extern "C" const int STREAM_URR_PTABLE;
extern "C" const int STREAM_VOLUME;
extern "C" const int STREAM_PHOTON;
constexpr int64_t DEFAULT_SEED = 1;
constexpr int N_STREAMS {6};
constexpr int STREAM_TRACKING {0};
constexpr int STREAM_TALLIES {1};
constexpr int STREAM_SOURCE {2};
constexpr int STREAM_URR_PTABLE {3};
constexpr int STREAM_VOLUME {4};
constexpr int STREAM_PHOTON {5};
constexpr int64_t DEFAULT_SEED {1};
//==============================================================================
//! Generate a pseudo-random number using a linear congruential generator.
//! @param prn_seeds Pseudorandom number seed array
//! @param stream Pseudorandom number stream index
//! @param prn_seed Pseudorandom number seed pointer
//! @return A random number between 0 and 1
//==============================================================================
extern "C" double prn(uint64_t * seeds, int stream);
double prn(uint64_t* prn_seed);
//==============================================================================
//! Generate a random number which is 'n' times ahead from the current seed.
//!
//! The result of this function will be the same as the result from calling
//! `prn()` 'n' times.
//! `prn()` 'n' times, though without the side effect of altering the RNG
//! state.
//! @param n The number of RNG seeds to skip ahead by
//! @param prn_seeds Pseudorandom number seed array
//! @param stream Pseudorandom number stream index
//! @param prn_seed Pseudorandom number seed
//! @return A random number between 0 and 1
//==============================================================================
extern "C" double future_prn(int64_t n, uint64_t * prn_seeds, int stream);
double future_prn(int64_t n, uint64_t prn_seed);
//==============================================================================
//! Set the RNG seeds to unique values based on the ID of the particle.
//! Set a RNG seed to a unique value based on a unique particle ID by striding
//! the seed.
//! @param prn_seeds Pseudorandom number seed array
//! @param id The particle ID
//==============================================================================
extern "C" void set_particle_seed(int64_t id, uint64_t * prn_seeds );
void init_seed(int64_t id, uint64_t* prn_seeds, int offset );
//==============================================================================
//! Advance the random number seed 'n' times from the current seed.
//! Set the RNG seeds to unique values based on the ID of the particle. This
//! function initializes the seeds for all RNG streams of the particle via
//! striding.
//! @param prn_seeds Pseudorandom number seed array
//! @param stream Pseudorandom number stream index
//! @param id The particle ID
//==============================================================================
void init_particle_seeds(int64_t id, uint64_t* prn_seeds );
//==============================================================================
//! Advance the random number seed 'n' times from the current seed. This
//! differs from the future_prn() function in that this function does alter
//! the RNG state.
//! @param prn_seed Pseudorandom number seed pointer
//! @param n The number of RNG seeds to skip ahead by
//==============================================================================
extern "C" void advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream);
void advance_prn_seed(int64_t n, uint64_t* prn_seed);
//==============================================================================
//! Advance a random number seed 'n' times.
@ -68,7 +79,7 @@ extern "C" void advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream);
//! @param seed The starting to seed to advance from
//==============================================================================
uint64_t future_seed(uint64_t n, uint64_t seed);
uint64_t future_seed(uint64_t n, uint64_t prn_seed);
//==============================================================================
// API FUNCTIONS

6
include/openmc/reaction_product.h
View file

@ -42,10 +42,8 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_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;
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu, uint64_t* prn_seed) const;
Particle::Type particle_; //!< Particle type
EmissionMode emission_mode_; //!< Emission mode

17
include/openmc/scattdata.h
View file

@ -65,11 +65,9 @@ class ScattData {
//! @param gout Sampled outgoing energy group.
//! @param mu Sampled cosine of the change-in-angle.
//! @param wgt Weight of the particle to be adjusted.
//! @param prn_seeds Array of pseudorandom stream seeds
//! @param stream Pseudorandom stream index
//! @param prn_seed Pseudorandom number seed pointer
virtual void
sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds,
int stream) = 0;
sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed) = 0;
//! \brief Initializes the ScattData object from a given scatter and
//! multiplicity matrix.
@ -112,10 +110,9 @@ class ScattData {
//! @param gin Incoming energy group.
//! @param gout Sampled outgoing energy group.
//! @param i_gout Sampled outgoing energy group index.
//! @param prn_seeds Array of pseudorandom stream seeds
//! @param stream Pseudorandom stream index
//! @param prn_seed Pseudorandom number seed pointer
void
sample_energy(int gin, int& gout, int& i_gout, uint64_t * prn_seeds, int stream);
sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seed);
//! \brief Provides a cross section value given certain parameters
//!
@ -166,7 +163,7 @@ class ScattDataLegendre: public ScattData {
calc_f(int gin, int gout, double mu);
void
sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, int stream);
sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
size_t
get_order() {return dist[0][0].size() - 1;};
@ -202,7 +199,7 @@ class ScattDataHistogram: public ScattData {
calc_f(int gin, int gout, double mu);
void
sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, int stream);
sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
size_t
get_order() {return dist[0][0].size();};
@ -243,7 +240,7 @@ class ScattDataTabular: public ScattData {
calc_f(int gin, int gout, double mu);
void
sample(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds, int stream);
sample(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed);
size_t
get_order() {return dist[0][0].size();};

7
include/openmc/secondary_correlated.h
View file

@ -38,10 +38,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
// energy property
std::vector<double>& energy() { return energy_; }

7
include/openmc/secondary_kalbach.h
View file

@ -29,10 +29,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Pseudorandom stream seed array
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int sample) const override;
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
//! Outgoing energy/angle at a single incoming energy
struct KMTable {

7
include/openmc/secondary_nbody.h
View file

@ -24,10 +24,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
int n_bodies_; //!< Number of particles distributed
double mass_ratio_; //!< Total mass of particles [neutron mass]

35
include/openmc/secondary_thermal.h
View file

@ -31,10 +31,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
};
@ -54,10 +53,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom number seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
double debye_waller_;
};
@ -78,10 +76,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom number seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
const std::vector<double>& energy_; //!< Energies at which cosines are tabulated
xt::xtensor<double, 2> mu_out_; //!< Cosines for each incident energy
@ -103,10 +100,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom number seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
const std::vector<double>& energy_; //!< Incident energies
xt::xtensor<double, 2> energy_out_; //!< Outgoing energies for each incident energy
@ -129,10 +125,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom number seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
private:
//! Secondary energy/angle distribution
struct DistEnergySab {

7
include/openmc/secondary_uncorrelated.h
View file

@ -29,10 +29,9 @@ public:
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
void sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const override;
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const override;
// Accessors
AngleDistribution& angle() { return angle_; }

11
include/openmc/simulation.h
View file

@ -16,10 +16,6 @@ constexpr int STATUS_EXIT_NORMAL {0};
constexpr int STATUS_EXIT_MAX_BATCH {1};
constexpr int STATUS_EXIT_ON_TRIGGER {2};
extern Particle::Bank * shared_fission_bank;
extern int shared_fission_bank_length;
extern int shared_fission_bank_max;
//==============================================================================
// Global variable declarations
//==============================================================================
@ -42,9 +38,7 @@ extern "C" int restart_batch; //!< batch at which a restart job resumed
extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied?
extern "C" int total_gen; //!< total number of generations simulated
extern double total_weight; //!< Total source weight in a batch
extern int64_t thread_work_index;
extern int64_t work_per_rank; //!< number of particles per MPI rank
extern int64_t work_per_thread; //!< number of particles on a given tread
extern int64_t work_per_rank; //!< number of particles per MPI rank
extern const RegularMesh* entropy_mesh;
extern const RegularMesh* ufs_mesh;
@ -55,8 +49,7 @@ extern std::vector<int64_t> work_index;
// Threadprivate variables
extern "C" bool trace; //!< flag to show debug information
#pragma omp threadprivate(current_work, work_per_thread, trace)
#pragma omp threadprivate(current_work, trace)
} // namespace simulation

9
include/openmc/source.h
View file

@ -38,10 +38,9 @@ public:
explicit SourceDistribution(pugi::xml_node node);
//! Sample from the external source distribution
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
//! \param[inout] prn_seed Pseudorandom seed pointer
//! \return Sampled site
Particle::Bank sample(uint64_t * prn_seed, int stream) const;
Particle::Bank sample(uint64_t* prn_seed) const;
// Properties
double strength() const { return strength_; }
@ -62,9 +61,9 @@ extern "C" void initialize_source();
//! Sample a site from all external source distributions in proportion to their
//! source strength
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[inout] prn_seed Pseudorandom seed pointer
//! \return Sampled source site
Particle::Bank sample_external_source(uint64_t * prn_seeds);
Particle::Bank sample_external_source(uint64_t* prn_seed);
//! Fill source bank at end of generation for fixed source simulations
void fill_source_bank_fixedsource();

4
include/openmc/surface.h
View file

@ -116,8 +116,8 @@ public:
//! \return Outgoing direction of the ray
virtual Direction reflect(Position r, Direction u) const;
virtual Direction diffuse_reflect(Position r, Direction u, uint64_t * prn_seeds,
int stream) const;
virtual Direction diffuse_reflect(Position r, Direction u,
uint64_t* prn_seed) const;
//! Evaluate the equation describing the surface.
//!

10
include/openmc/thermal.h
View file

@ -64,10 +64,9 @@ public:
//! \param[in] E_in Incident neutron energy in [eV]
//! \param[out] E_out Outgoing neutron energy in [eV]
//! \param[out] mu Outgoing scattering angle cosine
//! \param[inout] prn_seeds Pseudorandom seed array
//! \param[in] stream Pseudorandom stream index
//! \param[inout] prn_seed Pseudorandom seed pointer
void sample(const NuclideMicroXS& micro_xs, double E_in,
double* E_out, double* mu, uint64_t * prn_seeds, int stream);
double* E_out, double* mu, uint64_t* prn_seed);
private:
struct Reaction {
// Default constructor
@ -102,10 +101,9 @@ public:
//! \param[out] i_temp corresponding temperature index
//! \param[out] elastic Thermal elastic scattering cross section
//! \param[out] inelastic Thermal inelastic scattering cross section
//! \param[inout] prn_seeds Array of pseudorandom seeds
//! \param[in] stream Pseudorandom stream index
//! \param[inout] prn_seed Pseudorandom seed pointer
void calculate_xs(double E, double sqrtkT, int* i_temp, double* elastic,
double* inelastic, uint64_t * prn_seeds, int stream) const;
double* inelastic, uint64_t* prn_seed) const;
//! Determine whether table applies to a particular nuclide
//!

73
openmc/lib/math.py
View file

@ -1,10 +1,12 @@
from ctypes import (c_int, c_double, POINTER, c_uint64)
from ctypes import c_int, c_double, POINTER, c_uint64
import numpy as np
from numpy.ctypeslib import ndpointer
from . import _dll
from random import randint
_dll.t_percentile.restype = c_double
_dll.t_percentile.argtypes = [c_double, c_int]
@ -26,19 +28,19 @@ _dll.calc_zn_rad.argtypes = [c_int, c_double, ndpointer(c_double)]
_dll.rotate_angle_c.restype = None
_dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
POINTER(c_double), ndpointer(c_uint64), c_int]
POINTER(c_double), POINTER(c_uint64)]
_dll.maxwell_spectrum.restype = c_double
_dll.maxwell_spectrum.argtypes = [c_double, ndpointer(c_uint64), c_int]
_dll.maxwell_spectrum.argtypes = [c_double, POINTER(c_uint64)]
_dll.watt_spectrum.restype = c_double
_dll.watt_spectrum.argtypes = [c_double, c_double, ndpointer(c_uint64), c_int]
_dll.watt_spectrum.argtypes = [c_double, c_double, POINTER(c_uint64)]
_dll.broaden_wmp_polynomials.restype = None
_dll.broaden_wmp_polynomials.argtypes = [c_double, c_double, c_int,
ndpointer(c_double)]
_dll.normal_variate.restype = c_double
_dll.normal_variate.argtypes = [c_double, c_double, ndpointer(c_uint64), c_int]
_dll.normal_variate.argtypes = [c_double, c_double, POINTER(c_uint64)]
def t_percentile(p, df):
""" Calculate the percentile of the Student's t distribution with a
@ -185,7 +187,7 @@ def calc_zn_rad(n, rho):
return zn_rad
def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
def rotate_angle(uvw0, mu, phi, prn_seed=None):
""" Rotates direction cosines through a polar angle whose cosine is
mu and through an azimuthal angle sampled uniformly.
@ -197,10 +199,8 @@ def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
Polar angle cosine to rotate
phi : float
Azimuthal angle; if None, one will be sampled uniformly
prn_seeds : iterable of int
PRNG seed array
stream : int
PRNG stream index
prn_seed : int
PRNG seed; if None, one will be generated randomly
Returns
-------
@ -209,19 +209,21 @@ def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
"""
if prn_seed is None:
prn_seed = randint(0,10000000)
uvw0_arr = np.array(uvw0, dtype=np.float64)
prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
if phi is None:
_dll.rotate_angle_c(uvw0_arr, mu, None, prn_seeds_arr, stream)
_dll.rotate_angle_c(uvw0_arr, mu, None, c_uint64(prn_seed))
else:
_dll.rotate_angle_c(uvw0_arr, mu, c_double(phi), prn_seeds_arr, stream)
_dll.rotate_angle_c(uvw0_arr, mu, c_double(phi), c_uint64(prn_seed))
uvw = uvw0_arr
return uvw
def maxwell_spectrum(T, prn_seeds, stream):
def maxwell_spectrum(T, prn_seed=None):
""" Samples an energy from the Maxwell fission distribution based
on a direct sampling scheme.
@ -229,10 +231,8 @@ def maxwell_spectrum(T, prn_seeds, stream):
----------
T : float
Spectrum parameter
prn_seeds : iterable of int
PRNG seed array
stream : int
PRNG stream index
prn_seed : int
PRNG seed; if None, one will be generated randomly
Returns
-------
@ -240,13 +240,14 @@ def maxwell_spectrum(T, prn_seeds, stream):
Sampled outgoing energy
"""
if prn_seed is None:
prn_seed = randint(0,10000000)
prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
return _dll.maxwell_spectrum(T, prn_seeds_arr, stream)
return _dll.maxwell_spectrum(T, c_uint64(prn_seed))
def watt_spectrum(a, b, prn_seeds, stream):
def watt_spectrum(a, b, prn_seed=None):
""" Samples an energy from the Watt energy-dependent fission spectrum.
Parameters
@ -255,10 +256,8 @@ def watt_spectrum(a, b, prn_seeds, stream):
Spectrum parameter a
b : float
Spectrum parameter b
prn_seeds : iterable of int
PRNG seed array
stream : int
PRNG stream index
prn_seed : int
PRNG seed; if None, one will be generated randomly
Returns
-------
@ -266,13 +265,14 @@ def watt_spectrum(a, b, prn_seeds, stream):
Sampled outgoing energy
"""
prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
if prn_seed is None:
prn_seed = randint(0,10000000)
return _dll.watt_spectrum(a, b, prn_seeds_arr, stream)
return _dll.watt_spectrum(a, b, c_uint64(prn_seed))
def normal_variate(mean_value, std_dev, prn_seeds, stream):
def normal_variate(mean_value, std_dev, prn_seed=None):
""" Samples an energy from the Normal distribution.
Parameters
@ -281,10 +281,8 @@ def normal_variate(mean_value, std_dev, prn_seeds, stream):
Mean of the Normal distribution
std_dev : float
Standard deviation of the normal distribution
prn_seeds : iterable of int
PRNG seed array
stream : int
PRNG stream index
prn_seed : int
PRNG seed; if None, one will be generated randomly
Returns
-------
@ -292,10 +290,11 @@ def normal_variate(mean_value, std_dev, prn_seeds, stream):
Sampled outgoing normally distributed value
"""
prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
if prn_seed is None:
prn_seed = randint(0,10000000)
return _dll.normal_variate(mean_value, std_dev, prn_seeds_arr, stream)
return _dll.normal_variate(mean_value, std_dev, c_uint64(prn_seed))
def broaden_wmp_polynomials(E, dopp, n):

6
src/bremsstrahlung.cpp
View file

@ -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(p.prn_seeds_, p.stream_);
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(p.prn_seeds_, p.stream_) <= 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(p.prn_seeds_, p.stream_)*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

36
src/distribution.cpp
View file

@ -35,11 +35,11 @@ Discrete::Discrete(const double* x, const double* p, int n)
normalize();
}
double Discrete::sample(uint64_t * prn_seeds, int stream) const
double Discrete::sample(uint64_t* prn_seed) const
{
int n = x_.size();
if (n > 1) {
double xi = prn(prn_seeds, stream);
double xi = prn(prn_seed);
double c = 0.0;
for (int i = 0; i < n; ++i) {
c += p_[i];
@ -74,9 +74,9 @@ Uniform::Uniform(pugi::xml_node node)
b_ = params.at(1);
}
double Uniform::sample(uint64_t * prn_seeds, int stream) const
double Uniform::sample(uint64_t* prn_seed) const
{
return a_ + prn(prn_seeds, stream)*(b_ - a_);
return a_ + prn(prn_seed)*(b_ - a_);
}
//==============================================================================
@ -88,9 +88,9 @@ Maxwell::Maxwell(pugi::xml_node node)
theta_ = std::stod(get_node_value(node, "parameters"));
}
double Maxwell::sample(uint64_t * prn_seeds, int stream) const
double Maxwell::sample(uint64_t* prn_seed) const
{
return maxwell_spectrum(theta_, prn_seeds, stream);
return maxwell_spectrum(theta_, prn_seed);
}
//==============================================================================
@ -108,15 +108,15 @@ Watt::Watt(pugi::xml_node node)
b_ = params.at(1);
}
double Watt::sample(uint64_t * prn_seeds, int stream) const
double Watt::sample(uint64_t* prn_seed) const
{
return watt_spectrum(a_, b_, prn_seeds, stream);
return watt_spectrum(a_, b_, prn_seed);
}
//==============================================================================
// Normal implementation
//==============================================================================
Normal::Normal(pugi::xml_node node)
Normal::Normal(pugi::xml_node node)
{
auto params = get_node_array<double>(node,"parameters");
if (params.size() != 2)
@ -127,15 +127,15 @@ Normal::Normal(pugi::xml_node node)
std_dev_ = params.at(1);
}
double Normal::sample(uint64_t * prn_seeds, int stream) const
double Normal::sample(uint64_t* prn_seed) const
{
return normal_variate(mean_value_, std_dev_, prn_seeds, stream);
return normal_variate(mean_value_, std_dev_, prn_seed);
}
//==============================================================================
// Muir implementation
//==============================================================================
Muir::Muir(pugi::xml_node node)
Muir::Muir(pugi::xml_node node)
{
auto params = get_node_array<double>(node,"parameters");
if (params.size() != 3)
@ -147,9 +147,9 @@ Muir::Muir(pugi::xml_node node)
kt_ = params.at(2);
}
double Muir::sample(uint64_t * prn_seeds, int stream) const
double Muir::sample(uint64_t* prn_seed) const
{
return muir_spectrum(e0_, m_rat_, kt_, prn_seeds, stream);
return muir_spectrum(e0_, m_rat_, kt_, prn_seed);
}
//==============================================================================
@ -220,10 +220,10 @@ void Tabular::init(const double* x, const double* p, std::size_t n, const double
}
}
double Tabular::sample(uint64_t * prn_seeds, int stream) const
double Tabular::sample(uint64_t* prn_seed) const
{
// Sample value of CDF
double c = prn(prn_seeds, stream);
double c = prn(prn_seed);
// Find first CDF bin which is above the sampled value
double c_i = c_[0];
@ -263,11 +263,11 @@ double Tabular::sample(uint64_t * prn_seeds, int stream) const
// Equiprobable implementation
//==============================================================================
double Equiprobable::sample(uint64_t * prn_seeds, int stream) const
double Equiprobable::sample(uint64_t* prn_seed) const
{
std::size_t n = x_.size();
double r = prn(prn_seeds, stream);
double r = prn(prn_seed);
int i = std::floor((n - 1)*r);
double xl = x_[i];

7
src/distribution_angle.cpp
View file

@ -62,7 +62,7 @@ AngleDistribution::AngleDistribution(hid_t group)
}
}
double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) const
double AngleDistribution::sample(double E, uint64_t* prn_seed) const
{
// Determine number of incoming energies
auto n = energy_.size();
@ -83,11 +83,10 @@ double AngleDistribution::sample(double E, uint64_t * prn_seeds, int stream) con
}
// Sample between the ith and (i+1)th bin
if (r > prn(prn_seeds, stream)) ++i;
if (r > prn(prn_seed)) ++i;
//assert(i >= 0 );
// Sample i-th distribution
double mu = distribution_[i]->sample(prn_seeds, stream);
double mu = distribution_[i]->sample(prn_seed);
// Make sure mu is in range [-1,1] and return
if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);

22
src/distribution_energy.cpp
View file

@ -25,7 +25,7 @@ DiscretePhoton::DiscretePhoton(hid_t group)
read_attribute(group, "atomic_weight_ratio", A_);
}
double DiscretePhoton::sample(double E, uint64_t * prn_seeds, int stream) const
double DiscretePhoton::sample(double E, uint64_t* prn_seed) const
{
if (primary_flag_ == 2) {
return energy_ + A_/(A_+ 1)*E;
@ -44,7 +44,7 @@ LevelInelastic::LevelInelastic(hid_t group)
read_attribute(group, "mass_ratio", mass_ratio_);
}
double LevelInelastic::sample(double E, uint64_t * prn_seeds, int stream ) const
double LevelInelastic::sample(double E, uint64_t* prn_seed ) const
{
return mass_ratio_*(E - threshold_);
}
@ -146,7 +146,7 @@ ContinuousTabular::ContinuousTabular(hid_t group)
} // incoming energies
}
double ContinuousTabular::sample(double E, uint64_t * prn_seeds, int stream) const
double ContinuousTabular::sample(double E, uint64_t* prn_seed) const
{
// Read number of interpolation regions and incoming energies
bool histogram_interp;
@ -177,7 +177,7 @@ double ContinuousTabular::sample(double E, uint64_t * prn_seeds, int stream) con
if (histogram_interp) {
l = i;
} else {
l = r > prn(prn_seeds, stream) ? i + 1 : i;
l = r > prn(prn_seed) ? i + 1 : i;
}
// Interpolation for energy E1 and EK
@ -197,7 +197,7 @@ double ContinuousTabular::sample(double E, uint64_t * prn_seeds, int stream) con
// Determine outgoing energy bin
n_energy_out = distribution_[l].e_out.size();
n_discrete = distribution_[l].n_discrete;
double r1 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
double c_k = distribution_[l].c[0];
int k = 0;
int end = n_energy_out - 2;
@ -275,14 +275,14 @@ MaxwellEnergy::MaxwellEnergy(hid_t group)
close_dataset(dset);
}
double MaxwellEnergy::sample(double E, uint64_t * prn_seeds, int stream) const
double MaxwellEnergy::sample(double E, uint64_t* prn_seed) const
{
// Get temperature corresponding to incoming energy
double theta = theta_(E);
while (true) {
// Sample maxwell fission spectrum
double E_out = maxwell_spectrum(theta, prn_seeds, stream);
double E_out = maxwell_spectrum(theta, prn_seed);
// Accept energy based on restriction energy
if (E_out <= E - u_) return E_out;
@ -301,7 +301,7 @@ Evaporation::Evaporation(hid_t group)
close_dataset(dset);
}
double Evaporation::sample(double E, uint64_t * prn_seeds, int stream) const
double Evaporation::sample(double E, uint64_t* prn_seed) const
{
// Get temperature corresponding to incoming energy
double theta = theta_(E);
@ -313,7 +313,7 @@ double Evaporation::sample(double E, uint64_t * prn_seeds, int stream) const
// density function
double x;
while (true) {
x = -std::log((1.0 - v*prn(prn_seeds, stream))*(1.0 - v*prn(prn_seeds, stream)));
x = -std::log((1.0 - v*prn(prn_seed))*(1.0 - v*prn(prn_seed)));
if (x <= y) break;
}
@ -338,7 +338,7 @@ WattEnergy::WattEnergy(hid_t group)
close_dataset(dset);
}
double WattEnergy::sample(double E, uint64_t * prn_seeds, int stream) const
double WattEnergy::sample(double E, uint64_t* prn_seed) const
{
// Determine Watt parameters at incident energy
double a = a_(E);
@ -346,7 +346,7 @@ double WattEnergy::sample(double E, uint64_t * prn_seeds, int stream) const
while (true) {
// Sample energy-dependent Watt fission spectrum
double E_out = watt_spectrum(a, b, prn_seeds, stream);
double E_out = watt_spectrum(a, b, prn_seed);
// Accept energy based on restriction energy
if (E_out <= E - u_) return E_out;

16
src/distribution_multi.cpp
View file

@ -53,31 +53,31 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
}
}
Direction PolarAzimuthal::sample(uint64_t * prn_seeds, int stream) const
Direction PolarAzimuthal::sample(uint64_t* prn_seed) const
{
// Sample cosine of polar angle
double mu = mu_->sample(prn_seeds, stream);
double mu = mu_->sample(prn_seed);
if (mu == 1.0) return u_ref_;
// Sample azimuthal angle
double phi = phi_->sample(prn_seeds, stream);
double phi = phi_->sample(prn_seed);
// If the reference direction is along the z-axis, rotate the aziumthal angle
// to match spherical coordinate conventions.
// TODO: apply this change directly to rotate_angle
if (u_ref_.x == 0 && u_ref_.y == 0) phi += 0.5*PI;
return rotate_angle(u_ref_, mu, &phi, prn_seeds, stream);
return rotate_angle(u_ref_, mu, &phi, prn_seed);
}
//==============================================================================
// Isotropic implementation
//==============================================================================
Direction Isotropic::sample(uint64_t * prn_seeds, int stream) const
Direction Isotropic::sample(uint64_t* prn_seed) const
{
double phi = 2.0*PI*prn(prn_seeds, stream);
double mu = 2.0*prn(prn_seeds, stream) - 1.0;
double phi = 2.0*PI*prn(prn_seed);
double mu = 2.0*prn(prn_seed) - 1.0;
return {mu, std::sqrt(1.0 - mu*mu) * std::cos(phi),
std::sqrt(1.0 - mu*mu) * std::sin(phi)};
}
@ -86,7 +86,7 @@ Direction Isotropic::sample(uint64_t * prn_seeds, int stream) const
// Monodirectional implementation
//==============================================================================
Direction Monodirectional::sample(uint64_t * prn_seeds, int stream) const
Direction Monodirectional::sample(uint64_t* prn_seed) const
{
return u_ref_;
}

18
src/distribution_spatial.cpp
View file

@ -46,9 +46,9 @@ CartesianIndependent::CartesianIndependent(pugi::xml_node node)
}
}
Position CartesianIndependent::sample(uint64_t * prn_seeds, int stream) const
Position CartesianIndependent::sample(uint64_t* prn_seed) const
{
return {x_->sample(prn_seeds, stream), y_->sample(prn_seeds, stream), z_->sample(prn_seeds, stream)};
return {x_->sample(prn_seed), y_->sample(prn_seed), z_->sample(prn_seed)};
}
//==============================================================================
@ -107,11 +107,11 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
}
Position SphericalIndependent::sample(uint64_t * prn_seeds, int stream) const
Position SphericalIndependent::sample(uint64_t* prn_seed) const
{
double r = r_->sample(prn_seeds, stream);
double theta = theta_->sample(prn_seeds, stream);
double phi = phi_->sample(prn_seeds, stream);
double r = r_->sample(prn_seed);
double theta = theta_->sample(prn_seed);
double phi = phi_->sample(prn_seed);
double x = r*sin(theta)*cos(phi) + origin_.x;
double y = r*sin(theta)*sin(phi) + origin_.y;
double z = r*cos(theta) + origin_.z;
@ -135,9 +135,9 @@ SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
upper_right_ = Position{params[3], params[4], params[5]};
}
Position SpatialBox::sample(uint64_t * prn_seeds, int stream) const
Position SpatialBox::sample(uint64_t* prn_seed) const
{
Position xi {prn(prn_seeds, stream), prn(prn_seeds, stream), prn(prn_seeds, stream)};
Position xi {prn(prn_seed), prn(prn_seed), prn(prn_seed)};
return lower_left_ + xi*(upper_right_ - lower_left_);
}
@ -157,7 +157,7 @@ SpatialPoint::SpatialPoint(pugi::xml_node node)
r_ = Position{params.data()};
}
Position SpatialPoint::sample(uint64_t * prn_seeds, int stream) const
Position SpatialPoint::sample(uint64_t* prn_seed) const
{
return r_;
}

13
src/eigenvalue.cpp
View file

@ -114,17 +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
// Create pseudorandom number seed
uint64_t prn_seed;
// 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.
//std::cout << "MPI rank " << mpi::rank << " start = " << start << std::endl;
set_particle_seed(simulation::total_gen + overall_generation(), prn_seeds);
advance_prn_seed(start, prn_seeds, stream);
int64_t id = simulation::total_gen + overall_generation();
init_seed(id, &prn_seed, STREAM_TRACKING);
advance_prn_seed(start, &prn_seed);
// Determine how many fission sites we need to sample from the source bank
// and the probability for selecting a site.
@ -159,7 +158,7 @@ void synchronize_bank()
}
// Randomly sample sites needed
if (prn(prn_seeds, stream) < p_sample) {
if (prn(&prn_seed) < p_sample) {
temp_sites[index_temp] = site;
++index_temp;
}

34
src/math_functions.cpp
View file

@ -630,22 +630,22 @@ void calc_zn_rad(int n, double rho, double zn_rad[]) {
}
void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t * prn_seeds, int stream) {
Direction u = rotate_angle({uvw}, mu, phi, prn_seeds, stream);
void rotate_angle_c(double uvw[3], double mu, const double* phi, uint64_t* prn_seed) {
Direction u = rotate_angle({uvw}, mu, phi, prn_seed);
uvw[0] = u.x;
uvw[1] = u.y;
uvw[2] = u.z;
}
Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t * prn_seeds, int stream)
Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* prn_seed)
{
// Sample azimuthal angle in [0,2pi) if none provided
double phi_;
if (phi != nullptr) {
phi_ = (*phi);
} else {
phi_ = 2.0*PI*prn(prn_seeds, stream);
phi_ = 2.0*PI*prn(prn_seed);
}
// Precompute factors to save flops
@ -675,11 +675,11 @@ Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t * prn
}
double maxwell_spectrum(double T, uint64_t * prn_seeds, int stream) {
double maxwell_spectrum(double T, uint64_t* prn_seed) {
// Set the random numbers
double r1 = prn(prn_seeds, stream);
double r2 = prn(prn_seeds, stream);
double r3 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
double r2 = prn(prn_seed);
double r3 = prn(prn_seed);
// determine cosine of pi/2*r
double c = std::cos(PI / 2. * r3);
@ -691,33 +691,33 @@ double maxwell_spectrum(double T, uint64_t * prn_seeds, int stream) {
}
double normal_variate(double mean, double standard_deviation, uint64_t * prn_seeds, int stream) {
double normal_variate(double mean, double standard_deviation, uint64_t* prn_seed) {
// perhaps there should be a limit to the number of resamples
while ( true ) {
double v1 = 2 * prn(prn_seeds, stream) - 1.;
double v2 = 2 * prn(prn_seeds, stream) - 1.;
double v1 = 2 * prn(prn_seed) - 1.;
double v2 = 2 * prn(prn_seed) - 1.;
double r = std::pow(v1, 2) + std::pow(v2, 2);
double r2 = std::pow(r, 2);
if (r2 < 1) {
double z = std::sqrt(-2.0 * std::log(r2)/r2);
z *= (prn(prn_seeds, stream) <= 0.5) ? v1 : v2;
z *= (prn(prn_seed) <= 0.5) ? v1 : v2;
return mean + standard_deviation*z;
}
}
}
double muir_spectrum(double e0, double m_rat, double kt, uint64_t * prn_seeds, int stream) {
double muir_spectrum(double e0, double m_rat, double kt, uint64_t* prn_seed) {
// note sigma here is a factor of 2 shy of equation
// 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
double sigma = std::sqrt(2.*e0*kt/m_rat);
return normal_variate(e0, sigma, prn_seeds, stream);
return normal_variate(e0, sigma, prn_seed);
}
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);
double watt_spectrum(double a, double b, uint64_t* prn_seed) {
double w = maxwell_spectrum(a, prn_seed);
double E_out = w + 0.25 * a * a * b + (2. * prn(prn_seed) - 1.) * std::sqrt(a * a * b * w);
return E_out;
}

11
src/mgxs.cpp
View file

@ -529,7 +529,7 @@ Mgxs::get_xs(int xstype, int gin, const int* gout, const double* mu,
//==============================================================================
void
Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t * prn_seeds, int stream)
Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t* prn_seed)
{
// This method assumes that the temperature and angle indices are set
#ifdef _OPENMP
@ -544,8 +544,8 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t * prn_seeds, i
double prob_prompt = xs_t->prompt_nu_fission(cache[tid].a, gin);
// sample random numbers
double xi_pd = prn(prn_seeds, stream) * nu_fission;
double xi_gout = prn(prn_seeds, stream);
double xi_pd = prn(prn_seed) * nu_fission;
double xi_gout = prn(prn_seed);
// Select whether the neutron is prompt or delayed
if (xi_pd <= prob_prompt) {
@ -585,8 +585,7 @@ Mgxs::sample_fission_energy(int gin, int& dg, int& gout, uint64_t * prn_seeds, i
//==============================================================================
void
Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t * prn_seeds,
int stream)
Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t* prn_seed)
{
// This method assumes that the temperature and angle indices are set
// Sample the data
@ -595,7 +594,7 @@ Mgxs::sample_scatter(int gin, int& gout, double& mu, double& wgt, uint64_t * prn
#else
int tid = 0;
#endif
xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, prn_seeds, stream);
xs[cache[tid].t].scatter[cache[tid].a]->sample(gin, gout, mu, wgt, prn_seed);
}
//==============================================================================

6
src/nuclide.cpp
View file

@ -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(p.prn_seeds_, p.stream_)) ++i_temp;
if (f > prn(p.prn_seeds_ + p.stream_)) ++i_temp;
break;
}
@ -720,7 +720,7 @@ void Nuclide::calculate_sab_xs(int i_sab, double sab_frac, Particle& p)
int i_temp;
double elastic;
double inelastic;
data::thermal_scatt[i_sab]->calculate_xs(p.E_, p.sqrtkT_, &i_temp, &elastic, &inelastic, p.prn_seeds_, p.stream_);
data::thermal_scatt[i_sab]->calculate_xs(p.E_, p.sqrtkT_, &i_temp, &elastic, &inelastic, p.prn_seeds_ + p.stream_);
// Store the S(a,b) cross sections.
micro.thermal = sab_frac * (elastic + inelastic);
@ -759,7 +759,7 @@ void Nuclide::calculate_urr_xs(int i_temp, Particle& p) const
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<int64_t>(i_nuclide_ + 1), p.prn_seeds_, p.stream_);
double r = future_prn(static_cast<int64_t>(i_nuclide_ + 1), p.prn_seeds_[p.stream_]);
p.stream_ = STREAM_TRACKING;
int i_low = 0;

10
src/output.cpp
View file

@ -67,15 +67,7 @@ void title()
" ############### %%%%%%%%%%%%%%%%\n" <<
" ############ %%%%%%%%%%%%%%%\n" <<
" ######## %%%%%%%%%%%%%%\n" <<
" %%%%%%%%%%%\n" <<
" ) ( (\n" <<
" ( /( * ) ( ( )\\ ) )\\ )\n" <<
" ( ( ( ( )\\())` ) /( ( )\\ )\\ (()/( ( (()/(\n" <<
" )\\ )\\ )\\ )\\ ((_)\\ ( )(_))___ )((_)((((_)( /(_)))\\ /(_))\n" <<
" ((_) ((_)((_)((_) _((_)(_(_())|___|((_)_ )\\ _ )\\ (_)) ((_)(_))_\n" <<
" | __|\\ \\ / / | __|| \\| ||_ _| | _ ) (_)_\\(_)/ __|| __|| \\\n" <<
" | _| \\ V / | _| | .` | | | | _ \\ / _ \\ \\__ \\| _| | |) |\n" <<
" |___| \\_/ |___||_|\\_| |_| |___/ /_/ \\_\\ |___/|___||___/\n\n";
" %%%%%%%%%%%\n\n";
// Write version information
std::cout <<

18
src/particle.cpp
View file

@ -27,14 +27,8 @@
#include "openmc/tallies/tally_scoring.h"
#include "openmc/track_output.h"
// Explicit template instantiation definition
template class std::vector<openmc::Particle>;
namespace openmc {
std::vector<Particle> particle_bank;
//==============================================================================
// LocalCoord implementation
//==============================================================================
@ -140,15 +134,14 @@ void
Particle::transport()
{
// Display message if high verbosity or trace is on
//if (settings::verbosity >= 9 || simulation::trace) {
// write_message("Simulating Particle " + std::to_string(id_));
//}
if (settings::verbosity >= 9 || simulation::trace) {
write_message("Simulating Particle " + std::to_string(id_));
}
// Initialize number of events to zero
int n_event = 0;
// Add paricle's starting weight to count for normalizing tallies later
/*
#pragma omp atomic
simulation::total_weight += wgt_;
@ -162,7 +155,6 @@ Particle::transport()
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty()) zero_flux_derivs();
*/
while (true) {
// Set the random number stream
@ -238,7 +230,7 @@ Particle::transport()
} else if (macro_xs_.total == 0.0) {
d_collision = INFINITY;
} else {
d_collision = -std::log(prn(prn_seeds_, stream_)) / macro_xs_.total;
d_collision = -std::log(prn(prn_seeds_ + stream_)) / macro_xs_.total;
}
// Select smaller of the two distances
@ -474,7 +466,7 @@ Particle::cross_surface()
Direction u = (surf->bc_ == BC_REFLECT) ?
surf->reflect(this->r(), this->u()) :
surf->diffuse_reflect(this->r(), this->u(), prn_seeds_, stream_);
surf->diffuse_reflect(this->r(), this->u(), prn_seeds_ + stream_);
// Make sure new particle direction is normalized
this->u() = u / u.norm();

2
src/particle_restart.cpp
View file

@ -98,7 +98,7 @@ void run_particle_restart()
throw std::runtime_error{"Unexpected run mode: " +
std::to_string(previous_run_mode)};
}
set_particle_seed(particle_seed, p.prn_seeds_);
init_particle_seeds(particle_seed, p.prn_seeds_);
// Transport neutron
p.transport();

70
src/photon.cpp
View file

@ -290,12 +290,12 @@ PhotonInteraction::PhotonInteraction(hid_t group, int i_element)
}
void PhotonInteraction::compton_scatter(double alpha, bool doppler,
double* alpha_out, double* mu, int* i_shell, uint64_t * prn_seeds, int stream) const
double* alpha_out, double* mu, int* i_shell, uint64_t* prn_seed) const
{
double form_factor_xmax = 0.0;
while (true) {
// Sample Klein-Nishina distribution for trial energy and angle
std::tie(*alpha_out, *mu) = klein_nishina(alpha, prn_seeds, stream);
std::tie(*alpha_out, *mu) = klein_nishina(alpha, prn_seed);
// Note that the parameter used here does not correspond exactly to the
// momentum transfer q in ENDF-102 Eq. (27.2). Rather, this is the
@ -309,10 +309,10 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler,
}
// Perform rejection on form factor
if (prn(prn_seeds, stream) < form_factor_x / form_factor_xmax) {
if (prn(prn_seed) < form_factor_x / form_factor_xmax) {
if (doppler) {
double E_out;
this->compton_doppler(alpha, *mu, &E_out, i_shell, prn_seeds, stream);
this->compton_doppler(alpha, *mu, &E_out, i_shell, prn_seed);
*alpha_out = E_out/MASS_ELECTRON_EV;
} else {
*i_shell = -1;
@ -323,14 +323,14 @@ void PhotonInteraction::compton_scatter(double alpha, bool doppler,
}
void PhotonInteraction::compton_doppler(double alpha, double mu,
double* E_out, int* i_shell, uint64_t * prn_seeds, int stream) const
double* E_out, int* i_shell, uint64_t* prn_seed) const
{
auto n = data::compton_profile_pz.size();
int shell; // index for shell
while (true) {
// Sample electron shell
double rn = prn(prn_seeds, stream);
double rn = prn(prn_seed);
double c = 0.0;
for (shell = 0; shell < electron_pdf_.size(); ++shell) {
c += electron_pdf_(shell);
@ -377,7 +377,7 @@ void PhotonInteraction::compton_doppler(double alpha, double mu,
}
// Sample value on bounded cdf
c = prn(prn_seeds, stream)*c_max;
c = prn(prn_seed)*c_max;
// Determine pz corresponding to sampled cdf value
auto cdf_shell = xt::view(profile_cdf_, shell, xt::all());
@ -418,7 +418,7 @@ void PhotonInteraction::compton_doppler(double alpha, double mu,
if (E_out1 > 0.0) {
if (E_out2 > 0.0) {
// If both are positive, pick one at random
*E_out = prn(prn_seeds, stream) < 0.5 ? E_out1 : E_out2;
*E_out = prn(prn_seed) < 0.5 ? E_out1 : E_out2;
} else {
*E_out = E_out1;
}
@ -496,7 +496,7 @@ void PhotonInteraction::calculate_xs(Particle& p) const
xs.last_E = p.E_;
}
double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t * prn_seeds, int stream) const
double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t* prn_seed) const
{
double mu;
while (true) {
@ -507,7 +507,7 @@ double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t * prn_seeds, i
double F_max = coherent_int_form_factor_(x2_max);
// Sample cumulative distribution
double F = prn(prn_seeds, stream)*F_max;
double F = prn(prn_seed)*F_max;
// Determine x^2 corresponding to F
const auto& x {coherent_int_form_factor_.x()};
@ -519,14 +519,14 @@ double PhotonInteraction::rayleigh_scatter(double alpha, uint64_t * prn_seeds, i
// Calculate mu
mu = 1.0 - 2.0*x2/x2_max;
if (prn(prn_seeds, stream) < 0.5*(1.0 + mu*mu)) break;
if (prn(prn_seed) < 0.5*(1.0 + mu*mu)) break;
}
return mu;
}
void PhotonInteraction::pair_production(double alpha, double* E_electron,
double* E_positron, double* mu_electron, double* mu_positron, uint64_t * prn_seeds,
int stream) const
double* E_positron, double* mu_electron, double* mu_positron,
uint64_t* prn_seed) const
{
constexpr double r[] {
122.81, 73.167, 69.228, 67.301, 64.696, 61.228,
@ -593,12 +593,12 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
double u2 = phi2_max;
double e;
while (true) {
double rn = prn(prn_seeds, stream);
double rn = prn(prn_seed);
// Sample the index i in (1, 2) using the point probabilities
// p(1) = u_1/(u_1 + u_2) and p(2) = u_2/(u_1 + u_2)
int i;
if (prn(prn_seeds, stream) < u1/(u1 + u2)) {
if (prn(prn_seed) < u1/(u1 + u2)) {
i = 1;
// Sample e from pi_1 using the inverse transform method
@ -619,10 +619,10 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
t3 = b*b*(4.0 - 4.0*t2 - 3.0*std::log(1.0 + 1.0/(b*b)));
if (i == 1) {
double phi1 = 7.0/3.0 - t1 - 6.0*t2 - t3 + t4;
if (prn(prn_seeds, stream) <= phi1/phi1_max) break;
if (prn(prn_seed) <= phi1/phi1_max) break;
} else {
double phi2 = 11.0/6.0 - t1 - 3.0*t2 + 0.5*t3 + t4;
if (prn(prn_seeds, stream) <= phi2/phi2_max) break;
if (prn(prn_seed) <= phi2/phi2_max) break;
}
}
@ -635,13 +635,13 @@ void PhotonInteraction::pair_production(double alpha, double* E_electron,
// p(mu) = C/(1 - beta*mu)^2 using the inverse transform method.
double beta = std::sqrt(*E_electron*(*E_electron + 2.0*MASS_ELECTRON_EV))
/ (*E_electron + MASS_ELECTRON_EV) ;
double rn = 2.0*prn(prn_seeds, stream) - 1.0;
double rn = 2.0*prn(prn_seed) - 1.0;
*mu_electron = (rn + beta)/(rn*beta + 1.0);
// Sample the scattering angle of the positron
beta = std::sqrt(*E_positron*(*E_positron + 2.0*MASS_ELECTRON_EV))
/ (*E_positron + MASS_ELECTRON_EV);
rn = 2.0*prn(prn_seeds, stream) - 1.0;
rn = 2.0*prn(prn_seed) - 1.0;
*mu_positron = (rn + beta)/(rn*beta + 1.0);
}
@ -649,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(p.prn_seeds_, p.stream_) - 1.0;
double phi = 2.0*PI*prn(p.prn_seeds_, p.stream_);
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);
@ -661,7 +661,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
}
// Sample transition
double rn = prn(p.prn_seeds_, p.stream_);
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) {
@ -674,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(p.prn_seeds_, p.stream_) - 1.0;
double phi = 2.0*PI*prn(p.prn_seeds_, p.stream_);
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);
@ -711,7 +711,7 @@ void PhotonInteraction::atomic_relaxation(const ElectronSubshell& shell, Particl
// Non-member functions
//==============================================================================
std::pair<double, double> klein_nishina(double alpha, uint64_t * prn_seeds, int stream)
std::pair<double, double> klein_nishina(double alpha, uint64_t* prn_seed)
{
double alpha_out, mu;
double beta = 1.0 + 2.0*alpha;
@ -720,19 +720,19 @@ std::pair<double, double> klein_nishina(double alpha, uint64_t * prn_seeds, int
double t = beta/(beta + 8.0);
double x;
while (true) {
if (prn(prn_seeds, stream) < t) {
if (prn(prn_seed) < t) {
// Left branch of flow chart
double r = 2.0*prn(prn_seeds, stream);
double r = 2.0*prn(prn_seed);
x = 1.0 + alpha*r;
if (prn(prn_seeds, stream) < 4.0/x*(1.0 - 1.0/x)) {
if (prn(prn_seed) < 4.0/x*(1.0 - 1.0/x)) {
mu = 1 - r;
break;
}
} else {
// Right branch of flow chart
x = beta/(1.0 + 2.0*alpha*prn(prn_seeds, stream));
x = beta/(1.0 + 2.0*alpha*prn(prn_seed));
mu = 1.0 + (1.0 - x)/alpha;
if (prn(prn_seeds, stream) < 0.5*(mu*mu + 1.0/x)) break;
if (prn(prn_seed) < 0.5*(mu*mu + 1.0/x)) break;
}
}
alpha_out = alpha/x;
@ -740,24 +740,24 @@ std::pair<double, double> klein_nishina(double alpha, uint64_t * prn_seeds, int
} else {
// Koblinger's direct method
double gamma = 1.0 - std::pow(beta, -2);
double s = prn(prn_seeds, stream)*(4.0/alpha + 0.5*gamma +
double s = prn(prn_seed)*(4.0/alpha + 0.5*gamma +
(1.0 - (1.0 + beta)/(alpha*alpha))*std::log(beta));
if (s <= 2.0/alpha) {
// For first term, x = 1 + 2ar
// Therefore, a' = a/(1 + 2ar)
alpha_out = alpha/(1.0 + 2.0*alpha*prn(prn_seeds, stream));
alpha_out = alpha/(1.0 + 2.0*alpha*prn(prn_seed));
} else if (s <= 4.0/alpha) {
// For third term, x = beta/(1 + 2ar)
// Therefore, a' = a(1 + 2ar)/beta
alpha_out = alpha*(1.0 + 2.0*alpha*prn(prn_seeds, stream))/beta;
alpha_out = alpha*(1.0 + 2.0*alpha*prn(prn_seed))/beta;
} else if (s <= 4.0/alpha + 0.5*gamma) {
// For fourth term, x = 1/sqrt(1 - gamma*r)
// Therefore, a' = a*sqrt(1 - gamma*r)
alpha_out = alpha*std::sqrt(1.0 - gamma*prn(prn_seeds, stream));
alpha_out = alpha*std::sqrt(1.0 - gamma*prn(prn_seed));
} else {
// For third term, x = beta^r
// Therefore, a' = a/beta^r
alpha_out = alpha/std::pow(beta, prn(prn_seeds, stream));
alpha_out = alpha/std::pow(beta, prn(prn_seed));
}
// Calculate cosine of scattering angle based on basic relation

153
src/physics.cpp
View file

@ -79,7 +79,6 @@ void collision(Particle* p)
void sample_neutron_reaction(Particle* p)
{
//std::cout << "particle energy e = " << p->E_ << std::endl;
// Sample a nuclide within the material
int i_nuclide = sample_nuclide(p);
@ -130,18 +129,15 @@ void sample_neutron_reaction(Particle* p)
}
if (!p->alive_) return;
//std::cout << "before scatter particle energy e = " << p->E_ << " with nuclide = " << i_nuclide << std::endl;
// Sample a scattering reaction and determine the secondary energy of the
// exiting neutron
scatter(p, i_nuclide);
//std::cout << "after scatter particle energy e = " << p->E_ << std::endl;
// Advance URR seed stream 'N' times after energy changes
if (p->E_ != p->E_last_) {
p->stream_ = STREAM_URR_PTABLE;
advance_prn_seed(data::nuclides.size(), p->prn_seeds_, p->stream_);
p->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
@ -165,7 +161,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn(p->prn_seeds_, p->stream_) <= (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
@ -177,16 +173,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
p->fission_ = true;
for (int i = 0; i < nu; ++i) {
/*
// Create new bank site and get reference to last element
bank.emplace_back();
auto& site {bank.back()};
// Bank source neutrons by copying the particle data
site.r = p->r();
site.particle = Particle::Type::neutron;
site.wgt = 1. / weight;
*/
int idx;
#pragma omp atomic capture
idx = shared_fission_bank_length++;
@ -197,12 +184,10 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
site->parent_id = p->id_;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, rx, p->E_, site, p->prn_seeds_, p->stream_);
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;
p->delayed_group_ = site->delayed_group;
p->delayed_group_ = site.delayed_group;
// Increment the number of neutrons born delayed
if (p->delayed_group_ > 0) {
@ -241,13 +226,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(p->prn_seeds_, p->stream_) * 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->prn_seeds_, p->stream_);
p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_, p->stream_);
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;
@ -258,7 +243,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, p->prn_seeds_, p->stream_);
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.
@ -270,13 +255,13 @@ void sample_photon_reaction(Particle* p)
}
// Create Compton electron
double phi = 2.0*PI*prn(p->prn_seeds_, p->stream_);
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<int>(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, p->prn_seeds_, p->stream_);
Direction u = rotate_angle(p->u(), mu_electron, &phi, p->prn_seeds_ + p->stream_);
p->create_secondary(u, E_electron, Particle::Type::electron);
}
@ -290,7 +275,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->prn_seeds_, p->stream_);
p->u() = rotate_angle(p->u(), mu, &phi, p->prn_seeds_ + p->stream_);
p->event_ = EVENT_SCATTER;
p->event_mt_ = INCOHERENT;
return;
@ -322,8 +307,8 @@ void sample_photon_reaction(Particle* p)
// model in Serpent 2" by Toni Kaltiaisenaho
double mu;
while (true) {
double r = prn(p->prn_seeds_, p->stream_);
if (4.0*(1.0 - r)*r >= prn(p->prn_seeds_, p->stream_)) {
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);
@ -331,7 +316,7 @@ void sample_photon_reaction(Particle* p)
}
}
double phi = 2.0*PI*prn(p->prn_seeds_, p->stream_);
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);
@ -359,14 +344,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, p->prn_seeds_, p->stream_);
&mu_electron, &mu_positron, p->prn_seeds_ + p->stream_);
// Create secondary electron
Direction u = rotate_angle(p->u(), mu_electron, nullptr, p->prn_seeds_, p->stream_);
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, p->prn_seeds_, p->stream_);
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;
@ -400,8 +385,8 @@ void sample_positron_reaction(Particle* p)
}
// Sample angle isotropically
double mu = 2.0*prn(p->prn_seeds_, p->stream_) - 1.0;
double phi = 2.0*PI*prn(p->prn_seeds_, p->stream_);
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);
@ -419,7 +404,7 @@ void sample_positron_reaction(Particle* p)
int sample_nuclide(Particle* p)
{
// Sample cumulative distribution function
double cutoff = prn(p->prn_seeds_, p->stream_) * 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_]};
@ -444,7 +429,7 @@ int sample_nuclide(Particle* p)
int sample_element(Particle* p)
{
// Sample cumulative distribution function
double cutoff = prn(p->prn_seeds_, p->stream_) * 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_]};
@ -497,7 +482,7 @@ Reaction* sample_fission(int i_nuclide, 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->prn_seeds_, p->stream_) * 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
@ -524,7 +509,7 @@ void sample_photon_product(int i_nuclide, Particle* p, int* i_rx, int* i_product
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->prn_seeds_, p->stream_) * 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
@ -572,7 +557,7 @@ void absorption(Particle* p, int i_nuclide)
} else {
// See if disappearance reaction happens
if (p->neutron_xs_[i_nuclide].absorption >
prn(p->prn_seeds_, p->stream_) * 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_[
@ -600,7 +585,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(p->prn_seeds_, p->stream_) * (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
@ -616,7 +601,6 @@ void scatter(Particle* p, int i_nuclide)
// Determine temperature
double kT = nuc->multipole_ ? p->sqrtkT_*p->sqrtkT_ : nuc->kTs_[i_temp];
//std::cout << "we are doing an elastic scatter" << std::endl;
// Perform collision physics for elastic scattering
elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
@ -629,7 +613,6 @@ void scatter(Particle* p, int i_nuclide)
// =======================================================================
// S(A,B) SCATTERING
//std::cout << "we are doing an SAB scatter" << std::endl;
sab_scatter(i_nuclide, micro.index_sab, p);
p->event_mt_ = ELASTIC;
@ -637,7 +620,6 @@ void scatter(Particle* p, int i_nuclide)
}
if (!sampled) {
//std::cout << "we are doing an Inelastic scatter" << std::endl;
// =======================================================================
// INELASTIC SCATTERING
@ -664,9 +646,7 @@ void scatter(Particle* p, int i_nuclide)
// Perform collision physics for inelastic scattering
const auto& rx {nuc->reactions_[i]};
//std::cout << "Energy before Inelastic scatter = " << p->E_ << std::endl;
inelastic_scatter(nuc.get(), rx.get(), p);
//std::cout << "Energy after Inelastic scatter = " << p->E_ << std::endl;
p->event_mt_ = rx->mt_;
}
@ -679,8 +659,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(p->prn_seeds_, p->stream_) - 1.0;
double phi = 2.0*PI*prn(p->prn_seeds_, p->stream_);
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;
@ -707,7 +687,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->prn_seeds_, p->stream_);
p->neutron_xs_[i_nuclide].elastic, kT, p->prn_seeds_ + p->stream_);
}
// Velocity of center-of-mass
@ -725,9 +705,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
auto& d = rx.products_[0].distribution_[0];
auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
if (d_) {
mu_cm = d_->angle().sample(p->E_, p->prn_seeds_, p->stream_);
mu_cm = d_->angle().sample(p->E_, p->prn_seeds_ + p->stream_);
} else {
mu_cm = 2.0*prn(p->prn_seeds_, p->stream_) - 1.0;
mu_cm = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
}
// Determine direction cosines in CM
@ -736,7 +716,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, p->prn_seeds_, p->stream_);
v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->prn_seeds_ + p->stream_);
// Transform back to LAB frame
v_n += v_cm;
@ -765,15 +745,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_, p->prn_seeds_, p->stream_);
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->prn_seeds_, p->stream_);
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, uint64_t * prn_seeds, int stream)
Direction v_neut, double xs_eff, double kT, uint64_t* prn_seed)
{
// check if nuclide is a resonant scatterer
ResScatMethod sampling_method;
@ -805,7 +785,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
case ResScatMethod::cxs:
// sample target velocity with the constant cross section (cxs) approx.
return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream);
return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed);
case ResScatMethod::dbrc:
case ResScatMethod::rvs: {
@ -839,7 +819,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
if (i_E_up == i_E_low) {
// Handle degenerate case -- if the upper/lower bounds occur for the same
// index, then using cxs is probably a good approximation
return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream);
return sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed);
}
if (sampling_method == ResScatMethod::dbrc) {
@ -863,7 +843,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
Direction v_target;
while (true) {
// sample target velocity with the constant cross section (cxs) approx.
v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seeds, stream);
v_target = sample_cxs_target_velocity(nuc->awr_, E, u, kT, prn_seed);
Direction v_rel = v_neut - v_target;
E_rel = v_rel.dot(v_rel);
if (E_rel < E_up) break;
@ -872,7 +852,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
// perform Doppler broadening rejection correction (dbrc)
double xs_0K = nuc->elastic_xs_0K(E_rel);
double R = xs_0K / xs_max;
if (prn(prn_seeds, stream) < R) return v_target;
if (prn(prn_seed) < R) return v_target;
}
} else if (sampling_method == ResScatMethod::rvs) {
@ -892,10 +872,10 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
while (true) {
// directly sample Maxwellian
double E_t = -kT * std::log(prn(prn_seeds, stream));
double E_t = -kT * std::log(prn(prn_seed));
// sample a relative energy using the xs cdf
double cdf_rel = cdf_low + prn(prn_seeds, stream)*(cdf_up - cdf_low);
double cdf_rel = cdf_low + prn(prn_seed)*(cdf_up - cdf_low);
int i_E_rel = lower_bound_index(&nuc->xs_cdf_[i_E_low-1],
&nuc->xs_cdf_[i_E_up+1], cdf_rel);
double E_rel = nuc->energy_0K_[i_E_low + i_E_rel];
@ -913,7 +893,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
if (std::abs(mu) < 1.0) {
// set and accept target velocity
E_t /= nuc->awr_;
return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, prn_seeds, stream);
return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, prn_seed);
}
}
}
@ -924,7 +904,7 @@ Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
}
Direction
sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t * prn_seeds, int stream)
sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_t* prn_seed)
{
double beta_vn = std::sqrt(awr * E / kT);
double alpha = 1.0/(1.0 + std::sqrt(PI)*beta_vn/2.0);
@ -933,10 +913,10 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
double mu;
while (true) {
// Sample two random numbers
double r1 = prn(prn_seeds, stream);
double r2 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
double r2 = prn(prn_seed);
if (prn(prn_seeds, stream) < alpha) {
if (prn(prn_seed) < alpha) {
// With probability alpha, we sample the distribution p(y) =
// y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
// Carlo sampler
@ -948,7 +928,7 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
// e^(-y^2). This can be done with sampling scheme C61 from the Monte
// Carlo sampler
double c = std::cos(PI/2.0 * prn(prn_seeds, stream));
double c = std::cos(PI/2.0 * prn(prn_seed));
beta_vt_sq = -std::log(r1) - std::log(r2)*c*c;
}
@ -956,14 +936,14 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
double beta_vt = std::sqrt(beta_vt_sq);
// Sample cosine of angle between neutron and target velocity
mu = 2.0*prn(prn_seeds, stream) - 1.0;
mu = 2.0*prn(prn_seed) - 1.0;
// Determine rejection probability
double accept_prob = std::sqrt(beta_vn*beta_vn + beta_vt_sq -
2*beta_vn*beta_vt*mu) / (beta_vn + beta_vt);
// Perform rejection sampling on vt and mu
if (prn(prn_seeds, stream) < accept_prob) break;
if (prn(prn_seed) < accept_prob) break;
}
// Determine speed of target nucleus
@ -971,19 +951,19 @@ sample_cxs_target_velocity(double awr, double E, Direction u, double kT, uint64_
// Determine velocity vector of target nucleus based on neutron's velocity
// and the sampled angle between them
return vt * rotate_angle(u, mu, nullptr, prn_seeds, stream);
return vt * rotate_angle(u, mu, nullptr, prn_seed);
}
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t * prn_seeds, int stream)
void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Particle::Bank* site, uint64_t* prn_seed)
{
// Sample cosine of angle -- fission neutrons are always emitted
// isotropically. Sometimes in ACE data, fission reactions actually have
// an angular distribution listed, but for those that do, it's simply just
// a uniform distribution in mu
double mu = 2.0 * prn(prn_seeds, stream) - 1.0;
double mu = 2.0 * prn(prn_seed) - 1.0;
// Sample azimuthal angle uniformly in [0,2*pi)
double phi = 2.0*PI*prn(prn_seeds, stream);
double phi = 2.0*PI*prn(prn_seed);
site->u.x = mu;
site->u.y = std::sqrt(1.0 - mu*mu) * std::cos(phi);
site->u.z = std::sqrt(1.0 - mu*mu) * std::sin(phi);
@ -994,12 +974,12 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed);
double beta = nu_d / nu_t;
if (prn(prn_seeds, stream) < beta) {
if (prn(prn_seed) < beta) {
// ====================================================================
// DELAYED NEUTRON SAMPLED
// sampled delayed precursor group
double xi = prn(prn_seeds, stream)*nu_d;
double xi = prn(prn_seed)*nu_d;
double prob = 0.0;
int group;
for (group = 1; group < nuc->n_precursor_; ++group) {
@ -1023,7 +1003,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
while (true) {
// sample from energy/angle distribution -- note that mu has already been
// sampled above and doesn't need to be resampled
rx->products_[group].sample(E_in, site->E, mu, prn_seeds, stream);
rx->products_[group].sample(E_in, site->E, mu, prn_seed);
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(Particle::Type::neutron);
@ -1048,7 +1028,7 @@ void sample_fission_neutron(int i_nuclide, const Reaction* rx, double E_in, Part
// sample from prompt neutron energy distribution
int n_sample = 0;
while (true) {
rx->products_[0].sample(E_in, site->E, mu, prn_seeds, stream);
rx->products_[0].sample(E_in, site->E, mu, prn_seed);
// resample if energy is greater than maximum neutron energy
constexpr int neutron = static_cast<int>(Particle::Type::neutron);
@ -1073,40 +1053,33 @@ 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, p->prn_seeds_, p->stream_);
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
if (rx->scatter_in_cm_) {
double E_cm = E;
//std::cout << "E_cm = " << E_cm << std::endl;
// determine outgoing energy in lab
double A = nuc->awr_;
// std::cout << "A = " << nuc->awr_ << std::endl;
E = E_cm + (E_in + 2.0*mu*(A + 1.0) * std::sqrt(E_in*E_cm))
/ ((A + 1.0)*(A + 1.0));
// here's where it goes wrong
//std::cout << "E = " << E << std::endl;
// determine outgoing angle in lab
mu = mu*std::sqrt(E_cm/E) + 1.0/(A+1.0) * std::sqrt(E_in/E);
// std::cout << "mu = " << mu << std::endl;
}
// Because of floating-point roundoff, it may be possible for mu to be
// outside of the range [-1,1). In these cases, we just set mu to exactly -1
// or 1
if (std::abs(mu) > 1.0) mu = std::copysign(1.0, mu);
// std::cout << "mu = " << mu << std::endl;
// Set outgoing energy and scattering angle
// std::cout << "FINAL E = " << E << " mu = " << mu << std::endl;
p->E_ = E;
p->mu_ = mu;
// change direction of particle
p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_, p->stream_);
p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_ + p->stream_);
// evaluate yield
double yield = (*rx->products_[0].yield_)(E_in);
@ -1127,7 +1100,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<int>(y_t);
if (prn(p->prn_seeds_, p->stream_) <= 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) {
@ -1140,10 +1113,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, p->prn_seeds_, p->stream_);
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, p->prn_seeds_, p->stream_);
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);

2
src/physics_common.cpp
View file

@ -12,7 +12,7 @@ namespace openmc {
void russian_roulette(Particle* p)
{
if (p->wgt_ < settings::weight_cutoff) {
if (prn(p->prn_seeds_, p->stream_) < 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 {

16
src/physics_mg.cpp
View file

@ -79,10 +79,10 @@ void
scatter(Particle* p)
{
data::mg.macro_xs_[p->material_].sample_scatter(p->g_last_, p->g_, p->mu_,
p->wgt_, p->prn_seeds_, p->stream_);
p->wgt_, p->prn_seeds_ + p->stream_);
// Rotate the angle
p->u() = rotate_angle(p->u(), p->mu_, nullptr, p->prn_seeds_, p->stream_);
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_];
@ -104,7 +104,7 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn(p->prn_seeds_, p->stream_) <= (nu_t - int(nu_t))) {
if (prn(p->prn_seeds_ + p->stream_) <= (nu_t - int(nu_t))) {
nu++;
}
@ -129,10 +129,10 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
// Sample the cosine of the angle, assuming fission neutrons are emitted
// isotropically
double mu = 2.*prn(p->prn_seeds_, p->stream_) - 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(p->prn_seeds_, p->stream_ );
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);
@ -140,8 +140,8 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& 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, p->prn_seeds_,
p->stream_);
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
@ -181,7 +181,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->prn_seeds_, p->stream_) * 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;

9
src/plot.cpp
View file

@ -80,8 +80,7 @@ namespace model {
std::vector<Plot> plots;
std::unordered_map<int, int> plot_map;
uint64_t plotter_prn_seeds[N_STREAMS] = {1, 2, 3, 4, 5, 6};
int plotter_stream = STREAM_TRACKING;
uint64_t plotter_prn_seed = 1;
} // namespace model
@ -961,9 +960,9 @@ voxel_finalize(hid_t dspace, hid_t dset, hid_t memspace)
}
RGBColor random_color(void) {
return {int(prn(model::plotter_prn_seeds, model::plotter_stream)*255),
int(prn(model::plotter_prn_seeds, model::plotter_stream)*255),
int(prn(model::plotter_prn_seeds, model::plotter_stream)*255)};
return {int(prn(&model::plotter_prn_seed)*255),
int(prn(&model::plotter_prn_seed)*255),
int(prn(&model::plotter_prn_seed)*255)};
}
extern "C" int openmc_id_map(const void* plot, int32_t* data_out)

48
src/random_lcg.cpp
View file

@ -6,15 +6,6 @@
namespace openmc {
// Constants
extern "C" const int STREAM_TRACKING {0};
extern "C" const int STREAM_TALLIES {1};
extern "C" const int STREAM_SOURCE {2};
extern "C" const int STREAM_URR_PTABLE {3};
extern "C" const int STREAM_VOLUME {4};
extern "C" const int STREAM_PHOTON {5};
// Starting seed
int64_t master_seed {1};
@ -32,34 +23,40 @@ constexpr double prn_norm {1.0 / prn_mod}; // 2^-63
// PRN
//==============================================================================
extern "C" double
prn(uint64_t * prn_seeds, int stream)
double prn(uint64_t* prn_seed)
{
// This algorithm uses bit-masking to find the next integer(8) value to be
// used to calculate the random number.
prn_seeds[stream] = (prn_mult*prn_seeds[stream] + prn_add) & prn_mask;
*prn_seed = (prn_mult * (*prn_seed) + prn_add) & prn_mask;
// Once the integer is calculated, we just need to divide by 2**m,
// represented here as multiplying by a pre-calculated factor
return prn_seeds[stream] * prn_norm;
return (*prn_seed) * prn_norm;
}
//==============================================================================
// FUTURE_PRN
//==============================================================================
extern "C" double
future_prn(int64_t n, uint64_t * prn_seeds, int stream)
double future_prn(int64_t n, uint64_t prn_seed)
{
return future_seed(static_cast<uint64_t>(n), prn_seeds[stream]) * prn_norm;
return future_seed(static_cast<uint64_t>(n), prn_seed) * prn_norm;
}
//==============================================================================
// SET_PARTICLE_SEED
// INIT_SEED
//==============================================================================
extern "C" void
set_particle_seed(int64_t id, uint64_t * prn_seeds)
void init_seed(int64_t id, uint64_t* prn_seed, int offset)
{
*prn_seed = future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
}
//==============================================================================
// INIT_PARTICLE_SEEDS
//==============================================================================
void init_particle_seeds(int64_t id, uint64_t* prn_seeds)
{
//std::cout << "Master seed = " << master_seed << " id = " << id << std::endl;
for (int i = 0; i < N_STREAMS; i++) {
@ -71,18 +68,16 @@ set_particle_seed(int64_t id, uint64_t * prn_seeds)
// ADVANCE_PRN_SEED
//==============================================================================
extern "C" void
advance_prn_seed(int64_t n, uint64_t * prn_seeds, int stream)
void advance_prn_seed(int64_t n, uint64_t* prn_seed)
{
prn_seeds[stream] = future_seed(static_cast<uint64_t>(n), prn_seeds[stream]);
*prn_seed = future_seed(static_cast<uint64_t>(n), *prn_seed);
}
//==============================================================================
// FUTURE_SEED
//==============================================================================
uint64_t
future_seed(uint64_t n, uint64_t seed)
uint64_t future_seed(uint64_t n, uint64_t prn_seed)
{
// Make sure nskip is less than 2^M.
n &= prn_mask;
@ -113,7 +108,7 @@ future_seed(uint64_t n, uint64_t seed)
}
// With G and C, we can now find the new seed.
return (g_new * seed + c_new) & prn_mask;
return (g_new * prn_seed + c_new) & prn_mask;
}
//==============================================================================
@ -122,8 +117,7 @@ future_seed(uint64_t n, uint64_t seed)
extern "C" int64_t openmc_get_seed() {return master_seed;}
extern "C" void
openmc_set_seed(int64_t new_seed)
extern "C" void openmc_set_seed(int64_t new_seed)
{
master_seed = new_seed;
}

8
src/reaction_product.cpp
View file

@ -77,25 +77,25 @@ ReactionProduct::ReactionProduct(hid_t group)
}
void ReactionProduct::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
auto n = applicability_.size();
if (n > 1) {
double prob = 0.0;
double c = prn(prn_seeds, stream);
double c = prn(prn_seed);
for (int i = 0; i < n; ++i) {
// Determine probability that i-th energy distribution is sampled
prob += applicability_[i](E_in);
// If i-th distribution is sampled, sample energy from the distribution
if (c <= prob) {
distribution_[i]->sample(E_in, E_out, mu, prn_seeds, stream);
distribution_[i]->sample(E_in, E_out, mu, prn_seed);
break;
}
}
} else {
// If only one distribution is present, go ahead and sample it
distribution_[0]->sample(E_in, E_out, mu, prn_seeds, stream);
distribution_[0]->sample(E_in, E_out, mu, prn_seed);
}
}

26
src/scattdata.cpp
View file

@ -167,10 +167,10 @@ ScattData::base_combine(size_t max_order,
//==============================================================================
void
ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t * prn_seeds, int stream)
ScattData::sample_energy(int gin, int& gout, int& i_gout, uint64_t* prn_seed)
{
// Sample the outgoing group
double xi = prn(prn_seeds, stream);
double xi = prn(prn_seed);
double prob = 0.;
i_gout = 0;
for (gout = gmin[gin]; gout < gmax[gin]; ++gout) {
@ -348,21 +348,21 @@ ScattDataLegendre::calc_f(int gin, int gout, double mu)
void
ScattDataLegendre::sample(int gin, int& gout, double& mu, double& wgt,
uint64_t * prn_seeds, int stream)
uint64_t* prn_seed)
{
// Sample the outgoing energy using the base-class method
int i_gout;
sample_energy(gin, gout, i_gout, prn_seeds, stream);
sample_energy(gin, gout, i_gout, prn_seed);
// Now we can sample mu using the scattering kernel using rejection
// sampling from a rectangular bounding box
double M = max_val[gin][i_gout];
int samples;
for (samples = 0; samples < MAX_SAMPLE; ++samples) {
mu = 2. * prn(prn_seeds, stream) - 1.;
mu = 2. * prn(prn_seed) - 1.;
double f = calc_f(gin, gout, mu);
if (f > 0.) {
double u = prn(prn_seeds, stream) * M;
double u = prn(prn_seed) * M;
if (u <= f) break;
}
}
@ -537,14 +537,14 @@ ScattDataHistogram::calc_f(int gin, int gout, double mu)
void
ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt,
uint64_t * prn_seeds, int stream)
uint64_t* prn_seed)
{
// Sample the outgoing energy using the base-class method
int i_gout;
sample_energy(gin, gout, i_gout, prn_seeds, stream);
sample_energy(gin, gout, i_gout, prn_seed);
// Determine the outgoing cosine bin
double xi = prn(prn_seeds, stream);
double xi = prn(prn_seed);
int imu;
if (xi < dist[gin][i_gout][0]) {
@ -556,7 +556,7 @@ ScattDataHistogram::sample(int gin, int& gout, double& mu, double& wgt,
}
// Randomly select mu within the imu bin
mu = prn(prn_seeds, stream) * dmu + this->mu[imu];
mu = prn(prn_seed) * dmu + this->mu[imu];
if (mu < -1.) {
mu = -1.;
@ -741,15 +741,15 @@ ScattDataTabular::calc_f(int gin, int gout, double mu)
void
ScattDataTabular::sample(int gin, int& gout, double& mu, double& wgt,
uint64_t * prn_seeds, int stream)
uint64_t* prn_seed)
{
// Sample the outgoing energy using the base-class method
int i_gout;
sample_energy(gin, gout, i_gout, prn_seeds, stream);
sample_energy(gin, gout, i_gout, prn_seed);
// Determine the outgoing cosine bin
int NP = this->mu.shape()[0];
double xi = prn(prn_seeds, stream);
double xi = prn(prn_seed);
double c_k = dist[gin][i_gout][0];
int k;

12
src/secondary_correlated.cpp
View file

@ -153,14 +153,14 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
}
void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
// Before the secondary distribution refactor, an isotropic polar cosine was
// always sampled but then overwritten with the polar cosine sampled from the
// correlated distribution. To preserve the random number stream, we keep
// this dummy sampling here but can remove it later (will change answers)
mu = 2.0*prn(prn_seeds, stream) - 1.0;
mu = 2.0*prn(prn_seed) - 1.0;
// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
// Find energy bin and calculate interpolation factor -- if the energy is
@ -180,7 +180,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
}
// Sample between the ith and [i+1]th bin
int l = r > prn(prn_seeds, stream) ? i + 1 : i;
int l = r > prn(prn_seed) ? i + 1 : i;
// Interpolation for energy E1 and EK
int n_energy_out = distribution_[i].e_out.size();
@ -199,7 +199,7 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
// Determine outgoing energy bin
n_energy_out = distribution_[l].e_out.size();
n_discrete = distribution_[l].n_discrete;
double r1 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
double c_k = distribution_[l].c[0];
int k = 0;
int end = n_energy_out - 2;
@ -260,9 +260,9 @@ void CorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
// Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1) {
mu = distribution_[l].angle[k]->sample(prn_seeds, stream);
mu = distribution_[l].angle[k]->sample(prn_seed);
} else {
mu = distribution_[l].angle[k + 1]->sample(prn_seeds, stream);
mu = distribution_[l].angle[k + 1]->sample(prn_seed);
}
}

14
src/secondary_kalbach.cpp
View file

@ -113,14 +113,14 @@ KalbachMann::KalbachMann(hid_t group)
} // incoming energies
}
void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds, int stream) const
void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t* prn_seed) const
{
// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
// Before the secondary distribution refactor, an isotropic polar cosine was
// always sampled but then overwritten with the polar cosine sampled from the
// correlated distribution. To preserve the random number stream, we keep
// this dummy sampling here but can remove it later (will change answers)
mu = 2.0*prn(prn_seeds, stream) - 1.0;
mu = 2.0*prn(prn_seed) - 1.0;
// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
// Find energy bin and calculate interpolation factor -- if the energy is
@ -140,7 +140,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t * prn_
}
// Sample between the ith and [i+1]th bin
int l = r > prn(prn_seeds, stream) ? i + 1 : i;
int l = r > prn(prn_seed) ? i + 1 : i;
// Interpolation for energy E1 and EK
int n_energy_out = distribution_[i].e_out.size();
@ -159,7 +159,7 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t * prn_
// Determine outgoing energy bin
n_energy_out = distribution_[l].e_out.size();
n_discrete = distribution_[l].n_discrete;
double r1 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
double c_k = distribution_[l].c[0];
int k = 0;
int end = n_energy_out - 2;
@ -229,11 +229,11 @@ void KalbachMann::sample(double E_in, double& E_out, double& mu, uint64_t * prn_
}
// Sampled correlated angle from Kalbach-Mann parameters
if (prn(prn_seeds, stream) > km_r) {
double T = (2.0*prn(prn_seeds, stream) - 1.0) * std::sinh(km_a);
if (prn(prn_seed) > km_r) {
double T = (2.0*prn(prn_seed) - 1.0) * std::sinh(km_a);
mu = std::log(T + std::sqrt(T*T + 1.0))/km_a;
} else {
double r1 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
mu = std::log(r1*std::exp(km_a) + (1.0 - r1)*std::exp(-km_a))/km_a;
}
}

26
src/secondary_nbody.cpp
View file

@ -22,38 +22,38 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
}
void NBodyPhaseSpace::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// By definition, the distribution of the angle is isotropic for an N-body
// phase space distribution
mu = 2.0*prn(prn_seeds, stream) - 1.0;
mu = 2.0*prn(prn_seed) - 1.0;
// Determine E_max parameter
double Ap = mass_ratio_;
double E_max = (Ap - 1.0)/Ap * (A_/(A_ + 1.0)*E_in + Q_);
// x is essentially a Maxwellian distribution
double x = maxwell_spectrum(1.0, prn_seeds, stream);
double x = maxwell_spectrum(1.0, prn_seed);
double y;
double r1, r2, r3, r4, r5, r6;
switch (n_bodies_) {
case 3:
y = maxwell_spectrum(1.0, prn_seeds, stream);
y = maxwell_spectrum(1.0, prn_seed);
break;
case 4:
r1 = prn(prn_seeds, stream);
r2 = prn(prn_seeds, stream);
r3 = prn(prn_seeds, stream);
r1 = prn(prn_seed);
r2 = prn(prn_seed);
r3 = prn(prn_seed);
y = -std::log(r1*r2*r3);
break;
case 5:
r1 = prn(prn_seeds, stream);
r2 = prn(prn_seeds, stream);
r3 = prn(prn_seeds, stream);
r4 = prn(prn_seeds, stream);
r5 = prn(prn_seeds, stream);
r6 = prn(prn_seeds, stream);
r1 = prn(prn_seed);
r2 = prn(prn_seed);
r3 = prn(prn_seed);
r4 = prn(prn_seed);
r5 = prn(prn_seed);
r6 = prn(prn_seed);
y = -std::log(r1*r2*r3*r4) - std::log(r5) * std::pow(std::cos(PI/2.0*r6), 2);
break;
default:

32
src/secondary_thermal.cpp
View file

@ -32,8 +32,8 @@ CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs)
{ }
void
CoherentElasticAE::sample(double E_in, double& E_out, double& mu, uint64_t * prn_seeds,
int stream) const
CoherentElasticAE::sample(double E_in, double& E_out, double& mu,
uint64_t* prn_seed) const
{
// Get index and interpolation factor for elastic grid
int i;
@ -43,7 +43,7 @@ CoherentElasticAE::sample(double E_in, double& E_out, double& mu, uint64_t * prn
// Sample a Bragg edge between 1 and i
const auto& factors = xs_.factors();
double prob = prn(prn_seeds, stream) * factors[i+1];
double prob = prn(prn_seed) * factors[i+1];
int k = 0;
if (prob >= factors.front()) {
k = lower_bound_index(factors.begin(), factors.begin() + (i+1), prob);
@ -67,11 +67,11 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group)
void
IncoherentElasticAE::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
double c = 2 * E_in * debye_waller_;
mu = std::log(1.0 + prn(prn_seeds, stream)*(std::exp(2.0*c) - 1))/c - 1.0;
mu = std::log(1.0 + prn(prn_seed)*(std::exp(2.0*c) - 1))/c - 1.0;
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
E_out = E_in;
@ -90,7 +90,7 @@ IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(hid_t group,
void
IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// Get index and interpolation factor for elastic grid
int i;
@ -101,7 +101,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
// incoming energies.
// Sample outgoing cosine bin
int k = prn(prn_seeds, stream) * mu_out_.shape()[1];
int k = prn(prn_seed) * mu_out_.shape()[1];
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ik = mu_out_(i, k);
@ -129,7 +129,7 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(hid_t group,
void
IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// Get index and interpolation factor for inelastic grid
int i;
@ -147,10 +147,10 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
int n = energy_out_.shape()[1];
if (!skewed_) {
// All bins equally likely
j = prn(prn_seeds, stream) * n;
j = prn(prn_seed) * n;
} else {
// Distribution skewed away from edge points
double r = prn(prn_seeds, stream) * (n - 3);
double r = prn(prn_seed) * (n - 3);
if (r > 1.0) {
// equally likely N-4 middle bins
j = r + 1;
@ -178,7 +178,7 @@ IncoherentInelasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
// Sample outgoing cosine bin
int m = mu_out_.shape()[2];
int k = prn(prn_seeds, stream) * m;
int k = prn(prn_seed) * m;
// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
double mu_ijk = mu_out_(i, j, k);
@ -233,7 +233,7 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
void
IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// Get index and interpolation factor for inelastic grid
int i;
@ -241,7 +241,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
get_energy_index(energy_, E_in, i, f);
// Sample between ith and [i+1]th bin
int l = f > prn(prn_seeds, stream) ? i + 1 : i;
int l = f > prn(prn_seed) ? i + 1 : i;
// Determine endpoints on grid i
auto n = distribution_[i].e_out.size();
@ -259,7 +259,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
// Determine outgoing energy bin
// (First reset n_energy_out to the right value)
n = distribution_[l].n_e_out;
double r1 = prn(prn_seeds, stream);
double r1 = prn(prn_seed);
double c_j = distribution_[l].e_out_cdf[0];
double c_j1;
std::size_t j;
@ -298,7 +298,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
// Sample outgoing cosine bin
int n_mu = distribution_[l].mu.shape()[1];
std::size_t k = prn(prn_seeds, stream) * n_mu;
std::size_t k = prn(prn_seed) * n_mu;
// Rather than use the sampled discrete mu directly, it is smeared over
// a bin of width min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
@ -326,7 +326,7 @@ IncoherentInelasticAE::sample(double E_in, double& E_out, double& mu,
}
// Smear angle
mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seeds, stream) - 0.5);
mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5);
}
} // namespace openmc

8
src/secondary_uncorrelated.cpp
View file

@ -53,7 +53,7 @@ UncorrelatedAngleEnergy::UncorrelatedAngleEnergy(hid_t group)
void
UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// Sample cosine of scattering angle
if (fission_) {
@ -62,14 +62,14 @@ UncorrelatedAngleEnergy::sample(double E_in, double& E_out, double& mu,
mu = 1.0;
// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<< REMOVE THIS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<
} else if (!angle_.empty()) {
mu = angle_.sample(E_in, prn_seeds, stream);
mu = angle_.sample(E_in, prn_seed);
} else {
// no angle distribution given => assume isotropic for all energies
mu = 2.0*prn(prn_seeds, stream) - 1.0;
mu = 2.0*prn(prn_seed) - 1.0;
}
// Sample outgoing energy
E_out = energy_->sample(E_in, prn_seeds, stream);
E_out = energy_->sample(E_in, prn_seed);
}
} // namespace openmc

106
src/simulation.cpp
View file

@ -2,32 +2,23 @@
#include "openmc/bank.h"
#include "openmc/capi.h"
#include "openmc/cell.h"
#include "openmc/container_util.h"
#include "openmc/eigenvalue.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/material.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/output.h"
#include "openmc/particle.h"
#include "openmc/photon.h"
#include "openmc/physics.h"
#include "openmc/physics_mg.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/source.h"
#include "openmc/state_point.h"
#include "openmc/thermal.h"
#include "openmc/timer.h"
#include "openmc/tallies/derivative.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/tallies/trigger.h"
#include "openmc/track_output.h"
#ifdef _OPENMP
#include <omp.h>
@ -42,6 +33,7 @@
#include <string>
#include <chrono>
namespace openmc {
/*
@ -791,12 +783,10 @@ int openmc_run()
int err = 0;
int status = 0;
while (status == 0 && err == 0) {
err = openmc_next_batch(&status);
}
openmc_simulation_finalize();
return err;
}
@ -807,16 +797,9 @@ int openmc_simulation_init()
// Skip if simulation has already been initialized
if (simulation::initialized) return 0;
// Determine how much work each process should do
calculate_work();
// We assume that the number of fission sites will be at most
// 10x the number of neutrons run on this rank
int shared_fission_bank_size = 10 * simulation::work_per_rank;
// Initializes the shared fission bank
init_shared_fission_bank(shared_fission_bank_size);
// Allocate array for matching filter bins
#pragma omp parallel
@ -829,11 +812,10 @@ int openmc_simulation_init()
allocate_banks();
// Allocate tally results arrays if they're not allocated yet
//std::cout << "Trying to Initialize Results..." << std::endl;
for (auto& t : model::tallies) {
t->init_results();
}
//std::cout << "Success in Initializing Results!" << std::endl;
// Set up material nuclide index mapping
for (auto& mat : model::materials) {
@ -878,8 +860,6 @@ int openmc_simulation_finalize()
// Skip if simulation was never run
if (!simulation::initialized) return 0;
free_shared_fission_bank();
// Stop active batch timer and start finalization timer
simulation::time_active.stop();
@ -947,29 +927,9 @@ int openmc_next_batch(int* status)
// Start timer for transport
simulation::time_transport.start();
// ====================================================================
// LOOP OVER PARTICLES
simulation::current_work = 1;
simulation::current_work = 1;
//#pragma omp parallel
{
transport();
}
/*
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
simulation::current_work = i_work;
// grab source particle from bank
Particle p;
initialize_history(&p, simulation::current_work);
// transport particle
p.transport();
}
*/
transport();
// Accumulate time for transport
simulation::time_transport.stop();
@ -1026,9 +986,7 @@ int restart_batch;
bool satisfy_triggers {false};
int total_gen {0};
double total_weight;
int64_t thread_work_index;
int64_t work_per_rank;
int64_t work_per_thread;
const RegularMesh* entropy_mesh {nullptr};
const RegularMesh* ufs_mesh {nullptr};
@ -1252,30 +1210,14 @@ void finalize_generation()
}
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
//#ifdef _OPENMP
// Join the fission bank from each thread into one global fission bank
//join_bank_from_threads();
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
// We need to move all the stuff from the shared_fission_bank into the real one.
//std::vector<Particle::Bank> shared_fission_bank_vector(shared_fission_bank, shared_fission_bank + shared_fission_bank_length);
//simulation::fission_bank = shared_fission_bank_vector;
//std::cout << "Fission bank length = " << shared_fission_bank_length << std::endl;
for( int i = 0; i < shared_fission_bank_length; i++ )
simulation::fission_bank.push_back(shared_fission_bank[i]);
shared_fission_bank_length = 0;
//#endif
// Sorts the fission bank so as to allow for reproducibility
//std::sort(simulation::fission_bank.begin(), simulation::fission_bank.end(), bank_site_comparator());
//std::sort(simulation::fission_bank.begin(), simulation::fission_bank.end());
std::stable_sort(simulation::fission_bank.begin(), simulation::fission_bank.end());
/*
std::cout << "Fission bank on rank " << mpi::rank << " is of length " << simulation::fission_bank.size() << std::endl;
for (Particle::Bank p : simulation::fission_bank )
{
std::cout << "E = " << p.E << std::endl;
}
*/
// Sorts the fission bank so as to allow for reproducibility
std::stable_sort(simulation::fission_bank.begin(), simulation::fission_bank.end());
// Distribute fission bank across processors evenly
synchronize_bank();
@ -1302,10 +1244,8 @@ void finalize_generation()
void initialize_history(Particle* p, int64_t index_source)
{
//assert(index_source - 1 >= 0 );
// set defaults
p->from_source(&simulation::source_bank[index_source - 1]);
//p->from_source(&simulation::source_bank[index_source]);
// set identifier for particle
p->id_ = simulation::work_index[mpi::rank] + index_source;
@ -1313,8 +1253,7 @@ void initialize_history(Particle* p, int64_t index_source)
// set random number seed
int64_t particle_seed = (simulation::total_gen + overall_generation() - 1)
* settings::n_particles + p->id_;
//std::cout << "MPI rank " << mpi::rank << " is initializing particle " << index_source << " with ID/seed " << particle_seed << std::endl;
set_particle_seed(particle_seed, p->prn_seeds_);
init_particle_seeds(particle_seed, p->prn_seeds_);
// set particle trace
simulation::trace = false;
@ -1337,6 +1276,7 @@ void initialize_history(Particle* p, int64_t index_source)
}
}
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || simulation::trace) {
write_message("Simulating Particle " + std::to_string(p->id_));
@ -1396,32 +1336,6 @@ void calculate_work()
i_bank += work_i;
simulation::work_index[i + 1] = i_bank;
}
#ifdef _OPENMP
// Determine work per thread
int remainder_thread = simulation::work_per_rank % omp_get_max_threads();
#pragma omp parallel
{
simulation::work_per_thread = simulation::work_per_rank / omp_get_num_threads();
if (omp_get_thread_num() < remainder_thread) {
++simulation::work_per_thread;
}
}
int64_t work_i = 0;
#pragma omp parallel for ordered
for (int i = 0; i < omp_get_num_threads(); ++i) {
#pragma omp ordered
{
simulation::thread_work_index = work_i;
work_i += simulation::work_per_thread;
/*
std::cout << "Thread " << omp_get_thread_num() << ": " << simulation::work_per_thread
<< std::endl;
*/
}
}
#endif
}
#ifdef OPENMC_MPI

29
src/source.cpp
View file

@ -142,7 +142,7 @@ SourceDistribution::SourceDistribution(pugi::xml_node node)
}
Particle::Bank SourceDistribution::sample(uint64_t * prn_seeds, int stream) const
Particle::Bank SourceDistribution::sample(uint64_t* prn_seed) const
{
Particle::Bank site;
@ -158,7 +158,7 @@ Particle::Bank SourceDistribution::sample(uint64_t * prn_seeds, int stream) cons
site.particle = particle_;
// Sample spatial distribution
site.r = space_->sample(prn_seeds, stream);
site.r = space_->sample(prn_seed);
double xyz[] {site.r.x, site.r.y, site.r.z};
// Now search to see if location exists in geometry
@ -200,7 +200,7 @@ Particle::Bank SourceDistribution::sample(uint64_t * prn_seeds, int stream) cons
++n_accept;
// Sample angle
site.u = angle_->sample(prn_seeds, stream);
site.u = angle_->sample(prn_seed);
// Check for monoenergetic source above maximum particle energy
auto p = static_cast<int>(particle_);
@ -218,7 +218,7 @@ Particle::Bank SourceDistribution::sample(uint64_t * prn_seeds, int stream) cons
while (true) {
// Sample energy spectrum
site.E = energy_->sample(prn_seeds, stream);
site.E = energy_->sample(prn_seed);
// Resample if energy falls outside minimum or maximum particle energy
if (site.E < data::energy_max[p] && site.E > data::energy_min[p]) break;
@ -270,11 +270,11 @@ void initialize_source()
// initialize random number seed
int64_t id = simulation::total_gen*settings::n_particles +
simulation::work_index[mpi::rank] + i + 1;
uint64_t prn_seeds[N_STREAMS];
set_particle_seed(id, prn_seeds);
uint64_t prn_seed;
init_seed(id, &prn_seed, STREAM_SOURCE);
// sample external source distribution
simulation::source_bank[i] = sample_external_source(prn_seeds);
simulation::source_bank[i] = sample_external_source(&prn_seed);
}
}
@ -288,11 +288,8 @@ void initialize_source()
}
}
Particle::Bank sample_external_source(uint64_t * prn_seeds)
Particle::Bank sample_external_source(uint64_t* prn_seed)
{
// Set the random number generator to the source stream.
int stream = STREAM_SOURCE;
// Determine total source strength
double total_strength = 0.0;
for (auto& s : model::external_sources)
@ -301,7 +298,7 @@ Particle::Bank sample_external_source(uint64_t * prn_seeds)
// Sample from among multiple source distributions
int i = 0;
if (model::external_sources.size() > 1) {
double xi = prn(prn_seeds, stream)*total_strength;
double xi = prn(prn_seed)*total_strength;
double c = 0.0;
for (; i < model::external_sources.size(); ++i) {
c += model::external_sources[i].strength();
@ -310,7 +307,7 @@ Particle::Bank sample_external_source(uint64_t * prn_seeds)
}
// Sample source site from i-th source distribution
Particle::Bank site {model::external_sources[i].sample(prn_seeds, stream)};
Particle::Bank site {model::external_sources[i].sample(prn_seed)};
// If running in MG, convert site.E to group
if (!settings::run_CE) {
@ -334,11 +331,11 @@ void fill_source_bank_fixedsource()
// initialize random number seed
int64_t id = (simulation::total_gen + overall_generation()) *
settings::n_particles + simulation::work_index[mpi::rank] + i + 1;
uint64_t prn_seeds[N_STREAMS];
set_particle_seed(id, prn_seeds);
uint64_t prn_seed;
init_seed(id, &prn_seed, STREAM_SOURCE);
// sample external source distribution
simulation::source_bank[i] = sample_external_source(prn_seeds);
simulation::source_bank[i] = sample_external_source(&prn_seed);
}
}
}

26
src/surface.cpp
View file

@ -197,24 +197,24 @@ Surface::reflect(Position r, Direction u) const
}
Direction
Surface::diffuse_reflect(Position r, Direction u, uint64_t * prn_seeds, int stream) const
Surface::diffuse_reflect(Position r, Direction u, uint64_t* prn_seed) const
{
// Diffuse reflect direction according to the normal.
// cosine distribution
// cosine distribution
Direction n = this->normal(r);
n /= n.norm();
const double projection = n.dot(u);
// sample from inverse function, u=sqrt(rand) since p(u)=2u, so F(u)=u^2
const double mu = (projection>=0.0) ?
-std::sqrt(prn(prn_seeds, stream)) : std::sqrt(prn(prn_seeds, stream));
// sample azimuthal distribution uniformly
u = rotate_angle(n, mu, nullptr, prn_seeds, stream);
// normalize the direction
return u/u.norm();
// sample from inverse function, u=sqrt(rand) since p(u)=2u, so F(u)=u^2
const double mu = (projection>=0.0) ?
-std::sqrt(prn(prn_seed)) : std::sqrt(prn(prn_seed));
// sample azimuthal distribution uniformly
u = rotate_angle(n, mu, nullptr, prn_seed);
// normalize the direction
return u/u.norm();
}
CSGSurface::CSGSurface() : Surface{} {};

1
src/tallies/tally.cpp
View file

@ -804,7 +804,6 @@ Tally::init_triggers(pugi::xml_node node)
void Tally::init_results()
{
int n_scores = scores_.size() * nuclides_.size();
//std::cout << "TRYING TO ALLOCATE n filter bins, n scores = " << n_filter_bins_ << " " << n_scores << std::endl;
results_ = xt::empty<double>({n_filter_bins_, n_scores, 3});
}

12
src/thermal.cpp
View file

@ -151,7 +151,7 @@ ThermalScattering::ThermalScattering(hid_t group, const std::vector<double>& tem
void
ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
double* elastic, double* inelastic,
uint64_t * prn_seeds, int stream) const
uint64_t* prn_seed) const
{
// Determine temperature for S(a,b) table
double kT = sqrtkT*sqrtkT;
@ -173,7 +173,7 @@ ThermalScattering::calculate_xs(double E, double sqrtkT, int* i_temp,
// Randomly sample between temperature i and i+1
double f = (kT - kTs_[i]) / (kTs_[i+1] - kTs_[i]);
if (f > prn(prn_seeds, stream)) ++i;
if (f > prn(prn_seed)) ++i;
}
// Set temperature index
@ -266,13 +266,13 @@ ThermalData::calculate_xs(double E, double* elastic, double* inelastic) const
void
ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
double* E_out, double* mu, uint64_t * prn_seeds, int stream)
double* E_out, double* mu, uint64_t* prn_seed)
{
// Determine whether inelastic or elastic scattering will occur
if (prn(prn_seeds, stream) < micro_xs.thermal_elastic / micro_xs.thermal) {
elastic_.distribution->sample(E, *E_out, *mu, prn_seeds, stream);
if (prn(prn_seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
elastic_.distribution->sample(E, *E_out, *mu, prn_seed);
} else {
inelastic_.distribution->sample(E, *E_out, *mu, prn_seeds, stream);
inelastic_.distribution->sample(E, *E_out, *mu, prn_seed);
}
// Because of floating-point roundoff, it may be possible for mu to be

9
src/volume_calc.cpp
View file

@ -125,16 +125,15 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
std::vector<std::vector<int>> hits(n);
Particle p;
uint64_t prn_seeds[N_STREAMS];
int stream = STREAM_VOLUME;
// Sample locations and count hits
#pragma omp for
for (size_t i = i_start; i < i_end; i++) {
set_particle_seed(iterations * n_samples_ + i, prn_seeds);
uint64_t prn_seed;
int64_t id = iterations * n_samples_ + i;
init_seed(id, &prn_seed, STREAM_VOLUME);
p.n_coord_ = 1;
Position xi {prn(prn_seeds, stream), prn(prn_seeds, stream), prn(prn_seeds, stream)};
Position xi {prn(&prn_seed), prn(&prn_seed), prn(&prn_seed)};
p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
p.u() = {0.5, 0.5, 0.5};

29
tests/unit_tests/test_math.py
View file

@ -156,13 +156,11 @@ def test_rotate_angle():
uvw0 = np.array([1., 0., 0.])
phi = 0.
mu = 0.
prn_seeds = [1, 2, 3, 4, 5, 6]
stream = 0
# reference: mu of 0 pulls the vector the bottom, so:
ref_uvw = np.array([0., 0., -1.])
test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream)
test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi)
assert np.array_equal(ref_uvw, test_uvw)
@ -170,59 +168,56 @@ def test_rotate_angle():
mu = 1.
ref_uvw = np.array([1., 0., 0.])
test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream)
test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi)
assert np.array_equal(ref_uvw, test_uvw)
# Now to test phi is None
mu = 0.9
phi = None
prn_seed = 1
# When seed = 1, phi will be sampled as 1.9116495709698769
# The resultant reference is from hand-calculations given the above
ref_uvw = [0.9, 0.410813051297112, 0.1457142302040]
test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seeds, stream)
test_uvw = openmc.lib.math.rotate_angle(uvw0, mu, phi, prn_seed)
assert np.allclose(ref_uvw, test_uvw)
def test_maxwell_spectrum():
prn_seeds = [1, 2, 3, 4, 5, 6]
stream = 0
prn_seed = 1
T = 0.5
ref_val = 0.6129982175261098
test_val = openmc.lib.math.maxwell_spectrum(T, prn_seeds, stream)
test_val = openmc.lib.math.maxwell_spectrum(T, prn_seed)
assert ref_val == test_val
def test_watt_spectrum():
prn_seeds = [1, 2, 3, 4, 5, 6]
stream = 0
prn_seed = 1
a = 0.5
b = 0.75
ref_val = 0.6247242713640233
test_val = openmc.lib.math.watt_spectrum(a, b, prn_seeds, stream)
test_val = openmc.lib.math.watt_spectrum(a, b, prn_seed)
assert ref_val == test_val
def test_normal_dist():
prn_seeds = [1, 2, 3, 4, 5, 6]
stream = 0
prn_seed = 1
a = 14.08
b = 0.0
ref_val = 14.08
test_val = openmc.lib.math.normal_variate(a, b, prn_seeds, stream)
test_val = openmc.lib.math.normal_variate(a, b, prn_seed)
assert ref_val == pytest.approx(test_val)
prn_seeds = [1, 2, 3, 4, 5, 6]
stream = 0
prn_seed = 1
a = 14.08
b = 1.0
ref_val = 16.436645416691427
test_val = openmc.lib.math.normal_variate(a, b, prn_seeds, stream)
test_val = openmc.lib.math.normal_variate(a, b, prn_seed)
assert ref_val == pytest.approx(test_val)