diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 6a2e67049fd..405e37e9f04 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -39,6 +39,30 @@ class Particle : public ParticleData { double speed() const; + //! get the distance for particle movement + // + //! sample a distance for next collision based on the xs + //! and find the distance to the boundary along the particle's + //! direction and return the above two's minimum + double get_destination_distance(); + + // Advance particle in space and time + // Short-term solution until the surface source is revised and we can use + // this->move_distance(distance) + void is_hit_time_boundary(double distance, bool& hit_time_boundary); + + void set_particle_weight_to_zero_if_it_hit_time_boundary(bool hit_time_boundary); + + void score_non_mesh_track_length_tallies(double distance); + + void event_advance_deprecated(); + + //! score tallies corresponding particle movement + // + //! score tracklength tally, keff tally and differential tallies + //! \param distance the distance the particle is moved + void score_the_tallies(double distance); + //! moves the particle by the distance length to its next location //! \param length the distance the particle is moved void move_distance(double length); diff --git a/src/event.cpp b/src/event.cpp index ae914c8be34..e34a8df40d7 100644 --- a/src/event.cpp +++ b/src/event.cpp @@ -103,15 +103,32 @@ void process_calculate_xs_events(SharedArray& queue) simulation::time_event_calculate_xs.stop(); } +void openmc_event_advance_wrapper() +{ +#pragma omp parallel for schedule(runtime) + for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) { + int64_t buffer_idx = simulation::advance_particle_queue[i].idx; + Particle& p = simulation::particles[buffer_idx]; + p.event_advance(); + } +} + +void send_particles_to_other_queues(); + void process_advance_particle_events() { simulation::time_event_advance_particle.start(); + openmc_event_advance_wrapper(); + send_particles_to_other_queues(); + simulation::time_event_advance_particle.stop(); +} +void send_particles_to_other_queues() +{ #pragma omp parallel for schedule(runtime) for (int64_t i = 0; i < simulation::advance_particle_queue.size(); i++) { int64_t buffer_idx = simulation::advance_particle_queue[i].idx; Particle& p = simulation::particles[buffer_idx]; - p.event_advance(); if (!p.alive()) continue; if (p.collision_distance() > p.boundary().distance) { @@ -122,8 +139,6 @@ void process_advance_particle_events() } simulation::advance_particle_queue.resize(0); - - simulation::time_event_advance_particle.stop(); } void process_surface_crossing_events() diff --git a/src/particle.cpp b/src/particle.cpp index 64c50c9438f..2be24a72628 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -209,7 +209,89 @@ void Particle::event_calculate_xs() } } + +double Particle::get_destination_distance() +{ + // Find the distance to the nearest boundary + this->boundary() = distance_to_boundary(*this); + + // Sample a distance to collision + if (this->type() == ParticleType::electron || this->type() == ParticleType::positron) { + this->collision_distance() = 0.0; + } else if (this->macro_xs().total == 0.0) { + this->collision_distance() = INFINITY; + } else { + this->collision_distance() = -log(prn(this->current_seed())) / this->macro_xs().total; + } + + // Select smaller of the two distances + double distance = std::min(this->boundary().distance, this->collision_distance()); + return distance; +} + +void Particle::is_hit_time_boundary(double distance, bool& hit_time_boundary) +{ + for (int j = 0; j < this->n_coord(); ++j) { + this->coord(j).r += distance * this->coord(j).u; + } + this->time() += distance / this->speed(); + + // Kill particle if its time exceeds the cutoff + hit_time_boundary = false; + double time_cutoff = settings::time_cutoff[static_cast(this->type())]; + if (this->time() > time_cutoff) { + double dt = this->time() - time_cutoff; + this->time() = time_cutoff; + + double push_back_distance = this->speed() * dt; + this->move_distance(-push_back_distance); + hit_time_boundary = true; + } +} + +void Particle::set_particle_weight_to_zero_if_it_hit_time_boundary(bool hit_time_boundary) +{ + // Set particle weight to zero if it hit the time boundary + if (hit_time_boundary) { + this->wgt() = 0.0; + } +} + +void Particle::score_the_tallies(double distance) +{ + if (!model::active_tracklength_tallies.empty()) { + score_tracklength_tally(*this, distance); + } + // Score track-length estimate of k-eff + score_non_mesh_track_length_tallies(distance); +} + +void Particle::score_non_mesh_track_length_tallies(double distance) +{ + if (settings::run_mode == RunMode::EIGENVALUE && + this->type() == ParticleType::neutron) { + this->keff_tally_tracklength() += this->wgt() * distance * this->macro_xs().nu_fission; + } + + // Score flux derivative accumulators for differential tallies. + if (!model::active_tallies.empty()) { + score_track_derivative(*this, distance); + } +} + + void Particle::event_advance() +{ + double distance = get_destination_distance(); + bool hit_time_boundary = false; + is_hit_time_boundary(distance, hit_time_boundary); + + // Score track-length tallies + score_the_tallies(distance); + set_particle_weight_to_zero_if_it_hit_time_boundary(hit_time_boundary); +} + +void Particle::event_advance_deprecated() { // Find the distance to the nearest boundary boundary() = distance_to_boundary(*this); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e798161ec2f..551a69f7d8f 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2338,14 +2338,77 @@ void score_analog_tally_mg(Particle& p) match.bins_present_ = false; } -void score_tracklength_tally(Particle& p, double distance) + +double get_atom_density( + Particle& p, int i_log_union, const Tally& tally, int i_nuclide) { - // Determine the tracklength estimate of the flux - double flux = p.wgt() * distance; + double atom_density = 0.; + if (i_nuclide >= 0) { + if (p.material() != MATERIAL_VOID) { + const auto& mat = model::materials[p.material()]; + auto j = mat->mat_nuclide_index_[i_nuclide]; + if (j == C_NONE) { + // Determine log union grid index + if (i_log_union == C_NONE) { + int neutron = static_cast(ParticleType::neutron); + i_log_union = log(p.E() / data::energy_min[neutron]) / + simulation::log_spacing; + } - // Set 'none' value for log union grid index - int i_log_union = C_NONE; + // Update micro xs cache + if (!tally.multiply_density()) { + p.update_neutron_xs(i_nuclide, i_log_union); + atom_density = 1.0; + } + } else { + atom_density = + tally.multiply_density() ? mat->atom_density_(j) : 1.0; + } + } + } + return atom_density; +} + +void loop_over_filters(Particle& p, double flux, int i_log_union, + int i_tally, const Tally& tally) +{ + // Initialize an iterator over valid filter bin combinations. If there are + // no valid combinations, use a continue statement to ensure we skip the + // assume_separate break below. + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true, &p.filter_matches()); + if (filter_iter == end) + return; + // Loop over filter bins. + for (; filter_iter != end; ++filter_iter) { + auto filter_index = filter_iter.index_; + auto filter_weight = filter_iter.weight_; + + // Loop over nuclide bins. + for (auto i = 0; i < tally.nuclides_.size(); ++i) { + auto i_nuclide = tally.nuclides_[i]; + + double atom_density = + get_atom_density(p, i_log_union, tally, i_nuclide); + + // TODO: consider replacing this "if" with pointers or templates + // handles all estimators except tracklength estimator + // tracklength estimator tallies are done in the FilterBinIterator call above + if (settings::run_CE) { + score_general_ce_nonanalog(p, i_tally, i * tally.scores_.size(), + filter_index, filter_weight, i_nuclide, atom_density, flux); + } else { + score_general_mg(p, i_tally, i * tally.scores_.size(), filter_index, + filter_weight, i_nuclide, atom_density, flux); + } + } + } +} + +// compute all the active tracklength tallies +void loop_over_all_tallies(Particle& p, double flux, int i_log_union) +{ for (auto i_tally : model::active_tracklength_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2357,50 +2420,7 @@ void score_tracklength_tally(Particle& p, double distance) if (filter_iter == end) continue; - // Loop over filter bins. - for (; filter_iter != end; ++filter_iter) { - auto filter_index = filter_iter.index_; - auto filter_weight = filter_iter.weight_; - - // Loop over nuclide bins. - for (auto i = 0; i < tally.nuclides_.size(); ++i) { - auto i_nuclide = tally.nuclides_[i]; - - double atom_density = 0.; - if (i_nuclide >= 0) { - if (p.material() != MATERIAL_VOID) { - const auto& mat = model::materials[p.material()]; - auto j = mat->mat_nuclide_index_[i_nuclide]; - if (j == C_NONE) { - // Determine log union grid index - if (i_log_union == C_NONE) { - int neutron = static_cast(ParticleType::neutron); - i_log_union = std::log(p.E() / data::energy_min[neutron]) / - simulation::log_spacing; - } - - // Update micro xs cache - if (!tally.multiply_density()) { - p.update_neutron_xs(i_nuclide, i_log_union); - atom_density = 1.0; - } - } else { - atom_density = - tally.multiply_density() ? mat->atom_density_(j) : 1.0; - } - } - } - - // TODO: consider replacing this "if" with pointers or templates - if (settings::run_CE) { - score_general_ce_nonanalog(p, i_tally, i * tally.scores_.size(), - filter_index, filter_weight, i_nuclide, atom_density, flux); - } else { - score_general_mg(p, i_tally, i * tally.scores_.size(), filter_index, - filter_weight, i_nuclide, atom_density, flux); - } - } - } + loop_over_filters(p, flux, i_log_union, i_tally, tally); // If the user has specified that we can assume all tallies are spatially // separate, this implies that once a tally has been scored to, we needn't @@ -2409,6 +2429,16 @@ void score_tracklength_tally(Particle& p, double distance) if (settings::assume_separate) break; } +} + +void score_tracklength_tally(Particle& p, double distance) +{ + // Determine the tracklength estimate of the flux + double flux = p.wgt() * distance; + + // Set 'none' value for log union grid index + int i_log_union = C_NONE; + loop_over_all_tallies(p, flux, i_log_union); // Reset all the filter matches for the next tally event. for (auto& match : p.filter_matches())