From 4fb438cbead271c03fa7451fd506918255c4948a Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 12 Dec 2019 21:09:55 +0000 Subject: [PATCH] moved filter_matches arrays to be particle owned rather than global thread private --- include/openmc/geometry.h | 1 + include/openmc/particle.h | 4 + include/openmc/tallies/filter.h | 12 --- include/openmc/tallies/tally_scoring.h | 11 ++- src/output.cpp | 11 ++- src/simulation.cpp | 14 +--- src/tallies/filter.cpp | 4 - src/tallies/tally_scoring.cpp | 106 ++++++++++++------------- 8 files changed, 76 insertions(+), 87 deletions(-) diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index ef226db66..a7ac53ea3 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -7,6 +7,7 @@ #include #include "openmc/particle.h" +#include "openmc/tallies/filter.h" namespace openmc { diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 6d1e05a49..251781f57 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -14,8 +14,10 @@ #include "openmc/constants.h" #include "openmc/position.h" #include "openmc/random_lcg.h" +//#include "openmc/tallies/filter.h" namespace openmc { +class FilterMatch; //============================================================================== // Constants @@ -315,6 +317,8 @@ public: int64_t current_work_; // current work index std::vector flux_derivs_; // Derivatives of the current particle's weight + + std::vector filter_matches_; }; } // namespace openmc diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 7c46a24c2..dbe5252ea 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -32,11 +32,6 @@ public: } // namespace openmc -// Without an explicit instantiation of vector, the Intel compiler -// will complain about the threadprivate directive on filter_matches. Note that -// this has to happen *outside* of the openmc namespace -extern template class std::vector; - namespace openmc { //============================================================================== @@ -127,13 +122,6 @@ private: // Global variables //============================================================================== -namespace simulation { - -extern std::vector filter_matches; -#pragma omp threadprivate(filter_matches) - -} // namespace simulation - namespace model { extern "C" int32_t n_filters; extern std::vector> tally_filters; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 2c5a29c17..c73b22e7b 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -3,6 +3,7 @@ #include "openmc/particle.h" #include "openmc/tallies/tally.h" +#include "openmc/tallies/filter.h" namespace openmc { @@ -22,12 +23,12 @@ class FilterBinIter public: //! Construct an iterator over bins that match a given particle's state. - FilterBinIter(const Tally& tally, const Particle* p); + FilterBinIter(const Tally& tally, const Particle* p, std::vector * particle_filter_matches); //! Construct an iterator over all filter bin combinations. // //! \param end if true, the returned iterator indicates the end of a loop. - FilterBinIter(const Tally& tally, bool end); + FilterBinIter(const Tally& tally, bool end, std::vector * particle_filter_matches); bool operator==(const FilterBinIter& other) const {return index_ == other.index_;} @@ -39,6 +40,8 @@ public: int index_ {1}; double weight_ {1.}; + + std::vector & filter_matches_; private: void compute_index_weight(); @@ -73,7 +76,7 @@ void score_analog_tally_ce(Particle* p); //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_mg(const Particle* p); +void score_analog_tally_mg(Particle* p); //! Score tallies using a tracklength estimate of the flux. // @@ -89,7 +92,7 @@ void score_tracklength_tally(Particle* p, double distance); // //! \param p The particle being tracked //! \param tallies A vector of tallies to score to -void score_surface_tally(const Particle* p, const std::vector& tallies); +void score_surface_tally(Particle* p, const std::vector& tallies); } // namespace openmc diff --git a/src/output.cpp b/src/output.cpp index ff506410c..5a9d7de19 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -679,9 +679,14 @@ write_tallies() } } + // Initialize Filter Matches Object + std::vector filter_matches; + // Allocate space for tally filter matches + filter_matches.resize(model::tally_filters.size()); + // Loop over all filter bin combinations. - auto filter_iter = FilterBinIter(tally, false); - auto end = FilterBinIter(tally, true); + auto filter_iter = FilterBinIter(tally, false, &filter_matches); + auto end = FilterBinIter(tally, true, &filter_matches); for (; filter_iter != end; ++filter_iter) { auto filter_index = filter_iter.index_; @@ -692,7 +697,7 @@ write_tallies() if (filter_index % tally.strides(i) == 0) { auto i_filt = tally.filters(i); const auto& filt {*model::tally_filters[i_filt]}; - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches[i_filt]}; tallies_out << std::string(indent+1, ' ') << filt.text_label(match.i_bin_) << "\n"; } diff --git a/src/simulation.cpp b/src/simulation.cpp index 1b65231de..a8a3e5f28 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -821,12 +821,6 @@ int openmc_simulation_init() // Determine how much work each process should do calculate_work(); - // Allocate array for matching filter bins - #pragma omp parallel - { - simulation::filter_matches.resize(model::tally_filters.size()); - } - // Allocate source bank, and for eigenvalue simulations also allocate the // fission bank allocate_banks(); @@ -901,11 +895,6 @@ int openmc_simulation_finalize() // Write tally results to tallies.out if (settings::output_tallies && mpi::master) write_tallies(); - #pragma omp parallel - { - simulation::filter_matches.clear(); - } - // Deactivate all tallies for (auto& t : model::tallies) { t->active_ = false; @@ -1303,6 +1292,9 @@ void initialize_history(Particle* p, int64_t index_source) 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()); } int overall_generation() diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index d8fe58cc7..8d4092c1b 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -38,10 +38,6 @@ namespace openmc { // Global variables //============================================================================== -namespace simulation { - std::vector filter_matches; -} - namespace model { std::vector> tally_filters; std::unordered_map filter_map; diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 97ea37921..b3f1b3676 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -26,13 +26,13 @@ namespace openmc { // FilterBinIter implementation //============================================================================== -FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p) - : tally_{tally} +FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p, std::vector * particle_filter_matches) + : tally_{tally}, filter_matches_{*particle_filter_matches} { // Find all valid bins in each relevant filter if they have not already been // found for this event. for (auto i_filt : tally_.filters()) { - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; if (!match.bins_present_) { match.bins_.clear(); match.weights_.clear(); @@ -55,8 +55,8 @@ FilterBinIter::FilterBinIter(const Tally& tally, const Particle* p) this->compute_index_weight(); } -FilterBinIter::FilterBinIter(const Tally& tally, bool end) - : tally_{tally} +FilterBinIter::FilterBinIter(const Tally& tally, bool end, std::vector * particle_filter_matches) + : tally_{tally}, filter_matches_{*particle_filter_matches} { // Handle the special case for an iterator that points to the end. if (end) { @@ -65,7 +65,7 @@ FilterBinIter::FilterBinIter(const Tally& tally, bool end) } for (auto i_filt : tally_.filters()) { - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; if (!match.bins_present_) { match.bins_.clear(); match.weights_.clear(); @@ -97,7 +97,7 @@ FilterBinIter::operator++() bool done_looping = true; for (int i = tally_.filters().size()-1; i >= 0; --i) { auto i_filt = tally_.filters(i); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; if (match.i_bin_ < match.bins_.size()-1) { // The bin for this filter can be incremented. Increment it and do not // touch any of the remaining filters. @@ -130,7 +130,7 @@ FilterBinIter::compute_index_weight() weight_ = 1.; for (auto i = 0; i < tally_.filters().size(); ++i) { auto i_filt = tally_.filters(i); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches_[i_filt]}; auto i_bin = match.i_bin_; index_ += match.bins_[i_bin] * tally_.strides(i); weight_ *= match.weights_[i_bin]; @@ -144,12 +144,12 @@ FilterBinIter::compute_index_weight() //! Helper function used to increment tallies with a delayed group filter. void -score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) +score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index, std::vector & filter_matches) { // Save the original delayed group bin auto& tally {*model::tallies[i_tally]}; auto i_filt = tally.filters(tally.delayedgroup_filter_); - auto& dg_match {simulation::filter_matches[i_filt]}; + auto& dg_match {filter_matches[i_filt]}; auto i_bin = dg_match.i_bin_; auto original_bin = dg_match.bins_[i_bin]; dg_match.bins_[i_bin] = d_bin; @@ -158,7 +158,7 @@ score_fission_delayed_dg(int i_tally, int d_bin, double score, int score_index) auto filter_index = 0; for (auto i = 0; i < tally.filters().size(); ++i) { auto i_filt = tally.filters(i); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {filter_matches[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(i); } @@ -314,12 +314,12 @@ double score_neutron_heating(const Particle* p, const Tally& tally, double flux, //! neutrons produced with different energies. void -score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) +score_fission_eout(Particle* p, int i_tally, int i_score, int score_bin) { auto& tally {*model::tallies[i_tally]}; auto i_eout_filt = tally.filters()[tally.energyout_filter_]; - auto i_bin = simulation::filter_matches[i_eout_filt].i_bin_; - auto bin_energyout = simulation::filter_matches[i_eout_filt].bins_[i_bin]; + auto i_bin = p->filter_matches_[i_eout_filt].i_bin_; + auto bin_energyout = p->filter_matches_[i_eout_filt].bins_[i_bin]; const EnergyoutFilter& eo_filt {*dynamic_cast(model::tally_filters[i_eout_filt].get())}; @@ -356,7 +356,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) g_out = eo_filt.n_bins() - g_out - 1; // change outgoing energy bin - simulation::filter_matches[i_eout_filt].bins_[i_bin] = g_out; + p->filter_matches_[i_eout_filt].bins_[i_bin] = g_out; } else { @@ -373,7 +373,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) } else { auto i_match = lower_bound_index(eo_filt.bins().begin(), eo_filt.bins().end(), E_out); - simulation::filter_matches[i_eout_filt].bins_[i_bin] = i_match; + p->filter_matches_[i_eout_filt].bins_[i_bin] = i_match; } } @@ -387,7 +387,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) int filter_index = 0; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {p->filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); } @@ -415,13 +415,13 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {p->filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_weight *= match.weights_[i_bin]; } score_fission_delayed_dg(i_tally, d_bin, score*filter_weight, - i_score); + i_score, p->filter_matches_); } } @@ -433,7 +433,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) double filter_weight = 1.; for (auto j = 0; j < tally.filters().size(); ++j) { auto i_filt = tally.filters(j); - auto& match {simulation::filter_matches[i_filt]}; + auto& match {p->filter_matches_[i_filt]}; auto i_bin = match.i_bin_; filter_index += match.bins_[i_bin] * tally.strides(j); filter_weight *= match.weights_[i_bin]; @@ -447,7 +447,7 @@ score_fission_eout(const Particle* p, int i_tally, int i_score, int score_bin) } // Reset outgoing energy bin and score index - simulation::filter_matches[i_eout_filt].bins_[i_bin] = bin_energyout; + p->filter_matches_[i_eout_filt].bins_[i_bin] = bin_energyout; } //! Update tally results for continuous-energy tallies with any estimator. @@ -777,7 +777,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, * p->neutron_xs_[p->event_nuclide_].fission / p->neutron_xs_[p->event_nuclide_].absorption * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } continue; } else { @@ -811,7 +811,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto d = filt.groups()[d_bin]; score = simulation::keff * p->wgt_bank_ / p->n_bank_ * p->n_delayed_bank_[d-1] * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -836,7 +836,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, ->nu(E, ReactionProduct::EmissionMode::delayed, d); score = p->neutron_xs_[i_nuclide].fission * yield * atom_density * flux; - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -866,7 +866,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = p->neutron_xs_[j_nuclide].fission * yield * atom_density * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } } } @@ -917,7 +917,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, / p->neutron_xs_[p->event_nuclide_].absorption * rate * flux; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } continue; } else { @@ -971,7 +971,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto d = filt.groups()[d_bin]; if (d == g) score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } score = 0.; } @@ -996,7 +996,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, auto rate = rxn.products_[d].decay_rate_; score = p->neutron_xs_[i_nuclide].fission * yield * flux * atom_density * rate; - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1037,7 +1037,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, score = p->neutron_xs_[j_nuclide].fission * yield * flux * atom_density * rate; score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } } } @@ -1338,7 +1338,7 @@ score_general_ce(Particle* p, int i_tally, int start_index, //! argument is really just used for filter weights. void -score_general_mg(const Particle* p, int i_tally, int start_index, +score_general_mg(Particle* p, int i_tally, int start_index, int filter_index, int i_nuclide, double atom_density, double flux) { auto& tally {*model::tallies[i_tally]}; @@ -1719,7 +1719,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1760,7 +1760,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score *= atom_density * nuc_xs.get_xs(MG_GET_XS_FISSION, p_g) / macro_xs.get_xs(MG_GET_XS_FISSION, p_g); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1791,7 +1791,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, score = flux * macro_xs.get_xs( MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1834,7 +1834,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d) / abs_xs; } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -1895,7 +1895,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, auto dg = filt.groups()[d_bin]; if (dg == d + 1) score_fission_delayed_dg(i_tally, d_bin, score, - score_index); + score_index, p->filter_matches_); } score = 0.; } @@ -1925,7 +1925,7 @@ score_general_mg(const Particle* p, int i_tally, int start_index, * macro_xs.get_xs(MG_GET_XS_DELAYED_NU_FISSION, p_g, nullptr, nullptr, &d); } - score_fission_delayed_dg(i_tally, d_bin, score, score_index); + score_fission_delayed_dg(i_tally, d_bin, score, score_index, p->filter_matches_); } continue; } else { @@ -2041,8 +2041,8 @@ void score_analog_tally_ce(Particle* p) // 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); + auto filter_iter = FilterBinIter(tally, p, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2086,11 +2086,11 @@ void score_analog_tally_ce(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } -void score_analog_tally_mg(const Particle* p) +void score_analog_tally_mg(Particle* p) { for (auto i_tally : model::active_analog_tallies) { const Tally& tally {*model::tallies[i_tally]}; @@ -2098,8 +2098,8 @@ void score_analog_tally_mg(const Particle* p) // 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); + auto filter_iter = FilterBinIter(tally, p, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2131,7 +2131,7 @@ void score_analog_tally_mg(const Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } @@ -2147,8 +2147,8 @@ score_tracklength_tally(Particle* p, double distance) // 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); + auto filter_iter = FilterBinIter(tally, p, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2195,7 +2195,7 @@ score_tracklength_tally(Particle* p, double distance) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } @@ -2215,8 +2215,8 @@ void score_collision_tally(Particle* p) // 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); + auto filter_iter = FilterBinIter(tally, p, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2259,12 +2259,12 @@ void score_collision_tally(Particle* p) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; } void -score_surface_tally(const Particle* p, const std::vector& tallies) +score_surface_tally(Particle* p, const std::vector& tallies) { // No collision, so no weight change when survival biasing double flux = p->wgt_; @@ -2275,8 +2275,8 @@ score_surface_tally(const Particle* p, const std::vector& tallies) // 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); + auto filter_iter = FilterBinIter(tally, p, &p->filter_matches_); + auto end = FilterBinIter(tally, true, &p->filter_matches_); if (filter_iter == end) continue; // Loop over filter bins. @@ -2303,7 +2303,7 @@ score_surface_tally(const Particle* p, const std::vector& tallies) } // Reset all the filter matches for the next tally event. - for (auto& match : simulation::filter_matches) + for (auto& match : p->filter_matches_) match.bins_present_ = false; }