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

incorporated code review changes

Browse source
This commit is contained in:
John Tramm 2019-12-05 19:32:12 +00:00
parent 08ce3f4614
commit 428bab4fd2
14 changed files with 82 additions and 79 deletions
Download patch file Download diff file Expand all files Collapse all files

4
include/openmc/particle.h
View file

@ -218,6 +218,10 @@ public:
//! create a particle restart HDF5 file
void write_restart() const;
//! Gets the pointer to the particle's current PRN seed
uint64_t* current_seed();
const uint64_t* current_seed() const;
//==========================================================================
// Data members

8
include/openmc/random_lcg.h
View file

@ -43,11 +43,13 @@ double future_prn(int64_t n, uint64_t prn_seed);
//==============================================================================
//! 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
//! @param offset The offset from the master seed to be used (e.g., for creating
//! different streams)
//! @return The initialized seed value
//==============================================================================
void init_seed(int64_t id, uint64_t* prn_seeds, int offset );
uint64_t init_seed(int64_t id, int offset);
//==============================================================================
//! Set the RNG seeds to unique values based on the ID of the particle. This
@ -57,7 +59,7 @@ void init_seed(int64_t id, uint64_t* prn_seeds, int offset );
//! @param id The particle ID
//==============================================================================
void init_particle_seeds(int64_t id, uint64_t* prn_seeds );
void init_particle_seeds(int64_t id, uint64_t* prn_seeds);
//==============================================================================
//! Advance the random number seed 'n' times from the current seed. This

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.current_seed());
*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.current_seed()) <= 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.current_seed())*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

9
src/eigenvalue.cpp
View file

@ -114,16 +114,13 @@ void synchronize_bank()
fatal_error("No fission sites banked on MPI rank " + std::to_string(mpi::rank));
}
// 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.
int64_t id = simulation::total_gen + overall_generation();
init_seed(id, &prn_seed, STREAM_TRACKING);
advance_prn_seed(start, &prn_seed);
uint64_t seed = init_seed(id, STREAM_TRACKING);
advance_prn_seed(start, &seed);
// Determine how many fission sites we need to sample from the source bank
// and the probability for selecting a site.
@ -158,7 +155,7 @@ void synchronize_bank()
}
// Randomly sample sites needed
if (prn(&prn_seed) < p_sample) {
if (prn(&seed) < p_sample) {
temp_sites[index_temp] = site;
++index_temp;
}

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.current_seed())) ++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.current_seed());
// 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.current_seed());
p.stream_ = STREAM_TRACKING;
int i_low = 0;

7
src/particle.cpp
View file

@ -228,7 +228,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(this->current_seed())) / macro_xs_.total;
}
// Select smaller of the two distances
@ -461,7 +461,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(), this->current_seed());
// Make sure new particle direction is normalized
this->u() = u / u.norm();
@ -678,4 +678,7 @@ Particle::write_restart() const
} // #pragma omp critical
}
uint64_t* Particle::current_seed() {return prn_seeds_ + stream_;}
const uint64_t* Particle::current_seed() const {return prn_seeds_ + stream_;}
} // namespace openmc

10
src/photon.cpp
View file

@ -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.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
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.current_seed());
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.current_seed()) - 1.0;
double phi = 2.0*PI*prn(p.current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);

74
src/physics.cpp
View file

@ -134,7 +134,7 @@ void sample_neutron_reaction(Particle* p)
// 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_);
advance_prn_seed(data::nuclides.size(), p->current_seed());
p->stream_ = STREAM_TRACKING;
}
@ -159,7 +159,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
// Sample the number of neutrons produced
int nu = static_cast<int>(nu_t);
if (prn(p->prn_seeds_ + p->stream_) <= (nu_t - nu)) ++nu;
if (prn(p->current_seed()) <= (nu_t - nu)) ++nu;
// Begin banking the source neutrons
// First, if our bank is full then don't continue
@ -181,7 +181,7 @@ create_fission_sites(Particle* p, int i_nuclide, const Reaction* rx,
site.wgt = 1. / weight;
// Sample delayed group and angle/energy for fission reaction
sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->prn_seeds_ + p->stream_);
sample_fission_neutron(i_nuclide, rx, p->E_, &site, p->current_seed());
// Set the delayed group on the particle as well
p->delayed_group_ = site.delayed_group;
@ -223,13 +223,13 @@ void sample_photon_reaction(Particle* p)
// For tallying purposes, this routine might be called directly. In that
// case, we need to sample a reaction via the cutoff variable
double prob = 0.0;
double cutoff = prn(p->prn_seeds_ + p->stream_) * micro.total;
double cutoff = prn(p->current_seed()) * 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->current_seed());
p->u() = rotate_angle(p->u(), mu, nullptr, p->current_seed());
p->event_ = EVENT_SCATTER;
p->event_mt_ = COHERENT;
return;
@ -240,7 +240,7 @@ void sample_photon_reaction(Particle* p)
if (prob > cutoff) {
double alpha_out, mu;
int i_shell;
element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p->prn_seeds_ + p->stream_);
element.compton_scatter(alpha, true, &alpha_out, &mu, &i_shell, p->current_seed());
// Determine binding energy of shell. The binding energy is 0.0 if
// doppler broadening is not used.
@ -252,13 +252,13 @@ void sample_photon_reaction(Particle* p)
}
// Create Compton electron
double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
double phi = 2.0*PI*prn(p->current_seed());
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->current_seed());
p->create_secondary(u, E_electron, Particle::Type::electron);
}
@ -272,7 +272,7 @@ void sample_photon_reaction(Particle* p)
phi += PI;
p->E_ = alpha_out*MASS_ELECTRON_EV;
p->u() = rotate_angle(p->u(), mu, &phi, p->prn_seeds_ + p->stream_);
p->u() = rotate_angle(p->u(), mu, &phi, p->current_seed());
p->event_ = EVENT_SCATTER;
p->event_mt_ = INCOHERENT;
return;
@ -304,8 +304,8 @@ void sample_photon_reaction(Particle* p)
// model in Serpent 2" by Toni Kaltiaisenaho
double mu;
while (true) {
double r = prn(p->prn_seeds_ + p->stream_);
if (4.0*(1.0 - r)*r >= prn(p->prn_seeds_ + p->stream_)) {
double r = prn(p->current_seed());
if (4.0*(1.0 - r)*r >= prn(p->current_seed())) {
double rel_vel = std::sqrt(E_electron * (E_electron +
2.0*MASS_ELECTRON_EV)) / (E_electron + MASS_ELECTRON_EV);
mu = (2.0*r + rel_vel - 1.0) / (2.0*rel_vel*r - rel_vel + 1.0);
@ -313,7 +313,7 @@ void sample_photon_reaction(Particle* p)
}
}
double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
double phi = 2.0*PI*prn(p->current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
@ -341,14 +341,14 @@ void sample_photon_reaction(Particle* p)
double E_electron, E_positron;
double mu_electron, mu_positron;
element.pair_production(alpha, &E_electron, &E_positron,
&mu_electron, &mu_positron, p->prn_seeds_ + p->stream_);
&mu_electron, &mu_positron, p->current_seed());
// 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->current_seed());
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->current_seed());
p->create_secondary(u, E_positron, Particle::Type::positron);
p->event_ = EVENT_ABSORB;
@ -382,8 +382,8 @@ void sample_positron_reaction(Particle* p)
}
// Sample angle isotropically
double mu = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
double mu = 2.0*prn(p->current_seed()) - 1.0;
double phi = 2.0*PI*prn(p->current_seed());
Direction u;
u.x = mu;
u.y = std::sqrt(1.0 - mu*mu)*std::cos(phi);
@ -401,7 +401,7 @@ void sample_positron_reaction(Particle* p)
int sample_nuclide(Particle* p)
{
// Sample cumulative distribution function
double cutoff = prn(p->prn_seeds_ + p->stream_) * p->macro_xs_.total;
double cutoff = prn(p->current_seed()) * p->macro_xs_.total;
// Get pointers to nuclide/density arrays
const auto& mat {model::materials[p->material_]};
@ -426,7 +426,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->current_seed()) * p->macro_xs_.total;
// Get pointers to elements, densities
const auto& mat {model::materials[p->material_]};
@ -479,7 +479,7 @@ Reaction* sample_fission(int i_nuclide, 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->current_seed()) * p->neutron_xs_[i_nuclide].fission;
double prob = 0.0;
// Loop through each partial fission reaction type
@ -506,7 +506,7 @@ void sample_photon_product(int i_nuclide, 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->current_seed()) * p->neutron_xs_[i_nuclide].photon_prod;
double prob = 0.0;
// Loop through each reaction type
@ -554,7 +554,7 @@ void absorption(Particle* p, int i_nuclide)
} else {
// See if disappearance reaction happens
if (p->neutron_xs_[i_nuclide].absorption >
prn(p->prn_seeds_ + p->stream_) * p->neutron_xs_[i_nuclide].total) {
prn(p->current_seed()) * p->neutron_xs_[i_nuclide].total) {
// Score absorption estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE) {
global_tally_absorption += p->wgt_ * p->neutron_xs_[
@ -582,7 +582,7 @@ void scatter(Particle* p, int i_nuclide)
// For tallying purposes, this routine might be called directly. In that
// case, we need to sample a reaction via the cutoff variable
double cutoff = prn(p->prn_seeds_ + p->stream_) * (micro.total - micro.absorption);
double cutoff = prn(p->current_seed()) * (micro.total - micro.absorption);
bool sampled = false;
// Calculate elastic cross section if it wasn't precalculated
@ -656,8 +656,8 @@ void scatter(Particle* p, int i_nuclide)
int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
if (mat->p0_[i_nuc_mat]) {
// Sample isotropic-in-lab outgoing direction
double mu = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
double phi = 2.0*PI*prn(p->prn_seeds_ + p->stream_);
double mu = 2.0*prn(p->current_seed()) - 1.0;
double phi = 2.0*PI*prn(p->current_seed());
// Change direction of particle
p->u().x = mu;
@ -684,7 +684,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
Direction v_t {};
if (!p->neutron_xs_[i_nuclide].use_ptable) {
v_t = sample_target_velocity(nuc.get(), p->E_, p->u(), v_n,
p->neutron_xs_[i_nuclide].elastic, kT, p->prn_seeds_ + p->stream_);
p->neutron_xs_[i_nuclide].elastic, kT, p->current_seed());
}
// Velocity of center-of-mass
@ -702,9 +702,9 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
auto& d = rx.products_[0].distribution_[0];
auto d_ = dynamic_cast<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->current_seed());
} else {
mu_cm = 2.0*prn(p->prn_seeds_ + p->stream_) - 1.0;
mu_cm = 2.0*prn(p->current_seed()) - 1.0;
}
// Determine direction cosines in CM
@ -713,7 +713,7 @@ void elastic_scatter(int i_nuclide, const Reaction& rx, double kT,
// Rotate neutron velocity vector to new angle -- note that the speed of the
// neutron in CM does not change in elastic scattering. However, the speed
// will change when we convert back to LAB
v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->prn_seeds_ + p->stream_);
v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p->current_seed());
// Transform back to LAB frame
v_n += v_cm;
@ -742,11 +742,11 @@ 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->current_seed());
// 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->current_seed());
}
Direction sample_target_velocity(const Nuclide* nuc, double E, Direction u,
@ -1050,7 +1050,7 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
// sample outgoing energy and scattering cosine
double E;
double mu;
rx->products_[0].sample(E_in, E, mu, p->prn_seeds_ + p->stream_);
rx->products_[0].sample(E_in, E, mu, p->current_seed());
// if scattering system is in center-of-mass, transfer cosine of scattering
// angle and outgoing energy from CM to LAB
@ -1076,7 +1076,7 @@ void inelastic_scatter(const Nuclide* nuc, const Reaction* rx, Particle* p)
p->mu_ = mu;
// change direction of particle
p->u() = rotate_angle(p->u(), mu, nullptr, p->prn_seeds_ + p->stream_);
p->u() = rotate_angle(p->u(), mu, nullptr, p->current_seed());
// evaluate yield
double yield = (*rx->products_[0].yield_)(E_in);
@ -1097,7 +1097,7 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
double y_t = p->wgt_ * p->neutron_xs_[i_nuclide].photon_prod /
p->neutron_xs_[i_nuclide].total;
int y = static_cast<int>(y_t);
if (prn(p->prn_seeds_ + p->stream_) <= y_t - y) ++y;
if (prn(p->current_seed()) <= y_t - y) ++y;
// Sample each secondary photon
for (int i = 0; i < y; ++i) {
@ -1110,10 +1110,10 @@ void sample_secondary_photons(Particle* p, int i_nuclide)
auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx];
double E;
double mu;
rx->products_[i_product].sample(p->E_, E, mu, p->prn_seeds_ + p->stream_);
rx->products_[i_product].sample(p->E_, E, mu, p->current_seed());
// 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->current_seed());
// 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->current_seed()) < p->wgt_ / settings::weight_survive) {
p->wgt_ = settings::weight_survive;
p->wgt_last_ = p->wgt_;
} else {

14
src/physics_mg.cpp
View file

@ -78,10 +78,10 @@ void
scatter(Particle* p)
{
data::mg.macro_xs_[p->material_].sample_scatter(p->g_last_, p->g_, p->mu_,
p->wgt_, p->prn_seeds_ + p->stream_);
p->wgt_, p->current_seed());
// 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->current_seed());
// Update energy value for downstream compatability (in tallying)
p->E_ = data::mg.energy_bin_avg_[p->g_];
@ -103,7 +103,7 @@ create_fission_sites(Particle* p, std::vector<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->current_seed()) <= (nu_t - int(nu_t))) {
nu++;
}
@ -128,10 +128,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->current_seed()) - 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->current_seed() );
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,7 +140,7 @@ create_fission_sites(Particle* p, std::vector<Particle::Bank>& bank)
int dg;
int gout;
data::mg.macro_xs_[p->material_].sample_fission_energy(p->g_, dg, gout,
p->prn_seeds_ + p->stream_);
p->current_seed());
// 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
@ -180,7 +180,7 @@ absorption(Particle* p)
global_tally_absorption += p->wgt_absorb_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;
} else {
if (p->macro_xs_.absorption > prn(p->prn_seeds_ + p->stream_) * p->macro_xs_.total) {
if (p->macro_xs_.absorption > prn(p->current_seed()) * p->macro_xs_.total) {
#pragma omp atomic
global_tally_absorption += p->wgt_ * p->macro_xs_.nu_fission /
p->macro_xs_.absorption;

4
src/random_lcg.cpp
View file

@ -46,9 +46,9 @@ double future_prn(int64_t n, uint64_t prn_seed)
// INIT_SEED
//==============================================================================
void init_seed(int64_t id, uint64_t* prn_seed, int offset)
uint64_t init_seed(int64_t id, int offset)
{
*prn_seed = future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
return future_seed(static_cast<uint64_t>(id) * prn_stride, master_seed + offset);
}
//==============================================================================

2
src/secondary_thermal.cpp
View file

@ -125,7 +125,7 @@ IncoherentElasticAEDiscrete::sample(double E_in, double& E_out, double& mu,
mu_out_(i, k+1) + f*(mu_out_(i+1, k+1) - mu_out_(i, k+1));
// Smear cosine
mu += std::min(mu - mu_left, mu_right - mu)*(prn() - 0.5);
mu += std::min(mu - mu_left, mu_right - mu)*(prn(prn_seed) - 0.5);
// Energy doesn't change in elastic scattering
E_out = E_in;

10
src/source.cpp
View file

@ -270,11 +270,10 @@ 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_seed;
init_seed(id, &prn_seed, STREAM_SOURCE);
uint64_t seed = init_seed(id, STREAM_SOURCE);
// sample external source distribution
simulation::source_bank[i] = sample_external_source(&prn_seed);
simulation::source_bank[i] = sample_external_source(&seed);
}
}
@ -331,11 +330,10 @@ 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_seed;
init_seed(id, &prn_seed, STREAM_SOURCE);
uint64_t seed = init_seed(id, STREAM_SOURCE);
// sample external source distribution
simulation::source_bank[i] = sample_external_source(&prn_seed);
simulation::source_bank[i] = sample_external_source(&seed);
}
}
}

5
src/volume_calc.cpp
View file

@ -128,12 +128,11 @@ std::vector<VolumeCalculation::Result> VolumeCalculation::execute() const
// Sample locations and count hits
#pragma omp for
for (size_t i = i_start; i < i_end; i++) {
uint64_t prn_seed;
int64_t id = iterations * n_samples_ + i;
init_seed(id, &prn_seed, STREAM_VOLUME);
uint64_t seed = init_seed(id, STREAM_VOLUME);
p.n_coord_ = 1;
Position xi {prn(&prn_seed), prn(&prn_seed), prn(&prn_seed)};
Position xi {prn(&seed), prn(&seed), prn(&seed)};
p.r() = lower_left_ + xi*(upper_right_ - lower_left_);
p.u() = {0.5, 0.5, 0.5};