From 8d3f9d3c69d06e9233d6d064b6fd5ae714b5a5bf Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 27 Oct 2018 15:11:49 -0400 Subject: [PATCH] Move most EnergyFilter functionality to C++ --- CMakeLists.txt | 1 + include/openmc/constants.h | 1 + include/openmc/tallies/filter_energy.h | 55 +++++++++++ src/geometry.cpp | 2 - src/tallies/filter.cpp | 5 + src/tallies/filter_energy.cpp | 116 +++++++++++++++++++++++ src/tallies/tally_filter_energy.F90 | 124 +------------------------ 7 files changed, 182 insertions(+), 122 deletions(-) create mode 100644 include/openmc/tallies/filter_energy.h create mode 100644 src/tallies/filter_energy.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 4ac9d9c5d7..636c5fb532 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -438,6 +438,7 @@ add_library(libopenmc SHARED src/tallies/filter_cell.cpp src/tallies/filter_distribcell.cpp src/tallies/filter_energyfunc.cpp + src/tallies/filter_energy.cpp src/tallies/filter_legendre.cpp src/tallies/filter_material.cpp src/tallies/filter_mesh.cpp diff --git a/include/openmc/constants.h b/include/openmc/constants.h index d249e4cb5b..75da02c0a6 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -420,6 +420,7 @@ constexpr int DIFF_NUCLIDE_DENSITY {2}; constexpr int DIFF_TEMPERATURE {3}; constexpr int C_NONE {-1}; +constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE // Interpolation rules enum class Interpolation { diff --git a/include/openmc/tallies/filter_energy.h b/include/openmc/tallies/filter_energy.h new file mode 100644 index 0000000000..9f613e33d4 --- /dev/null +++ b/include/openmc/tallies/filter_energy.h @@ -0,0 +1,55 @@ +#ifndef OPENMC_TALLIES_FILTER_ENERGY_H +#define OPENMC_TALLIES_FILTER_ENERGY_H + +//#include +//#include +#include + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Bins the incident neutron energy. +//============================================================================== + +class EnergyFilter : public Filter +{ +public: + ~EnergyFilter() = default; + + std::string type() const override {return "energy";} + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + std::vector bins_; + bool matches_transport_groups_ {false}; +}; + +//============================================================================== +//! Bins the outgoing neutron energy. +//! +//! Only scattering events use the get_all_bins functionality. Nu-fission +//! tallies manually iterate over the filter bins. +//============================================================================== + +class EnergyoutFilter : public EnergyFilter +{ +public: + std::string type() const override {return "energyout";} + + void get_all_bins(const Particle* p, int estimator, FilterMatch& match) + const override; + + std::string text_label(int bin) const override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_ENERGY_H diff --git a/src/geometry.cpp b/src/geometry.cpp index 28885a5db3..cff5657a56 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -16,8 +16,6 @@ namespace openmc { std::vector overlap_check_count; -constexpr int F90_NONE {0}; //TODO: replace usage of this with C_NONE - //============================================================================== extern "C" bool diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 44758f3799..cd79d44ae9 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -11,6 +11,7 @@ #include "openmc/tallies/filter_cellfrom.h" #include "openmc/tallies/filter_distribcell.h" #include "openmc/tallies/filter_energyfunc.h" +#include "openmc/tallies/filter_energy.h" #include "openmc/tallies/filter_legendre.h" #include "openmc/tallies/filter_material.h" #include "openmc/tallies/filter_mesh.h" @@ -105,6 +106,10 @@ extern "C" { tally_filters.push_back(new DistribcellFilter()); } else if (type_ == "energyfunction") { tally_filters.push_back(new EnergyFunctionFilter()); + } else if (type_ == "energy") { + tally_filters.push_back(new EnergyFilter()); + } else if (type_ == "energyout") { + tally_filters.push_back(new EnergyoutFilter()); } else if (type_ == "legendre") { tally_filters.push_back(new LegendreFilter()); } else if (type_ == "material") { diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp new file mode 100644 index 0000000000..35d34bf0a4 --- /dev/null +++ b/src/tallies/filter_energy.cpp @@ -0,0 +1,116 @@ +#include "openmc/tallies/filter_energy.h" + +#include "openmc/constants.h" // For F90_NONE +#include "openmc/mgxs_interface.h" +#include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// EnergyFilter implementation +//============================================================================== + +// Used to grab the rev_energy_bins array defined in Fortran. +extern "C" double* rev_energy_bins_ptr(); + +void +EnergyFilter::from_xml(pugi::xml_node node) +{ + bins_ = get_node_array(node, "bins"); + n_bins_ = bins_.size() - 1; + + // In MG mode, check if the filter bins match the transport bins. + // We can save tallying time if we know that the tally bins match the energy + // group structure. In that case, the matching bin index is simply the group + // (after flipping for the different ordering of the library and tallying + // systems). + if (!settings::run_CE) { + if (n_bins_ == num_energy_groups) { + matches_transport_groups_ = true; + double* rev_energy_bins = rev_energy_bins_ptr(); + for (auto i = 0; i < n_bins_ + 1; i++) { + if (rev_energy_bins[i] != bins_[i]) { + matches_transport_groups_ = false; + break; + } + } + } + } +} + +void +EnergyFilter::get_all_bins(const Particle* p, int estimator, FilterMatch& match) +const +{ + if (p->g != F90_NONE && matches_transport_groups_) { + if (estimator == ESTIMATOR_TRACKLENGTH) { + //TODO: off-by-one + match.bins_.push_back(num_energy_groups - p->g + 1); + } else { + //TODO: off-by-one + match.bins_.push_back(num_energy_groups - p->last_g + 1); + } + match.weights_.push_back(1.0); + + } else { + // Get the pre-collision energy of the particle. + auto E = p->last_E; + + // Bin the energy. + if (E >= bins_.front() && E <= bins_.back()) { + //TODO: off-by-one + auto bin = lower_bound_index(bins_.begin(), bins_.end(), E) + 1; + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } + } +} + +void +EnergyFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); +} + +std::string +EnergyFilter::text_label(int bin) const +{ + std::stringstream out; + out << "Incoming Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); +} + +//============================================================================== +// EnergyoutFilter implementation +//============================================================================== + +void +EnergyoutFilter::get_all_bins(const Particle* p, int estimator, + FilterMatch& match) const +{ + if (p->g != F90_NONE && matches_transport_groups_) { + match.bins_.push_back(num_energy_groups - p->g + 1); + match.weights_.push_back(1.0); + + } else { + if (p->E >= bins_.front() && p->E <= bins_.back()) { + //TODO: off-by-one + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->E) + 1; + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } + } +} + +std::string +EnergyoutFilter::text_label(int bin) const +{ + std::stringstream out; + out << "Outgoing Energy [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); +} + +}// namespace openmc diff --git a/src/tallies/tally_filter_energy.F90 b/src/tallies/tally_filter_energy.F90 index 70142156e2..338384f9d9 100644 --- a/src/tallies/tally_filter_energy.F90 +++ b/src/tallies/tally_filter_energy.F90 @@ -10,7 +10,7 @@ module tally_filter_energy use particle_header, only: Particle use settings, only: run_CE use string, only: to_str - use tally_filter_header + use tally_filter_cpp use xml_interface implicit none @@ -22,16 +22,13 @@ module tally_filter_energy ! ENERGYFILTER bins the incident neutron energy. !=============================================================================== - type, public, extends(TallyFilter) :: EnergyFilter + type, public, extends(CppTallyFilter) :: EnergyFilter real(8), allocatable :: bins(:) ! True if transport group number can be used directly to get bin number logical :: matches_transport_groups = .false. contains procedure :: from_xml => from_xml_energy - procedure :: get_all_bins => get_all_bins_energy - procedure :: to_statepoint => to_statepoint_energy - procedure :: text_label => text_label_energy end type EnergyFilter !=============================================================================== @@ -41,11 +38,6 @@ module tally_filter_energy !=============================================================================== type, public, extends(EnergyFilter) :: EnergyoutFilter - contains - ! Inherit from_xml from EnergyFilter - procedure :: get_all_bins => get_all_bins_energyout - procedure :: to_statepoint => to_statepoint_energyout - procedure :: text_label => text_label_energyout end type EnergyoutFilter contains @@ -60,6 +52,8 @@ contains integer :: n + call this % from_xml_cpp_inner(node) + n = node_word_count(node, "bins") ! Allocate and store bins @@ -81,116 +75,6 @@ contains end if end subroutine from_xml_energy - subroutine get_all_bins_energy(this, p, estimator, match) - class(EnergyFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: E - - n = this % n_bins - - if (p % g /= NONE .and. this % matches_transport_groups) then - if (estimator == ESTIMATOR_TRACKLENGTH) then - call match % bins_push_back(num_energy_groups - p % g + 1) - call match % weights_push_back(ONE) - else - call match % bins_push_back(num_energy_groups - p % last_g + 1) - call match % weights_push_back(ONE) - end if - - else - ! Pre-collision energy of particle - E = p % last_E - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, E) - if (bin /= NO_BIN_FOUND) then - call match % bins_push_back(bin) - call match % weights_push_back(ONE) - end if - end if - end subroutine get_all_bins_energy - - subroutine to_statepoint_energy(this, filter_group) - class(EnergyFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "energy") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_energy - - function text_label_energy(this, bin) result(label) - class(EnergyFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Incoming Energy [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_energy - -!=============================================================================== -! EnergyoutFilter methods -!=============================================================================== - - subroutine get_all_bins_energyout(this, p, estimator, match) - class(EnergyoutFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - - n = this % n_bins - - if (p % g /= NONE .and. this % matches_transport_groups) then - ! Tallies are ordered in increasing groups, group indices - ! however are the opposite, so switch - call match % bins_push_back(num_energy_groups - p % g + 1) - call match % weights_push_back(ONE) - - else - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, p % E) - if (bin /= NO_BIN_FOUND) then - call match % bins_push_back(bin) - call match % weights_push_back(ONE) - end if - end if - end subroutine get_all_bins_energyout - - subroutine to_statepoint_energyout(this, filter_group) - class(EnergyoutFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "energyout") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_energyout - - function text_label_energyout(this, bin) result(label) - class(EnergyoutFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - real(8) :: E0, E1 - - E0 = this % bins(bin) - E1 = this % bins(bin + 1) - label = "Outgoing Energy [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_energyout - !=============================================================================== ! C API FUNCTIONS !===============================================================================