Merge pull request #1037 from hanzhuoran/zernike

Added a ZernikeRadialFilter

docs/source/pythonapi/base.rst

@@ -122,6 +122,7 @@ Constructing Tallies
    openmc.SpatialLegendreFilter
    openmc.SphericalHarmonicsFilter
    openmc.ZernikeFilter
+   openmc.ZernikeRadialFilter
    openmc.ParticleFilter
    openmc.Mesh
    openmc.Trigger

openmc/capi/filter.py

@@ -20,7 +20,7 @@ __all__ = ['Filter', 'AzimuthalFilter', 'CellFilter',
            'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MeshFilter',
            'MeshSurfaceFilter', 'MuFilter', 'PolarFilter', 'SphericalHarmonicsFilter',
            'SpatialLegendreFilter', 'SurfaceFilter',
-           'UniverseFilter', 'ZernikeFilter', 'filters']
+           'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters']
 
 # Tally functions
 _dll.openmc_cell_filter_get_bins.argtypes = [
@@ -360,6 +360,10 @@ class ZernikeFilter(Filter):
         _dll.openmc_zernike_filter_set_order(self._index, order)
 
 
+class ZernikeRadialFilter(ZernikeFilter):
+    filter_type = 'zernikeradial'
+
+
 _FILTER_TYPE_MAP = {
     'azimuthal': AzimuthalFilter,
     'cell': CellFilter,
@@ -380,7 +384,8 @@ _FILTER_TYPE_MAP = {
     'spatiallegendre': SpatialLegendreFilter,
     'surface': SurfaceFilter,
     'universe': UniverseFilter,
-    'zernike': ZernikeFilter
+    'zernike': ZernikeFilter,
+    'zernikeradial': ZernikeRadialFilter
 }
 
 

openmc/capi/math.py

@@ -21,6 +21,9 @@ _dll.calc_rn_c.argtypes = [c_int, ndpointer(c_double), ndpointer(c_double)]
 _dll.calc_zn_c.restype = None
 _dll.calc_zn_c.argtypes = [c_int, c_double, c_double, ndpointer(c_double)]
 
+_dll.calc_zn_rad_c.restype = None
+_dll.calc_zn_rad_c.argtypes = [c_int, c_double, ndpointer(c_double)]
+
 _dll.rotate_angle_c.restype = None
 _dll.rotate_angle_c.argtypes = [ndpointer(c_double), c_double,
                                 POINTER(c_double)]
@@ -155,6 +158,32 @@ def calc_zn(n, rho, phi):
     return zn
 
 
+def calc_zn_rad(n, rho):
+    """ Calculate the even orders in n-th order modified Zernike polynomial
+    moment with no azimuthal dependency (m=0) for a given radial location in
+    the unit disk. The normalization of the polynomials is such that the
+    integral of Z_pq*Z_pq over the unit disk is exactly pi.
+
+    Parameters
+    ----------
+    n : int
+        Maximum order
+    rho : float
+        Radial location in the unit disk
+
+    Returns
+    -------
+    numpy.ndarray
+        Corresponding resulting list of coefficients
+
+    """
+
+    num_bins = n // 2 + 1
+    zn_rad = np.zeros(num_bins, dtype=np.float64)
+    _dll.calc_zn_rad_c(n, rho, zn_rad)
+    return zn_rad
+    
+
 def rotate_angle(uvw0, mu, phi=None):
     """ Rotates direction cosines through a polar angle whose cosine is
     mu and through an azimuthal angle sampled uniformly.

openmc/filter.py

@@ -22,7 +22,7 @@ _FILTER_TYPES = (
     'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy',
     'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup',
     'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre',
-    'sphericalharmonics', 'zernike', 'particle'
+    'sphericalharmonics', 'zernike', 'zernikeradial', 'particle'
 )
 
 _CURRENT_NAMES = (

openmc/filter_expansion.py

@@ -234,8 +234,8 @@ class SphericalHarmonicsFilter(ExpansionFilter):
     r"""Score spherical harmonic expansion moments up to specified order.
 
     This filter allows you to obtain real spherical harmonic moments of either
-    the particle's direction or the cosine of the scattering angle. Specifying a
-    filter with order :math:`\ell` tallies moments for all orders from 0 to
+    the particle's direction or the cosine of the scattering angle. Specifying
+    a filter with order :math:`\ell` tallies moments for all orders from 0 to
     :math:`\ell`.
 
     Parameters
@@ -342,11 +342,11 @@ class ZernikeFilter(ExpansionFilter):
         \frac{n+m}{2} - k)! (\frac{n-m}{2} - k)!} \rho^{n-2k}.
 
     With this definition, the integral of :math:`(Z_n^m)^2` over the unit disk
-    is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where :math:`\epsilon_m` is
-    2 if :math:`m` equals 0 and 1 otherwise.
+    is :math:`\frac{\epsilon_m\pi}{2n+2}` for each polynomial where 
+    :math:`\epsilon_m` is 2 if :math:`m` equals 0 and 1 otherwise.
 
-    Specifying a filter with order N tallies moments for all :math:`n` from 0 to
-    N and each value of :math:`m`. The ordering of the Zernike polynomial
+    Specifying a filter with order N tallies moments for all :math:`n` from 0
+    to N and each value of :math:`m`. The ordering of the Zernike polynomial
     moments follows the ANSI Z80.28 standard, where the one-dimensional index
     :math:`j` corresponds to the :math:`n` and :math:`m` by
 
@@ -463,3 +463,63 @@ class ZernikeFilter(ExpansionFilter):
         subelement.text = str(self.r)
 
         return element
+
+
+class ZernikeRadialFilter(ZernikeFilter):
+    r"""Score the :math:`m = 0` (radial variation only) Zernike moments up to
+    specified order.
+
+    The Zernike polynomials are defined the same as in :class:`ZernikeFilter`.
+
+    .. math::
+
+        Z_n^{0}(\rho, \theta) = R_n^{0}(\rho)
+
+    where the radial polynomials are
+
+    .. math::
+        R_n^{0}(\rho) = \sum\limits_{k=0}^{n/2} \frac{(-1)^k (n-k)!}{k! ((
+        \frac{n}{2} - k)!)^{2}} \rho^{n-2k}.
+
+    With this definition, the integral of :math:`(Z_n^0)^2` over the unit disk
+    is :math:`\frac{\pi}{n+1}`.
+
+    If there is only radial dependency, the polynomials are integrated over
+    the azimuthal angles. The only terms left are :math:`Z_n^{0}(\rho, \theta) 
+    = R_n^{0}(\rho)`. Note that :math:`n` could only be even orders.
+    Therefore, for a radial Zernike polynomials up to order of :math:`n`,
+    there are :math:`\frac{n}{2} + 1` terms in total. The indexing is from the
+    lowest even order (0) to highest even order.
+
+    Parameters
+    ----------
+    order : int
+        Maximum radial Zernike polynomial order
+    x : float
+        x-coordinate of center of circle for normalization
+    y : float
+        y-coordinate of center of circle for normalization
+    r : int or None
+        Radius of circle for normalization
+
+    Attributes
+    ----------
+    order : int
+        Maximum radial Zernike polynomial order
+    x : float
+        x-coordinate of center of circle for normalization
+    y : float
+        y-coordinate of center of circle for normalization
+    r : int or None
+        Radius of circle for normalization
+    id : int
+        Unique identifier for the filter
+    num_bins : int
+        The number of filter bins
+
+    """
+
+    @ExpansionFilter.order.setter
+    def order(self, order):
+        ExpansionFilter.order.__set__(self, order)
+        self.bins = ['Z{},0'.format(n) for n in range(0, order+1, 2)]

src/constants.F90

@@ -357,7 +357,7 @@ module constants
   integer, parameter :: NO_BIN_FOUND = -1
 
   ! Tally filter and map types
-  integer, parameter :: N_FILTER_TYPES = 21
+  integer, parameter :: N_FILTER_TYPES = 22
   integer, parameter :: &
        FILTER_UNIVERSE       = 1,  &
        FILTER_MATERIAL       = 2,  &
@@ -379,7 +379,9 @@ module constants
        FILTER_SPH_HARMONICS  = 18, &
        FILTER_SPTL_LEGENDRE  = 19, &
        FILTER_ZERNIKE        = 20, &
-       FILTER_PARTICLE       = 21
+       FILTER_ZERNIKE_RADIAL = 21, &
+       FILTER_PARTICLE       = 22
+
 
   ! Mesh types
   integer, parameter :: &

src/math.F90

@@ -11,6 +11,7 @@ module math
   public :: calc_pn
   public :: calc_rn
   public :: calc_zn
+  public :: calc_zn_rad
   public :: evaluate_legendre
   public :: rotate_angle
   public :: maxwell_spectrum
@@ -68,6 +69,14 @@ module math
       real(C_DOUBLE), intent(out) :: zn(((n + 1) * (n + 2)) / 2)
     end subroutine calc_zn
 
+    pure subroutine calc_zn_rad(n, rho, zn_rad) bind(C, name='calc_zn_rad_c')
+      use ISO_C_BINDING
+      implicit none
+      integer(C_INT), value, intent(in) :: n
+      real(C_DOUBLE), value, intent(in) :: rho
+      real(C_DOUBLE), intent(out) :: zn_rad((n / 2) + 1)
+    end subroutine calc_zn_rad
+
     subroutine rotate_angle_c_intfc(uvw, mu, phi) bind(C, name='rotate_angle_c')
       use ISO_C_BINDING
       implicit none

src/math_functions.cpp

@@ -587,6 +587,33 @@ void calc_zn_c(int n, double rho, double phi, double zn[]) {
 
 }
 
+void calc_zn_rad_c(int n, double rho, double zn_rad[]) {
+  // Calculate R_p0(rho) as Zn_p0(rho)
+  // Set up the array of the coefficients
+
+  double q = 0;
+
+  // R_00 is always 1
+  zn_rad[0] = 1;
+
+  // Fill in the rest of the array (Eq 3.8 and Eq 3.10 in Chong)
+  for (int p = 2; p <= n; p += 2) {
+    int index = int(p/2);
+    if (p == 2) {
+    // Setting up R_22 to calculate R_20 (Eq 3.10 in Chong)
+      double R_22 = rho * rho;
+      zn_rad[index] = 2 * R_22 - zn_rad[0];
+    } else {
+      double k1 = ((p + q) * (p - q) * (p - 2)) / 2.;
+      double k2 = 2 * p * (p - 1) * (p - 2);
+      double k3 = -q * q * (p - 1) - p * (p - 1) * (p - 2);
+      double k4 = (-p * (p + q - 2) * (p - q - 2)) / 2.;
+      zn_rad[index] = 
+        ((k2 * rho * rho + k3) * zn_rad[index-1] + k4 * zn_rad[index-2]) / k1;
+    }
+  }
+}
+
 
 void rotate_angle_c(double uvw[3], double mu, double* phi) {
   // Copy original directional cosines

src/math_functions.h

@@ -91,6 +91,24 @@ extern "C" void calc_rn_c(int n, const double uvw[3], double rn[]);
 
 extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
 
+//==============================================================================
+//! Calculate only the even order components of n-th order modified Zernike 
+//! polynomial moment with azimuthal dependency m = 0 for a given radial (rho)
+//! location on the unit disk.
+//!
+//! Since m = 0, n could only be even orders. Z_q0 = R_q0
+//!
+//! See calc_zn_c for methodology.
+//!
+//! @param n       The maximum order requested
+//! @param rho     The radial parameter to specify location on the unit disk
+//! @param phi     The angle parameter to specify location on the unit disk
+//! @param zn_rad  The requested moments of order 0 to n (inclusive)
+//!   evaluated at rho and phi when m = 0.
+//==============================================================================
+
+extern "C" void calc_zn_rad_c(int n, double rho, double zn_rad[]);
+
 //==============================================================================
 //! Rotate the direction cosines through a polar angle whose cosine is mu and
 //! through an azimuthal angle sampled uniformly.

src/tallies/tally_filter.F90

@@ -91,6 +91,8 @@ contains
         type_ = 'zernike'
       type is (ParticleFilter)
         type_ = 'particle'
+      type is (ZernikeRadialFilter)
+        type_ = 'zernikeradial'
       end select
 
       ! Convert Fortran string to null-terminated C string. We assume the
@@ -164,6 +166,8 @@ contains
           allocate(UniverseFilter :: filters(index) % obj)
         case ('zernike')
           allocate(ZernikeFilter :: filters(index) % obj)
+        case ('zernikeradial')
+          allocate(ZernikeRadialFilter :: filters(index) % obj)
         case default
           err = E_UNASSIGNED
           call set_errmsg("Unknown filter type: " // trim(type_))

src/tallies/tally_filter_zernike.F90

@@ -5,7 +5,7 @@ module tally_filter_zernike
   use constants
   use error
   use hdf5_interface
-  use math,                only: calc_zn
+  use math,                only: calc_zn, calc_zn_rad
   use particle_header,     only: Particle
   use string,              only: to_str
   use tally_filter_header
@@ -13,6 +13,10 @@ module tally_filter_zernike
 
   implicit none
   private
+  public :: openmc_zernike_filter_get_order
+  public :: openmc_zernike_filter_get_params
+  public :: openmc_zernike_filter_set_order
+  public :: openmc_zernike_filter_set_params
 
 !===============================================================================
 ! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position
@@ -24,19 +28,43 @@ module tally_filter_zernike
     real(8) :: y
     real(8) :: r
   contains
-    procedure :: from_xml
-    procedure :: get_all_bins
-    procedure :: to_statepoint
-    procedure :: text_label
+    procedure :: calc_n_bins => calc_n_bins_zn
+    procedure :: from_xml => from_xml_zn
+    procedure :: get_all_bins => get_all_bins_zn
+    procedure :: to_statepoint => to_statepoint_zn
+    procedure :: text_label => text_label_zn
   end type ZernikeFilter
 
+!===============================================================================
+! ZERNIKERADIALFILTER gives even order radial Zernike polynomial moments of a
+! particle's position
+!===============================================================================
+
+  type, public, extends(ZernikeFilter) :: ZernikeRadialFilter
+  contains
+    procedure :: calc_n_bins => calc_n_bins_zn_rad
+    ! Inherit from_xml from ZernikeFilter
+    procedure :: get_all_bins => get_all_bins_zn_rad
+    procedure :: to_statepoint => to_statepoint_zn_rad
+    procedure :: text_label => text_label_zn_rad
+  end type ZernikeRadialFilter
+
 contains
 
 !===============================================================================
 ! ZernikeFilter methods
 !===============================================================================
 
-  subroutine from_xml(this, node)
+  function calc_n_bins_zn(this) result(n_bins)
+    class(ZernikeFilter), intent(in) :: this
+    integer :: n
+    integer :: n_bins
+
+    n = this % order
+    n_bins = ((n+1) * (n+2))/2
+  end function calc_n_bins_zn
+
+  subroutine from_xml_zn(this, node)
     class(ZernikeFilter), intent(inout) :: this
     type(XMLNode), intent(in) :: node
 
@@ -50,10 +78,10 @@ contains
     ! Get specified order
     call get_node_value(node, "order", n)
     this % order = n
-    this % n_bins = ((n + 1)*(n + 2))/2
-  end subroutine from_xml
+    this % n_bins = this % calc_n_bins()
+  end subroutine from_xml_zn
 
-  subroutine get_all_bins(this, p, estimator, match)
+  subroutine get_all_bins_zn(this, p, estimator, match)
     class(ZernikeFilter), intent(in)  :: this
     type(Particle),      intent(in)  :: p
     integer,             intent(in)  :: estimator
@@ -67,18 +95,20 @@ contains
     x = p % coord(1) % xyz(1) - this % x
     y = p % coord(1) % xyz(2) - this % y
     r = sqrt(x*x + y*y)/this % r
-    theta = atan2(y, x)
+    if (r <= ONE) then
+      theta = atan2(y, x)
+      
+      ! Get moments for Zernike polynomial orders 0..n
+      call calc_zn(this % order, r, theta, zn)
 
-    ! Get moments for Zernike polynomial orders 0..n
-    call calc_zn(this % order, r, theta, zn)
+      do i = 1, this % n_bins
+        call match % bins % push_back(i)
+        call match % weights % push_back(zn(i))
+      end do
+    endif
+  end subroutine get_all_bins_zn
 
-    do i = 1, this % n_bins
-      call match % bins % push_back(i)
-      call match % weights % push_back(zn(i))
-    end do
-  end subroutine get_all_bins
-
-  subroutine to_statepoint(this, filter_group)
+  subroutine to_statepoint_zn(this, filter_group)
     class(ZernikeFilter), intent(in) :: this
     integer(HID_T),      intent(in) :: filter_group
 
@@ -88,9 +118,9 @@ contains
     call write_dataset(filter_group, "x", this % x)
     call write_dataset(filter_group, "y", this % y)
     call write_dataset(filter_group, "r", this % r)
-  end subroutine to_statepoint
+  end subroutine to_statepoint_zn
 
-  function text_label(this, bin) result(label)
+  function text_label_zn(this, bin) result(label)
     class(ZernikeFilter), intent(in) :: this
     integer,             intent(in) :: bin
     character(MAX_LINE_LEN)         :: label
@@ -108,7 +138,66 @@ contains
         exit
       end if
     end do
-  end function text_label
+  end function text_label_zn
+
+!===============================================================================
+! ZernikeRadialFilter methods
+!===============================================================================
+
+  function calc_n_bins_zn_rad(this) result(num_n_bins)
+    class(ZernikeRadialFilter), intent(in) :: this
+    integer :: n
+    integer :: num_n_bins
+
+    n = this % order
+    num_n_bins = n/2 + 1
+  end function calc_n_bins_zn_rad
+
+  subroutine get_all_bins_zn_rad(this, p, estimator, match)
+    class(ZernikeRadialFilter), intent(in)  :: this
+    type(Particle),      intent(in)  :: p
+    integer,             intent(in)  :: estimator
+    type(TallyFilterMatch),   intent(inout) :: match
+
+    integer :: i
+    real(8) :: x, y, r
+    real(C_DOUBLE) :: zn_rad(this % n_bins)
+
+    ! Determine normalized (r,theta) positions
+    x = p % coord(1) % xyz(1) - this % x
+    y = p % coord(1) % xyz(2) - this % y
+    r = sqrt(x*x + y*y)/this % r
+    if (r <= ONE) then
+
+      ! Get moments for even order Zernike polynomial orders 0..n
+      call calc_zn_rad(this % order, r, zn_rad)
+
+      do i = 1, this % n_bins
+        call match % bins % push_back(i)
+        call match % weights % push_back(zn_rad(i))
+      end do
+    endif
+  end subroutine get_all_bins_zn_rad
+
+  subroutine to_statepoint_zn_rad(this, filter_group)
+    class(ZernikeRadialFilter), intent(in) :: this
+    integer(HID_T),      intent(in) :: filter_group
+
+    call write_dataset(filter_group, "type", "zernikeradial")
+    call write_dataset(filter_group, "n_bins", this % n_bins)
+    call write_dataset(filter_group, "order", this % order)
+    call write_dataset(filter_group, "x", this % x)
+    call write_dataset(filter_group, "y", this % y)
+    call write_dataset(filter_group, "r", this % r)
+  end subroutine to_statepoint_zn_rad
+
+  function text_label_zn_rad(this, bin) result(label)
+    class(ZernikeRadialFilter), intent(in) :: this
+    integer,             intent(in) :: bin
+    character(MAX_LINE_LEN)         :: label
+
+    label = "Zernike expansion, Z" // trim(to_str(2*(bin-1))) // ",0"
+  end function text_label_zn_rad
 
 !===============================================================================
 !                               C API FUNCTIONS
@@ -125,6 +214,8 @@ contains
       select type (f => filters(index) % obj)
       type is (ZernikeFilter)
         order = f % order
+      type is (ZernikeRadialFilter)
+        order = f % order
       class default
         err = E_INVALID_TYPE
         call set_errmsg("Not a Zernike filter.")
@@ -144,7 +235,7 @@ contains
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
-      type is (ZernikeFilter)
+      class is (ZernikeFilter)
         x = f % x
         y = f % y
         r = f % r
@@ -168,6 +259,9 @@ contains
       type is (ZernikeFilter)
         f % order = order
         f % n_bins = ((order + 1)*(order + 2))/2
+      type is (ZernikeRadialFilter)
+        f % order = order
+        f % n_bins = order/2 + 1
       class default
         err = E_INVALID_TYPE
         call set_errmsg("Not a Zernike filter.")
@@ -187,7 +281,7 @@ contains
     err = verify_filter(index)
     if (err == 0) then
       select type (f => filters(index) % obj)
-      type is (ZernikeFilter)
+      class is (ZernikeFilter)
         if (present(x)) f % x = x
         if (present(y)) f % y = y
         if (present(r)) f % r = r

src/tallies/tally_header.F90

@@ -340,6 +340,9 @@ contains
       type is (ZernikeFilter)
         j = FILTER_ZERNIKE
         this % estimator = ESTIMATOR_COLLISION
+      type is (ZernikeRadialFilter)
+        j = FILTER_ZERNIKE_RADIAL
+        this % estimator = ESTIMATOR_COLLISION
       type is (ParticleFilter)
         j = FILTER_PARTICLE
       end select

tests/unit_tests/test_filters.py

@@ -100,6 +100,23 @@ def test_zernike():
     assert elem.attrib['type'] == 'zernike'
     assert elem.find('order').text == str(n)
 
+def test_zernike_radial():
+    n = 4
+    f = openmc.ZernikeRadialFilter(n, 0., 0., 1.)
+    assert f.order == n
+    assert f.bins[0] == 'Z0,0'
+    assert f.bins[-1] == 'Z{},0'.format(n)
+    assert len(f.bins) == n//2 + 1
+
+    # Make sure __repr__ works
+    repr(f)
+
+    # to_xml_element()
+    elem = f.to_xml_element()
+    assert elem.tag == 'filter'
+    assert elem.attrib['type'] == 'zernikeradial'
+    assert elem.find('order').text == str(n)
+
 
 def test_first_moment(run_in_tmpdir, box_model):
     plain_tally = openmc.Tally()