From 8bb0fa7791d995f02a94de19f64db683c527c951 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 11 Feb 2019 13:58:12 -0500 Subject: [PATCH] Move write_tallies to C++ --- include/openmc/tallies/tally.h | 40 +++++ src/input_xml.F90 | 8 + src/output.F90 | 235 +---------------------------- src/output.cpp | 182 ++++++++++++++++++++++- src/simulation.cpp | 10 +- src/tallies/tally.cpp | 264 +++++++++++++++++---------------- 6 files changed, 366 insertions(+), 373 deletions(-) diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 2091ece0e..74adc7ed6 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -2,6 +2,7 @@ #define OPENMC_TALLIES_TALLY_H #include "openmc/constants.h" +#include "openmc/tallies/filter.h" #include "openmc/tallies/trigger.h" #include "pugixml.hpp" @@ -21,6 +22,8 @@ class Tally { public: Tally() {} + void init_from_xml(pugi::xml_node node); + void set_scores(pugi::xml_node node); void set_scores(std::vector scores); @@ -50,6 +53,8 @@ public: int id_; //!< user-defined identifier + std::string name_; //!< user-defined name + int type_ {TALLY_VOLUME}; //!< volume, surface current //! Event type that contributes to this tally @@ -93,6 +98,41 @@ private: int32_t n_filter_bins_ {0}; }; +//============================================================================== +//! An iterator over all combinations of a tally's matching filter bins. +// +//! This iterator handles two distinct tasks. First, it maps the N-dimensional +//! space created by the indices of N filters onto a 1D sequence. In other +//! words, it provides a single number that uniquely identifies a combination of +//! bins for many filters. Second, it handles the task of finding each all +//! valid combinations of filter bins given that each filter can have 1 or 2 or +//! many bins that are valid for the current tally event. +//============================================================================== + +class FilterBinIter +{ +public: + FilterBinIter(const Tally& tally, Particle* p); + + FilterBinIter(const Tally& tally, bool end); + + bool operator==(const FilterBinIter& other) const + {return index_ == other.index_;} + + bool operator!=(const FilterBinIter& other) const + {return !(*this == other);} + + FilterBinIter& operator++(); + + int index_ {1}; + double weight_ {1.}; + +private: + void compute_index_weight(); + + const Tally& tally_; +}; + //============================================================================== // Global variable declarations //============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index aa38f7db6..ec4440a12 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -473,6 +473,12 @@ contains type(TallyDerivative), pointer :: deriv interface + subroutine tally_init_from_xml(tally_ptr, xml_node) bind(C) + import C_PTR + type(C_PTR), value :: tally_ptr + type(C_PTR) :: xml_node + end subroutine + subroutine tally_set_scores(tally_ptr, xml_node) bind(C) import C_PTR type(C_PTR), value :: tally_ptr @@ -630,6 +636,8 @@ contains ! Get pointer to tally xml node node_tal = node_tal_list(i) + call tally_init_from_xml(t % ptr, node_tal % ptr) + ! Copy and set tally id if (check_for_node(node_tal, "id")) then call get_node_value(node_tal, "id", tally_id) diff --git a/src/output.F90 b/src/output.F90 index 02f5ec316..aa66ae017 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -39,6 +39,9 @@ module output import Particle type(Particle), intent(in) :: p end subroutine + + subroutine write_tallies() bind(C) + end subroutine end interface contains @@ -441,238 +444,6 @@ contains end subroutine print_results -!=============================================================================== -! WRITE_TALLIES creates an output file and writes out the mean values of all -! tallies and their standard deviations -!=============================================================================== - - subroutine write_tallies() bind(C) - - integer :: i ! index in tallies array - integer :: j ! level in tally hierarchy - integer :: k ! loop index for scoring bins - integer :: n ! loop index for nuclides - integer :: h ! loop index for tally filters - integer :: indent ! number of spaces to preceed output - integer :: filter_index ! index in results array for filters - integer :: score_index ! scoring bin index - integer :: i_nuclide ! index in nuclides array - integer :: i_filt - integer :: unit_tally ! tallies.out file unit - integer :: nr ! number of realizations - real(8) :: t_value ! t-values for confidence intervals - real(8) :: alpha ! significance level for CI - real(8) :: x(2) ! mean and standard deviation - character(MAX_FILE_LEN) :: filename ! name of output file - character(36) :: score_names(N_SCORE_TYPES) ! names of scoring function - character(36) :: score_name ! names of scoring function - ! to be applied at write-time - integer, allocatable :: filter_bins(:) - character(MAX_WORD_LEN) :: temp_name - type(TallyDerivative), pointer :: deriv - - ! Skip if there are no tallies - if (n_tallies == 0) return - - allocate(filter_bins(n_filters)) - - ! Initialize names for scores - score_names(abs(SCORE_FLUX)) = "Flux" - score_names(abs(SCORE_TOTAL)) = "Total Reaction Rate" - score_names(abs(SCORE_SCATTER)) = "Scattering Rate" - score_names(abs(SCORE_NU_SCATTER)) = "Scattering Production Rate" - score_names(abs(SCORE_ABSORPTION)) = "Absorption Rate" - score_names(abs(SCORE_FISSION)) = "Fission Rate" - score_names(abs(SCORE_NU_FISSION)) = "Nu-Fission Rate" - score_names(abs(SCORE_KAPPA_FISSION)) = "Kappa-Fission Rate" - score_names(abs(SCORE_EVENTS)) = "Events" - score_names(abs(SCORE_DECAY_RATE)) = "Decay Rate" - score_names(abs(SCORE_DELAYED_NU_FISSION)) = "Delayed-Nu-Fission Rate" - score_names(abs(SCORE_PROMPT_NU_FISSION)) = "Prompt-Nu-Fission Rate" - score_names(abs(SCORE_INVERSE_VELOCITY)) = "Flux-Weighted Inverse Velocity" - score_names(abs(SCORE_FISS_Q_PROMPT)) = "Prompt fission power" - score_names(abs(SCORE_FISS_Q_RECOV)) = "Recoverable fission power" - score_names(abs(SCORE_CURRENT)) = "Current" - - ! Create filename for tally output - filename = trim(path_output) // "tallies.out" - - ! Open tally file for writing - open(FILE=filename, NEWUNIT=unit_tally, STATUS='replace', ACTION='write') - - ! Calculate t-value for confidence intervals - if (confidence_intervals) then - alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, n_realizations - 1) - else - t_value = ONE - end if - - TALLY_LOOP: do i = 1, n_tallies - associate (t => tallies(i) % obj) - nr = t % n_realizations - - if (confidence_intervals) then - ! Calculate t-value for confidence intervals - alpha = ONE - CONFIDENCE_LEVEL - t_value = t_percentile(ONE - alpha/TWO, nr - 1) - else - t_value = ONE - end if - - ! Write header block - if (t % name == "") then - call header("TALLY " // trim(to_str(t % id())), 1, unit=unit_tally) - else - call header("TALLY " // trim(to_str(t % id())) // ": " & - // trim(t % name), 1, unit=unit_tally) - endif - - ! Write derivative information. - if (t % deriv() /= C_NONE) then - !associate(deriv => tally_derivs(t % deriv())) - deriv => tally_deriv_c(t % deriv()) - select case (deriv % variable) - case (DIFF_DENSITY) - write(unit=unit_tally, fmt="(' Density derivative Material ',A)") & - to_str(deriv % diff_material) - case (DIFF_NUCLIDE_DENSITY) - write(unit=unit_tally, fmt="(' Nuclide density derivative & - &Material ',A,' Nuclide ',A)") & - trim(to_str(deriv % diff_material)), & - trim(nuclides(deriv % diff_nuclide) % name) - case (DIFF_TEMPERATURE) - write(unit=unit_tally, fmt="(' Temperature derivative Material ',& - &A)") to_str(deriv % diff_material) - case default - call fatal_error("Differential tally dependent variable for tally "& - // trim(to_str(t % id())) // " not defined in output.F90.") - end select - !end associate - end if - - ! WARNING: Admittedly, the logic for moving for printing results is - ! extremely confusing and took quite a bit of time to get correct. The - ! logic is structured this way since it is not practical to have a do - ! loop for each filter variable (given that only a few filters are likely - ! to be used for a given tally. - - ! Initialize bins, filter level, and indentation - do h = 1, t % n_filters() - filter_bins(t % filter(h) + 1) = 0 - end do - j = 1 - indent = 0 - - print_bin: do - find_bin: do - ! Check for no filters - if (t % n_filters() == 0) exit find_bin - - ! Increment bin combination - filter_bins(t % filter(j) + 1) = filter_bins(t % filter(j) + 1) + 1 - - ! ================================================================= - ! REACHED END OF BINS FOR THIS FILTER, MOVE TO NEXT FILTER - - if (filter_bins(t % filter(j) + 1) > & - filters(t % filter(j) + 1) % obj % n_bins) then - ! If this is the first filter, then exit - if (j == 1) exit print_bin - - filter_bins(t % filter(j) + 1) = 0 - j = j - 1 - indent = indent - 2 - - ! ================================================================= - ! VALID BIN -- WRITE FILTER INFORMATION OR EXIT TO WRITE RESULTS - - else - ! Check if this is last filter - if (j == t % n_filters()) exit find_bin - - ! Print current filter information - i_filt = t % filter(j) + 1 - write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), & - trim(filters(i_filt) % obj % & - text_label(filter_bins(i_filt))) - indent = indent + 2 - j = j + 1 - end if - - end do find_bin - - ! Print filter information - if (t % n_filters() > 0) then - i_filt = t % filter(j) + 1 - write(UNIT=unit_tally, FMT='(1X,2A)') repeat(" ", indent), & - trim(filters(i_filt) % obj % & - text_label(filter_bins(i_filt))) - end if - - ! Determine scoring index for this bin combination -- note that unlike - ! in the score_tally subroutine, we have to use max(bins,1) since all - ! bins below the lowest filter level will be zeros - - filter_index = 1 - do h = 1, t % n_filters() - filter_index = filter_index & - + (max(filter_bins(t % filter(h)+1) ,1) - 1) * t % stride(h) - end do - - ! Write results for this filter bin combination - score_index = 0 - if (t % n_filters() > 0) indent = indent + 2 - do n = 1, t % n_nuclide_bins() - ! Write label for nuclide - i_nuclide = t % nuclide_bins(n) - if (i_nuclide == -1) then - write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - "Total Material" - else - if (run_CE) then - write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(nuclides(i_nuclide+1) % name) - else - call get_name_c(i_nuclide+1, len(temp_name), temp_name) - write(UNIT=unit_tally, FMT='(1X,2A,1X,A)') repeat(" ", indent), & - trim(temp_name) - end if - end if - - indent = indent + 2 - do k = 1, t % n_score_bins() - score_index = score_index + 1 - - associate(r => t % results(RESULT_SUM:RESULT_SUM_SQ, :, :)) - - if (t % score_bins(k) > 0) then - score_name = reaction_name(t % score_bins(k)) - else - score_name = score_names(abs(t % score_bins(k))) - end if - x(:) = mean_stdev(r(:, score_index, filter_index), nr) - write(UNIT=unit_tally, FMT='(1X,2A,1X,A,"+/- ",A)') & - repeat(" ", indent), score_name, & - to_str(x(1)), trim(to_str(t_value * x(2))) - end associate - end do - indent = indent - 2 - - end do - indent = indent - 2 - - if (t % n_filters() == 0) exit print_bin - - end do print_bin - - end associate - end do TALLY_LOOP - - close(UNIT=unit_tally) - - end subroutine write_tallies - !=============================================================================== ! MEAN_STDEV computes the sample mean and standard deviation of the mean of a ! single tally score diff --git a/src/output.cpp b/src/output.cpp index 16b96d4be..9a64fd41b 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -2,27 +2,37 @@ #include // for std::transform #include // for strlen -#include // for setw -#include -#include #include +#include // for setw, setprecision +#include // for left +#include +#include +#include +#include #include "openmc/capi.h" #include "openmc/cell.h" #include "openmc/constants.h" +#include "openmc/error.h" #include "openmc/geometry.h" #include "openmc/lattice.h" +#include "openmc/math_functions.h" #include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" +#include "openmc/nuclide.h" #include "openmc/plot.h" +#include "openmc/reaction.h" #include "openmc/settings.h" #include "openmc/surface.h" +#include "openmc/tallies/derivative.h" +#include "openmc/tallies/tally.h" namespace openmc { //============================================================================== -void -header(const char* msg, int level) { +std::string +header(const char* msg) { // Determine how many times to repeat the '=' character. int n_prefix = (63 - strlen(msg)) / 2; int n_suffix = n_prefix; @@ -39,10 +49,18 @@ header(const char* msg, int level) { out << "> " << upper << " <"; for (int i = 0; i < n_suffix; i++) out << '='; + return out.str(); +} + +std::string header(const std::string& msg) {return header(msg.c_str());} + +void +header(const char* msg, int level) { + auto out = header(msg); + // Print header based on verbosity level. - if (settings::verbosity >= level) { - std::cout << out.str() << "\n\n"; - } + if (settings::verbosity >= level) + std::cout << out << "\n\n"; } //============================================================================== @@ -230,4 +248,152 @@ print_overlap_check() } } +//============================================================================== + +std::pair +mean_stdev(double sum, double sum_sq, int n) +{ + double mean, std_dev; + mean = sum / n; + if (n > 1) { + std_dev = std::sqrt((sum_sq / n - mean*mean) / (n - 1)); + } else { + std_dev = 0; + } + return {mean, std_dev}; +} + +const std::unordered_map score_names = { + {SCORE_FLUX, "Flux"}, + {SCORE_TOTAL, "Total Reaction Rate"}, + {SCORE_SCATTER, "Scattering Rate"}, + {SCORE_NU_SCATTER, "Scattering Production Rate"}, + {SCORE_ABSORPTION, "Absorption Rate"}, + {SCORE_FISSION, "Fission Rate"}, + {SCORE_NU_FISSION, "Nu-Fission Rate"}, + {SCORE_KAPPA_FISSION, "Kappa-Fission Rate"}, + {SCORE_EVENTS, "Events"}, + {SCORE_DECAY_RATE, "Decay Rate"}, + {SCORE_DELAYED_NU_FISSION, "Delayed-Nu-Fission Rate"}, + {SCORE_PROMPT_NU_FISSION, "Prompt-Nu-Fission Rate"}, + {SCORE_INVERSE_VELOCITY, "Flux-Weighted Inverse Velocity"}, + {SCORE_FISS_Q_PROMPT, "Prompt fission power"}, + {SCORE_FISS_Q_RECOV, "Recoverable fission power"}, + {SCORE_CURRENT, "Current"}, +}; + +//! Create an ASCII output file showing all tally results. + +extern "C" void +write_tallies() +{ + if (model::tallies.empty()) return; + + // Open the tallies.out file. + std::ofstream tallies_out; + tallies_out.open("tallies.out", std::ios::out | std::ios::trunc); + tallies_out << std::setprecision(6); + + // Loop over each tally. + for (auto i_tally = 0; i_tally < model::tallies.size(); ++i_tally) { + const auto& tally {*model::tallies[i_tally]}; + auto results = tally_results(i_tally+1); + // TODO: get this directly from the tally object when it's been translated + int32_t n_realizations; + auto err = openmc_tally_get_n_realizations(i_tally+1, &n_realizations); + + // Calculate t-value for confidence intervals + double t_value = 1; + if (settings::confidence_intervals) { + auto alpha = 1 - CONFIDENCE_LEVEL; + t_value = t_percentile_c(1 - alpha*0.5, n_realizations - 1); + } + + // Write header block. + std::string tally_header("TALLY " + std::to_string(tally.id_)); + if (!tally.name_.empty()) tally_header += ": " + tally.name_; + tallies_out << "\n" << header(tally_header) << "\n\n"; + + // Write derivative information. + if (tally.deriv_ != C_NONE) { + const auto& deriv {model::tally_derivs[tally.deriv_]}; + switch (deriv.variable) { + case DIFF_DENSITY: + tallies_out << " Density derivative Material " + << std::to_string(deriv.diff_material) << "\n"; + break; + case DIFF_NUCLIDE_DENSITY: + tallies_out << " Nuclide density derivative Material " + << std::to_string(deriv.diff_material) << " Nuclide " + // TODO: off-by-one + << data::nuclides[deriv.diff_nuclide-1]->name_ << "\n"; + break; + case DIFF_TEMPERATURE: + tallies_out << " Temperature derivative Material " + << std::to_string(deriv.diff_material) << "\n"; + break; + default: + fatal_error("Differential tally dependent variable for tally " + + std::to_string(tally.id_) + " not defined in output.cpp"); + } + } + + // Loop over all filter bin combinations. + auto filter_iter = FilterBinIter(tally, false); + auto end = FilterBinIter(tally, true); + for (; filter_iter != end; ++filter_iter) { + auto filter_index = filter_iter.index_; + + // Print info about this combination of filter bins. The stride check + // prevents redundant output. + int indent = 0; + for (auto i = 0; i < tally.filters().size(); ++i) { + if ((filter_index-1) % 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]}; + tallies_out << std::string(indent+1, ' ') + << filt.text_label(match.i_bin_) << "\n"; + } + indent += 2; + } + + // Loop over all nuclide and score combinations. + int score_index = 0; + for (auto i_nuclide : tally.nuclides_) { + // Write label for this nuclide bin. + if (i_nuclide == -1) { + tallies_out << std::string(indent+1, ' ') << "Total Material\n"; + } else { + if (settings::run_CE) { + tallies_out << std::string(indent+1, ' ') + << data::nuclides[i_nuclide]->name_ << "\n"; + } else { + tallies_out << std::string(indent+1, ' ') + << data::nuclides_MG[i_nuclide].name << "\n"; + } + } + + // Write the score, mean, and uncertainty. + indent += 2; + for (auto score : tally.scores_) { + std::string score_name = score > 0 ? reaction_name(score) + : score_names.at(score); + double mean, stdev; + //TODO: off-by-one + std::tie(mean, stdev) = mean_stdev( + results(filter_index-1, score_index, RESULT_SUM), + results(filter_index-1, score_index, RESULT_SUM_SQ), + n_realizations); + tallies_out << std::string(indent+1, ' ') << std::left + << std::setw(36) << score_name << " " << mean << " +/- " + << t_value * stdev << "\n"; + score_index += 1; + } + indent -= 2; + } + } + } +} + } // namespace openmc diff --git a/src/simulation.cpp b/src/simulation.cpp index fbbcca9da..978f57e0f 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -144,11 +144,6 @@ int openmc_simulation_finalize() simulation::time_active.stop(); simulation::time_finalize.start(); -#pragma omp parallel - { - simulation::filter_matches.clear(); - } - // Deallocate Fortran variables, set tallies to inactive for (auto& mat : model::materials) { mat->mat_nuclide_index_.clear(); @@ -165,6 +160,11 @@ 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 (int i = 1; i <= n_tallies; ++i) { openmc_tally_set_active(i, false); diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 704bedaa5..f1d402979 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -63,124 +63,6 @@ double global_tally_collision; double global_tally_tracklength; double global_tally_leakage; -//============================================================================== -//! An iterator over all combinations of a tally's matching filter bins. -// -//! This iterator handles two distinct tasks. First, it maps the N-dimensional -//! space created by the indices of N filters onto a 1D sequence. In other -//! words, it provides a single number that uniquely identifies a combination of -//! bins for many filters. Second, it handles the task of finding each all -//! valid combinations of filter bins given that each filter can have 1 or 2 or -//! many bins that are valid for the current tally event. -//============================================================================== - -class FilterBinIter -{ -public: - FilterBinIter(const Tally& tally, Particle* p, bool end) - : tally_{tally} - { - // Handle the special case for an iterator that points to the end. - if (end) { - index_ = -1; - return; - } - - // 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]}; - if (!match.bins_present_) { - match.bins_.clear(); - match.weights_.clear(); - model::tally_filters[i_filt]->get_all_bins(p, tally_.estimator_, match); - match.bins_present_ = true; - } - - // If there are no valid bins for this filter, then there are no valid - // filter bin combinations so all iterators are end iterators. - if (match.bins_.size() == 0) { - index_ = -1; - return; - } - - // Set the index of the bin used in the first filter combination - match.i_bin_ = 1; - } - - // Compute the initial index and weight. - compute_index_weight(); - } - - bool - operator==(const FilterBinIter& other) - { - return index_ == other.index_; - } - - bool - operator!=(const FilterBinIter& other) - { - return !(*this == other); - } - - FilterBinIter& - operator++() - { - // Find the next valid combination of filter bins. To do this, we search - // backwards through the filters until we find the first filter whose bins - // can be incremented. - 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]}; - if (match.i_bin_< match.bins_.size()) { - // The bin for this filter can be incremented. Increment it and do not - // touch any of the remaining filters. - ++match.i_bin_; - done_looping = false; - break; - } else { - // This bin cannot be incremented so reset it and continue to the next - // filter. - match.i_bin_ = 1; - } - } - - if (done_looping) { - // We have visited every valid combination. All done! - index_ = -1; - } else { - // The loop found a new valid combination. Compute the corresponding - // index and weight. - compute_index_weight(); - } - - return *this; - } - - int index_ {1}; - double weight_ {1.}; - -private: - void - compute_index_weight() - { - index_ = 1; - 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 i_bin = match.i_bin_; - //TODO: off-by-one - index_ += (match.bins_[i_bin-1] - 1) * tally_.strides(i); - weight_ *= match.weights_[i_bin-1]; - } - } - - const Tally& tally_; -}; - int score_str_to_int(std::string score_str) { @@ -341,6 +223,12 @@ score_str_to_int(std::string score_str) // Tally object implementation //============================================================================== +void +Tally::init_from_xml(pugi::xml_node node) +{ + if (check_for_node(node, "name")) name_ = get_node_value(node, "name"); +} + void Tally::set_filters(const int32_t filter_indices[], int n) { @@ -632,6 +520,123 @@ Tally::init_triggers(pugi::xml_node node, int i_tally) } } +//============================================================================== +// FilterBinIter implementation +//============================================================================== + +FilterBinIter::FilterBinIter(const Tally& tally, Particle* p) + : tally_{tally} +{ + // 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]}; + if (!match.bins_present_) { + match.bins_.clear(); + match.weights_.clear(); + model::tally_filters[i_filt]->get_all_bins(p, tally_.estimator_, match); + match.bins_present_ = true; + } + + // If there are no valid bins for this filter, then there are no valid + // filter bin combinations so all iterators are end iterators. + if (match.bins_.size() == 0) { + index_ = -1; + return; + } + + // Set the index of the bin used in the first filter combination + match.i_bin_ = 1; + } + + // Compute the initial index and weight. + compute_index_weight(); +} + +FilterBinIter::FilterBinIter(const Tally& tally, bool end) + : tally_{tally} +{ + // Handle the special case for an iterator that points to the end. + if (end) { + index_ = -1; + return; + } + + for (auto i_filt : tally_.filters()) { + auto& match {simulation::filter_matches[i_filt]}; + if (!match.bins_present_) { + match.bins_.clear(); + match.weights_.clear(); + for (auto i = 0; i < model::tally_filters[i_filt]->n_bins_; ++i) { + // TODO: off-by-one + match.bins_.push_back(i+1); + match.weights_.push_back(1.0); + } + match.bins_present_ = true; + } + + if (match.bins_.size() == 0) { + index_ = -1; + return; + } + + match.i_bin_ = 1; + } + + // Compute the initial index and weight. + compute_index_weight(); +} + +FilterBinIter& +FilterBinIter::operator++() +{ + // Find the next valid combination of filter bins. To do this, we search + // backwards through the filters until we find the first filter whose bins + // can be incremented. + 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]}; + if (match.i_bin_< match.bins_.size()) { + // The bin for this filter can be incremented. Increment it and do not + // touch any of the remaining filters. + ++match.i_bin_; + done_looping = false; + break; + } else { + // This bin cannot be incremented so reset it and continue to the next + // filter. + match.i_bin_ = 1; + } + } + + if (done_looping) { + // We have visited every valid combination. All done! + index_ = -1; + } else { + // The loop found a new valid combination. Compute the corresponding + // index and weight. + compute_index_weight(); + } + + return *this; +} + +void +FilterBinIter::compute_index_weight() +{ + index_ = 1; + 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 i_bin = match.i_bin_; + //TODO: off-by-one + index_ += (match.bins_[i_bin-1] - 1) * tally_.strides(i); + weight_ *= match.weights_[i_bin-1]; + } +} + //============================================================================== // Non-member functions //============================================================================== @@ -2511,8 +2516,8 @@ 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, false); - auto end = FilterBinIter(tally, nullptr, true); + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true); if (filter_iter == end) continue; // Loop over filter bins. @@ -2576,8 +2581,8 @@ score_analog_tally_mg(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, false); - auto end = FilterBinIter(tally, nullptr, true); + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true); if (filter_iter == end) continue; // Loop over filter bins. @@ -2636,8 +2641,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, false); - auto end = FilterBinIter(tally, nullptr, true); + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true); if (filter_iter == end) continue; // Loop over filter bins. @@ -2718,8 +2723,8 @@ 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, false); - auto end = FilterBinIter(tally, nullptr, true); + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true); if (filter_iter == end) continue; // Loop over filter bins. @@ -2786,8 +2791,8 @@ score_surface_tally_inner(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, false); - auto end = FilterBinIter(tally, nullptr, true); + auto filter_iter = FilterBinIter(tally, p); + auto end = FilterBinIter(tally, true); if (filter_iter == end) continue; // Loop over filter bins. @@ -3172,6 +3177,9 @@ extern "C" { int active_surface_tallies_size() {return model::active_surface_tallies.size();} + void tally_init_from_xml(Tally* tally, pugi::xml_node* node) + {tally->init_from_xml(*node);} + int tally_get_id_c(Tally* tally) {return tally->id_;} void tally_set_id_c(Tally* tally, int id) {tally->id_ = id;}