crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Fixed one issue with event-based mode when converting to unified functions. Needed to be attempting to revive after BOTH surface crossing and collions events. Had previously only been doing it for one.

Browse source
This commit is contained in:
John Tramm 2020-01-06 22:47:38 +00:00
parent e937e9c951
commit e59a70ebb7
1 changed files with 6 additions and 371 deletions
Download patch file Download diff file Expand all files Collapse all files

377
src/simulation.cpp
View file

@ -201,86 +201,18 @@ void process_calculate_xs_events(QueueItem * queue, int n)
{
// Sort queue by energy
std::sort(queue, queue+n);
//std::sort(queue, queue+n, by_energy);
// Then, stable sort by material (so as to preserve energy ordering)
//std::stable_sort(queue, queue+n, by_material);
// Save last_ members, find grid index
int lost_particles = 0;
#pragma omp parallel for reduction(+:lost_particles) schedule(dynamic,DYNAMIC_SIZE)
#pragma omp parallel for schedule(dynamic,DYNAMIC_SIZE)
for (int i = 0; i < n; i++) {
Particle *p = particles + queue[i].idx;
//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();
// Reset event variables
p->event_ = EVENT_KILL;
p->event_nuclide_ = NUCLIDE_NONE;
p->event_mt_ = REACTION_NONE;
// 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;
lost_particles += 1;
continue;
}
// 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);
*/
p->event_calculate_xs_I();
}
if( lost_particles > 0 )
exit(1);
#pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
for( int i = 0; i < n; i++ )
{
Particle * p = particles + queue[i].idx;
/*
// Calculate microscopic and macroscopic cross sections
if (p->material_ != MATERIAL_VOID) {
if (settings::run_CE) {
if (p->material_ != p->material_last_ || p->sqrtkT_ != p->sqrtkT_last_) {
// 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 MG not supported
} 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;
}
*/
p->event_calculate_xs_II();
}
@ -300,75 +232,9 @@ void process_calculate_xs_events(QueueItem * queue, int n)
void process_advance_particle_events()
{
//for (auto& p : advance_particle_queue) {
#pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
for (int i = 0; i < advance_particle_queue_length; i++) {
Particle * p = particles + advance_particle_queue[i].idx;
//p->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->current_seed())) / 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].idx = advance_particle_queue[i].idx;
surface_crossing_queue[idx].E = p->E_;
surface_crossing_queue[idx].material = p->material_;
surface_crossing_queue[idx].type = p->type_;
distance = p->boundary_.distance;
} else {
#pragma omp atomic capture
idx = collision_queue_length++;
collision_queue[idx].idx = advance_particle_queue[i].idx;
collision_queue[idx].E = p->E_;
collision_queue[idx].material = p->material_;
collision_queue[idx].type = p->type_;
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) {
p->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);
}
*/
p->event_advance();
if( p->collision_distance_ > p->boundary_.distance )
{
@ -396,48 +262,13 @@ void process_advance_particle_events()
void process_surface_crossing_events()
{
//for (auto& p : surface_crossing_queue) {
#pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
for (int i = 0; i < surface_crossing_queue_length; i++) {
Particle * p = particles + surface_crossing_queue[i].idx;
/*
// 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);
}
*/
p->event_cross_surface();
p->event_revive_from_secondary();
if (p->alive_)
{
dispatch_xs_event(surface_crossing_queue[i].idx);
}
}
surface_crossing_queue_length = 0;
@ -445,166 +276,18 @@ void process_surface_crossing_events()
void process_collision_events()
{
//for (auto& p : collision_queue) {
#pragma omp parallel for schedule(dynamic,DYNAMIC_SIZE)
for (int i = 0; i < collision_queue_length; i++) {
Particle * p = particles + collision_queue[i].idx;
/*
//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) {
p->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->n_bank_second_ = 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 = u.rotate(m);
} 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);
*/
p->event_collide();
//if (!p->alive_ && !simulation::secondary_bank.empty()) {
/*
if (!p->alive_ && !p->secondary_bank_.empty()) {
revive_particle_from_secondary(p);
}
*/
p->event_revive_from_secondary();
if (p->alive_)
{
dispatch_xs_event(collision_queue[i].idx);
//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);
}
}
}
*/
double get_time()
{
@ -648,35 +331,17 @@ void transport()
while (remaining_work > 0) {
start = get_time();
// 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
#pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
for (int i = 0; i < n_particles; i++) {
initialize_history(particles + i, source_offset + i + 1);
}
/*
for( int p = 0; p < mpi::n_procs; p++)
{
MPI_Barrier(mpi::intracomm);
if( p == mpi::rank )
{
*/
#pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
for (int i = 0; i < n_particles; i++) {
initialize_history(particles + i, source_offset + i + 1);
}
/*
}
}
*/
//std::cout << "Enqueing particles for XS Lookups..." << std::endl;
// Add all particles to advance particle queue
// TODO: Parallelize
#pragma omp parallel for schedule(dynamic, DYNAMIC_SIZE)
for (int i = 0; i < n_particles; i++) {
dispatch_xs_event(i);
@ -688,62 +353,32 @@ void transport()
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;
*/
/*
Particle * p = particles + 1;
std::cout << "E = " << p->E_ << " and Position {" <<
p->r().x << ", " <<
p->r().y << ", " <<
p->r().z << "}" << 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;
start = get_time();
process_calculate_xs_events(calculate_fuel_xs_queue, calculate_fuel_xs_queue_length);
stop = get_time();
time_fuel_xs += (stop-start);
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;
start = get_time();
process_calculate_xs_events(calculate_nonfuel_xs_queue, calculate_nonfuel_xs_queue_length);
stop = get_time();
time_nonfuel_xs += (stop-start);
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;
start = get_time();
process_advance_particle_events();
stop = get_time();
time_advance += (stop-start);
} 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;
start = get_time();
process_surface_crossing_events();
stop = get_time();
time_surf += (stop-start);
} 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;
start = get_time();
process_collision_events();
stop = get_time();