crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

updating branch with develop.

Browse source 
Merge branch 'develop' of https://github.com/mit-crpg/openmc into new-update
This commit is contained in:
Jose Salcedo Perez 2018-04-12 16:35:52 +00:00
commit 17b2cd6c65
27 changed files with 2556 additions and 797 deletions
Download patch file Download diff file Expand all files Collapse all files

4
CMakeLists.txt
View file

@ -420,13 +420,17 @@ set(LIBOPENMC_FORTRAN_SRC
src/tallies/tally_filter_distribcell.F90
src/tallies/tally_filter_energy.F90
src/tallies/tally_filter_energyfunc.F90
src/tallies/tally_filter_legendre.F90
src/tallies/tally_filter_material.F90
src/tallies/tally_filter_mesh.F90
src/tallies/tally_filter_meshsurface.F90
src/tallies/tally_filter_mu.F90
src/tallies/tally_filter_polar.F90
src/tallies/tally_filter_sph_harm.F90
src/tallies/tally_filter_sptl_legendre.F90
src/tallies/tally_filter_surface.F90
src/tallies/tally_filter_universe.F90
src/tallies/tally_filter_zernike.F90
src/tallies/tally_header.F90
src/tallies/trigger.F90
src/tallies/trigger_header.F90

5
docs/source/capi/index.rst
View file

@ -247,11 +247,12 @@ Functions
:return: Return status (negative if an error occurs)
:rtype: int
.. c:function:: int openmc_get_nuclide_index(char name[], int* index)
.. c:function:: int openmc_get_nuclide_index(const char name[], int* index)
Get the index in the nuclides array for a nuclide with a given name
:param char[] name: Name of the nuclide
:param name: Name of the nuclide
:type name: const char[]
:param int* index: Index in the nuclides array
:return: Return status (negative if an error occurs)
:rtype: int

13
docs/source/examples/expansion-filters.rst Normal file
View file

@ -0,0 +1,13 @@
.. _notebook_expansion:
=====================
Functional Expansions
=====================
.. only:: html
.. notebook:: ../../../examples/jupyter/expansion-filters.ipynb
.. only:: latex
IPython notebooks must be viewed in the online HTML documentation.

14
docs/source/examples/index.rst
View file

@ -1,13 +1,12 @@
.. _examples:
=================
Example Notebooks
=================
========
Examples
========
The following series of Jupyter_ Notebooks provide examples for usage of OpenMC
features via the :ref:`pythonapi`.
.. _Jupyter: https://jupyter.org/
The following series of `Jupyter <https://jupyter.org/>`_ Notebooks provide
examples for how to use various features of OpenMC by leveraging the
:ref:`pythonapi`.
-----------
Basic Usage
@ -20,6 +19,7 @@ Basic Usage
post-processing
pandas-dataframes
tally-arithmetic
expansion-filters
search
triso
candu

4
docs/source/pythonapi/base.rst
View file

@ -118,6 +118,10 @@ Constructing Tallies
openmc.DistribcellFilter
openmc.DelayedGroupFilter
openmc.EnergyFunctionFilter
openmc.LegendreFilter
openmc.SpatialLegendreFilter
openmc.SphericalHarmonicsFilter
openmc.ZernikeFilter
openmc.Mesh
openmc.Trigger
openmc.TallyDerivative

463
examples/jupyter/expansion-filters.ipynb Normal file
View file

File diff suppressed because one or more lines are too long

18
include/openmc.h
View file

@ -38,10 +38,12 @@ extern "C" {
void openmc_get_filter_next_id(int32_t* id);
int openmc_get_keff(double k_combined[]);
int openmc_get_material_index(int32_t id, int32_t* index);
int openmc_get_nuclide_index(char name[], int* index);
int openmc_get_nuclide_index(const char name[], int* index);
int openmc_get_tally_index(int32_t id, int32_t* index);
void openmc_hard_reset();
void openmc_init(const int* intracomm);
int openmc_legendre_filter_get_order(int32_t index, int* order);
int openmc_legendre_filter_set_order(int32_t index, int order);
int openmc_load_nuclide(char name[]);
int openmc_material_add_nuclide(int32_t index, const char name[], double density);
int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n);
@ -62,6 +64,15 @@ extern "C" {
void openmc_simulation_init();
int openmc_source_bank(struct Bank** ptr, int64_t* n);
int openmc_source_set_strength(int32_t index, double strength);
int openmc_spatial_legendre_filter_get_order(int32_t index, int* order);
int openmc_spatial_legendre_filter_get_params(int32_t index, int* axis, double* min, double* max);
int openmc_spatial_legendre_filter_set_order(int32_t index, int order);
int openmc_spatial_legendre_filter_set_params(int32_t index, const int* axis,
const double* min, const double* max);
int openmc_sphharm_filter_get_order(int32_t index, int* order);
int openmc_sphharm_filter_get_cosine(int32_t index, char cosine[]);
int openmc_sphharm_filter_set_order(int32_t index, int order);
int openmc_sphharm_filter_set_cosine(int32_t index, const char cosine[]);
void openmc_statepoint_write(const char filename[]);
int openmc_tally_get_id(int32_t index, int32_t* id);
int openmc_tally_get_filters(int32_t index, int32_t** indices, int* n);
@ -73,6 +84,11 @@ extern "C" {
int openmc_tally_set_id(int32_t index, int32_t id);
int openmc_tally_set_nuclides(int32_t index, int n, const char** nuclides);
int openmc_tally_set_scores(int32_t index, int n, const int* scores);
int openmc_zernike_filter_get_order(int32_t index, int* order);
int openmc_zernike_filter_get_params(int32_t index, double* x, double* y, double* r);
int openmc_zernike_filter_set_order(int32_t index, int order);
int openmc_zernike_filter_set_params(int32_t index, const double* x,
const double* y, const double* r);
// Error codes
extern int E_UNASSIGNED;

1
openmc/__init__.py
View file

@ -15,6 +15,7 @@ from openmc.surface import *
from openmc.universe import *
from openmc.lattice import *
from openmc.filter import *
from openmc.filter_expansion import *
from openmc.trigger import *
from openmc.tally_derivative import *
from openmc.tallies import *

4
openmc/checkvalue.py
View file

@ -196,7 +196,7 @@ def check_less_than(name, value, maximum, equality=False):
maximum : object
Maximum value to check against
equality : bool, optional
Whether equality is allowed. Defaluts to False.
Whether equality is allowed. Defaults to False.
"""
@ -223,7 +223,7 @@ def check_greater_than(name, value, minimum, equality=False):
minimum : object
Minimum value to check against
equality : bool, optional
Whether equality is allowed. Defaluts to False.
Whether equality is allowed. Defaults to False.
"""

785
openmc/filter.py
View file

File diff suppressed because it is too large Load diff

457
openmc/filter_expansion.py Normal file
View file

@ -0,0 +1,457 @@
from numbers import Integral, Real
from xml.etree import ElementTree as ET
import numpy as np
import pandas as pd
import openmc.checkvalue as cv
from . import Filter
class ExpansionFilter(Filter):
"""Abstract filter class for functional expansions."""
def __init__(self, order, filter_id=None):
self.order = order
self.id = filter_id
@property
def order(self):
return self._order
@order.setter
def order(self, order):
cv.check_type('expansion order', order, Integral)
cv.check_greater_than('expansion order', order, 0, equality=True)
self._order = order
def to_xml_element(self):
"""Return XML Element representing the filter.
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing Legendre filter data
"""
element = ET.Element('filter')
element.set('id', str(self.id))
element.set('type', self.short_name.lower())
subelement = ET.SubElement(element, 'order')
subelement.text = str(self.order)
return element
class LegendreFilter(ExpansionFilter):
r"""Score Legendre expansion moments up to specified order.
This filter allows scores to be multiplied by Legendre polynomials of the
change in particle angle ($\mu$) up to a user-specified order.
Parameters
----------
order : int
Maximum Legendre polynomial order
filter_id : int or None
Unique identifier for the filter
Attributes
----------
order : int
Maximum Legendre polynomial order
id : int
Unique identifier for the filter
num_bins : int
The number of filter bins
"""
def __hash__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
return hash(string)
def __repr__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tID', self.id)
return string
@ExpansionFilter.order.setter
def order(self, order):
ExpansionFilter.order.__set__(self, order)
self.bins = ['P{}'.format(i) for i in range(order + 1)]
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(group['order'].value, filter_id)
return out
class SpatialLegendreFilter(ExpansionFilter):
r"""Score Legendre expansion moments in space up to specified order.
This filter allows scores to be multiplied by Legendre polynomials of the
the particle's position along a particular axis, normalized to a given
range, up to a user-specified order.
Parameters
----------
order : int
Maximum Legendre polynomial order
axis : {'x', 'y', 'z'}
Axis along which to take the expansion
minimum : float
Minimum value along selected axis
maximum : float
Maximum value along selected axis
filter_id : int or None
Unique identifier for the filter
Attributes
----------
order : int
Maximum Legendre polynomial order
axis : {'x', 'y', 'z'}
Axis along which to take the expansion
minimum : float
Minimum value along selected axis
maximum : float
Maximum value along selected axis
id : int
Unique identifier for the filter
num_bins : int
The number of filter bins
"""
def __init__(self, order, axis, minimum, maximum, filter_id=None):
super().__init__(order, filter_id)
self.axis = axis
self.minimum = minimum
self.maximum = maximum
def __hash__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tAxis', self.axis)
string += '{: <16}=\t{}\n'.format('\tMin', self.minimum)
string += '{: <16}=\t{}\n'.format('\tMax', self.maximum)
return hash(string)
def __repr__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tAxis', self.axis)
string += '{: <16}=\t{}\n'.format('\tMin', self.minimum)
string += '{: <16}=\t{}\n'.format('\tMax', self.maximum)
string += '{: <16}=\t{}\n'.format('\tID', self.id)
return string
@ExpansionFilter.order.setter
def order(self, order):
ExpansionFilter.order.__set__(self, order)
self.bins = ['P{}'.format(i) for i in range(order + 1)]
@property
def axis(self):
return self._axis
@axis.setter
def axis(self, axis):
cv.check_value('axis', axis, ('x', 'y', 'z'))
self._axis = axis
@property
def minimum(self):
return self._minimum
@minimum.setter
def minimum(self, minimum):
cv.check_type('minimum', minimum, Real)
self._minimum = minimum
@property
def maximum(self):
return self._maximum
@maximum.setter
def maximum(self, maximum):
cv.check_type('maximum', maximum, Real)
self._maximum = maximum
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
order = group['order'].value
axis = group['axis'].value.decode()
min_, max_ = group['min'].value, group['max'].value
return cls(order, axis, min_, max_, filter_id)
def to_xml_element(self):
"""Return XML Element representing the filter.
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing Legendre filter data
"""
element = super().to_xml_element()
subelement = ET.SubElement(element, 'axis')
subelement.text = self.axis
subelement = ET.SubElement(element, 'min')
subelement.text = str(self.minimum)
subelement = ET.SubElement(element, 'max')
subelement.text = str(self.maximum)
return element
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
:math:`\ell`.
Parameters
----------
order : int
Maximum spherical harmonics order, :math:`\ell`
filter_id : int or None
Unique identifier for the filter
Attributes
----------
order : int
Maximum spherical harmonics order, :math:`\ell`
id : int
Unique identifier for the filter
cosine : {'scatter', 'particle'}
How to handle the cosine term.
num_bins : int
The number of filter bins
"""
def __init__(self, order, filter_id=None):
super().__init__(order, filter_id)
self._cosine = 'particle'
def __hash__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tCosine', self.cosine)
return hash(string)
def __repr__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tCosine', self.cosine)
string += '{: <16}=\t{}\n'.format('\tID', self.id)
return string
@ExpansionFilter.order.setter
def order(self, order):
ExpansionFilter.order.__set__(self, order)
self.bins = ['Y{},{}'.format(n, m)
for n in range(order + 1)
for m in range(-n, n + 1)]
@property
def cosine(self):
return self._cosine
@cosine.setter
def cosine(self, cosine):
cv.check_value('Spherical harmonics cosine treatment', cosine,
('scatter', 'particle'))
self._cosine = cosine
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
out = cls(group['order'].value, filter_id)
out.cosine = group['cosine'].value.decode()
return out
def to_xml_element(self):
"""Return XML Element representing the filter.
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing spherical harmonics filter data
"""
element = super().to_xml_element()
element.set('cosine', self.cosine)
return element
class ZernikeFilter(ExpansionFilter):
r"""Score Zernike expansion moments in space up to specified order.
This filter allows scores to be multiplied by Zernike polynomials of the
particle's position normalized to a given unit circle, up to a
user-specified order. The Zernike polynomials follow the definition by `Noll
<https://doi.org/10.1364/JOSA.66.000207>`_ and are defined as
.. math::
Z_n^m(\rho, \theta) = \sqrt{2n + 2} R_n^m(\rho) \cos (m\theta), \quad m > 0
Z_n^{m}(\rho, \theta) = \sqrt{2n + 2} R_n^{m}(\rho) \sin (m\theta), \quad m < 0
Z_n^{m}(\rho, \theta) = \sqrt{n + 1} R_n^{m}(\rho), \quad m = 0
where the radial polynomials are
.. math::
R_n^m(\rho) = \sum\limits_{k=0}^{(n-m)/2} \frac{(-1)^k (n-k)!}{k! (
\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 exactly :math:`\pi` for each 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
.. math::
j = \frac{n(n + 2) + m}{2}.
Parameters
----------
order : int
Maximum 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 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
"""
def __init__(self, order, x=0.0, y=0.0, r=1.0, filter_id=None):
super().__init__(order, filter_id)
self.x = x
self.y = y
self.r = r
def __hash__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
return hash(string)
def __repr__(self):
string = type(self).__name__ + '\n'
string += '{: <16}=\t{}\n'.format('\tOrder', self.order)
string += '{: <16}=\t{}\n'.format('\tID', self.id)
return string
@ExpansionFilter.order.setter
def order(self, order):
ExpansionFilter.order.__set__(self, order)
self.bins = ['Z{},{}'.format(n, m)
for n in range(order + 1)
for m in range(-n, n + 1, 2)]
@property
def x(self):
return self._x
@x.setter
def x(self, x):
cv.check_type('x', x, Real)
self._x = x
@property
def y(self):
return self._y
@y.setter
def y(self, y):
cv.check_type('y', y, Real)
self._y = y
@property
def r(self):
return self._r
@r.setter
def r(self, r):
cv.check_type('r', r, Real)
self._r = r
@classmethod
def from_hdf5(cls, group, **kwargs):
if group['type'].value.decode() != cls.short_name.lower():
raise ValueError("Expected HDF5 data for filter type '"
+ cls.short_name.lower() + "' but got '"
+ group['type'].value.decode() + " instead")
filter_id = int(group.name.split('/')[-1].lstrip('filter '))
order = group['order'].value
x, y, r = group['x'].value, group['y'].value, group['r'].value
return cls(order, x, y, r, filter_id)
def to_xml_element(self):
"""Return XML Element representing the filter.
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing Zernike filter data
"""
element = super().to_xml_element()
subelement = ET.SubElement(element, 'x')
subelement.text = str(self.x)
subelement = ET.SubElement(element, 'y')
subelement.text = str(self.y)
subelement = ET.SubElement(element, 'r')
subelement.text = str(self.r)
return element

3
openmc/mesh.py
View file

@ -38,6 +38,9 @@ class Mesh(IDManagerMixin):
are given, it is assumed that the mesh is an x-y mesh.
width : Iterable of float
The width of mesh cells in each direction.
indices : list of tuple
A list of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1,
1), ...]
"""

12
openmc/mgxs/mgxs.py
View file

@ -1164,12 +1164,14 @@ class MGXS(metaclass=ABCMeta):
if not isinstance(tally_filter, (openmc.EnergyFilter,
openmc.EnergyoutFilter)):
continue
elif len(tally_filter.bins) != len(fine_edges):
elif len(tally_filter.bins) != len(fine_edges) - 1:
continue
elif not np.allclose(tally_filter.bins, fine_edges):
elif not np.allclose(tally_filter.bins[:, 0], fine_edges[:-1]):
continue
else:
tally_filter.bins = coarse_groups.group_edges
cedge = coarse_groups.group_edges
tally_filter.values = cedge
tally_filter.bins = np.vstack((cedge[:-1], cedge[1:])).T
mean = np.add.reduceat(mean, energy_indices, axis=i)
std_dev = np.add.reduceat(std_dev**2, energy_indices,
axis=i)
@ -2738,7 +2740,7 @@ class TransportXS(MGXS):
if self._rxn_rate_tally is None:
# Switch EnergyoutFilter to EnergyFilter.
old_filt = self.tallies['scatter-1'].filters[-1]
new_filt = openmc.EnergyFilter(old_filt.bins)
new_filt = openmc.EnergyFilter(old_filt.values)
self.tallies['scatter-1'].filters[-1] = new_filt
self._rxn_rate_tally = \
@ -2757,7 +2759,7 @@ class TransportXS(MGXS):
# Switch EnergyoutFilter to EnergyFilter.
old_filt = self.tallies['scatter-1'].filters[-1]
new_filt = openmc.EnergyFilter(old_filt.bins)
new_filt = openmc.EnergyFilter(old_filt.values)
self.tallies['scatter-1'].filters[-1] = new_filt
# Compute total cross section

203
openmc/tallies.py
View file

@ -218,10 +218,9 @@ class Tally(IDManagerMixin):
f = h5py.File(self._sp_filename, 'r')
# Extract Tally data from the file
data = f['tallies/tally {0}/results'.format(
self.id)].value
sum = data[:,:,0]
sum_sq = data[:,:,1]
data = f['tallies/tally {0}/results'.format(self.id)].value
sum = data[:, :, 0]
sum_sq = data[:, :, 1]
# Reshape the results arrays
sum = np.reshape(sum, self.shape)
@ -273,8 +272,8 @@ class Tally(IDManagerMixin):
# Convert NumPy array to SciPy sparse LIL matrix
if self.sparse:
self._mean = \
sps.lil_matrix(self._mean.flatten(), self._mean.shape)
self._mean = sps.lil_matrix(self._mean.flatten(),
self._mean.shape)
if self.sparse:
return np.reshape(self._mean.toarray(), self.shape)
@ -295,8 +294,8 @@ class Tally(IDManagerMixin):
# Convert NumPy array to SciPy sparse LIL matrix
if self.sparse:
self._std_dev = \
sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
self._std_dev = sps.lil_matrix(self._std_dev.flatten(),
self._std_dev.shape)
self.with_batch_statistics = True
@ -436,17 +435,16 @@ class Tally(IDManagerMixin):
# Convert NumPy arrays to SciPy sparse LIL matrices
if sparse and not self.sparse:
if self._sum is not None:
self._sum = \
sps.lil_matrix(self._sum.flatten(), self._sum.shape)
self._sum = sps.lil_matrix(self._sum.flatten(), self._sum.shape)
if self._sum_sq is not None:
self._sum_sq = \
sps.lil_matrix(self._sum_sq.flatten(), self._sum_sq.shape)
self._sum_sq = sps.lil_matrix(self._sum_sq.flatten(),
self._sum_sq.shape)
if self._mean is not None:
self._mean = \
sps.lil_matrix(self._mean.flatten(), self._mean.shape)
self._mean = sps.lil_matrix(self._mean.flatten(),
self._mean.shape)
if self._std_dev is not None:
self._std_dev = \
sps.lil_matrix(self._std_dev.flatten(), self._std_dev.shape)
self._std_dev = sps.lil_matrix(self._std_dev.flatten(),
self._std_dev.shape)
self._sparse = True
@ -776,11 +774,11 @@ class Tally(IDManagerMixin):
other_sum = other_copy.get_reshaped_data(value='sum')
if join_right:
merged_sum = \
np.concatenate((self_sum, other_sum), axis=merge_axis)
merged_sum = np.concatenate((self_sum, other_sum),
axis=merge_axis)
else:
merged_sum = \
np.concatenate((other_sum, self_sum), axis=merge_axis)
merged_sum = np.concatenate((other_sum, self_sum),
axis=merge_axis)
merged_tally._sum = np.reshape(merged_sum, merged_tally.shape)
@ -790,11 +788,11 @@ class Tally(IDManagerMixin):
other_sum_sq = other_copy.get_reshaped_data(value='sum_sq')
if join_right:
merged_sum_sq = \
np.concatenate((self_sum_sq, other_sum_sq), axis=merge_axis)
merged_sum_sq = np.concatenate((self_sum_sq, other_sum_sq),
axis=merge_axis)
else:
merged_sum_sq = \
np.concatenate((other_sum_sq, self_sum_sq), axis=merge_axis)
merged_sum_sq = np.concatenate((other_sum_sq, self_sum_sq),
axis=merge_axis)
merged_tally._sum_sq = np.reshape(merged_sum_sq, merged_tally.shape)
@ -804,11 +802,11 @@ class Tally(IDManagerMixin):
other_mean = other_copy.get_reshaped_data(value='mean')
if join_right:
merged_mean = \
np.concatenate((self_mean, other_mean), axis=merge_axis)
merged_mean = np.concatenate((self_mean, other_mean),
axis=merge_axis)
else:
merged_mean = \
np.concatenate((other_mean, self_mean), axis=merge_axis)
merged_mean = np.concatenate((other_mean, self_mean),
axis=merge_axis)
merged_tally._mean = np.reshape(merged_mean, merged_tally.shape)
@ -818,11 +816,11 @@ class Tally(IDManagerMixin):
other_std_dev = other_copy.get_reshaped_data(value='std_dev')
if join_right:
merged_std_dev = \
np.concatenate((self_std_dev, other_std_dev), axis=merge_axis)
merged_std_dev = np.concatenate((self_std_dev, other_std_dev),
axis=merge_axis)
else:
merged_std_dev = \
np.concatenate((other_std_dev, self_std_dev), axis=merge_axis)
merged_std_dev = np.concatenate((other_std_dev, self_std_dev),
axis=merge_axis)
merged_tally._std_dev = np.reshape(merged_std_dev, merged_tally.shape)
@ -1003,35 +1001,6 @@ class Tally(IDManagerMixin):
return filter_found
def get_filter_index(self, filter_type, filter_bin):
"""Returns the index in the Tally's results array for a Filter bin
Parameters
----------
filter_type : openmc.FilterMeta
Type of the filter, e.g. MeshFilter
filter_bin : int or tuple
The bin is an integer ID for 'material', 'surface', 'cell',
'cellborn', and 'universe' Filters. The bin is an integer for the
cell instance ID for 'distribcell' Filters. The bin is a 2-tuple of
floats for 'energy' and 'energyout' filters corresponding to the
energy boundaries of the bin of interest. The bin is a (x,y,z)
3-tuple for 'mesh' filters corresponding to the mesh cell of
interest.
Returns
-------
The index in the Tally data array for this filter bin
"""
# Find the equivalent Filter in this Tally's list of Filters
filter_found = self.find_filter(filter_type)
# Get the index for the requested bin from the Filter and return it
filter_index = filter_found.get_bin_index(filter_bin)
return filter_index
def get_nuclide_index(self, nuclide):
"""Returns the index in the Tally's results array for a Nuclide bin
@ -1151,48 +1120,28 @@ class Tally(IDManagerMixin):
# Loop over all of the Tally's Filters
for i, self_filter in enumerate(self.filters):
user_filter = False
# If a user-requested Filter, get the user-requested bins
for j, test_filter in enumerate(filters):
if type(self_filter) is test_filter:
bins = filter_bins[j]
user_filter = True
break
else:
# If not a user-requested Filter, get all bins
if isinstance(self_filter, openmc.DistribcellFilter):
# Create list of cell instance IDs for distribcell Filters
bins = list(range(self_filter.num_bins))
# If not a user-requested Filter, get all bins
if not user_filter:
# Create list of 2- or 3-tuples tuples for mesh cell bins
if isinstance(self_filter, openmc.MeshFilter):
bins = list(self_filter.mesh.indices)
# Create list of 2-tuples for energy boundary bins
elif isinstance(self_filter, (openmc.EnergyFilter,
openmc.EnergyoutFilter, openmc.MuFilter,
openmc.PolarFilter, openmc.AzimuthalFilter)):
bins = []
for k in range(self_filter.num_bins):
bins.append((self_filter.bins[k], self_filter.bins[k+1]))
# Create list of cell instance IDs for distribcell Filters
elif isinstance(self_filter, openmc.DistribcellFilter):
bins = [b for b in range(self_filter.num_bins)]
# EnergyFunctionFilters don't have bins so just add a None
elif isinstance(self_filter, openmc.EnergyFunctionFilter):
# EnergyFunctionFilters don't have bins so just add a None
bins = [None]
# Create list of IDs for bins for all other filter types
else:
# Create list of IDs for bins for all other filter types
bins = self_filter.bins
# Initialize a NumPy array for the Filter bin indices
filter_indices.append(np.zeros(len(bins), dtype=np.int))
# Add indices for each bin in this Filter to the list
for j, bin in enumerate(bins):
filter_index = self.get_filter_index(type(self_filter), bin)
filter_indices[i][j] = filter_index
indices = np.array([self_filter.get_bin_index(b) for b in bins])
filter_indices.append(indices)
# Account for stride in each of the previous filters
for indices in filter_indices[:i]:
@ -1956,14 +1905,14 @@ class Tally(IDManagerMixin):
elif isinstance(filter1, openmc.EnergyFunctionFilter):
filter1_bins = [None]
else:
filter1_bins = [filter1.get_bin(i) for i in range(filter1.num_bins)]
filter1_bins = filter1.bins
if isinstance(filter2, openmc.DistribcellFilter):
filter2_bins = [b for b in range(filter2.num_bins)]
elif isinstance(filter2, openmc.EnergyFunctionFilter):
filter2_bins = [None]
else:
filter2_bins = [filter2.get_bin(i) for i in range(filter2.num_bins)]
filter2_bins = filter2.bins
# Create variables to store views of data in the misaligned structure
mean = {}
@ -2604,7 +2553,8 @@ class Tally(IDManagerMixin):
new_tally = self * -1
return new_tally
def get_slice(self, scores=[], filters=[], filter_bins=[], nuclides=[]):
def get_slice(self, scores=[], filters=[], filter_bins=[], nuclides=[],
squeeze=False):
"""Build a sliced tally for the specified filters, scores and nuclides.
This method constructs a new tally to encapsulate a subset of the data
@ -2615,26 +2565,26 @@ class Tally(IDManagerMixin):
Parameters
----------
scores : list of str
A list of one or more score strings
(e.g., ['absorption', 'nu-fission']; default is [])
A list of one or more score strings (e.g., ['absorption',
'nu-fission']
filters : Iterable of openmc.FilterMeta
An iterable of filter types
(e.g., [MeshFilter, EnergyFilter]; default is [])
An iterable of filter types (e.g., [MeshFilter, EnergyFilter])
filter_bins : list of Iterables
A list of tuples of filter bins corresponding to the filter_types
parameter (e.g., [(1,), ((0., 0.625e-6),)]; default is []). Each
tuple contains bins to slice for the corresponding filter type in
the filters parameter. Each bins is the integer ID for 'material',
A list of iterables of filter bins corresponding to the specified
filter types (e.g., [(1,), ((0., 0.625e-6),)]). Each iterable
contains bins to slice for the corresponding filter type in the
filters parameter. Each bin is the integer ID for 'material',
'surface', 'cell', 'cellborn', and 'universe' Filters. Each bin is
an integer for the cell instance ID for 'distribcell' Filters. Each
bin is a 2-tuple of floats for 'energy' and 'energyout' filters
corresponding to the energy boundaries of the bin of interest. The
bin is an (x,y,z) 3-tuple for 'mesh' filters corresponding to the
mesh cell of interest. The order of the bins in the list must
correspond to the filter_types parameter.
correspond to the `filters` argument.
nuclides : list of str
A list of nuclide name strings
(e.g., ['U235', 'U238']; default is [])
A list of nuclide name strings (e.g., ['U235', 'U238'])
squeeze : bool
Whether to remove filters with only a single bin in the sliced tally
Returns
-------
@ -2714,32 +2664,29 @@ class Tally(IDManagerMixin):
# Determine the filter indices from any of the requested filters
for i, filter_type in enumerate(filters):
find_filter = new_tally.find_filter(filter_type)
f = new_tally.find_filter(filter_type)
# Remove filters with only a single bin if requested
if squeeze:
if len(filter_bins[i]) == 1:
new_tally.filters.remove(f)
continue
else:
raise RuntimeError('Cannot remove sliced filter with '
'more than one bin.')
# Remove and/or reorder filter bins to user specifications
bin_indices = []
bin_indices = [f.get_bin_index(b)
for b in filter_bins[i]]
bin_indices = np.unique(bin_indices)
for filter_bin in filter_bins[i]:
bin_index = find_filter.get_bin_index(filter_bin)
if issubclass(filter_type, openmc.RealFilter):
bin_indices.extend([bin_index, bin_index+1])
else:
bin_indices.append(bin_index)
# Set bins for mesh/distribcell filters apart from others
if filter_type is openmc.MeshFilter:
bins = find_filter.mesh
elif filter_type is openmc.DistribcellFilter:
bins = find_filter.bins
else:
bins = np.unique(find_filter.bins[bin_indices])
# Create new filter
new_filter = filter_type(bins)
# Set bins for sliced filter
new_filter = copy.copy(f)
new_filter.bins = [f.bins[i] for i in bin_indices]
# Set number of bins manually for mesh/distribcell filters
if filter_type is openmc.DistribcellFilter:
new_filter._num_bins = find_filter._num_bins
new_filter._num_bins = f._num_bins
# Replace existing filter with new one
for j, test_filter in enumerate(new_tally.filters):
@ -2816,9 +2763,7 @@ class Tally(IDManagerMixin):
elif isinstance(find_filter, openmc.EnergyFunctionFilter):
filter_bins = [None]
else:
num_bins = find_filter.num_bins
filter_bins = \
[(find_filter.get_bin(i)) for i in range(num_bins)]
filter_bins = find_filter.bins
# Only sum across bins specified by the user
else:
@ -2970,9 +2915,7 @@ class Tally(IDManagerMixin):
elif isinstance(find_filter, openmc.EnergyFunctionFilter):
filter_bins = [None]
else:
num_bins = find_filter.num_bins
filter_bins = \
[(find_filter.get_bin(i)) for i in range(num_bins)]
filter_bins = find_filter.bins
# Only average across bins specified by the user
else:

8
src/constants.F90
View file

@ -358,7 +358,7 @@ module constants
integer, parameter :: NO_BIN_FOUND = -1
! Tally filter and map types
integer, parameter :: N_FILTER_TYPES = 16
integer, parameter :: N_FILTER_TYPES = 20
integer, parameter :: &
FILTER_UNIVERSE = 1, &
FILTER_MATERIAL = 2, &
@ -375,7 +375,11 @@ module constants
FILTER_DELAYEDGROUP = 13, &
FILTER_ENERGYFUNCTION = 14, &
FILTER_CELLFROM = 15, &
FILTER_MESHSURFACE = 16
FILTER_MESHSURFACE = 16, &
FILTER_LEGENDRE = 17, &
FILTER_SPH_HARMONICS = 18, &
FILTER_SPTL_LEGENDRE = 19, &
FILTER_ZERNIKE = 20
! Mesh types
integer, parameter :: &

109
src/math.F90
View file

@ -181,7 +181,7 @@ contains
end function calc_pn
!===============================================================================
! CALC_RN calculates the n-th order spherical harmonics for a given angle
! CALC_RN calculates the n-th order real spherical harmonics for a given angle
! (in terms of (u,v,w)). All Rn,m values are provided (where -n<=m<=n)
!===============================================================================
@ -574,6 +574,113 @@ contains
end function calc_rn
!===============================================================================
! CALC_ZN calculates the n-th order modified Zernike polynomial moment for a
! given angle (rho, theta) location in the unit disk. The normlization of the
! polynomials is such that the integral of Z_pq*Z_pq over the unit disk is
! exactly pi
!===============================================================================
subroutine calc_zn(n, rho, phi, zn)
! This procedure uses the modified Kintner's method for calculating Zernike
! polynomials as outlined in Chong, C. W., Raveendran, P., & Mukundan,
! R. (2003). A comparative analysis of algorithms for fast computation of
! Zernike moments. Pattern Recognition, 36(3), 731-742.
integer, intent(in) :: n ! Maximum order
real(8), intent(in) :: rho ! Radial location in the unit disk
real(8), intent(in) :: phi ! Theta (radians) location in the unit disk
real(8), intent(out) :: zn(:) ! The resulting list of coefficients
real(8) :: sin_phi, cos_phi ! Sine and Cosine of phi
real(8) :: sin_phi_vec(n+1) ! Contains sin(n*phi)
real(8) :: cos_phi_vec(n+1) ! Contains cos(n*phi)
real(8) :: zn_mat(n+1, n+1) ! Matrix form of the coefficients which is
! easier to work with
real(8) :: k1, k2, k3, k4 ! Variables for R_m_n calculation
real(8) :: sqrt_norm ! normalization for radial moments
integer :: i,p,q ! Loop counters
real(8), parameter :: SQRT_N_1(0:10) = [&
sqrt(1.0_8), sqrt(2.0_8), sqrt(3.0_8), sqrt(4.0_8), &
sqrt(5.0_8), sqrt(6.0_8), sqrt(7.0_8), sqrt(8.0_8), &
sqrt(9.0_8), sqrt(10.0_8), sqrt(11.0_8)]
real(8), parameter :: SQRT_2N_2(0:10) = SQRT_N_1*sqrt(2.0_8)
! n == radial degree
! m == azimuthal frequency
! ==========================================================================
! Determine vector of sin(n*phi) and cos(n*phi). This takes advantage of the
! following recurrence relations so that only a single sin/cos have to be
! evaluated (http://mathworld.wolfram.com/Multiple-AngleFormulas.html)
!
! sin(nx) = 2 cos(x) sin((n-1)x) - sin((n-2)x)
! cos(nx) = 2 cos(x) cos((n-1)x) - cos((n-2)x)
sin_phi = sin(phi)
cos_phi = cos(phi)
sin_phi_vec(1) = 1.0_8
cos_phi_vec(1) = 1.0_8
sin_phi_vec(2) = 2.0_8 * cos_phi
cos_phi_vec(2) = cos_phi
do i = 3, n+1
sin_phi_vec(i) = 2.0_8 * cos_phi * sin_phi_vec(i-1) - sin_phi_vec(i-2)
cos_phi_vec(i) = 2.0_8 * cos_phi * cos_phi_vec(i-1) - cos_phi_vec(i-2)
end do
do i = 1, n+1
sin_phi_vec(i) = sin_phi_vec(i) * sin_phi
end do
! ==========================================================================
! Calculate R_pq(rho)
! Fill the main diagonal first (Eq. 3.9 in Chong)
do p = 0, n
zn_mat(p+1, p+1) = rho**p
end do
! Fill in the second diagonal (Eq. 3.10 in Chong)
do q = 0, n-2
zn_mat(q+2+1, q+1) = (q+2) * zn_mat(q+2+1, q+2+1) - (q+1) * zn_mat(q+1, q+1)
end do
! Fill in the rest of the values using the original results (Eq. 3.8 in Chong)
do p = 4, n
k2 = 2 * p * (p - 1) * (p - 2)
do q = p-4, 0, -2
k1 = (p + q) * (p - q) * (p - 2) / 2
k3 = -q**2*(p - 1) - p * (p - 1) * (p - 2)
k4 = -p * (p + q - 2) * (p - q - 2) / 2
zn_mat(p+1, q+1) = ((k2 * rho**2 + k3) * zn_mat(p-2+1, q+1) + k4 * zn_mat(p-4+1, q+1)) / k1
end do
end do
! Roll into a single vector for easier computation later
! The vector is ordered (0,0), (1,-1), (1,1), (2,-2), (2,0),
! (2, 2), .... in (n,m) indices
! Note that the cos and sin vectors are offset by one
! sin_phi_vec = [sin(x), sin(2x), sin(3x) ...]
! cos_phi_vec = [1.0, cos(x), cos(2x)... ]
i = 1
do p = 0, n
do q = -p, p, 2
if (q < 0) then
zn(i) = zn_mat(p+1, abs(q)+1) * sin_phi_vec(abs(q)) * SQRT_2N_2(p)
else if (q == 0) then
zn(i) = zn_mat(p+1, q+1) * SQRT_N_1(p)
else
zn(i) = zn_mat(p+1, q+1) * cos_phi_vec(abs(q)+1) * SQRT_2N_2(p)
end if
i = i + 1
end do
end do
end subroutine calc_zn
!===============================================================================
! EXPAND_HARMONIC expands a given series of real spherical harmonics
!===============================================================================

116
src/tallies/tally_filter.F90
View file

@ -17,13 +17,17 @@ module tally_filter
use tally_filter_distribcell
use tally_filter_energy
use tally_filter_energyfunc
use tally_filter_legendre
use tally_filter_material
use tally_filter_mesh
use tally_filter_meshsurface
use tally_filter_mu
use tally_filter_polar
use tally_filter_sph_harm
use tally_filter_sptl_legendre
use tally_filter_surface
use tally_filter_universe
use tally_filter_zernike
implicit none
@ -42,60 +46,60 @@ contains
integer :: i
character(20) :: type_
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
! Get type as a Fortran string
select type (f => filters(index) % obj)
type is (AzimuthalFilter)
type_ = 'azimuthal'
type is (CellFilter)
type_ = 'cell'
type is (CellbornFilter)
type_ = 'cellborn'
type is (CellfromFilter)
type_ = 'cellfrom'
type is (DelayedGroupFilter)
type_ = 'delayedgroup'
type is (DistribcellFilter)
type_ = 'distribcell'
type is (EnergyFilter)
type_ = 'energy'
type is (EnergyoutFilter)
type_ = 'energyout'
type is (EnergyFunctionFilter)
type_ = 'energyfunction'
type is (MaterialFilter)
type_ = 'material'
type is (MeshFilter)
type_ = 'mesh'
type is (MeshSurfaceFilter)
type_ = 'meshsurface'
type is (MuFilter)
type_ = 'mu'
type is (PolarFilter)
type_ = 'polar'
type is (SurfaceFilter)
type_ = 'surface'
type is (UniverseFilter)
type_ = 'universe'
end select
err = verify_filter(index)
if (err == 0) then
! Get type as a Fortran string
select type (f => filters(index) % obj)
type is (AzimuthalFilter)
type_ = 'azimuthal'
type is (CellFilter)
type_ = 'cell'
type is (CellbornFilter)
type_ = 'cellborn'
type is (CellfromFilter)
type_ = 'cellfrom'
type is (DelayedGroupFilter)
type_ = 'delayedgroup'
type is (DistribcellFilter)
type_ = 'distribcell'
type is (EnergyFilter)
type_ = 'energy'
type is (EnergyoutFilter)
type_ = 'energyout'
type is (EnergyFunctionFilter)
type_ = 'energyfunction'
type is (LegendreFilter)
type_ = 'legendre'
type is (MaterialFilter)
type_ = 'material'
type is (MeshFilter)
type_ = 'mesh'
type is (MeshSurfaceFilter)
type_ = 'meshsurface'
type is (MuFilter)
type_ = 'mu'
type is (PolarFilter)
type_ = 'polar'
type is (SphericalHarmonicsFilter)
type_ = 'sphericalharmonics'
type is (SpatialLegendreFilter)
type_ = 'spatiallegendre'
type is (SurfaceFilter)
type_ = 'surface'
type is (UniverseFilter)
type_ = 'universe'
type is (ZernikeFilter)
type_ = 'zernike'
end select
! Convert Fortran string to null-terminated C string. We assume the
! caller has allocated a char array buffer
do i = 1, len_trim(type_)
type(i) = type_(i:i)
end do
type(len_trim(type_) + 1) = C_NULL_CHAR
err = 0
else
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array is out of bounds.")
! Convert Fortran string to null-terminated C string. We assume the
! caller has allocated a char array buffer
do i = 1, len_trim(type_)
type(i) = type_(i:i)
end do
type(len_trim(type_) + 1) = C_NULL_CHAR
end if
end function openmc_filter_get_type
@ -135,6 +139,8 @@ contains
allocate(EnergyoutFilter :: filters(index) % obj)
case ('energyfunction')
allocate(EnergyFunctionFilter :: filters(index) % obj)
case ('legendre')
allocate(LegendreFilter :: filters(index) % obj)
case ('material')
allocate(MaterialFilter :: filters(index) % obj)
case ('mesh')
@ -145,10 +151,16 @@ contains
allocate(MuFilter :: filters(index) % obj)
case ('polar')
allocate(PolarFilter :: filters(index) % obj)
case ('sphericalharmonics')
allocate(SphericalHarmonicsFilter :: filters(index) % obj)
case ('spatiallegendre')
allocate(SpatialLegendreFilter :: filters(index) % obj)
case ('surface')
allocate(SurfaceFilter :: filters(index) % obj)
case ('universe')
allocate(UniverseFilter :: filters(index) % obj)
case ('zernike')
allocate(ZernikeFilter :: filters(index) % obj)
case default
err = E_UNASSIGNED
call set_errmsg("Unknown filter type: " // trim(type_))

75
src/tallies/tally_filter_energy.F90
View file

@ -204,28 +204,21 @@ contains
integer(C_INT32_T), intent(out) :: n
integer(C_INT) :: err
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
select type (f => filters(index) % obj)
type is (EnergyFilter)
energies = C_LOC(f % bins)
n = size(f % bins)
err = 0
type is (EnergyoutFilter)
energies = C_LOC(f % bins)
n = size(f % bins)
err = 0
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (EnergyFilter)
energies = C_LOC(f % bins)
n = size(f % bins)
err = 0
type is (EnergyoutFilter)
energies = C_LOC(f % bins)
n = size(f % bins)
err = 0
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to get energy bins on a non-energy filter.")
end select
else
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
err = E_INVALID_TYPE
call set_errmsg("Tried to get energy bins on a non-energy filter.")
end select
end if
end function openmc_energy_filter_get_bins
@ -237,31 +230,23 @@ contains
real(C_DOUBLE), intent(in) :: energies(n)
integer(C_INT) :: err
err = 0
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
select type (f => filters(index) % obj)
type is (EnergyFilter)
f % n_bins = n - 1
if (allocated(f % bins)) deallocate(f % bins)
allocate(f % bins(n))
f % bins(:) = energies
type is (EnergyoutFilter)
f % n_bins = n - 1
if (allocated(f % bins)) deallocate(f % bins)
allocate(f % bins(n))
f % bins(:) = energies
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (EnergyFilter)
f % n_bins = n - 1
if (allocated(f % bins)) deallocate(f % bins)
allocate(f % bins(n))
f % bins(:) = energies
type is (EnergyoutFilter)
f % n_bins = n - 1
if (allocated(f % bins)) deallocate(f % bins)
allocate(f % bins(n))
f % bins(:) = energies
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to get energy bins on a non-energy filter.")
end select
else
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
err = E_INVALID_TYPE
call set_errmsg("Tried to get energy bins on a non-energy filter.")
end select
end if
end function openmc_energy_filter_set_bins

22
src/tallies/tally_filter_header.F90
View file

@ -15,6 +15,7 @@ module tally_filter_header
implicit none
private
public :: free_memory_tally_filter
public :: verify_filter
public :: openmc_extend_filters
public :: openmc_filter_get_id
public :: openmc_filter_set_id
@ -148,6 +149,27 @@ contains
largest_filter_id = 0
end subroutine free_memory_tally_filter
!===============================================================================
! VERIFY_FILTER makes sure that given a filter index, the size of the filters
! array is sufficient and a filter object has already been allocated.
!===============================================================================
function verify_filter(index) result(err)
integer(C_INT32_T), intent(in) :: index
integer(C_INT) :: err
err = 0
if (index >= 1 .and. index <= n_filters) then
if (.not. allocated(filters(index) % obj)) then
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
end if
end function verify_filter
!===============================================================================
! C API FUNCTIONS
!===============================================================================

125
src/tallies/tally_filter_legendre.F90 Normal file
View file

@ -0,0 +1,125 @@
module tally_filter_legendre
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
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
public :: openmc_legendre_filter_get_order
public :: openmc_legendre_filter_set_order
!===============================================================================
! 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
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(8) :: wgt
! TODO: Use recursive formula to calculate higher orders
do i = 0, this % order
wgt = calc_pn(i, p % mu)
call match % bins % push_back(i + 1)
call match % weights % push_back(wgt)
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
!===============================================================================
function openmc_legendre_filter_get_order(index, order) result(err) bind(C)
! Get the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), intent(out) :: order
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (LegendreFilter)
order = f % order
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to get order on a non-expansion filter.")
end select
end if
end function openmc_legendre_filter_get_order
function openmc_legendre_filter_set_order(index, order) result(err) bind(C)
! Set the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), value :: order
integer(C_INT) :: err
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
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to set order on a non-expansion filter.")
end select
end if
end function openmc_legendre_filter_set_order
end module tally_filter_legendre

71
src/tallies/tally_filter_material.F90
View file

@ -126,25 +126,18 @@ contains
integer(C_INT32_T), intent(out) :: n
integer(C_INT) :: err
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
select type (f => filters(index) % obj)
type is (MaterialFilter)
bins = C_LOC(f % materials)
n = size(f % materials)
err = 0
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (MaterialFilter)
bins = C_LOC(f % materials)
n = size(f % materials)
err = 0
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to get material filter bins on a &
&non-material filter.")
end select
else
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
err = E_INVALID_TYPE
call set_errmsg("Tried to get material filter bins on a &
&non-material filter.")
end select
end if
end function openmc_material_filter_get_bins
@ -158,34 +151,26 @@ contains
integer :: i
err = 0
if (index >= 1 .and. index <= n_filters) then
if (allocated(filters(index) % obj)) then
select type (f => filters(index) % obj)
type is (MaterialFilter)
f % n_bins = n
if (allocated(f % materials)) deallocate(f % materials)
allocate(f % materials(n))
f % materials(:) = bins
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (MaterialFilter)
f % n_bins = n
if (allocated(f % materials)) deallocate(f % materials)
allocate(f % materials(n))
f % materials(:) = bins
! Generate mapping from material indices to filter bins.
call f % map % clear()
do i = 1, n
call f % map % set(f % materials(i), i)
end do
! Generate mapping from material indices to filter bins.
call f % map % clear()
do i = 1, n
call f % map % set(f % materials(i), i)
end do
class default
err = E_INVALID_TYPE
call set_errmsg("Tried to set material filter bins on a &
&non-material filter.")
end select
else
err = E_ALLOCATE
call set_errmsg("Filter type has not been set yet.")
end if
else
err = E_OUT_OF_BOUNDS
call set_errmsg("Index in filters array out of bounds.")
err = E_INVALID_TYPE
call set_errmsg("Tried to set material filter bins on a &
&non-material filter.")
end select
end if
end function openmc_material_filter_set_bins

242
src/tallies/tally_filter_sph_harm.F90 Normal file
View file

@ -0,0 +1,242 @@
module tally_filter_sph_harm
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
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
public :: openmc_sphharm_filter_get_order
public :: openmc_sphharm_filter_get_cosine
public :: openmc_sphharm_filter_set_order
public :: openmc_sphharm_filter_set_cosine
integer, public, parameter :: COSINE_SCATTER = 1
integer, public, parameter :: COSINE_PARTICLE = 2
!===============================================================================
! SPHERICALHARMONICSFILTER gives spherical harmonics expansion moments of a
! tally score
!===============================================================================
type, public, extends(TallyFilter) :: SphericalHarmonicsFilter
integer :: order
integer :: cosine
contains
procedure :: from_xml
procedure :: get_all_bins
procedure :: to_statepoint
procedure :: text_label
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(8) :: wgt
real(8) :: rn(2*this % order + 1)
! TODO: Use recursive formula to calculate higher orders
j = 0
do n = 0, this % order
! Determine cosine term for scatter expansion if necessary
if (this % cosine == COSINE_SCATTER) then
wgt = calc_pn(n, p % mu)
else
wgt = ONE
end if
! Calculate n-th order spherical harmonics for (u,v,w)
num_nm = 2*n + 1
rn(1:num_nm) = calc_rn(n, p % last_uvw)
! 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 * rn(i))
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
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_sphharm_filter_get_order(index, order) result(err) bind(C)
! Get the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), intent(out) :: order
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (SphericalHarmonicsFilter)
order = f % order
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spherical harmonics filter.")
end select
end if
end function openmc_sphharm_filter_get_order
function openmc_sphharm_filter_get_cosine(index, cosine) result(err) bind(C)
! Get the order of an expansion filter
integer(C_INT32_T), value :: index
character(kind=C_CHAR), intent(out) :: cosine(*)
integer(C_INT) :: err
integer :: i
character(10) :: cosine_
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (SphericalHarmonicsFilter)
select case (f % cosine)
case (COSINE_SCATTER)
cosine_ = 'scatter'
case (COSINE_PARTICLE)
cosine_ = 'particle'
end select
! Convert to C string
do i = 1, len_trim(cosine_)
cosine(i) = cosine_(i:i)
end do
cosine(len_trim(cosine_) + 1) = C_NULL_CHAR
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spherical harmonics filter.")
end select
end if
end function openmc_sphharm_filter_get_cosine
function openmc_sphharm_filter_set_order(index, order) result(err) bind(C)
! Set the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), value :: order
integer(C_INT) :: err
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
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spherical harmonics filter.")
end select
end if
end function openmc_sphharm_filter_set_order
function openmc_sphharm_filter_set_cosine(index, cosine) result(err) bind(C)
! Set the cosine parameter
integer(C_INT32_T), value :: index
character(kind=C_CHAR), intent(in) :: cosine(*)
integer(C_INT) :: err
character(:), allocatable :: cosine_
! Convert C string to Fortran string
cosine_ = to_f_string(cosine)
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (SphericalHarmonicsFilter)
select case (cosine_)
case ('scatter')
f % cosine = COSINE_SCATTER
case ('particle')
f % cosine = COSINE_PARTICLE
end select
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spherical harmonics filter.")
end select
end if
end function openmc_sphharm_filter_set_cosine
end module tally_filter_sph_harm

228
src/tallies/tally_filter_sptl_legendre.F90 Normal file
View file

@ -0,0 +1,228 @@
module tally_filter_sptl_legendre
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use error
use hdf5_interface
use math, only: calc_pn
use particle_header, only: Particle
use string, only: to_str, to_lower
use tally_filter_header
use xml_interface
implicit none
private
public :: openmc_spatial_legendre_filter_get_order
public :: openmc_spatial_legendre_filter_get_params
public :: openmc_spatial_legendre_filter_set_order
public :: openmc_spatial_legendre_filter_set_params
integer, parameter :: AXIS_X = 1
integer, parameter :: AXIS_Y = 2
integer, parameter :: AXIS_Z = 3
!===============================================================================
! SPATIALLEGENDREFILTER gives Legendre moments of the particle's normalized
! position along an axis
!===============================================================================
type, public, extends(TallyFilter) :: SpatialLegendreFilter
integer(C_INT) :: order
integer(C_INT) :: axis
real(C_DOUBLE) :: min
real(C_DOUBLE) :: max
contains
procedure :: from_xml
procedure :: get_all_bins
procedure :: to_statepoint
procedure :: text_label
end type SpatialLegendreFilter
contains
!===============================================================================
! SpatialLegendreFilter methods
!===============================================================================
subroutine from_xml(this, node)
class(SpatialLegendreFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
character(MAX_WORD_LEN) :: axis
! Get attributes from XML
call get_node_value(node, "order", this % order)
call get_node_value(node, "axis", axis)
select case (to_lower(axis))
case ('x')
this % axis = AXIS_X
case ('y')
this % axis = AXIS_Y
case ('z')
this % axis = AXIS_Z
end select
call get_node_value(node, "min", this % min)
call get_node_value(node, "max", this % max)
this % n_bins = this % order + 1
end subroutine from_xml
subroutine get_all_bins(this, p, estimator, match)
class(SpatialLegendreFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: i
real(8) :: wgt
real(8) :: x ! Position on specified axis
real(8) :: x_norm ! Normalized position
x = p % coord(1) % xyz(this % axis)
if (this % min <= x .and. x <= this % max) then
! Calculate normalized position between min and max
x_norm = TWO*(x - this % min)/(this % max - this % min) - ONE
! TODO: Use recursive formula to calculate higher orders
do i = 0, this % order
wgt = calc_pn(i, x_norm)
call match % bins % push_back(i + 1)
call match % weights % push_back(wgt)
end do
end if
end subroutine get_all_bins
subroutine to_statepoint(this, filter_group)
class(SpatialLegendreFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
character(kind=C_CHAR) :: axis
call write_dataset(filter_group, "type", "spatiallegendre")
call write_dataset(filter_group, "n_bins", this % n_bins)
call write_dataset(filter_group, "order", this % order)
select case (this % axis)
case (AXIS_X)
axis = 'x'
case (AXIS_Y)
axis = 'y'
case (AXIS_Z)
axis = 'z'
end select
call write_dataset(filter_group, 'axis', axis)
call write_dataset(filter_group, 'min', this % min)
call write_dataset(filter_group, 'max', this % max)
end subroutine to_statepoint
function text_label(this, bin) result(label)
class(SpatialLegendreFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
character(1) :: axis
select case (this % axis)
case (AXIS_X)
axis = 'x'
case (AXIS_Y)
axis = 'y'
case (AXIS_Z)
axis = 'z'
end select
label = "Legendre expansion, " // axis // " axis, P" // &
trim(to_str(bin - 1))
end function text_label
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_spatial_legendre_filter_get_order(index, order) result(err) bind(C)
! Get the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), intent(out) :: order
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (SpatialLegendreFilter)
order = f % order
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spatial Legendre filter.")
end select
end if
end function openmc_spatial_legendre_filter_get_order
function openmc_spatial_legendre_filter_set_order(index, order) result(err) bind(C)
! Set the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), value :: order
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (SpatialLegendreFilter)
f % order = order
f % n_bins = order + 1
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spatial Legendre filter.")
end select
end if
end function openmc_spatial_legendre_filter_set_order
function openmc_spatial_legendre_filter_get_params(index, axis, min, max) &
result(err) bind(C)
! Get the parameters for a spatial Legendre filter
integer(C_INT32_T), value :: index
integer(C_INT), intent(out) :: axis
real(C_DOUBLE), intent(out) :: min
real(C_DOUBLE), intent(out) :: max
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type(f => filters(index) % obj)
type is (SpatialLegendreFilter)
axis = f % axis
min = f % min
max = f % max
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spatial Legendre filter.")
end select
end if
end function openmc_spatial_legendre_filter_get_params
function openmc_spatial_legendre_filter_set_params(index, axis, min, max) &
result(err) bind(C)
! Set the parameters for a spatial Legendre filter
integer(C_INT32_T), value :: index
integer(C_INT), intent(in), optional :: axis
real(C_DOUBLE), intent(in), optional :: min
real(C_DOUBLE), intent(in), optional :: max
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type(f => filters(index) % obj)
type is (SpatialLegendreFilter)
if (present(axis)) f % axis = axis
if (present(min)) f % min = min
if (present(max)) f % max = max
class default
err = E_INVALID_TYPE
call set_errmsg("Not a spatial Legendre filter.")
end select
end if
end function openmc_spatial_legendre_filter_set_params
end module tally_filter_sptl_legendre

203
src/tallies/tally_filter_zernike.F90 Normal file
View file

@ -0,0 +1,203 @@
module tally_filter_zernike
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T
use constants
use error
use hdf5_interface
use math, only: calc_zn
use particle_header, only: Particle
use string, only: to_str
use tally_filter_header
use xml_interface
implicit none
private
!===============================================================================
! ZERNIKEFILTER gives Zernike polynomial moments of a particle's position
!===============================================================================
type, public, extends(TallyFilter) :: ZernikeFilter
integer :: order
real(8) :: x
real(8) :: y
real(8) :: r
contains
procedure :: from_xml
procedure :: get_all_bins
procedure :: to_statepoint
procedure :: text_label
end type ZernikeFilter
contains
!===============================================================================
! ZernikeFilter methods
!===============================================================================
subroutine from_xml(this, node)
class(ZernikeFilter), intent(inout) :: this
type(XMLNode), intent(in) :: node
integer :: n
! Get center of cylinder and radius
call get_node_value(node, "x", this % x)
call get_node_value(node, "y", this % y)
call get_node_value(node, "r", this % r)
! 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
subroutine get_all_bins(this, p, estimator, match)
class(ZernikeFilter), intent(in) :: this
type(Particle), intent(in) :: p
integer, intent(in) :: estimator
type(TallyFilterMatch), intent(inout) :: match
integer :: i
real(8) :: x, y, r, theta
real(8) :: zn(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
theta = atan2(y, x)
! 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
end subroutine get_all_bins
subroutine to_statepoint(this, filter_group)
class(ZernikeFilter), intent(in) :: this
integer(HID_T), intent(in) :: filter_group
call write_dataset(filter_group, "type", "zernike")
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
function text_label(this, bin) result(label)
class(ZernikeFilter), intent(in) :: this
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
integer :: n, m
integer :: first, last
do n = 0, this % order
last = (n + 1)*(n + 2)/2
if (bin <= last) then
first = last - n
m = -n + (bin - first)*2
label = "Zernike expansion, Z" // trim(to_str(n)) // "," &
// trim(to_str(m))
exit
end if
end do
end function text_label
!===============================================================================
! C API FUNCTIONS
!===============================================================================
function openmc_zernike_filter_get_order(index, order) result(err) bind(C)
! Get the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), intent(out) :: order
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (ZernikeFilter)
order = f % order
class default
err = E_INVALID_TYPE
call set_errmsg("Not a Zernike filter.")
end select
end if
end function openmc_zernike_filter_get_order
function openmc_zernike_filter_get_params(index, x, y, r) result(err) bind(C)
! Get the Zernike filter parameters
integer(C_INT32_T), value :: index
real(C_DOUBLE), intent(out) :: x
real(C_DOUBLE), intent(out) :: y
real(C_DOUBLE), intent(out) :: r
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (ZernikeFilter)
x = f % x
y = f % y
r = f % r
class default
err = E_INVALID_TYPE
call set_errmsg("Not a Zernike filter.")
end select
end if
end function openmc_zernike_filter_get_params
function openmc_zernike_filter_set_order(index, order) result(err) bind(C)
! Set the order of an expansion filter
integer(C_INT32_T), value :: index
integer(C_INT), value :: order
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (ZernikeFilter)
f % order = order
f % n_bins = ((order + 1)*(order + 2))/2
class default
err = E_INVALID_TYPE
call set_errmsg("Not a Zernike filter.")
end select
end if
end function openmc_zernike_filter_set_order
function openmc_zernike_filter_set_params(index, x, y, r) result(err) bind(C)
! Set the Zernike filter parameters
integer(C_INT32_T), value :: index
real(C_DOUBLE), intent(in), optional :: x
real(C_DOUBLE), intent(in), optional :: y
real(C_DOUBLE), intent(in), optional :: r
integer(C_INT) :: err
err = verify_filter(index)
if (err == 0) then
select type (f => filters(index) % obj)
type is (ZernikeFilter)
if (present(x)) f % x = x
if (present(y)) f % y = y
if (present(r)) f % r = r
class default
err = E_INVALID_TYPE
call set_errmsg("Not a Zernike filter.")
end select
end if
end function openmc_zernike_filter_set_params
end module tally_filter_zernike

14
src/tallies/tally_header.F90
View file

@ -354,6 +354,20 @@ contains
j = FILTER_AZIMUTHAL
type is (EnergyFunctionFilter)
j = FILTER_ENERGYFUNCTION
type is (LegendreFilter)
j = FILTER_LEGENDRE
this % estimator = ESTIMATOR_ANALOG
type is (SphericalHarmonicsFilter)
j = FILTER_SPH_HARMONICS
if (filt % cosine == COSINE_SCATTER) then
this % estimator = ESTIMATOR_ANALOG
end if
type is (SpatialLegendreFilter)
j = FILTER_SPTL_LEGENDRE
this % estimator = ESTIMATOR_COLLISION
type is (ZernikeFilter)
j = FILTER_ZERNIKE
this % estimator = ESTIMATOR_COLLISION
end select
this % find_filter(j) = i
end do

6
tests/regression_tests/tally_slice_merge/test.py
View file

@ -87,12 +87,14 @@ class TallySliceMergeTestHarness(PyAPITestHarness):
# Slice the tallies by cell filter bins
cell_filter_prod = itertools.product(tallies, self.cell_filters)
tallies = map(lambda tf: tf[0].get_slice(filters=[type(tf[1])],
filter_bins=[tf[1].get_bin(0)]), cell_filter_prod)
filter_bins=[(tf[1].bins[0],)]),
cell_filter_prod)
# Slice the tallies by energy filter bins
energy_filter_prod = itertools.product(tallies, self.energy_filters)
tallies = map(lambda tf: tf[0].get_slice(filters=[type(tf[1])],
filter_bins=[(tf[1].get_bin(0),)]), energy_filter_prod)
filter_bins=[(tf[1].bins[0],)]),
energy_filter_prod)
# Slice the tallies by nuclide
nuclide_prod = itertools.product(tallies, self.nuclides)

148
tests/unit_tests/test_filters.py Normal file
View file

@ -0,0 +1,148 @@
from math import sqrt, pi
import openmc
from pytest import fixture, approx
@fixture(scope='module')
def box_model():
model = openmc.model.Model()
m = openmc.Material()
m.add_nuclide('U235', 1.0)
m.set_density('g/cm3', 1.0)
box = openmc.model.get_rectangular_prism(10., 10., boundary_type='vacuum')
c = openmc.Cell(fill=m, region=box)
model.geometry.root_universe = openmc.Universe(cells=[c])
model.settings.particles = 100
model.settings.batches = 10
model.settings.inactive = 0
model.settings.source = openmc.Source(space=openmc.stats.Point())
return model
def test_legendre():
n = 5
f = openmc.LegendreFilter(n)
assert f.order == n
assert f.bins[0] == 'P0'
assert f.bins[-1] == 'P5'
assert len(f.bins) == n + 1
# Make sure __repr__ works
repr(f)
# to_xml_element()
elem = f.to_xml_element()
assert elem.tag == 'filter'
assert elem.attrib['type'] == 'legendre'
assert elem.find('order').text == str(n)
def test_spatial_legendre():
n = 5
axis = 'x'
f = openmc.SpatialLegendreFilter(n, axis, -10., 10.)
assert f.order == n
assert f.axis == axis
assert f.minimum == -10.
assert f.maximum == 10.
assert f.bins[0] == 'P0'
assert f.bins[-1] == 'P5'
assert len(f.bins) == n + 1
# Make sure __repr__ works
repr(f)
# to_xml_element()
elem = f.to_xml_element()
assert elem.tag == 'filter'
assert elem.attrib['type'] == 'spatiallegendre'
assert elem.find('order').text == str(n)
assert elem.find('axis').text == str(axis)
def test_spherical_harmonics():
n = 3
f = openmc.SphericalHarmonicsFilter(n)
f.cosine = 'particle'
assert f.order == n
assert f.bins[0] == 'Y0,0'
assert f.bins[-1] == 'Y{0},{0}'.format(n, n)
assert len(f.bins) == (n + 1)**2
# Make sure __repr__ works
repr(f)
# to_xml_element()
elem = f.to_xml_element()
assert elem.tag == 'filter'
assert elem.attrib['type'] == 'sphericalharmonics'
assert elem.attrib['cosine'] == f.cosine
assert elem.find('order').text == str(n)
def test_zernike():
n = 4
f = openmc.ZernikeFilter(n, 0., 0., 1.)
assert f.order == n
assert f.bins[0] == 'Z0,0'
assert f.bins[-1] == 'Z{0},{0}'.format(n)
assert len(f.bins) == (n + 1)*(n + 2)//2
# Make sure __repr__ works
repr(f)
# to_xml_element()
elem = f.to_xml_element()
assert elem.tag == 'filter'
assert elem.attrib['type'] == 'zernike'
assert elem.find('order').text == str(n)
def test_first_moment(run_in_tmpdir, box_model):
plain_tally = openmc.Tally()
plain_tally.scores = ['flux', 'scatter']
# Create tallies with expansion filters
leg_tally = openmc.Tally()
leg_tally.filters = [openmc.LegendreFilter(3)]
leg_tally.scores = ['scatter']
leg_sptl_tally = openmc.Tally()
leg_sptl_tally.filters = [openmc.SpatialLegendreFilter(3, 'x', -5., 5.)]
leg_sptl_tally.scores = ['scatter']
sph_scat_filter = openmc.SphericalHarmonicsFilter(5)
sph_scat_filter.cosine = 'scatter'
sph_scat_tally = openmc.Tally()
sph_scat_tally.filters = [sph_scat_filter]
sph_scat_tally.scores = ['scatter']
sph_flux_filter = openmc.SphericalHarmonicsFilter(5)
sph_flux_filter.cosine = 'particle'
sph_flux_tally = openmc.Tally()
sph_flux_tally.filters = [sph_flux_filter]
sph_flux_tally.scores = ['flux']
zernike_tally = openmc.Tally()
zernike_tally.filters = [openmc.ZernikeFilter(3, r=10.)]
zernike_tally.scores = ['scatter']
# Add tallies to model and ensure they all use the same estimator
box_model.tallies = [plain_tally, leg_tally, leg_sptl_tally,
sph_scat_tally, sph_flux_tally, zernike_tally]
for t in box_model.tallies:
t.estimator = 'analog'
box_model.run()
# Check that first moment matches the score from the plain tally
with openmc.StatePoint('statepoint.10.h5') as sp:
# Get scores from tally without expansion filters
flux, scatter = sp.tallies[plain_tally.id].mean.ravel()
# Check that first moment matches
first_score = lambda t: sp.tallies[t.id].mean.ravel()[0]
assert first_score(leg_tally) == scatter
assert first_score(leg_sptl_tally) == scatter
assert first_score(sph_scat_tally) == scatter
assert first_score(sph_flux_tally) == approx(flux)
assert first_score(zernike_tally) == approx(scatter)