From 00cb9e93ae1c062e792703b4c65cf677a1797a22 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sat, 13 Oct 2018 19:53:57 -0400 Subject: [PATCH] Move azimuthal, mu, polar filters to C++ --- .../openmc/tallies/tally_filter_azimuthal.h | 86 ++++++++++++++ include/openmc/tallies/tally_filter_mu.h | 79 +++++++++++++ include/openmc/tallies/tally_filter_polar.h | 86 ++++++++++++++ src/tallies/tally_filter.cpp | 11 +- src/tallies/tally_filter_azimuthal.F90 | 107 +----------------- src/tallies/tally_filter_distribcell.F90 | 1 - src/tallies/tally_filter_mu.F90 | 97 +--------------- src/tallies/tally_filter_polar.F90 | 105 +---------------- 8 files changed, 264 insertions(+), 308 deletions(-) create mode 100644 include/openmc/tallies/tally_filter_azimuthal.h create mode 100644 include/openmc/tallies/tally_filter_mu.h create mode 100644 include/openmc/tallies/tally_filter_polar.h diff --git a/include/openmc/tallies/tally_filter_azimuthal.h b/include/openmc/tallies/tally_filter_azimuthal.h new file mode 100644 index 0000000000..182e32bc07 --- /dev/null +++ b/include/openmc/tallies/tally_filter_azimuthal.h @@ -0,0 +1,86 @@ +#ifndef OPENMC_TALLY_FILTER_AZIMUTHAL_H +#define OPENMC_TALLY_FILTER_AZIMUTHAL_H + +#include +#include +#include + +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class AzimuthalFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "azimuthal";} + + virtual ~AzimuthalFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + auto bins = get_node_array(node, "bins"); + + if (bins.size() > 1) { + bins_ = bins; + + } else { + // Allow a user to input a lone number which will mean that you subdivide + // [-pi,pi) evenly with the input being the number of bins + + int n_angle = bins[0]; + + if (n_angle <= 1) fatal_error("Number of bins for azimuthal filter must " + "be greater than 1."); + + double d_angle = 2.0 * PI / n_angle; + bins_.resize(n_angle + 1); + for (int i = 0; i < n_angle; i++) bins_[i] = -PI + i * d_angle; + bins_[n_angle] = PI; + } + + n_bins_ = bins_.size() - 1; + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + double phi; + if (estimator == ESTIMATOR_TRACKLENGTH) { + phi = std::atan2(p->coord[0].uvw[1], p->coord[0].uvw[0]); + } else { + phi = std::atan2(p->last_uvw[1], p->last_uvw[0]); + } + + if (phi >= bins_[0] && phi <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), phi) + 1; + match.bins.push_back(bin); + match.weights.push_back(1); + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); + } + + virtual std::string + text_label(int bin) const override + { + std::stringstream out; + out << "Azimuthal Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); + } + +protected: + std::vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_AZIMUTHAL_H diff --git a/include/openmc/tallies/tally_filter_mu.h b/include/openmc/tallies/tally_filter_mu.h new file mode 100644 index 0000000000..e76d56e43b --- /dev/null +++ b/include/openmc/tallies/tally_filter_mu.h @@ -0,0 +1,79 @@ +#ifndef OPENMC_TALLY_FILTER_MU_H +#define OPENMC_TALLY_FILTER_MU_H + +#include +#include +#include + +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class MuFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "mu";} + + virtual ~MuFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + auto bins = get_node_array(node, "bins"); + + if (bins.size() > 1) { + bins_ = bins; + + } else { + // Allow a user to input a lone number which will mean that you subdivide + // [-1,1) evenly with the input being the number of bins + + int n_angle = bins[0]; + + if (n_angle <= 1) fatal_error("Number of bins for mu filter must " + "be greater than 1."); + + double d_angle = 2.0 / n_angle; + bins_.resize(n_angle + 1); + for (int i = 0; i < n_angle; i++) bins_[i] = -1 + i * d_angle; + bins_[n_angle] = PI; + } + + n_bins_ = bins_.size() - 1; + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + if (p->mu >= bins_[0] && p->mu <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), p->mu) + 1; + match.bins.push_back(bin); + match.weights.push_back(1); + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); + } + + virtual std::string + text_label(int bin) const override + { + std::stringstream out; + out << "Change-in-Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); + } + +protected: + std::vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_MU_H diff --git a/include/openmc/tallies/tally_filter_polar.h b/include/openmc/tallies/tally_filter_polar.h new file mode 100644 index 0000000000..d9f003330b --- /dev/null +++ b/include/openmc/tallies/tally_filter_polar.h @@ -0,0 +1,86 @@ +#ifndef OPENMC_TALLY_FILTER_POLAR_H +#define OPENMC_TALLY_FILTER_POLAR_H + +#include +#include +#include + +#include "openmc/error.h" +#include "openmc/search.h" +#include "openmc/tallies/tally_filter.h" + + +namespace openmc { + +class PolarFilter : public TallyFilter +{ +public: + virtual std::string type() const override {return "polar";} + + virtual ~PolarFilter() override = default; + + virtual void + from_xml(pugi::xml_node node) override + { + auto bins = get_node_array(node, "bins"); + + if (bins.size() > 1) { + bins_ = bins; + + } else { + // Allow a user to input a lone number which will mean that you subdivide + // [0,pi] evenly with the input being the number of bins + + int n_angle = bins[0]; + + if (n_angle <= 1) fatal_error("Number of bins for polar filter must " + "be greater than 1."); + + double d_angle = PI / n_angle; + bins_.resize(n_angle + 1); + for (int i = 0; i < n_angle; i++) bins_[i] = i * d_angle; + bins_[n_angle] = PI; + } + + n_bins_ = bins_.size() - 1; + } + + virtual void + get_all_bins(Particle* p, int estimator, TallyFilterMatch& match) + const override + { + double theta; + if (estimator == ESTIMATOR_TRACKLENGTH) { + theta = std::acos(p->coord[0].uvw[2]); + } else { + theta = std::acos(p->last_uvw[2]); + } + + if (theta >= bins_[0] && theta <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), theta) + 1; + match.bins.push_back(bin); + match.weights.push_back(1); + } + } + + virtual void + to_statepoint(hid_t filter_group) const override + { + TallyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); + } + + virtual std::string + text_label(int bin) const override + { + std::stringstream out; + out << "Polar Angle [" << bins_[bin-1] << ", " << bins_[bin] << ")"; + return out.str(); + } + +protected: + std::vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLY_FILTER_POLAR_H diff --git a/src/tallies/tally_filter.cpp b/src/tallies/tally_filter.cpp index b0f54dd1ee..05de1b46c1 100644 --- a/src/tallies/tally_filter.cpp +++ b/src/tallies/tally_filter.cpp @@ -3,12 +3,15 @@ #include #include "openmc/constants.h" // for MAX_LINE_LEN; +#include "openmc/tallies/tally_filter_azimuthal.h" #include "openmc/tallies/tally_filter_cell.h" #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_mesh.h" #include "openmc/tallies/tally_filter_meshsurface.h" +#include "openmc/tallies/tally_filter_mu.h" +#include "openmc/tallies/tally_filter_polar.h" namespace openmc { @@ -80,7 +83,9 @@ extern "C" { allocate_filter(const char* type) { std::string type_ {type}; - if (type_ == "cell") { + if (type_ == "azimuthal") { + tally_filters.push_back(new AzimuthalFilter()); + } else if (type_ == "cell") { tally_filters.push_back(new CellFilter()); } else if (type_ == "cellborn") { tally_filters.push_back(new CellbornFilter()); @@ -92,6 +97,10 @@ extern "C" { tally_filters.push_back(new MeshFilter()); } else if (type_ == "meshsurface") { tally_filters.push_back(new MeshSurfaceFilter()); + } else if (type_ == "mu") { + tally_filters.push_back(new MuFilter()); + } else if (type_ == "polar") { + tally_filters.push_back(new PolarFilter()); } else { return nullptr; } diff --git a/src/tallies/tally_filter_azimuthal.F90 b/src/tallies/tally_filter_azimuthal.F90 index 5dc5f9e2ce..69513c26d9 100644 --- a/src/tallies/tally_filter_azimuthal.F90 +++ b/src/tallies/tally_filter_azimuthal.F90 @@ -1,120 +1,15 @@ module tally_filter_azimuthal - 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 string, only: to_str use tally_filter_header - use xml_interface implicit none - private !=============================================================================== ! AZIMUTHALFILTER bins the incident neutron azimuthal angle (relative to the ! global xy-plane). !=============================================================================== - type, public, extends(TallyFilter) :: AzimuthalFilter - real(8), allocatable :: bins(:) - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_azimuthal - procedure :: to_statepoint => to_statepoint_azimuthal - procedure :: text_label => text_label_azimuthal + type, extends(CppTallyFilter) :: AzimuthalFilter end type AzimuthalFilter -contains - - subroutine from_xml(this, node) - class(AzimuthalFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: n - integer :: n_angle - real(8) :: d_angle - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - 1 - allocate(this % bins(n)) - call get_node_array(node, "bins", this % bins) - - ! Allow a user to input a lone number which will mean that you - ! subdivide [-pi,pi) evenly with the input being the number of - ! bins - if (n == 1) then - n_angle = int(this % bins(1)) - if (n_angle > 1) then - this % n_bins = n_angle - d_angle = TWO * PI / n_angle - deallocate(this % bins) - allocate(this % bins(n_angle + 1)) - do i = 1, n_angle - this % bins(i) = -PI + (i - 1) * d_angle - end do - this % bins(n_angle + 1) = PI - else - call fatal_error("Number of bins for azimuthal filter must be& - & greater than 1.") - end if - end if - end subroutine from_xml - - subroutine get_all_bins_azimuthal(this, p, estimator, match) - class(AzimuthalFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: phi - - n = this % n_bins - - ! Make sure the correct direction vector is used. - if (estimator == ESTIMATOR_TRACKLENGTH) then - phi = atan2(p % coord(1) % uvw(2), p % coord(1) % uvw(1)) - else - phi = atan2(p % last_uvw(2), p % last_uvw(1)) - end if - - ! Search to find azimuthal angle bin. - bin = binary_search(this % bins, n + 1, phi) - if (bin /= NO_BIN_FOUND) then - call match % bins_push_back(bin) - call match % weights_push_back(ONE) - end if - - end subroutine get_all_bins_azimuthal - - subroutine to_statepoint_azimuthal(this, filter_group) - class(AzimuthalFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "azimuthal") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_azimuthal - - function text_label_azimuthal(this, bin) result(label) - class(AzimuthalFilter), 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 = "Azimuthal Angle [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_azimuthal - end module tally_filter_azimuthal diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index 4a4b286852..8b25e27471 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -11,7 +11,6 @@ module tally_filter_distribcell !=============================================================================== type, public, extends(CppTallyFilter) :: DistribcellFilter - integer :: cell contains procedure :: initialize => initialize_distribcell end type DistribcellFilter diff --git a/src/tallies/tally_filter_mu.F90 b/src/tallies/tally_filter_mu.F90 index ecd0a75af5..8f8d14a7d5 100644 --- a/src/tallies/tally_filter_mu.F90 +++ b/src/tallies/tally_filter_mu.F90 @@ -1,110 +1,15 @@ module tally_filter_mu - use, intrinsic :: ISO_C_BINDING - - use algorithm, only: binary_search - use constants, only: ONE, TWO, MAX_LINE_LEN, NO_BIN_FOUND - use error, only: fatal_error - use hdf5_interface - use particle_header, only: Particle - use string, only: to_str use tally_filter_header - use xml_interface implicit none - private !=============================================================================== ! MUFILTER bins the incoming-outgoing direction cosine. This is only used for ! scatter reactions. !=============================================================================== - type, public, extends(TallyFilter) :: MuFilter - real(8), allocatable :: bins(:) - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_mu - procedure :: to_statepoint => to_statepoint_mu - procedure :: text_label => text_label_mu + type, extends(CppTallyFilter) :: MuFilter end type MuFilter -contains - - subroutine from_xml(this, node) - class(MuFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: n_angle - integer :: n - real(8) :: d_angle - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - 1 - allocate(this % bins(n)) - call get_node_array(node, "bins", this % bins) - - ! Allow a user to input a lone number which will mean that you - ! subdivide [-1,1] evenly with the input being the number of bins - if (n == 1) then - n_angle = int(this % bins(1)) - if (n_angle > 1) then - this % n_bins = n_angle - d_angle = TWO / n_angle - deallocate(this % bins) - allocate(this % bins(n_angle + 1)) - do i = 1, n_angle - this % bins(i) = -ONE + (i - 1) * d_angle - end do - this % bins(n_angle + 1) = ONE - else - call fatal_error("Number of bins for mu filter must be& - & greater than 1.") - end if - end if - end subroutine from_xml - - subroutine get_all_bins_mu(this, p, estimator, match) - class(MuFilter), 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 - - ! Search to find incoming energy bin. - bin = binary_search(this % bins, n + 1, p % mu) - if (bin /= NO_BIN_FOUND) then - call match % bins_push_back(bin) - call match % weights_push_back(ONE) - end if - end subroutine get_all_bins_mu - - subroutine to_statepoint_mu(this, filter_group) - class(MuFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "mu") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_mu - - function text_label_mu(this, bin) result(label) - class(MuFilter), 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 = "Change-in-Angle [" // trim(to_str(E0)) // ", " & - // trim(to_str(E1)) // ")" - end function text_label_mu - end module tally_filter_mu diff --git a/src/tallies/tally_filter_polar.F90 b/src/tallies/tally_filter_polar.F90 index 46627230ac..957f0f51ce 100644 --- a/src/tallies/tally_filter_polar.F90 +++ b/src/tallies/tally_filter_polar.F90 @@ -1,118 +1,15 @@ module tally_filter_polar - 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 string, only: to_str use tally_filter_header - use xml_interface implicit none - private !=============================================================================== ! POLARFILTER bins the incident neutron polar angle (relative to the global ! z-axis). !=============================================================================== - type, public, extends(TallyFilter) :: PolarFilter - real(8), allocatable :: bins(:) - contains - procedure :: from_xml - procedure :: get_all_bins => get_all_bins_polar - procedure :: to_statepoint => to_statepoint_polar - procedure :: text_label => text_label_polar + type, extends(CppTallyFilter) :: PolarFilter end type PolarFilter -contains - - subroutine from_xml(this, node) - class(PolarFilter), intent(inout) :: this - type(XMLNode), intent(in) :: node - - integer :: i - integer :: n_angle - integer :: n - real(8) :: d_angle - - n = node_word_count(node, "bins") - - ! Allocate and store bins - this % n_bins = n - 1 - allocate(this % bins(n)) - call get_node_array(node, "bins", this % bins) - - ! Allow a user to input a lone number which will mean that you - ! subdivide [0,pi] evenly with the input being the number of bins - if (n == 1) then - n_angle = int(this % bins(1)) - if (n_angle > 1) then - this % n_bins = n_angle - d_angle = PI / real(n_angle,8) - deallocate(this % bins) - allocate(this % bins(n_angle + 1)) - do i = 1, n_angle - this % bins(i) = (i - 1) * d_angle - end do - this % bins(n_angle + 1) = PI - else - call fatal_error("Number of bins for polar filter must be& - & greater than 1.") - end if - end if - end subroutine from_xml - - subroutine get_all_bins_polar(this, p, estimator, match) - class(PolarFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: n - integer :: bin - real(8) :: theta - - n = this % n_bins - - ! Make sure the correct direction vector is used. - if (estimator == ESTIMATOR_TRACKLENGTH) then - theta = acos(p % coord(1) % uvw(3)) - else - theta = acos(p % last_uvw(3)) - end if - - ! Search to find polar angle bin. - bin = binary_search(this % bins, n + 1, theta) - if (bin /= NO_BIN_FOUND) then - call match % bins_push_back(bin) - call match % weights_push_back(ONE) - end if - end subroutine get_all_bins_polar - - subroutine to_statepoint_polar(this, filter_group) - class(PolarFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "polar") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", this % bins) - end subroutine to_statepoint_polar - - function text_label_polar(this, bin) result(label) - class(PolarFilter), 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 = "Polar Angle [" // trim(to_str(E0)) // ", " // trim(to_str(E1)) & - // ")" - end function text_label_polar - end module tally_filter_polar