diff --git a/include/openmc/tallies/tally_filter_energyfunc.h b/include/openmc/tallies/tally_filter_energyfunc.h new file mode 100644 index 0000000000..f180a57a33 --- /dev/null +++ b/include/openmc/tallies/tally_filter_energyfunc.h @@ -0,0 +1,89 @@ +#ifndef OPENMC_TALLY_FILTER_ENERGYFUNC_H +#define OPENMC_TALLY_FILTER_ENERGYFUNC_H + +#include // for setprecision +#include // for scientific +#include +#include + +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/settings.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class EnergyFunctionFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "energyfunction";} + + EnergyFunctionFilter() + : TallyFilter {} + { + n_bins_ = 1; + } + + virtual ~EnergyFunctionFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + if (!settings::run_CE) + fatal_error("EnergyFunction filters are only supported for " + "continuous-energy transport calculations"); + + if (!check_for_node(node, "energy")) + fatal_error("Energy grid not specified for EnergyFunction filter."); + + energy_ = get_node_array(node, "energy"); + + if (!check_for_node(node, "y")) + fatal_error("y values not specified for EnergyFunction filter."); + + y_ = get_node_array(node, "y"); + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + if (p->last_E >= energy_.front() && p->last_E <= energy_.back()) { + // Search for the incoming energy bin. + auto i = lower_bound_index(energy_.begin(), energy_.end(), p->last_E); + + // Compute the interpolation factor between the nearest bins. + double f = (p->last_E - energy_[i]) / (energy_[i+1] - energy_[i]); + + // Interpolate on the lin-lin grid. + match.bins.push_back(1); + match.weights.push_back((1-f) * y_[i] + f * y_[i+1]); + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "energy", energy_); + write_dataset(filter_group, "y", y_); + } + + virtual std::string + text_label(int bin) const override + { + std::stringstream out; + out << std::scientific << std::setprecision(1) + << "Energy Function f" + << "([ " << energy_.front() << ", ..., " << energy_.back() << "]) = " + << "[" << y_.front() << ", ..., " << y_.back() << "]"; + return out.str(); + } + + std::vector energy_; + std::vector y_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_ENERGYFUNC_H diff --git a/src/tallies/tally_filter.cpp b/src/tallies/tally_filter.cpp index c4ce32e7e3..cfc54a6a24 100644 --- a/src/tallies/tally_filter.cpp +++ b/src/tallies/tally_filter.cpp @@ -8,6 +8,7 @@ #include "openmc/tallies/tally_filter_cellborn.h" #include "openmc/tallies/tally_filter_cellfrom.h" #include "openmc/tallies/tally_filter_distribcell.h" +#include "openmc/tallies/tally_filter_energyfunc.h" #include "openmc/tallies/tally_filter_legendre.h" #include "openmc/tallies/tally_filter_material.h" #include "openmc/tallies/tally_filter_mesh.h" @@ -98,6 +99,8 @@ extern "C" { tally_filters.push_back(new CellFromFilter()); } else if (type_ == "distribcell") { tally_filters.push_back(new DistribcellFilter()); + } else if (type_ == "energyfunction") { + tally_filters.push_back(new EnergyFunctionFilter()); } else if (type_ == "legendre") { tally_filters.push_back(new LegendreFilter()); } else if (type_ == "material") { diff --git a/src/tallies/tally_filter_energyfunc.F90 b/src/tallies/tally_filter_energyfunc.F90 index 0581fdacc8..de41b0c59b 100644 --- a/src/tallies/tally_filter_energyfunc.F90 +++ b/src/tallies/tally_filter_energyfunc.F90 @@ -1,124 +1,15 @@ module tally_filter_energyfunc - use, intrinsic :: ISO_C_BINDING - - use algorithm, only: binary_search - use constants - use error, only: fatal_error - use hdf5_interface - use particle_header, only: Particle - use settings, only: run_CE - use string, only: to_str use tally_filter_header - use xml_interface implicit none - private !=============================================================================== ! EnergyFunctionFilter multiplies tally scores by an arbitrary function of ! incident energy described by a piecewise linear-linear interpolation. !=============================================================================== - type, public, extends(TallyFilter) :: EnergyFunctionFilter - real(8), allocatable :: energy(:) - real(8), allocatable :: y(:) - - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_energyfunction - procedure :: to_statepoint => to_statepoint_energyfunction - procedure :: text_label => text_label_energyfunction + type, extends(CppTallyFilter) :: EnergyFunctionFilter end type EnergyFunctionFilter -contains - - subroutine from_xml(this, node) - class(EnergyFunctionFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: n - - this % n_bins = 1 - ! Make sure this is continuous-energy mode. - if (.not. run_CE) then - call fatal_error("EnergyFunction filters are only supported for & - &continuous-energy transport calculations") - end if - - ! Allocate and store energy grid. - if (.not. check_for_node(node, "energy")) then - call fatal_error("Energy grid not specified for EnergyFunction & - &filter.") - end if - n = node_word_count(node, "energy") - allocate(this % energy(n)) - call get_node_array(node, "energy", this % energy) - - ! Allocate and store interpolant values. - if (.not. check_for_node(node, "y")) then - call fatal_error("y values not specified for EnergyFunction & - &filter.") - end if - n = node_word_count(node, "y") - allocate(this % y(n)) - call get_node_array(node, "y", this % y) - end subroutine from_xml - - subroutine get_all_bins_energyfunction(this, p, estimator, match) - class(EnergyFunctionFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n, indx - real(8) :: E, f, weight - - select type(this) - type is (EnergyFunctionFilter) - n = size(this % energy) - - ! Get pre-collision energy of particle - E = p % last_E - - ! Search to find incoming energy bin. - indx = binary_search(this % energy, n, E) - - ! Compute an interpolation factor between nearest bins. - f = (E - this % energy(indx)) & - / (this % energy(indx+1) - this % energy(indx)) - - ! Interpolate on the lin-lin grid. - call match % bins_push_back(1) - weight = (ONE - f) * this % y(indx) + f * this % y(indx+1) - call match % weights_push_back(weight) - end select - end subroutine get_all_bins_energyfunction - - subroutine to_statepoint_energyfunction(this, filter_group) - class(EnergyFunctionFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - select type(this) - type is (EnergyFunctionFilter) - call write_dataset(filter_group, "type", "energyfunction") - call write_dataset(filter_group, "energy", this % energy) - call write_dataset(filter_group, "y", this % y) - end select - end subroutine to_statepoint_energyfunction - - function text_label_energyfunction(this, bin) result(label) - class(EnergyFunctionFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - select type(this) - type is (EnergyFunctionFilter) - write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") & - "Energy Function f([", this % energy(1), ", ..., ", & - this % energy(size(this % energy)), "]) = [", this % y(1), & - ", ..., ", this % y(size(this % y)), "]" - end select - end function text_label_energyfunction - end module tally_filter_energyfunc