Move azimuthal, mu, polar filters to C++

include/openmc/tallies/tally_filter_azimuthal.h

@@ -0,0 +1,86 @@
+#ifndef OPENMC_TALLY_FILTER_AZIMUTHAL_H
+#define OPENMC_TALLY_FILTER_AZIMUTHAL_H
+
+#include <sstream>
+#include <cmath>
+#include <vector>
+
+#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<double>(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<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLY_FILTER_AZIMUTHAL_H

include/openmc/tallies/tally_filter_mu.h

@@ -0,0 +1,79 @@
+#ifndef OPENMC_TALLY_FILTER_MU_H
+#define OPENMC_TALLY_FILTER_MU_H
+
+#include <sstream>
+#include <cmath>
+#include <vector>
+
+#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<double>(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<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLY_FILTER_MU_H

include/openmc/tallies/tally_filter_polar.h

@@ -0,0 +1,86 @@
+#ifndef OPENMC_TALLY_FILTER_POLAR_H
+#define OPENMC_TALLY_FILTER_POLAR_H
+
+#include <sstream>
+#include <cmath>
+#include <vector>
+
+#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<double>(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<double> bins_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLY_FILTER_POLAR_H

src/tallies/tally_filter.cpp

@@ -3,12 +3,15 @@
 #include <string>
 
 #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;
     }

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

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

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

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