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

Removed old style history() function from particle object, and moved it into new functions in simulation file. Also had to fix particle_restart function to follow similar logic that was moved to initialize_history() function.

Browse source
This commit is contained in:
John Tramm 2020-01-07 23:03:13 +00:00
parent bee602e036
commit 8aaa97a0bd
5 changed files with 56 additions and 351 deletions
Download patch file Download diff file Expand all files Collapse all files

3
include/openmc/particle.h
View file

@ -228,9 +228,6 @@ public:
void event_collide();
void event_revive_from_secondary();
void event_death();
void transport_history_based();
//! Transport a particle from birth to death
void transport();
//! Cross a surface and handle boundary conditions
void cross_surface();

2
include/openmc/simulation.h
View file

@ -74,6 +74,8 @@ void initialize_generation();
void initialize_history(Particle* p, int64_t index_source);
void transport_history_based_inner(Particle& p);
//! Finalize a batch
//!
//! Handles synchronization and accumulation of tallies, calculation of Shannon

298
src/particle.cpp
View file

@ -393,304 +393,6 @@ Particle::event_death()
global_tally_leakage += tally_leakage_;
}
void
Particle::transport_history_based()
{
n_event_ = 0;
while(true)
{
this->event_calculate_xs_I();
this->event_calculate_xs_II();
this->event_advance();
if( collision_distance_ > boundary_.distance )
this->event_cross_surface();
else
this->event_collide();
this->event_revive_from_secondary();
if(!alive_)
break;
}
this->event_death();
}
void
Particle::transport()
{
// Initialize number of events to zero
int n_event = 0;
/*
// Display message if high verbosity or trace is on
if (settings::verbosity >= 9 || trace_) {
write_message("Simulating Particle " + std::to_string(id_));
}
// Add paricle's starting weight to count for normalizing tallies later
#pragma omp atomic
simulation::total_weight += wgt_;
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
for (auto& micro : neutron_xs_) micro.last_E = 0.0;
}
// Prepare to write out particle track.
if (write_track_) add_particle_track(*this);
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty())
flux_derivs_.resize(model::tally_derivs.size(), 0.0);
*/
while (true) {
// Set the random number stream
if (type_ == Particle::Type::neutron) {
stream_ = STREAM_TRACKING;
} else {
stream_ = STREAM_PHOTON;
}
// Store pre-collision particle properties
wgt_last_ = wgt_;
E_last_ = E_;
u_last_ = this->u();
r_last_ = this->r();
// Reset event variables
event_ = EVENT_KILL;
event_nuclide_ = NUCLIDE_NONE;
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 (coord_[n_coord_ - 1].cell == C_NONE) {
if (!find_cell(this, false)) {
this->mark_as_lost("Could not find the cell containing particle "
+ std::to_string(id_));
return;
}
// Set birth cell attribute
if (cell_born_ == C_NONE) cell_born_ = coord_[n_coord_ - 1].cell;
}
// Write particle track.
if (write_track_) write_particle_track(*this);
if (settings::check_overlaps) check_cell_overlap(this);
// Calculate microscopic and macroscopic cross sections
if (material_ != MATERIAL_VOID) {
if (settings::run_CE) {
if (material_ != material_last_ || sqrtkT_ != 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[material_]->calculate_xs(*this);
}
} else {
// Get the MG data; unlike the CE case above, we have to re-calculate
// cross sections for every collision since the cross sections may
// be angle-dependent
data::mg.macro_xs_[material_].calculate_xs(*this);
// Update the particle's group while we know we are multi-group
g_last_ = g_;
}
} else {
macro_xs_.total = 0.0;
macro_xs_.absorption = 0.0;
macro_xs_.fission = 0.0;
macro_xs_.nu_fission = 0.0;
}
// Find the distance to the nearest boundary
auto boundary = distance_to_boundary(this);
// Sample a distance to collision
double d_collision;
if (type_ == Particle::Type::electron ||
type_ == Particle::Type::positron) {
d_collision = 0.0;
} else if (macro_xs_.total == 0.0) {
d_collision = INFINITY;
} else {
d_collision = -std::log(prn(this->current_seed())) / macro_xs_.total;
}
// Select smaller of the two distances
double distance = std::min(boundary.distance, d_collision);
// Advance particle
for (int j = 0; j < n_coord_; ++j) {
coord_[j].r += distance * coord_[j].u;
}
// Score track-length tallies
if (!model::active_tracklength_tallies.empty()) {
score_tracklength_tally(this, distance);
}
// Score track-length estimate of k-eff
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
type_ == Particle::Type::neutron) {
tally_tracklength_ += wgt_ * distance * macro_xs_.nu_fission;
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) {
score_track_derivative(this, distance);
}
if (d_collision > boundary.distance) {
// ====================================================================
// PARTICLE CROSSES SURFACE
// Set surface that particle is on and adjust coordinate levels
surface_ = boundary.surface_index;
n_coord_ = boundary.coord_level;
// Saving previous cell data
for (int j = 0; j < n_coord_; ++j) {
cell_last_[j] = coord_[j].cell;
}
n_coord_last_ = n_coord_;
if (boundary.lattice_translation[0] != 0 ||
boundary.lattice_translation[1] != 0 ||
boundary.lattice_translation[2] != 0) {
// Particle crosses lattice boundary
cross_lattice(this, boundary);
event_ = EVENT_LATTICE;
} else {
// Particle crosses surface
this->cross_surface();
event_ = EVENT_SURFACE;
}
// Score cell to cell partial currents
if (!model::active_surface_tallies.empty()) {
score_surface_tally(this, model::active_surface_tallies);
}
} else {
// ====================================================================
// PARTICLE HAS COLLISION
// Score collision estimate of keff
if (settings::run_mode == RUN_MODE_EIGENVALUE &&
type_ == Particle::Type::neutron) {
tally_collision_ += wgt_ * macro_xs_.nu_fission
/ 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(this, model::active_meshsurf_tallies);
// Clear surface component
surface_ = 0;
if (settings::run_CE) {
collision(this);
} else {
collision_mg(this);
}
// 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(this);
if (!model::active_analog_tallies.empty()) {
if (settings::run_CE) {
score_analog_tally_ce(this);
} else {
score_analog_tally_mg(this);
}
}
// Reset banked weight during collision
n_bank_ = 0;
n_bank_second_ = 0;
wgt_bank_ = 0.0;
for (int& v : n_delayed_bank_) v = 0;
// Reset fission logical
fission_ = false;
// Save coordinates for tallying purposes
r_last_current_ = this->r();
// Set last material to none since cross sections will need to be
// re-evaluated
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 < n_coord_ - 1; ++j) {
if (coord_[j + 1].rotated) {
// If next level is rotated, apply rotation matrix
const auto& m {model::cells[coord_[j].cell]->rotation_};
const auto& u {coord_[j].u};
coord_[j + 1].u = u.rotate(m);
} else {
// Otherwise, copy this level's direction
coord_[j+1].u = coord_[j].u;
}
}
// Score flux derivative accumulators for differential tallies.
if (!model::active_tallies.empty()) score_collision_derivative(this);
}
// If particle has too many events, display warning and kill it
++n_event;
if (n_event == MAX_EVENTS) {
warning("Particle " + std::to_string(id_) +
" underwent maximum number of events.");
alive_ = false;
}
// 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 (secondary_bank_.empty()) break;
//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
if (write_track_) add_particle_track(*this);
}
}
#ifdef DAGMC
if (settings::dagmc) simulation::history.reset();
#endif
// Finish particle track output.
if (write_track_) {
write_particle_track(*this);
finalize_particle_track(*this);
}
// Contribute tally reduction variables to global accumulator
#pragma omp atomic
global_tally_absorption += tally_absorption_;
#pragma omp atomic
global_tally_collision += tally_collision_;
#pragma omp atomic
global_tally_tracklength += tally_tracklength_;
#pragma omp atomic
global_tally_leakage += tally_leakage_;
}
void
Particle::cross_surface()

22
src/particle_restart.cpp
View file

@ -11,6 +11,8 @@
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/derivative.h"
#include "openmc/track_output.h"
#include <algorithm> // for copy
#include <array>
@ -100,8 +102,26 @@ void run_particle_restart()
}
init_particle_seeds(particle_seed, p.seeds_);
// Force calculation of cross-sections by setting last energy to zero
if (settings::run_CE) {
for (auto& micro : p.neutron_xs_) micro.last_E = 0.0;
}
// Prepare to write out particle track.
if (p.write_track_) add_particle_track(p);
// Every particle starts with no accumulated flux derivative.
if (!model::active_tallies.empty())
{
p.flux_derivs_.resize(model::tally_derivs.size(), 0.0);
std::fill(p.flux_derivs_.begin(), p.flux_derivs_.end(), 0.0);
}
// Allocate space for tally filter matches
p.filter_matches_.resize(model::tally_filters.size());
// Transport neutron
p.transport();
transport_history_based_inner(p);
// Write output if particle made it
print_particle(&p);

82
src/simulation.cpp
View file

@ -147,19 +147,6 @@ void free_shared_fission_bank(void)
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(*p);
}
void dispatch_xs_event(int i)
{
Particle * p = particles + i;
@ -303,7 +290,34 @@ double get_time()
return (double) us_since_epoch / 1.0e6;
}
void transport()
void transport_history_based_inner(Particle& p)
{
while(true) {
p.event_calculate_xs_I();
p.event_calculate_xs_II();
p.event_advance();
if( p.collision_distance_ > p.boundary_.distance )
p.event_cross_surface();
else
p.event_collide();
p.event_revive_from_secondary();
if(!p.alive_)
break;
}
p.event_death();
}
void transport_history_based()
{
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
Particle p;
initialize_history(&p, i_work);
transport_history_based_inner(p);
}
}
void transport_event_based()
{
double stop, start;
double time_init = 0;
@ -401,15 +415,6 @@ void transport()
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;
*/
std::cout << "Event kernels retired: " << event_kernel_executions << std::endl;
}
@ -422,7 +427,6 @@ void transport()
std::cout << "Surface Time: " << time_surf << std::endl;
std::cout << "Collision Time: " << time_collision<< std::endl;
}
//shared_fission_bank_length = 0;
free_event_queues();
}
@ -579,32 +583,9 @@ int openmc_next_batch(int* status)
// Start timer for transport
simulation::time_transport.start();
/*
// ====================================================================
// LOOP OVER PARTICLES
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
// grab source particle from bank
Particle p;
initialize_history(&p, i_work);
transport_history_based();
// transport particle
p.transport();
}
*/
/*
#pragma omp parallel for schedule(runtime)
for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
// grab source particle from bank
Particle p;
initialize_history(&p, i_work);
// transport particle
p.transport_history_based();
}
*/
transport();
//transport_event_based();
// Accumulate time for transport
simulation::time_transport.stop();
@ -896,12 +877,15 @@ void initialize_history(Particle* p, int64_t index_source)
// set defaults
p->from_source(&simulation::source_bank[index_source - 1]);
p->current_work_ = index_source;
/*
p->n_event_ = 0;
// Initialize global tally variables
p->tally_absorption_ = 0.0;
p->tally_collision_ = 0.0;
p->tally_tracklength_ = 0.0;
p->tally_leakage_ = 0.0;
*/
// set identifier for particle
p->id_ = simulation::work_index[mpi::rank] + index_source;