From 8aaa97a0bda893016cafa719132a6ed5b1b2fdfc Mon Sep 17 00:00:00 2001 From: John Tramm Date: Tue, 7 Jan 2020 23:03:13 +0000 Subject: [PATCH] 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. --- include/openmc/particle.h | 3 - include/openmc/simulation.h | 2 + src/particle.cpp | 298 ------------------------------------ src/particle_restart.cpp | 22 ++- src/simulation.cpp | 82 ++++------ 5 files changed, 56 insertions(+), 351 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index b88b7eef5..daf574985 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -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(); diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 6b06b646b..a38db4a34 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -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 diff --git a/src/particle.cpp b/src/particle.cpp index 6e7984763..b2f32a8a2 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -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() diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 053193207..dc39a2b00 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -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 // for copy #include @@ -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); diff --git a/src/simulation.cpp b/src/simulation.cpp index 6b9fad9cb..1a38353c3 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -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;