mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-25 12:35:29 -04:00
made secondary particle bank particle owned instead of thread private
This commit is contained in:
John Tramm
parent
f2bd362b4b
commit
ff27ed2bbe
8 changed files with 628 additions and 16 deletions
|
|
@ -22,12 +22,13 @@ namespace simulation {
|
|||
|
||||
extern std::vector<Particle::Bank> source_bank;
|
||||
extern std::vector<Particle::Bank> fission_bank;
|
||||
extern std::vector<Particle::Bank> secondary_bank;
|
||||
//extern std::vector<Particle::Bank> secondary_bank;
|
||||
#ifdef _OPENMP
|
||||
extern std::vector<Particle::Bank> master_fission_bank;
|
||||
#endif
|
||||
|
||||
#pragma omp threadprivate(fission_bank, secondary_bank)
|
||||
//#pragma omp threadprivate(fission_bank, secondary_bank)
|
||||
#pragma omp threadprivate(fission_bank)
|
||||
|
||||
} // namespace simulation
|
||||
|
||||
|
|
|
|||
|
|
@ -294,6 +294,9 @@ public:
|
|||
// Current PRNG state
|
||||
uint64_t seeds_[N_STREAMS]; // current seeds
|
||||
int stream_; // current RNG stream
|
||||
|
||||
// Secondary bank
|
||||
std::vector<Particle::Bank> secondary_bank_;
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ namespace simulation {
|
|||
|
||||
std::vector<Particle::Bank> source_bank;
|
||||
std::vector<Particle::Bank> fission_bank;
|
||||
std::vector<Particle::Bank> secondary_bank;
|
||||
//std::vector<Particle::Bank> secondary_bank;
|
||||
#ifdef _OPENMP
|
||||
std::vector<Particle::Bank> master_fission_bank;
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -84,9 +84,11 @@ Particle::clear()
|
|||
void
|
||||
Particle::create_secondary(Direction u, double E, Type type)
|
||||
{
|
||||
simulation::secondary_bank.emplace_back();
|
||||
//simulation::secondary_bank.emplace_back();
|
||||
secondary_bank_.emplace_back();
|
||||
|
||||
auto& bank {simulation::secondary_bank.back()};
|
||||
//auto& bank {simulation::secondary_bank.back()};
|
||||
auto& bank {secondary_bank_.back()};
|
||||
bank.particle = type;
|
||||
bank.wgt = wgt_;
|
||||
bank.r = this->r();
|
||||
|
|
@ -368,10 +370,13 @@ Particle::transport()
|
|||
// Check for secondary particles if this particle is dead
|
||||
if (!alive_) {
|
||||
// If no secondary particles, break out of event loop
|
||||
if (simulation::secondary_bank.empty()) break;
|
||||
//if (simulation::secondary_bank.empty()) break;
|
||||
if (secondary_bank_.empty()) break;
|
||||
|
||||
this->from_source(&simulation::secondary_bank.back());
|
||||
simulation::secondary_bank.pop_back();
|
||||
//this->from_source(&simulation::secondary_bank.back());
|
||||
this->from_source(&secondary_bank_.back());
|
||||
//simulation::secondary_bank.pop_back();
|
||||
secondary_bank_.pop_back();
|
||||
n_event = 0;
|
||||
|
||||
// Enter new particle in particle track file
|
||||
|
|
|
|||
|
|
@ -98,11 +98,13 @@ void sample_neutron_reaction(Particle* p)
|
|||
create_fission_sites(p, i_nuclide, rx, simulation::fission_bank);
|
||||
} else if (settings::run_mode == RUN_MODE_FIXEDSOURCE &&
|
||||
settings::create_fission_neutrons) {
|
||||
create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank);
|
||||
//create_fission_sites(p, i_nuclide, rx, simulation::secondary_bank);
|
||||
create_fission_sites(p, i_nuclide, rx, p->secondary_bank_);
|
||||
|
||||
// Make sure particle population doesn't grow out of control for
|
||||
// subcritical multiplication problems.
|
||||
if (simulation::secondary_bank.size() >= 10000) {
|
||||
//if (simulation::secondary_bank.size() >= 10000) {
|
||||
if (p->secondary_bank_.size() >= 10000) {
|
||||
fatal_error("The secondary particle bank appears to be growing without "
|
||||
"bound. You are likely running a subcritical multiplication problem "
|
||||
"with k-effective close to or greater than one.");
|
||||
|
|
|
|||
|
|
@ -51,7 +51,8 @@ sample_reaction(Particle* p)
|
|||
create_fission_sites(p, simulation::fission_bank);
|
||||
} else if ((settings::run_mode == RUN_MODE_FIXEDSOURCE) &&
|
||||
(settings::create_fission_neutrons)) {
|
||||
create_fission_sites(p, simulation::secondary_bank);
|
||||
//create_fission_sites(p, simulation::secondary_bank);
|
||||
create_fission_sites(p, p->secondary_bank_);
|
||||
}
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -28,6 +28,600 @@
|
|||
#include <algorithm>
|
||||
#include <string>
|
||||
|
||||
namespace openmc {
|
||||
|
||||
/*
|
||||
extern std::vector<Particle*> calculate_fuel_xs_queue;
|
||||
extern std::vector<Particle*> calculate_nonfuel_xs_queue;
|
||||
extern std::vector<Particle*> advance_particle_queue;
|
||||
extern std::vector<Particle*> surface_crossing_queue;
|
||||
extern std::vector<Particle*> collision_queue;
|
||||
#pragma omp threadprivate(calculate_fuel_xs_queue, calculate_nonfuel_xs_queue, advance_particle_queue, surface_crossing_queue, collision_queue)
|
||||
|
||||
std::vector<Particle*> calculate_fuel_xs_queue;
|
||||
std::vector<Particle*> calculate_nonfuel_xs_queue;
|
||||
std::vector<Particle*> advance_particle_queue;
|
||||
std::vector<Particle*> surface_crossing_queue;
|
||||
std::vector<Particle*> collision_queue;
|
||||
*/
|
||||
|
||||
int * calculate_fuel_xs_queue;
|
||||
int * calculate_nonfuel_xs_queue;
|
||||
int * advance_particle_queue;
|
||||
int * surface_crossing_queue;
|
||||
int * collision_queue;
|
||||
Particle * particles;
|
||||
|
||||
int calculate_fuel_xs_queue_length = 0;
|
||||
int calculate_nonfuel_xs_queue_length = 0;
|
||||
int advance_particle_queue_length = 0;
|
||||
int surface_crossing_queue_length = 0;
|
||||
int collision_queue_length = 0;
|
||||
|
||||
const int MAX_PARTICLES_IN_FLIGHT = 100;
|
||||
|
||||
void init_event_queues(void)
|
||||
{
|
||||
int n_particles = MAX_PARTICLES_IN_FLIGHT;
|
||||
calculate_fuel_xs_queue = new int[n_particles];
|
||||
calculate_nonfuel_xs_queue = new int[n_particles];
|
||||
advance_particle_queue = new int[n_particles];
|
||||
surface_crossing_queue = new int[n_particles];
|
||||
collision_queue = new int[n_particles];
|
||||
particles = new Particle[n_particles];
|
||||
}
|
||||
|
||||
void free_event_queues(void)
|
||||
{
|
||||
delete calculate_fuel_xs_queue;
|
||||
delete calculate_nonfuel_xs_queue;
|
||||
delete advance_particle_queue;
|
||||
delete surface_crossing_queue;
|
||||
delete collision_queue;
|
||||
delete[] particles;
|
||||
}
|
||||
|
||||
constexpr size_t MAX_PARTICLES_PER_THREAD {100};
|
||||
|
||||
Particle::Bank * shared_fission_bank;
|
||||
int shared_fission_bank_length = 0;
|
||||
int shared_fission_bank_max;
|
||||
|
||||
void init_shared_fission_bank(int max)
|
||||
{
|
||||
shared_fission_bank_max = max;
|
||||
shared_fission_bank = new Particle::Bank[max];
|
||||
}
|
||||
|
||||
void free_shared_fission_bank(void)
|
||||
{
|
||||
delete[] shared_fission_bank;
|
||||
shared_fission_bank_length = 0;
|
||||
}
|
||||
|
||||
// TODO: What is going on here?
|
||||
void revive_particle_from_secondary(Particle* p)
|
||||
{
|
||||
//p->from_source(&simulation::secondary_bank.back());
|
||||
p->from_source(&p->secondary_bank_.back());
|
||||
//simulation::secondary_bank.pop_back();
|
||||
p->secondary_bank_.pop_back();
|
||||
// n_event = 0;
|
||||
|
||||
// Enter new particle in particle track file
|
||||
if (p->write_track_) add_particle_track();
|
||||
}
|
||||
|
||||
void dispatch_xs_event(int i)
|
||||
{
|
||||
Particle * p = particles + i;
|
||||
int idx;
|
||||
if (p->material_ == MATERIAL_VOID) {
|
||||
|
||||
#pragma omp atomic capture
|
||||
idx = calculate_nonfuel_xs_queue_length++;
|
||||
//std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl;
|
||||
calculate_nonfuel_xs_queue[idx] = i;
|
||||
} else {
|
||||
if (model::materials[p->material_]->fissionable_) {
|
||||
#pragma omp atomic capture
|
||||
idx = calculate_fuel_xs_queue_length++;
|
||||
//std::cout << "Dispatching particle to Fuel XS queue idx = " << idx << std::endl;
|
||||
calculate_fuel_xs_queue[idx] = i;
|
||||
} else {
|
||||
#pragma omp atomic capture
|
||||
idx = calculate_nonfuel_xs_queue_length++;
|
||||
//std::cout << "Dispatching particle to non Fuel XS queue idx = " << idx << std::endl;
|
||||
calculate_nonfuel_xs_queue[idx] = i;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void process_calculate_xs_events(int * queue, int n)
|
||||
{
|
||||
// Save last_ members, find grid index
|
||||
for (int i = 0; i < n; i++) {
|
||||
Particle *p = particles + queue[i];
|
||||
//std::cout << "particle offset = " << queue[i] << std::endl;
|
||||
// Set the random number stream
|
||||
// TODO: Move RNG seeds to particle storage
|
||||
if (p->type_ == Particle::Type::neutron) {
|
||||
p->stream_ = STREAM_TRACKING;
|
||||
} else {
|
||||
p->stream_ = STREAM_PHOTON;
|
||||
}
|
||||
|
||||
// Store pre-collision particle properties
|
||||
p->wgt_last_ = p->wgt_;
|
||||
p->E_last_ = p->E_;
|
||||
p->u_last_ = p->u();
|
||||
p->r_last_ = p->r();
|
||||
|
||||
// If the cell hasn't been determined based on the particle's location,
|
||||
// initiate a search for the current cell. This generally happens at the
|
||||
// beginning of the history and again for any secondary particles
|
||||
if (p->coord_[p->n_coord_ - 1].cell == C_NONE) {
|
||||
if (!find_cell(p, false)) {
|
||||
p->mark_as_lost("Could not find the cell containing particle "
|
||||
+ std::to_string(p->id_));
|
||||
return;
|
||||
}
|
||||
|
||||
// set birth cell attribute
|
||||
if (p->cell_born_ == C_NONE) p->cell_born_ = p->coord_[p->n_coord_ - 1].cell;
|
||||
}
|
||||
|
||||
// Write particle track.
|
||||
if (p->write_track_) write_particle_track(*p);
|
||||
|
||||
if (settings::check_overlaps) check_cell_overlap(p);
|
||||
|
||||
if (settings::run_CE) {
|
||||
if (p->material_ == p->material_last_ && p->sqrtkT_ != p->sqrtkT_last_) {
|
||||
// Remove particle from queue
|
||||
}
|
||||
}
|
||||
|
||||
// Find energy index on energy grid
|
||||
// TODO: Calculate this separately?
|
||||
int neutron = static_cast<int>(Particle::Type::neutron);
|
||||
p->macro_xs_.i_grid = std::log(p->E_/data::energy_min[neutron]) / simulation::log_spacing;
|
||||
}
|
||||
|
||||
// Calculate nuclide micros
|
||||
for (int i = 0; i < data::nuclides.size(); ++i) {
|
||||
// loop over particles
|
||||
for (int j = 0; j < n; j++) {
|
||||
//Particle * p = particles + queue[i];
|
||||
Particle * p = particles + queue[j];
|
||||
if (p->material_ == MATERIAL_VOID) continue;
|
||||
|
||||
// If material doesn't have this nuclide, skip it
|
||||
const auto& mat {model::materials[p->material_]};
|
||||
if (mat->mat_nuclide_index_[i] == -1) continue;
|
||||
|
||||
// ======================================================================
|
||||
// CHECK FOR S(A,B) TABLE
|
||||
|
||||
// Check if this nuclide matches one of the S(a,b) tables specified.
|
||||
// This relies on thermal_tables_ being sorted by .index_nuclide
|
||||
int i_sab = C_NONE;
|
||||
double sab_frac = 0.0;
|
||||
for (const auto& sab : mat->thermal_tables_) {
|
||||
if (i == sab.index_nuclide) {
|
||||
// Get index in sab_tables
|
||||
i_sab = sab.index_table;
|
||||
sab_frac = sab.fraction;
|
||||
|
||||
// If particle energy is greater than the highest energy for the
|
||||
// S(a,b) table, then don't use the S(a,b) table
|
||||
//if (p->E_ > data::thermal_scatt[i_sab]->threshold()) i_sab = C_NONE;
|
||||
if (p->E_ > data::thermal_scatt[i_sab]->energy_max_) i_sab = C_NONE;
|
||||
}
|
||||
}
|
||||
|
||||
// ======================================================================
|
||||
// CALCULATE MICROSCOPIC CROSS SECTION
|
||||
|
||||
// Calculate microscopic cross section for this nuclide
|
||||
const auto& micro {p->neutron_xs_[i]};
|
||||
if (p->E_ != micro.last_E
|
||||
|| p->sqrtkT_ != micro.last_sqrtkT
|
||||
|| i_sab != micro.index_sab
|
||||
|| sab_frac != micro.sab_frac) {
|
||||
data::nuclides[i]->calculate_xs(i_sab, p->macro_xs_.i_grid, sab_frac, *p);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
for (int i = 0; i < n; i++) {
|
||||
Particle * p = particles + queue[i];
|
||||
// Calculate microscopic and macroscopic cross sections
|
||||
if (p->material_ != MATERIAL_VOID) {
|
||||
// Only works for CE, no MG support
|
||||
//if (settings::run_CE) {
|
||||
// If the material is the same as the last material and the
|
||||
// temperature hasn't changed, we don't need to lookup cross
|
||||
// sections again.
|
||||
model::materials[p->material_]->calculate_xs(*p);
|
||||
}
|
||||
//else {
|
||||
//}
|
||||
else {
|
||||
p->macro_xs_.total = 0.0;
|
||||
p->macro_xs_.absorption = 0.0;
|
||||
p->macro_xs_.fission = 0.0;
|
||||
p->macro_xs_.nu_fission = 0.0;
|
||||
}
|
||||
|
||||
int idx;
|
||||
#pragma omp atomic capture
|
||||
idx = advance_particle_queue_length++;
|
||||
advance_particle_queue[idx] = queue[i];
|
||||
}
|
||||
}
|
||||
|
||||
void process_advance_particle_events()
|
||||
{
|
||||
//for (auto& p : advance_particle_queue) {
|
||||
for (int i = 0; i < advance_particle_queue_length; i++) {
|
||||
Particle * p = particles + advance_particle_queue[i];
|
||||
simulation::trace == (p->id_ == 0);
|
||||
|
||||
// Sample a distance to collision
|
||||
double d_collision;
|
||||
if (p->type_ == Particle::Type::electron ||
|
||||
p->type_ == Particle::Type::positron) {
|
||||
d_collision = 0.0;
|
||||
} else if (p->macro_xs_.total == 0.0) {
|
||||
d_collision = INFINITY;
|
||||
} else {
|
||||
d_collision = -std::log(prn(p->prn_seeds_, p->stream_)) / p->macro_xs_.total;
|
||||
}
|
||||
|
||||
// -------------- break here? -------------------
|
||||
|
||||
// Find the distance to the nearest boundary
|
||||
p->boundary_ = distance_to_boundary(p);
|
||||
|
||||
// Select smaller of the two distances
|
||||
double distance;
|
||||
int idx;
|
||||
if (p->boundary_.distance < d_collision) {
|
||||
#pragma omp atomic capture
|
||||
idx = surface_crossing_queue_length++;
|
||||
surface_crossing_queue[idx] = advance_particle_queue[i];
|
||||
distance = p->boundary_.distance;
|
||||
} else {
|
||||
#pragma omp atomic capture
|
||||
idx = collision_queue_length++;
|
||||
collision_queue[idx] = advance_particle_queue[i];
|
||||
distance = d_collision;
|
||||
}
|
||||
|
||||
// -------------- break here? -------------------
|
||||
|
||||
// Advance particle
|
||||
for (int j = 0; j < p->n_coord_; ++j) {
|
||||
p->coord_[j].r += distance * p->coord_[j].u;
|
||||
}
|
||||
|
||||
// -------------- break here? -------------------
|
||||
|
||||
// Score track-length tallies
|
||||
if (!model::active_tracklength_tallies.empty()) {
|
||||
score_tracklength_tally(p, distance);
|
||||
}
|
||||
|
||||
// Score track-length estimate of k-eff
|
||||
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
|
||||
p->type_ == Particle::Type::neutron) {
|
||||
global_tally_tracklength += p->wgt_ * distance * p->macro_xs_.nu_fission;
|
||||
}
|
||||
|
||||
// Score flux derivative accumulators for differential tallies.
|
||||
if (!model::active_tallies.empty()) {
|
||||
score_track_derivative(p, distance);
|
||||
}
|
||||
}
|
||||
|
||||
advance_particle_queue_length = 0;
|
||||
}
|
||||
|
||||
void process_surface_crossing_events()
|
||||
{
|
||||
//for (auto& p : surface_crossing_queue) {
|
||||
for (int i = 0; i < surface_crossing_queue_length; i++) {
|
||||
Particle * p = particles + surface_crossing_queue[i];
|
||||
// Set surface that particle is on and adjust coordinate levels
|
||||
p->surface_ = p->boundary_.surface_index;
|
||||
p->n_coord_ = p->boundary_.coord_level;
|
||||
|
||||
// Saving previous cell data
|
||||
for (int j = 0; j < p->n_coord_; ++j) {
|
||||
p->cell_last_[j] = p->coord_[j].cell;
|
||||
}
|
||||
p->n_coord_last_ = p->n_coord_;
|
||||
|
||||
if (p->boundary_.lattice_translation[0] != 0 ||
|
||||
p->boundary_.lattice_translation[1] != 0 ||
|
||||
p->boundary_.lattice_translation[2] != 0) {
|
||||
// Particle crosses lattice boundary
|
||||
cross_lattice(p, p->boundary_);
|
||||
p->event_ = EVENT_LATTICE;
|
||||
} else {
|
||||
// Particle crosses surface
|
||||
p->cross_surface();
|
||||
p->event_ = EVENT_SURFACE;
|
||||
}
|
||||
// Score cell to cell partial currents
|
||||
if (!model::active_surface_tallies.empty()) {
|
||||
score_surface_tally(p, model::active_surface_tallies);
|
||||
}
|
||||
|
||||
//if (!p->alive_ && !simulation::secondary_bank.empty()) {
|
||||
if (!p->alive_ && !p->secondary_bank_.empty()) {
|
||||
revive_particle_from_secondary(p);
|
||||
}
|
||||
|
||||
if (p->alive_)
|
||||
{
|
||||
dispatch_xs_event(surface_crossing_queue[i]);
|
||||
}
|
||||
}
|
||||
|
||||
surface_crossing_queue_length = 0;
|
||||
}
|
||||
|
||||
void process_collision_events()
|
||||
{
|
||||
//for (auto& p : collision_queue) {
|
||||
for (int i = 0; i < collision_queue_length; i++) {
|
||||
Particle * p = particles + collision_queue[i];
|
||||
//std::cout << "Beginning collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
|
||||
// Score collision estimate of keff
|
||||
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
|
||||
p->type_ == Particle::Type::neutron) {
|
||||
global_tally_collision += p->wgt_ * p->macro_xs_.nu_fission
|
||||
/ p->macro_xs_.total;
|
||||
}
|
||||
|
||||
// Score surface current tallies -- this has to be done before the collision
|
||||
// since the direction of the particle will change and we need to use the
|
||||
// pre-collision direction to figure out what mesh surfaces were crossed
|
||||
|
||||
if (!model::active_meshsurf_tallies.empty())
|
||||
score_surface_tally(p, model::active_meshsurf_tallies);
|
||||
|
||||
//std::cout << "After surface tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
|
||||
|
||||
// Clear surface component
|
||||
p->surface_ = 0;
|
||||
|
||||
if (settings::run_CE) {
|
||||
collision(p);
|
||||
} else {
|
||||
collision_mg(p);
|
||||
}
|
||||
|
||||
//std::cout << "After collision() of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
|
||||
|
||||
// Score collision estimator tallies -- this is done after a collision
|
||||
// has occurred rather than before because we need information on the
|
||||
// outgoing energy for any tallies with an outgoing energy filter
|
||||
if (!model::active_collision_tallies.empty()) score_collision_tally(p);
|
||||
if (!model::active_analog_tallies.empty()) {
|
||||
if (settings::run_CE) {
|
||||
score_analog_tally_ce(p);
|
||||
} else {
|
||||
score_analog_tally_mg(p);
|
||||
}
|
||||
}
|
||||
//std::cout << "After analog tally of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
|
||||
|
||||
// Reset banked weight during collision
|
||||
p->n_bank_ = 0;
|
||||
p->wgt_bank_ = 0.0;
|
||||
for (int& v : p->n_delayed_bank_) v = 0;
|
||||
|
||||
// Reset fission logical
|
||||
p->fission_ = false;
|
||||
|
||||
// Save coordinates for tallying purposes
|
||||
p->r_last_current_ = p->r();
|
||||
|
||||
// Set last material to none since cross sections will need to be
|
||||
// re-evaluated
|
||||
p->material_last_ = C_NONE;
|
||||
|
||||
// Set all directions to base level -- right now, after a collision, only
|
||||
// the base level directions are changed
|
||||
for (int j = 0; j < p->n_coord_ - 1; ++j) {
|
||||
if (p->coord_[j + 1].rotated) {
|
||||
// If next level is rotated, apply rotation matrix
|
||||
const auto& m {model::cells[p->coord_[j].cell]->rotation_};
|
||||
const auto& u {p->coord_[j].u};
|
||||
p->coord_[j + 1].u.x = m[3]*u.x + m[4]*u.y + m[5]*u.z;
|
||||
p->coord_[j + 1].u.y = m[6]*u.x + m[7]*u.y + m[8]*u.z;
|
||||
p->coord_[j + 1].u.z = m[9]*u.x + m[10]*u.y + m[11]*u.z;
|
||||
} else {
|
||||
// Otherwise, copy this level's direction
|
||||
p->coord_[j+1].u = p->coord_[j].u;
|
||||
}
|
||||
}
|
||||
|
||||
// Score flux derivative accumulators for differential tallies.
|
||||
if (!model::active_tallies.empty()) score_collision_derivative(p);
|
||||
|
||||
//if (!p->alive_ && !simulation::secondary_bank.empty()) {
|
||||
if (!p->alive_ && !p->secondary_bank_.empty()) {
|
||||
revive_particle_from_secondary(p);
|
||||
}
|
||||
|
||||
if (p->alive_)
|
||||
{
|
||||
dispatch_xs_event(collision_queue[i]);
|
||||
//std::cout << "Ended collision of particle id " << collision_queue[i] << " with energy E = " << p->E_ << std::endl;
|
||||
assert(std::isfinite(p->E_) );
|
||||
}
|
||||
}
|
||||
|
||||
collision_queue_length = 0;
|
||||
}
|
||||
|
||||
void check_energies(void)
|
||||
{
|
||||
int * Q;
|
||||
int n;
|
||||
|
||||
|
||||
Q = calculate_fuel_xs_queue;
|
||||
n = calculate_fuel_xs_queue_length;
|
||||
for( int i = 0; i < n; i++ )
|
||||
{
|
||||
if( !std::isfinite(particles[Q[i]].E_ ) )
|
||||
{
|
||||
std::cout << "NAN energy particle found at index xs FUEL " << Q[i] << std::endl;
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
Q = calculate_nonfuel_xs_queue;
|
||||
n = calculate_nonfuel_xs_queue_length ;
|
||||
for( int i = 0; i < n; i++ )
|
||||
{
|
||||
if( !std::isfinite(particles[Q[i]].E_ ) )
|
||||
{
|
||||
std::cout << "NAN energy particle found at index xs Non fuel " << Q[i] << std::endl;
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
Q = advance_particle_queue;
|
||||
n = advance_particle_queue_length ;
|
||||
for( int i = 0; i < n; i++ )
|
||||
{
|
||||
if( !std::isfinite(particles[Q[i]].E_ ) )
|
||||
{
|
||||
std::cout << "NAN energy particle found at index advance particle " << Q[i] << std::endl;
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
Q = surface_crossing_queue;
|
||||
n = surface_crossing_queue_length ;
|
||||
for( int i = 0; i < n; i++ )
|
||||
{
|
||||
if( !std::isfinite(particles[Q[i]].E_ ) )
|
||||
{
|
||||
std::cout << "NAN energy particle found at index surface crossing " << Q[i] << std::endl;
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
Q = collision_queue;
|
||||
n = collision_queue_length ;
|
||||
for( int i = 0; i < n; i++ )
|
||||
{
|
||||
if( !std::isfinite(particles[Q[i]].E_ ) )
|
||||
{
|
||||
std::cout << "NAN energy particle found at index collision " << Q[i] << std::endl;
|
||||
assert(0);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void transport()
|
||||
{
|
||||
int remaining_work = simulation::work_per_rank;
|
||||
int source_offset = 0;
|
||||
|
||||
// Subiterations to complete sets of particles
|
||||
while (remaining_work > 0) {
|
||||
|
||||
// TODO: Do this only once per PI
|
||||
init_event_queues();
|
||||
|
||||
// Figure out work for this subiteration
|
||||
int n_particles = MAX_PARTICLES_IN_FLIGHT;
|
||||
if( n_particles > remaining_work)
|
||||
n_particles = remaining_work;
|
||||
|
||||
//std::cout << "Initializing particle histories..." << std::endl;
|
||||
// Initialize all histories
|
||||
// TODO: Parallelize
|
||||
for (int i = 0; i < n_particles; i++) {
|
||||
initialize_history(particles + i, source_offset + i);
|
||||
}
|
||||
|
||||
//std::cout << "Enqueing particles for XS Lookups..." << std::endl;
|
||||
// Add all particles to advance particle queue
|
||||
// TODO: Parallelize
|
||||
for (int i = 0; i < n_particles; i++) {
|
||||
dispatch_xs_event(i);
|
||||
}
|
||||
|
||||
int event_kernel_executions = 0;
|
||||
while (true) {
|
||||
event_kernel_executions++;
|
||||
/*
|
||||
std::cout << "Fuel XS Lookups = " << calculate_fuel_xs_queue_length << std::endl;
|
||||
std::cout << "Non Fuel XS Lookups = " << calculate_nonfuel_xs_queue_length << std::endl;
|
||||
std::cout << "Advance Particles = " << advance_particle_queue_length << std::endl;
|
||||
std::cout << "Surface Crossings = " << surface_crossing_queue_length << std::endl;
|
||||
std::cout << "Collisions = " << collision_queue_length << std::endl;
|
||||
*/
|
||||
int max = std::max({calculate_fuel_xs_queue_length, calculate_nonfuel_xs_queue_length, advance_particle_queue_length, surface_crossing_queue_length, collision_queue_length});
|
||||
check_energies();
|
||||
if (max == 0) {
|
||||
break;
|
||||
} else if (max == calculate_fuel_xs_queue_length) {
|
||||
//std::cout << "pre fuel XS check..." << std::endl;
|
||||
//check_energies(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
|
||||
//std::cout << "Performing Fuel XS Lookups..." << std::endl;
|
||||
process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
|
||||
calculate_fuel_xs_queue_length = 0;
|
||||
} else if (max == calculate_nonfuel_xs_queue_length) {
|
||||
//std::cout << "pre non fuel XS check..." << std::endl;
|
||||
//check_energies(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
|
||||
// std::cout << "Performing Non Fuel XS Lookups..." << std::endl;
|
||||
process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
|
||||
calculate_nonfuel_xs_queue_length = 0;
|
||||
} else if (max == advance_particle_queue_length) {
|
||||
//std::cout << "pre advancing check..." << std::endl;
|
||||
//check_energies(advance_particle_queue, advance_particle_queue_length);
|
||||
//std::cout << "Advancing Particles..." << std::endl;
|
||||
process_advance_particle_events();
|
||||
} else if (max == surface_crossing_queue_length) {
|
||||
//std::cout << "pre surface crossing check..." << std::endl;
|
||||
//check_energies(surface_crossing_queue, surface_crossing_queue_length);
|
||||
//std::cout << "Surface Crossings..." << std::endl;
|
||||
process_surface_crossing_events();
|
||||
} else if (max == collision_queue_length) {
|
||||
//std::cout << "pre Colliding check..." << std::endl;
|
||||
//check_energies(collision_queue, collision_queue_length);
|
||||
//std::cout << "Colliding..." << std::endl;
|
||||
process_collision_events();
|
||||
}
|
||||
}
|
||||
|
||||
remaining_work -= n_particles;
|
||||
source_offset += n_particles;
|
||||
|
||||
// Should all be zero
|
||||
/*
|
||||
calculate_fuel_xs_queue_length = 0;
|
||||
calculate_nonfuel_xs_queue_length = 0;
|
||||
advance_particle_queue_length = 0;
|
||||
surface_crossing_queue_length = 0;
|
||||
collision_queue_length = 0;
|
||||
*/
|
||||
|
||||
// TODO: Do this only once per PI
|
||||
free_event_queues();
|
||||
//std::cout << "Event kernels retired: " << event_kernel_executions << std::endl;
|
||||
}
|
||||
//shared_fission_bank_length = 0;
|
||||
}
|
||||
|
||||
} // namespace openmc
|
||||
|
||||
//==============================================================================
|
||||
// C API functions
|
||||
//==============================================================================
|
||||
|
|
@ -437,10 +1031,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();
|
||||
#endif
|
||||
// We need to move all the stuff from the shared_fission_bank into the real one.
|
||||
for( int i = 0; i < shared_fission_bank_length; i++ )
|
||||
simulation::fission_bank.push_back(shared_fission_bank[i]);
|
||||
shared_fission_bank_length = 0;
|
||||
|
||||
//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();
|
||||
|
|
|
|||
|
|
@ -1225,8 +1225,10 @@ score_general_ce(Particle* p, int i_tally, int start_index,
|
|||
// We need to substract the energy of the secondary particles since
|
||||
// they will be transported individually later
|
||||
for (auto i = 0; i < p->n_bank_second_; ++i) {
|
||||
auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i;
|
||||
const auto& bank = simulation::secondary_bank[i_bank];
|
||||
//auto i_bank = simulation::secondary_bank.size() - p->n_bank_second_ + i;
|
||||
auto i_bank = p->secondary_bank_.size() - p->n_bank_second_ + i;
|
||||
//const auto& bank = simulation::secondary_bank[i_bank];
|
||||
const auto& bank = p->secondary_bank_[i_bank];
|
||||
if (bank.particle == Particle::Type::photon ||
|
||||
bank.particle == Particle::Type::neutron) {
|
||||
score -= bank.E;
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue