Move Legendre and sph harmonics filters to C++

include/openmc/tallies/tally_filter_legendre.h

@@ -0,0 +1,57 @@
+#ifndef OPENMC_TALLY_FILTER_LEGENDRE_H
+#define OPENMC_TALLY_FILTER_LEGENDRE_H
+
+#include <string>
+
+#include "openmc/cell.h"
+#include "openmc/error.h"
+#include "openmc/math_functions.h"
+#include "openmc/tallies/tally_filter.h"
+
+
+namespace openmc {
+
+class LegendreFilter : public TallyFilter
+{
+public:
+  virtual std::string type() const override {return "legendre";}
+
+  virtual ~LegendreFilter() override = default;
+
+  virtual void
+  from_xml(pugi::xml_node node) override
+  {
+    order_ = std::stoi(get_node_value(node, "order"));
+    n_bins_ = order_ + 1;
+  }
+
+  virtual void
+  get_all_bins(Particle* p, int estimator, TallyFilterMatch& match)
+  const override
+  {
+    double wgt[n_bins_];
+    calc_pn_c(order_, p->mu, wgt);
+    for (int i = 0; i < n_bins_; i++) {
+      match.bins.push_back(i + 1);
+      match.weights.push_back(wgt[i]);
+    }
+  }
+
+  virtual void
+  to_statepoint(hid_t filter_group) const override
+  {
+    TallyFilter::to_statepoint(filter_group);
+    write_dataset(filter_group, "order", order_);
+  }
+
+  virtual std::string
+  text_label(int bin) const override
+  {
+    return "Legendre expansion, P" + std::to_string(bin - 1);
+  }
+
+  int order_;
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLY_FILTER_LEGENDRE_H

include/openmc/tallies/tally_filter_sph_harm.h

@@ -0,0 +1,106 @@
+#ifndef OPENMC_TALLY_FILTER_SPHERICAL_HARMONICS_H
+#define OPENMC_TALLY_FILTER_SPHERICAL_HARMONICS_H
+
+#include <string>
+
+#include "openmc/error.h"
+#include "openmc/tallies/tally_filter.h"
+
+
+namespace openmc {
+
+//TODO: those integer values are not needed when Fortran interop is removed
+enum class SphericalHarmonicsCosine {
+  scatter = 1, particle = 2
+};
+
+class SphericalHarmonicsFilter : public TallyFilter
+{
+public:
+  virtual std::string type() const override {return "sphericalharmonics";}
+
+  virtual ~SphericalHarmonicsFilter() override = default;
+
+  virtual void
+  from_xml(pugi::xml_node node) override
+  {
+    order_ = std::stoi(get_node_value(node, "order"));
+    n_bins_ = (order_ + 1) * (order_ + 1);
+
+    if (check_for_node(node, "cosine")) {
+      auto cos = get_node_value(node, "cosine", true);
+      std::cout << cos << "\n";
+      std::cout << (cos == "scatter") << " " << (cos == "particle") << "\n";
+      if (cos == "scatter") {
+        cosine_ = SphericalHarmonicsCosine::scatter;
+      } else if (cos == "particle") {
+        cosine_ = SphericalHarmonicsCosine::particle;
+      } else {
+        std::stringstream err_msg;
+        err_msg << "Unrecognized cosine type, \"" << cos
+                << "\" in spherical harmonics filter";
+        fatal_error(err_msg);
+      }
+    }
+  }
+
+  virtual void
+  get_all_bins(Particle* p, int estimator, TallyFilterMatch& match)
+  const override
+  {
+    // Determine cosine term for scatter expansion if necessary
+    double wgt[order_ + 1];
+    if (cosine_ == SphericalHarmonicsCosine::scatter) {
+      calc_pn_c(order_, p->mu, wgt);
+    } else {
+      for (int i = 0; i < order_ + 1; i++) wgt[i] = 1;
+    }
+
+    // Find the Rn,m values
+    double rn[n_bins_];
+    calc_rn_c(order_, p->last_uvw, rn);
+
+    int j = 0;
+    for (int n = 0; n < order_ + 1; n++) {
+      // Calculate n-th order spherical harmonics for (u,v,w)
+      int num_nm = 2*n + 1;
+
+      // Append the matching (bin,weight) for each moment
+      for (int i = 0; i < num_nm; i++) {
+        match.weights.push_back(wgt[n] * rn[j]);
+        match.bins.push_back(++j);
+      }
+    }
+  }
+
+  virtual void
+  to_statepoint(hid_t filter_group) const override
+  {
+    TallyFilter::to_statepoint(filter_group);
+    write_dataset(filter_group, "order", order_);
+    if (cosine_ == SphericalHarmonicsCosine::scatter) {
+      write_dataset(filter_group, "cosine", "scatter");
+    } else {
+      write_dataset(filter_group, "cosine", "particle");
+    }
+  }
+
+  virtual std::string
+  text_label(int bin) const override
+  {
+    std::stringstream out;
+    for (int n = 0; n < order_ + 1; n++) {
+      if (bin <= (n + 1) * (n + 1)) {
+        int m = (bin - n*n - 1) - n;
+        out << "Spherical harmonic expansion, Y" << n << "," << m;
+        return out.str();
+      }
+    }
+  }
+
+  int order_;
+  SphericalHarmonicsCosine cosine_ {SphericalHarmonicsCosine::particle};
+};
+
+} // namespace openmc
+#endif // OPENMC_TALLY_FILTER_SPHERICAL_HARMONICS_H

src/tallies/tally_filter.cpp

@@ -8,11 +8,13 @@
 #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_legendre.h"
 #include "openmc/tallies/tally_filter_material.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"
+#include "openmc/tallies/tally_filter_sph_harm.h"
 #include "openmc/tallies/tally_filter_surface.h"
 #include "openmc/tallies/tally_filter_universe.h"
 
@@ -96,6 +98,8 @@ extern "C" {
       tally_filters.push_back(new CellFromFilter());
     } else if (type_ == "distribcell") {
       tally_filters.push_back(new DistribcellFilter());
+    } else if (type_ == "legendre") {
+      tally_filters.push_back(new LegendreFilter());
     } else if (type_ == "material") {
       tally_filters.push_back(new MaterialFilter());
     } else if (type_ == "mesh") {
@@ -108,6 +112,8 @@ extern "C" {
       tally_filters.push_back(new PolarFilter());
     } else if (type_ == "surface") {
       tally_filters.push_back(new SurfaceFilter());
+    } else if (type_ == "sphericalharmonics") {
+      tally_filters.push_back(new SphericalHarmonicsFilter());
     } else if (type_ == "universe") {
       tally_filters.push_back(new UniverseFilter());
     } else {
@@ -149,6 +155,12 @@ extern "C" {
     *n = filt->cells_.size();
   }
 
+  int legendre_filter_get_order(LegendreFilter* filt)
+  {return filt->order_;}
+
+  void legendre_filter_set_order(LegendreFilter* filt, int order)
+  {filt->order_ = order;}
+
   void
   material_filter_get_bins(MaterialFilter* filt, int32_t** bins, int32_t* n)
   {
@@ -164,7 +176,7 @@ extern "C" {
     for (int i = 0; i < n; i++) filt->materials_[i] = bins[i];
     filt->n_bins_ = filt->materials_.size();
     filt->map_.clear();
-    for (int i = 0; i < n; i++) filt->map_[bins[i]] = i;
+    for (int i = 0; i < n; i++) filt->map_[filt->materials_[i]] = i;
   }
 
   int mesh_filter_get_mesh(MeshFilter* filt) {return filt->mesh_;}
@@ -186,6 +198,22 @@ extern "C" {
     filt->n_bins_ = 4 * meshes[mesh]->n_dimension_;
     for (auto dim : meshes[mesh]->shape_) filt->n_bins_ *= dim;
   }
+
+  int sphharm_filter_get_order(SphericalHarmonicsFilter* filt)
+  {return filt->order_;}
+
+  int sphharm_filter_get_cosine(SphericalHarmonicsFilter* filt)
+  {return static_cast<int>(filt->cosine_);}
+
+  void
+  sphharm_filter_set_order(SphericalHarmonicsFilter* filt, int order)
+  {
+    filt->order_ = order;
+    filt->n_bins_ = (order + 1) * (order + 1);
+  }
+
+  void sphharm_filter_set_cosine(SphericalHarmonicsFilter* filt, int cosine)
+  {filt->cosine_ = static_cast<SphericalHarmonicsCosine>(cosine);}
 }
 
 } // namespace openmc

src/tallies/tally_filter_legendre.F90

@@ -2,14 +2,8 @@ module tally_filter_legendre
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
   use error
-  use hdf5_interface
-  use math,                only: calc_pn
-  use particle_header,     only: Particle
-  use string,              only: to_str
   use tally_filter_header
-  use xml_interface
 
   implicit none
   private
@@ -20,63 +14,11 @@ module tally_filter_legendre
 ! LEGENDREFILTER gives Legendre moments of the change in scattering angle
 !===============================================================================
 
-  type, public, extends(TallyFilter) :: LegendreFilter
-    integer(C_INT) :: order
-  contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+  type, public, extends(CppTallyFilter) :: LegendreFilter
   end type LegendreFilter
 
 contains
 
-!===============================================================================
-! LegendreFilter methods
-!===============================================================================
-
-  subroutine from_xml(this, node)
-    class(LegendreFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    ! Get specified order
-    call get_node_value(node, "order", this % order)
-    this % n_bins = this % order + 1
-  end subroutine from_xml
-
-  subroutine get_all_bins(this, p, estimator, match)
-    class(LegendreFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i
-    real(C_DOUBLE) :: wgt(this % n_bins)
-
-    call calc_pn(this % order, p % mu, wgt)
-    do i = 1, this % n_bins
-      call match % bins_push_back(i)
-      call match % weights_push_back(wgt(i))
-    end do
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
-    class(LegendreFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "legendre")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(LegendreFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    label = "Legendre expansion, P" // trim(to_str(bin - 1))
-  end function text_label
-
 !===============================================================================
 !                               C API FUNCTIONS
 !===============================================================================
@@ -87,11 +29,19 @@ contains
     integer(C_INT),     intent(out) :: order
     integer(C_INT) :: err
 
+    interface
+      function legendre_filter_get_order(filt) result(order) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: filt
+        integer(C_INT)     :: order
+      end function
+    end interface
+
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (LegendreFilter)
-        order = f % order
+        order = legendre_filter_get_order(f % ptr)
       class default
         err = E_INVALID_TYPE
         call set_errmsg("Tried to get order on a non-expansion filter.")
@@ -106,12 +56,20 @@ contains
     integer(C_INT),     value :: order
     integer(C_INT) :: err
 
+    interface
+      subroutine legendre_filter_set_order(filt, order) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),    value :: filt
+        integer(C_INT), value :: order
+      end subroutine legendre_filter_set_order
+    end interface
+
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (LegendreFilter)
-        f % order = order
-        f % n_bins = order + 1
+        call legendre_filter_set_order(f % ptr, order)
+        f % n_bins = f % n_bins_cpp()
       class default
         err = E_INVALID_TYPE
         call set_errmsg("Tried to set order on a non-expansion filter.")

src/tallies/tally_filter_sph_harm.F90

@@ -2,14 +2,9 @@ module tally_filter_sph_harm
 
   use, intrinsic :: ISO_C_BINDING
 
-  use constants
   use error
-  use hdf5_interface
-  use math,                only: calc_pn, calc_rn
-  use particle_header,     only: Particle
   use string,              only: to_str, to_lower, to_f_string
   use tally_filter_header
-  use xml_interface
 
   implicit none
   private
@@ -26,112 +21,25 @@ module tally_filter_sph_harm
 ! tally score
 !===============================================================================
 
-  type, public, extends(TallyFilter) :: SphericalHarmonicsFilter
-    integer :: order
-    integer :: cosine
+  type, public, extends(CppTallyFilter) :: SphericalHarmonicsFilter
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+    procedure :: cosine
   end type SphericalHarmonicsFilter
 
 contains
 
-!===============================================================================
-! SphericalHarmonicsFilter methods
-!===============================================================================
-
-  subroutine from_xml(this, node)
-    class(SphericalHarmonicsFilter), intent(inout) :: this
-    type(XMLNode), intent(in) :: node
-
-    character(MAX_WORD_LEN) :: temp_str
-
-    ! Get specified order
-    call get_node_value(node, "order", this % order)
-    this % n_bins = (this % order + 1)**2
-
-    ! Determine how cosine term is to be treated
-    if (check_for_node(node, "cosine")) then
-      call get_node_value(node, "cosine", temp_str)
-      select case (to_lower(temp_str))
-      case ('scatter')
-        this % cosine = COSINE_SCATTER
-      case ('particle')
-        this % cosine = COSINE_PARTICLE
-      end select
-    else
-      this % cosine = COSINE_PARTICLE
-    end if
-  end subroutine from_xml
-
-  subroutine get_all_bins(this, p, estimator, match)
-    class(SphericalHarmonicsFilter), intent(in)  :: this
-    type(Particle),      intent(in)  :: p
-    integer,             intent(in)  :: estimator
-    type(TallyFilterMatch),   intent(inout) :: match
-
-    integer :: i, j, n
-    integer :: num_nm
-    real(C_DOUBLE) :: wgt(this % order + 1)
-    real(C_DOUBLE) :: rn(this % n_bins)
-
-    ! Determine cosine term for scatter expansion if necessary
-    if (this % cosine == COSINE_SCATTER) then
-      call calc_pn(this % order, p % mu, wgt)
-    else
-      wgt = ONE
-    end if
-
-    ! Find the Rn,m values
-    call calc_rn(this % order, p % last_uvw, rn)
-
-    j = 0
-    do n = 0, this % order
-      ! Calculate n-th order spherical harmonics for (u,v,w)
-      num_nm = 2*n + 1
-
-      ! Append matching (bin,weight) for each moment
-      do i = 1, num_nm
-        j = j + 1
-        call match % bins_push_back(j)
-        call match % weights_push_back(wgt(n + 1) * rn(j))
-      end do
-    end do
-
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
-    class(SphericalHarmonicsFilter), intent(in) :: this
-    integer(HID_T),      intent(in) :: filter_group
-
-    call write_dataset(filter_group, "type", "sphericalharmonics")
-    call write_dataset(filter_group, "n_bins", this % n_bins)
-    call write_dataset(filter_group, "order", this % order)
-    if (this % cosine == COSINE_SCATTER) then
-      call write_dataset(filter_group, "cosine", "scatter")
-    else
-      call write_dataset(filter_group, "cosine", "particle")
-    end if
-  end subroutine to_statepoint
-
-  function text_label(this, bin) result(label)
-    class(SphericalHarmonicsFilter), intent(in) :: this
-    integer,             intent(in) :: bin
-    character(MAX_LINE_LEN)         :: label
-
-    integer :: n, m
-
-    do n = 0, this % order
-      if (bin <= (n + 1)**2) then
-        m = (bin - n**2 - 1) - n
-        label = "Spherical harmonic expansion, Y" // trim(to_str(n)) // &
-             "," // trim(to_str(m))
-        exit
-      end if
-    end do
-  end function text_label
+  function cosine(this) result(val)
+    class(SphericalHarmonicsFilter) :: this
+    integer                         :: val
+    interface
+      function sphharm_filter_get_cosine(filt) result(val) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: filt
+        integer(C_INT) :: val
+      end function
+    end interface
+    val = sphharm_filter_get_cosine(this % ptr)
+  end function cosine
 
 !===============================================================================
 !                               C API FUNCTIONS
@@ -143,11 +51,19 @@ contains
     integer(C_INT),     intent(out) :: order
     integer(C_INT) :: err
 
+    interface
+      function sphharm_filter_get_order(filt) result(order) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR), value :: filt
+        integer(C_INT)     :: order
+      end function sphharm_filter_get_order
+    end interface
+
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (SphericalHarmonicsFilter)
-        order = f % order
+        order = sphharm_filter_get_order(f % ptr)
       class default
         err = E_INVALID_TYPE
         call set_errmsg("Not a spherical harmonics filter.")
@@ -169,7 +85,7 @@ contains
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (SphericalHarmonicsFilter)
-        select case (f % cosine)
+        select case (f % cosine())
         case (COSINE_SCATTER)
           cosine_ = 'scatter'
         case (COSINE_PARTICLE)
@@ -196,12 +112,19 @@ contains
     integer(C_INT),     value :: order
     integer(C_INT) :: err
 
+    interface
+      subroutine sphharm_filter_set_order(filt, order) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),    value :: filt
+        integer(C_INT), value :: order
+      end subroutine
+    end interface
+
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
       type is (SphericalHarmonicsFilter)
-        f % order = order
-        f % n_bins = (order + 1)**2
+        call sphharm_filter_set_order(f % ptr, order)
       class default
         err = E_INVALID_TYPE
         call set_errmsg("Not a spherical harmonics filter.")
@@ -218,6 +141,14 @@ contains
 
     character(:), allocatable :: cosine_
 
+    interface
+      subroutine sphharm_filter_set_cosine(filt, cosine) bind(C)
+        import C_PTR, C_INT
+        type(C_PTR),    value :: filt
+        integer(C_INT), value :: cosine
+      end subroutine
+    end interface
+
     ! Convert C string to Fortran string
     cosine_ = to_f_string(cosine)
 
@@ -227,9 +158,9 @@ contains
       type is (SphericalHarmonicsFilter)
         select case (cosine_)
         case ('scatter')
-          f % cosine = COSINE_SCATTER
+          call sphharm_filter_set_cosine(f % ptr, COSINE_SCATTER)
         case ('particle')
-          f % cosine = COSINE_PARTICLE
+          call sphharm_filter_set_cosine(f % ptr, COSINE_PARTICLE)
         end select
 
       class default

src/tallies/tally_header.F90

@@ -344,7 +344,7 @@ contains
         this % estimator = ESTIMATOR_ANALOG
       type is (SphericalHarmonicsFilter)
         j = FILTER_SPH_HARMONICS
-        if (filt % cosine == COSINE_SCATTER) then
+        if (filt % cosine() == COSINE_SCATTER) then
           this % estimator = ESTIMATOR_ANALOG
         end if
       type is (SpatialLegendreFilter)