diff --git a/CMakeLists.txt b/CMakeLists.txt index 04a42817bb..ca42bb5ba1 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/docs/source/capi/index.rst b/docs/source/capi/index.rst index d63ed1f3ae..0e6a2536c6 100644 --- a/docs/source/capi/index.rst +++ b/docs/source/capi/index.rst @@ -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 diff --git a/docs/source/examples/expansion-filters.rst b/docs/source/examples/expansion-filters.rst new file mode 100644 index 0000000000..f06d2bde6b --- /dev/null +++ b/docs/source/examples/expansion-filters.rst @@ -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. diff --git a/docs/source/examples/index.rst b/docs/source/examples/index.rst index 5c1efb57ce..89a1f1fe0b 100644 --- a/docs/source/examples/index.rst +++ b/docs/source/examples/index.rst @@ -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 `_ 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 diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 070b5bd201..ef692755d5 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -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 diff --git a/examples/jupyter/expansion-filters.ipynb b/examples/jupyter/expansion-filters.ipynb new file mode 100644 index 0000000000..7870faffe5 --- /dev/null +++ b/examples/jupyter/expansion-filters.ipynb @@ -0,0 +1,463 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "OpenMC's general tally system accommodates a wide range of tally *filters*. While most filters are meant to identify regions of phase space that contribute to a tally, there are a special set of functional expansion filters that will multiply the tally by a set of orthogonal functions, e.g. Legendre polynomials, so that continuous functions of space or angle can be reconstructed from the tallied moments.\n", + "\n", + "In this example, we will determine the spatial dependence of the flux along the $z$ axis by making a Legendre polynomial expansion. Let us represent the flux along the z axis, $\\phi(z)$, by the function\n", + "\n", + "$$ \\phi(z') = \\sum\\limits_{n=0}^N a_n P_n(z') $$\n", + "\n", + "where $z'$ is the position normalized to the range [-1, 1]. Since $P_n(z')$ are known functions, our only task is to determine the expansion coefficients, $a_n$. By the orthogonality properties of the Legendre polynomials, one can deduce that the coefficients, $a_n$, are given by\n", + "\n", + "$$ a_n = \\frac{2n + 1}{2} \\int_{-1}^1 dz' P_n(z') \\phi(z').$$\n", + "\n", + "Thus, the problem reduces to finding the integral of the flux times each Legendre polynomial -- a problem which can be solved by using a Monte Carlo tally. By using a Legendre polynomial filter, we obtain stochastic estimates of these integrals for each polynomial order." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import openmc\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To begin, let us first create a simple model. The model will be a slab of fuel material with reflective boundaries conditions in the x- and y-directions and vacuum boundaries in the z-direction. However, to make the distribution slightly more interesting, we'll put some B4C in the middle of the slab." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# Define fuel and B4C materials\n", + "fuel = openmc.Material()\n", + "fuel.add_element('U', 1.0, enrichment=4.5)\n", + "fuel.add_nuclide('O16', 2.0)\n", + "fuel.set_density('g/cm3', 10.0)\n", + "\n", + "b4c = openmc.Material()\n", + "b4c.add_element('B', 4.0)\n", + "b4c.add_element('C', 1.0)\n", + "b4c.set_density('g/cm3', 2.5)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# Define surfaces used to construct regions\n", + "zmin, zmax = -10., 10.\n", + "box = openmc.model.get_rectangular_prism(10., 10., boundary_type='reflective')\n", + "bottom = openmc.ZPlane(z0=zmin, boundary_type='vacuum')\n", + "boron_lower = openmc.ZPlane(z0=-0.5)\n", + "boron_upper = openmc.ZPlane(z0=0.5)\n", + "top = openmc.ZPlane(z0=zmax, boundary_type='vacuum')\n", + "\n", + "# Create three cells and add them to geometry\n", + "fuel1 = openmc.Cell(fill=fuel, region=box & +bottom & -boron_lower)\n", + "absorber = openmc.Cell(fill=b4c, region=box & +boron_lower & -boron_upper)\n", + "fuel2 = openmc.Cell(fill=fuel, region=box & +boron_upper & -top)\n", + "geom = openmc.Geometry([fuel1, absorber, fuel2])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the starting source, we'll use a uniform distribution over the entire box geometry." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "settings = openmc.Settings()\n", + "spatial_dist = openmc.stats.Box(*geom.bounding_box)\n", + "settings.source = openmc.Source(space=spatial_dist)\n", + "settings.batches = 210\n", + "settings.inactive = 10\n", + "settings.particles = 1000" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Defining the tally is relatively straightforward. One simply needs to list 'flux' as a score and then add an expansion filter. For this case, we will want to use the `SpatialLegendreFilter` class which multiplies tally scores by Legendre polynomials evaluated on normalized spatial positions along an axis." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a flux tally\n", + "flux_tally = openmc.Tally()\n", + "flux_tally.scores = ['flux']\n", + "\n", + "# Create a Legendre polynomial expansion filter and add to tally\n", + "order = 8\n", + "expand_filter = openmc.SpatialLegendreFilter(order, 'z', zmin, zmax)\n", + "flux_tally.filters.append(expand_filter)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The last thing we need to do is create a `Tallies` collection and export the entire model, which we'll do using the `Model` convenience class." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "tallies = openmc.Tallies([flux_tally])\n", + "model = openmc.model.Model(geometry=geom, settings=settings, tallies=tallies)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Running a simulation is now as simple as calling the `run()` method of `Model`." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "0.3016877031715249+/-0.0006126949350699303" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model.run(output=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now that the run is finished, we need to load the results from the statepoint file." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "with openmc.StatePoint('statepoint.210.h5') as sp:\n", + " df = sp.tallies[flux_tally.id].get_pandas_dataframe()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We've used the `get_pandas_dataframe()` method that returns tally data as a Pandas dataframe. Let's see what the raw data looks like." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
spatiallegendrenuclidescoremeanstd. dev.
0P0totalflux36.4348280.076755
1P1totalflux0.0217970.043545
2P2totalflux-4.3808920.025739
3P3totalflux0.0014480.020740
4P4totalflux-0.2954390.014215
5P5totalflux0.0035140.012017
6P6totalflux0.1052130.010103
7P7totalflux0.0025950.009265
8P8totalflux-0.0961970.007513
\n", + "
" + ], + "text/plain": [ + " spatiallegendre nuclide score mean std. dev.\n", + "0 P0 total flux 3.64e+01 7.68e-02\n", + "1 P1 total flux 2.18e-02 4.35e-02\n", + "2 P2 total flux -4.38e+00 2.57e-02\n", + "3 P3 total flux 1.45e-03 2.07e-02\n", + "4 P4 total flux -2.95e-01 1.42e-02\n", + "5 P5 total flux 3.51e-03 1.20e-02\n", + "6 P6 total flux 1.05e-01 1.01e-02\n", + "7 P7 total flux 2.60e-03 9.27e-03\n", + "8 P8 total flux -9.62e-02 7.51e-03" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "df" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the expansion coefficients are given as\n", + "\n", + "$$ a_n = \\frac{2n + 1}{2} \\int_{-1}^1 dz' P_n(z') \\phi(z')$$\n", + "\n", + "we just need to multiply the Legendre moments by $(2n + 1)/2$." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "n = np.arange(order + 1)\n", + "a_n = (2*n + 1)/2 * df['mean']" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To plot the flux distribution, we can use the `numpy.polynomial.Legendre` class which represents a truncated Legendre polynomial series. Since we really want to plot $\\phi(z)$ and not $\\phi(z')$ we first need to perform a change of variables. Since\n", + "\n", + "$$ \\lvert \\phi(z) dz \\rvert = \\lvert \\phi(z') dz' \\rvert $$\n", + "\n", + "and, for this case, $z = 10z'$, it follows that\n", + "\n", + "$$ \\phi(z) = \\frac{\\phi(z')}{10} = \\sum_{n=0}^N \\frac{a_n}{10} P_n(z'). $$" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "phi = np.polynomial.Legendre(a_n/10, domain=(zmin, zmax))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's plot it and see how our flux looks!" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEKCAYAAAAB0GKPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8FfXZ9/HPlT2BEBIStiSQsO9r2BQUXFjUilZb0d4u\nbZVH617t/tRarb3vLrdPa+tGLXWp+1ZRcUGRxYUloOwCIWFJDCQQIASy53r+OAM90iwnkMmcnFzv\n1+u8cs7MnJlvJie5MvOb+f1EVTHGGGOaEuZ1AGOMMW2DFQxjjDEBsYJhjDEmIFYwjDHGBMQKhjHG\nmIBYwTDGGBMQKxjGGGMCYgXDGGNMQKxgGGOMCUiE1wFaUnJysmZkZHgdwxhj2ow1a9bsV9WUQJYN\nqYKRkZFBdna21zGMMabNEJFdgS5rp6SMMcYExAqGMcaYgFjBMMYYExArGMYYYwJiBcMYY0xArGAY\nY4wJiBUMY4wxAQmp+zCMCYSqUllTR2l5NaUVNZRWVFNaXk1FdS1VtUpNbR3VtXVU1SrVNXXU1n19\nGGOR+tcrfjOOD32sCoo6X7/+2n9Z//nHtyHHv4o4r4Uw4cTziHAhLiqc2KgI4iLDiYsKJy46gqS4\nKJLjo4iLsl9v07Jc+0SJSDrwNNAN3+/CPFX980nLfAf4Cb7fjSPATaq6zpm305lWC9SoapZbWU3o\nqaiuJW//UXKKytheVEbBwXL2lpaz93AFew9XcLSq1uuIrouLCie5YzTdO8XQu0scGckdyOjSgX5d\nO9I3pQMR4XaCwTSPm/+C1AB3qepaEYkH1ojIIlXd7LdMHnC2qh4UkVnAPGCC3/xpqrrfxYwmBKgq\nO4rLWLPrIGt3HWLN7oPkFpdx/MAgTKB7pxi6JcQwsHs8Zw1IIbljNJ1iI+kUE3Hia0xkONERYUSG\n+x4R4UJUeBjhYXLi6EH9Dg30axn8XwDif5Qg/z5aQE4coZz4ivzHsgrUnTjycL46z+vUl6OmVjlW\nXUt5VQ1HK2s5VlXLsaoaDhytYn9ZJfuP+L7uPVzB0m3FvLwm/0TEmMgwhvZMYHhqAuMykjijbxcS\nO0S16M/FhB7XCoaqFgKFzvMjIrIFSAU2+y3zqd9bVgBpbuUxoeVoZQ0f5+znoy+LWPxlEUVHKgHo\nHBfJ2F6JXDCsO/27xdOva0cykzsQExnuceLmC6eBc19+EpuxvqOVNew6cIyt+0rZkF/KhoJDvLh6\nD09+uhMRGNYzgbMHpHDB8B4M7hH/tVNsxgCIqja91OluRCQDWAYMU9XSBpa5Gxikqtc7r/OAg/j+\n2XpcVec18L65wFyAXr16jd21K+BuUUwbU1unfJyzn9fW5vPepr1UVNcRHx3BlAHJnD0ghayMJPok\nd7A/dM1QU1vHuvzDfJKzn4+372fN7oPU1il9kjtw0YgefCsrnfSkOK9jGheJyJpAT/m7XjBEpCOw\nFHhAVV9rYJlpwCPAZFU94ExLVdUCEekKLAJuVdVljW0rKytLrfPB0HOgrJJnV+7m2ZW72FdaSUJs\nJBeN6MGFI3owLiOJSDsX32IOlFXy7qa9vL2+kBW5B1Bg6oAUrp7Um6kDuhIWZsU41ARNwRCRSOAt\n4D1VfbCBZUYArwOzVHVbA8vcC5Sp6h8b254VjNCyp+QYjyzZwWtr86msqePsASlcOT6daYO6Eh3R\n9k4xtTVfHSrnhVW7eX71HoqPVDKgW0duOac/Fw7vQbgVjpARFAVDfOcFngJKVPWOBpbpBSwGrvFv\nzxCRDkCY0/bRAd8Rxn2q+m5j27SCERqKSiv4y+IcXli9GxHhsjFpfH9yBv26xnsdrV2qrq3j7fWF\n/PWjHHKKyuiT0oEfzxjIjKHd7fRfCAiWgjEZWA5sAOqcyT8HegGo6mMi8gRwGXC84aFGVbNEpA++\now7wNcw/p6oPNLVNKxhtW2VNLY8vzeWRJTnU1CpXjEvn1nP60z0hxutoBqirU97dtJc/fbCNbfvK\nmNgniXsuGsqQnp28jmZOQ1AUDC9YwWi7lm0r5lcLNpG3/ygXDO/OT2YOoneXDl7HMvWoqa3j+VW7\neXDRNg6XV3PdGZn8aMZAYqPsNGFb1JyCYbeCGk8dqajm129u5pU1+WQmd+Dp743nrAEBjRZpPBIR\nHsbVkzK4eGQqf3j/S+Z/ksfiL/fxu8tGMKFPF6/jGRfZ5SXGM6vySpj15+W8tjafW6b14907plix\naEMS4iL5zSXDee6GCdQpXDFvBf+9cAvVtXVNv9m0SVYwTKurq1P+/MF2rpj3GWEivHzjJO6eMdCu\nfGqjzuibzLt3TOE7E3rx+LJc5sxbQeHhcq9jGRdYwTCt6nB5NTc8nc3/+2Abl4xKZeHtUxjbO8nr\nWOY0xUVF8MClw3noytF8WVjKBX9ezsfbrVefUGMFw7SanKIjXPLwJyzdVsy93xjCg98eScdoa0YL\nJReP7MmCWyfTNT6Ga/+xiudW7vY6kmlBVjBMq1iVV8I3H/mUIxXVPHfDRK47M9Ou4Q9RfVM68spN\nk5jSP5mfv76BB97e/B9dxJu2yQqGcd3CDYX8199Xktwxmtd/cCbjM+0UVKiLj4nkiWuyuHZSb/62\nPI9bn19LVY01hrd1dj7AuOqZFbu4542NjE7vzBPXjiPJutBuNyLCw/j17GGkJcbxwMItlFdl8+h/\njW2TPQcbHzvCMK558pM8fvmvjZwzsCvP3TDRikU7dcNZffjtpcNZsq2Y6/6xirLKGq8jmVNkBcO4\n4u8f53Hvm5uZPqSb/VdpuGpCL/50xShW7zzIdfNXcazKikZbZAXDtLh/fJLH/W9tZtaw7jz8nTFE\nRdjHzMDsUak8NGc0a3cfZO7Ta6ioDv1hckON/SabFvXa2nx+/eZmZg7tzkNXjraxKszXXDiiB7+/\nfCQf5+znlufW2l3hbYz9NpsW88HmffzolfWc2a8Lf75ylBULU6/Lx6Zx/yXD+GBLEXe/vI46u+S2\nzbCrpEyLWJVXws3PrWVoz048fnWWdfNhGnX1xN6Ullfzh/e2kto5lh/PHOR1JBMAKxjmtO0oLuP6\np1aTmhjLP64bZ3dvm4D8YGpf8g+W88iSHaQlxnHVhF5eRzJNsN9sc1oOHavi+qeyiQgP46nvjqdL\nx2ivI5k2QkS4f/ZQCg+X88s3NtKjcwzTBnb1OpZphJ1kNqesuraOHzy7loKD5Tx+9VjSk+K8jmTa\nmIjwMP561RgGdovntuc+Z0dxmdeRTCOsYJhToqrcu2ATn+44wG+/OZxxGdbdhzk1HaMj+Nu1WURG\nhDH36WyOVFR7Hck0wAqGOSXPr9rDsyt3c+PZfbl8bJrXcUwbl9o5loevGsPOA8f44Ut25VSwsoJh\nmm1D/mHuXbCJswak8KMZA72OY0LEpL5d+MUFg1m0eR9/WZzjdRxTD9cKhoiki8hHIrJZRDaJyO31\nLCMi8pCI5IjIehEZ4zfvWhHZ7jyudSunaZ5Dx6q46dk1JHeM4k9XjCI8zLooNy3nu2dm8M0xqfzp\nw20s317sdRxzEjePMGqAu1R1CDARuFlEhpy0zCygv/OYCzwKICJJwK+ACcB44FcikuhiVhOAujrl\nrpfWsa+0goe/M8Y6EzQtTkR44JLh9O/akTtf/IKiIxVeRzJ+XCsYqlqoqmud50eALUDqSYvNBp5W\nnxVAZxHpAcwAFqlqiaoeBBYBM93KagLz+LJcPvyyiP974RBG97L6bdwRGxXOX68aQ1llDXe++IUN\nvhREWqUNQ0QygNHAypNmpQJ7/F7nO9Maml7fuueKSLaIZBcX2yGsW9btOcT/vr+VC4f34JpJvb2O\nY0LcgG7x3PuNoXySc4BHl1h7RrBwvWCISEfgVeAOVS1t6fWr6jxVzVLVrJSUlJZevQGOVtZwx4tf\n0DU+mt9eOtyGVjWt4opx6Vw8sicPLtpG9s4Sr+MYXC4YIhKJr1g8q6qv1bNIAZDu9zrNmdbQdOOB\n37y9mZ0HjvLgFaNIiIv0Oo5pJ0SEBy4dRmpiLD98aR1HbeAlz7l5lZQAfwe2qOqDDSy2ALjGuVpq\nInBYVQuB94DpIpLoNHZPd6aZVvbuxr08v2oPN53dl4l9ungdx7Qz8TGR/O+3RrHn4DEeWLjF6zjt\nnpt9SZ0JXA1sEJEvnGk/B3oBqOpjwELgAiAHOAZ815lXIiL3A6ud992nqnZM2sr2lVbw09fWMyIt\ngTvOG+B1HNNOjc9M4oYpfZi3LJfzh3Sz/qY8JKqhcwVCVlaWZmdnex0jJKgqNzydzfLt+1l4+xT6\npnT0OpJpxyqqa7n4rx9z6Fg17995Fp3j7JLuliIia1Q1K5Bl7U5vU683vviKD7YU8aMZA61YGM/F\nRIbz4LdHUXK0il++scnrOO2WFQzzH4qOVHDvm5sY06sz3z0z0+s4xgAwLDWB287tz5vrvuL9TXu9\njtMuWcEwX6Oq3POvTRyrquX3l4+0rj9MULlpal8GdY/nl29spNR6tW11VjDM17y9oZB3N+3lzvMG\n0K+rnYoywSUyPIzfXTaC4iOV/P7dL72O0+5YwTAnlByt4p43NjEyLYEbptipKBOcRqb7TpX+c8Vu\nVtsNfa3KCoY54bcLt1BaXs3vLx9JRLh9NEzwumv6ANISY/npq+upqK71Ok67YX8VDAArcg/wypp8\n5p7Vh4Hd472OY0yj4qIieODS4ewoPsojH1lfU63FCoahqqaO//uvjaQlxnLrOf29jmNMQM4ekMI3\nR6fy6NId5BQd8TpOu2AFw/C35bnkFJVx/+xhxEaFex3HmID9/MLBxEaG86sFmwilm5CDlRWMdm73\ngWM89OF2Zg3rzrRB1uWCaVuSO0Zz94yBfJJzgIUb7N4Mt1nBaMdUlXsWbCQiTLjnGycPhmhM2/Cd\nCb0Z2rMT97+12Xq0dZkVjHbsvU17WbK1mDvPH0CPhFiv4xhzSsLDhPtmD2NvaQUPLd7udZyQZgWj\nnaqoruX+t7YwqHs8152R4XUcY07L2N6JfDsrjb8vz7MGcBdZwWinHl+aS8Ghcu69eKjdc2FCwk9m\nDiIuKpx73rAGcLfYX4p2qOBQOY8uzeHC4T1sUCQTMrp0jOZHMwby6Q5rAHeLFYx26L8XbkEVfnbB\nIK+jGNOirprQm0Hd4/nvd7bYHeAusILRzqzMPcBb6wv5P2f3JS0xzus4xrSo8DDhlxcNIf9gOfM/\nyfM6TsixgtGO1NYpv35zMz0TYrjp7L5exzHGFWf2S+a8wV155KMdFB2p8DpOSLGC0Y68uHoPmwtL\n+dkFg+2ObhPSfn7BYCqqa3nw/W1eRwkpVjDaidKKav74/lbGZyZx0YgeXscxxlV9UjpyzaQMXsze\nw6avDnsdJ2S4VjBEZL6IFInIxgbm/0hEvnAeG0WkVkSSnHk7RWSDMy/brYztySMf7eDgsSruuWgI\nIjaKngl9t5/bn4TYSH7z1ha7zLaFuHmE8SQws6GZqvoHVR2lqqOAnwFLVdV/NJRpzvwsFzO2CwWH\nfA2Al45KZVhqgtdxjGkVCXGR3HneAD7LPcCizfu8jhMSXCsYqroMCHQ4rCuB593K0t798b2tANw1\nY6DHSYxpXVdN6EW/rh357cItVNfWeR2nzfO8DUNE4vAdibzqN1mB90VkjYjMbeL9c0UkW0Syi4uL\n3YzaJm0sOMzrnxfwvTMzSe1s/UWZ9iUyPIyfzRrEzgPHeGHVbq/jtHmeFwzgG8AnJ52OmqyqY4BZ\nwM0iclZDb1bVeaqapapZKSkpbmdtU1SV3y7cQmJcJD+YZpfRmvbpnEFdGZ+ZxJ8/3G692Z6mYCgY\nczjpdJSqFjhfi4DXgfEe5Grzlmwt5tMdB7j93P50ion0Oo4xnhARfjprEPvLqvjb8lyv47RpnhYM\nEUkAzgbe8JvWQUTijz8HpgP1XmllGlZTW8dvF24ho0scV03o7XUcYzw1plcis4Z1Z96yXIqPVHod\np81y87La54HPgIEiki8i3xeRG0XkRr/FLgXeV9WjftO6AR+LyDpgFfC2qr7rVs5Q9fKafLYXlfGT\nmYOIigiGA0ljvHX3jIFU1tTxFxsz45RFuLViVb0ygGWexHf5rf+0XGCkO6nah6OVNTy4aBtjeycy\nc1h3r+MYExT6pnRkzrh0nlu5m++emUlmcgevI7U59q9nCHpieR7FRyr5+QWD7SY9Y/zcfl5/IsPD\nTlxqbprHCkaIKTnqa9ibObQ7Y3sneh3HmKDSNT6GG6Zk8vaGQr7Yc8jrOG2OFYwQ88hHORyrquHu\nGQO8jmJMUJp7dl+6dIjif96xLkOaywpGCPnqUDlPr9jFN8ek0a9rvNdxjAlKHaMjuO3c/qzILWHJ\nNrvZtzmsYISQvyzeDgp3nNff6yjGBLUrx/ciPSmWP763lbo6O8oIlBWMEJFbXMZL2flcNaGXjaRn\nTBOiIsK4/dwBbPqqlHc32fjfgWr0sloRWRDAOkpU9bqWiWNO1YOLthEdEcbN0/p5HcWYNuHS0ak8\nuiSHBxdtY8bQ7oSH2RWFTWnqPozBwPWNzBfg4ZaLY07FxoLDvLW+kFum9SMlPtrrOMa0CeFhwl3T\nB/KDZ9fyr88LuGxsmteRgl5TBeMXqrq0sQVE5NctmMecgj++v5WE2EhuOKuP11GMaVNmDu3O0J6d\n+NOH2/jGyJ7WK0ITGt07qvpSUysIZBnjnlV5JSzZWsxNU/uSEGsdDBrTHGFhwt3TB7KnpJyXsvd4\nHSfoBVRORWSRiHT2e50oIu+5F8sEQlX5/btf0jU+mmsnZXgdx5g2aerAFMb2TuQvi7dTUV3rdZyg\nFujxV7KqnrgtUlUPAl3diWQC9dHWIrJ3HeS2c/sTGxXudRxj2iQR4UczBrKvtJJnPtvldZygFmjB\nqBORXsdfiEhvfKPiGY/U1Sl/eG8bvbvEccW4dK/jGNOmTezThSn9k3l06Q7KbJClBgVaMH6Br8vx\nZ0Tkn8Ay4GfuxTJNeXtDIVsKS/nh+QOIDLeGOmNO113TB1JytIr5H+d5HSVoNfmXRnzdnW4CxgAv\nAi8AY1XV2jA8Ulun/OmDbQzo1pFvjOjpdRxjQsKo9M6cP6Qbf1uWy6FjVV7HCUpNFgz19c61UFX3\nq+pbzmN/K2QzDXhz3VfsKD7KnecNIMxuNjKmxdw1fQBlVTU8vsyGcq1PoOcy1orIOFeTmIDU1Nbx\n5w+3M7hHJ2YMtcGRjGlJg7p34hsjevLUpzs5UGZDuZ4s0IIxAfhMRHaIyHoR2SAi690MZur3+ucF\n5O0/yp3n9bejC2NccNu5/amormWeHWX8h0CHaJ3hagoTkOraOh5avJ1hqZ04f0g3r+MYE5L6de3I\nxSN78vRnu7h+Sh/rbsdPoEcYEcBeVd0FZAKzgcOupTL1enVNPntKyvnh+QNs6FVjXHTbuf2prKnl\n8aU7vI4SVAItGK8CtSLSD5gHpAPPNfYGEZkvIkUisrGB+VNF5LCIfOE87vGbN1NEtopIjoj8NMCM\nIa2qpo6/LM5hVHpnpg20eyaNcVOflI5cMiqVf67cRdGRCq/jBI2Ab9xT1Rrgm8BfVPVHQI8m3vMk\nMLOJZZar6ijncR+AiITj6wF3FjAEuFJEhgSYM2S9mL2HgkN2dGFMa7n13P5U1yqPLbG2jOMCLRjV\nInIlcA3wljOt0Z7uVHUZUHIKmcYDOaqaq6pV+O77mH0K6wkZFdW1PLw4h6zeiUzpn+x1HGPahczk\nDlw6OpVnV+6iqNSOMiDwgvFdYBLwgKrmiUgm8EwLbH+SiKwTkXdEZKgzLRXw7zYy35lWLxGZKyLZ\nIpJdXBya4/O+sGo3e0sr7OjCmFZ26zn9qKlTHllibRnQRMEQkXkicimwR1VvU9XnAVQ1T1V/d5rb\nXgv0VtWRwF+Af53KSlR1nqpmqWpWSkrKaUYKPhXVtTy8ZAcT+yRxRj87ujCmNfXu0oHLxqTy3Krd\n7D1sRxlNHWH8HRgJLBSRD0XkJyIysiU2rKqlqlrmPF8IRIpIMlCAr1H9uDRnWrv0zxW7KD5SyZ3n\nDfA6ijHt0q3n9KeuTnl0SY7XUTzX1ABKK1X1XlWdAnwb2A3c5VzVNF9Evn2qGxaR7k4/VYjIeCfL\nAWA10F9EMkUkCpgDBDK2eMg5VlXDo0t2MLlfMhP6dPE6jjHtUnpSHJePTeP5VXsoPFzudRxPBdzN\nqaoeUNXnVfUaVR2F70qm/g0tLyLPA58BA0UkX0S+LyI3isiNziKXAxtFZB3wEDBHfWqAW4D3gC3A\nS6q66dS+vbbt6c92ceBoFXee3+BuNsa0gpun9aNOlUc+at9tGQHd6S0i0cBlQIb/e45fClsfVb2y\nsXWq6l+BvzYwbyGwMJBsoaqssobHl+7g7AEpjO2d5HUcY9q19KQ4vpWVzour93DT1L707BzrdSRP\nBHqE8Qa+S1trgKN+D+OSpz7dycFj1fzwfGu7MCYY3HJOPxTl4Y/ab1tGoH1JpalqUzfhmRZSVlnD\n35bnMm1gCiPTOzf9BmOM61I7x/LtrHReyvYdZaQlxnkdqdUFeoTxqYgMdzWJOeGfK3Zx6Fg1t9uV\nUcYElZun9UMQHm6nbRmBFozJwBqnfyfr3txFx6pq+NuyXM4ekMIoO7owJqj07BzLt8el8cqaPXx1\nqP1dMRVowZiF74qo6cA3gIucr6aFPbtiNweOVnHbuXZllDHB6Maz+6IKj7XDnmwDKhiququ+h9vh\n2pvyqloeX5bL5H7JjO2d6HUcY0w90hLjuGxMGi+s3tPu+phqqmuQtU2tIJBlTGCeX7Wb/WWVdnRh\nTJD7wbS+1NZpuxv7u6mrpAY30VYhQEIL5mm3KqpreWzpDib16cL4TLvvwphg1rtLB2aP6smzK3dx\n09S+JHdsH6PyNVUwBgWwjtqWCNLevbh6D0VHKvnznNFeRzHGBODmaf14/fMC/rY8l5/NGux1nFbR\naMGwdorWUVlTy6NLdjA+I4mJfezowpi2oG9KRy4a0ZNnPtvFjWf1JbFDlNeRXBdwX1LGPS9n57O3\ntILbz+tv410Y04bcMq0fx6pqmf9JntdRWoUVDI9V1dTx6JIdjO2dyBl9rUdaY9qSgd3jmTWsO09+\nspPD5dVex3FdQAWjvjG1RWRqi6dph15dm0/BoXJuO9eOLoxpi245px9HKmt48pOdXkdxXaBHGC85\ngyeJiMSKyF+A/3YzWHtQXVvHwx/lMDK9M2fZWN3GtElDeyZw3uCuzP8kjyMVoX2UEWjBmIBvFLxP\n8Q1w9BVwpluh2ovXPy8g/2A5d9jRhTFt2q3n9OdweTXPrAjt64QCLRjVQDkQC8QAeapa51qqdqC2\nTnlsyQ6G9uzE1IGhNxa5Me3JyPTOnDUghSeW53GsqsbrOK4JtGCsxlcwxgFTgCtF5GXXUrUD727c\nS+7+o77eL+3owpg277Zz+lFytIrnVu72OoprAi0Y31fVe1S1WlULVXU27XSc7Zag6huEpU9yB2YM\n7e51HGNMC8jKSOKMvl14bGkuFdWheT9zoAWjSER6+T+ApW4GC2VLtxWzubCUG6f2JTzMji6MCRW3\nTOvH/rJKXl2b73UUVwQ64t7bgOLrOyoGyAS2AkNdyhXSHvloBz0SYrhkVKrXUYwxLWhS3y6MTEvg\n8aW5XJGVTkR4aN3qFmj35sNVdYTztT8wHvissfeIyHwRKRKRjQ3M/47fYEyfishIv3k7nelfiEh2\nc76hYLd6ZwmrdpYw96w+REWE1ofJmPZORLhpal92lxzjnY17vY7T4k7pL5aqrsV3qW1jngQaGwc8\nDzhbVYcD9wPzTpo/TVVHqWrWqWQMVo98lENShyjmjOvldRRjjAumD+lOn5QOPLpkB6rqdZwWFdAp\nKRH5od/LMGAMvnsxGqSqy0Qko5H5n/q9XAGkBZKlLdv01WE+2lrM3dMHEBsV7nUcY4wLwsKEG8/q\ny49fXc/y7fs5a0DoXDYf6BFGvN8jGl+bxuwWzPF94B2/1wq8LyJrRGRuY28Ukbkiki0i2cXFxS0Y\nqeU9umQHHaMjuHpShtdRjDEumj26J907xfDoktAaxjWgIwxV/bVbAURkGr6CMdlv8mRVLRCRrsAi\nEflSVZc1kG0ezumsrKysoD3+yy0u4+0Nhfyfs/qSEBvpdRxjjIuiI8K5fkomv3l7C5/vPsjoXqEx\n5HJTQ7S+KSILGnqc7sZFZATwBDBbVQ8cn66qBc7XIuB1fI3sbdq8ZblEhofxvckZXkcxxrSCOeN7\nkRAbyWNLQ+coo6kjjD+6tWHnXo7XgKtVdZvf9A5AmKoecZ5PB+5zK0drKD5SyWtrC7g8K42u8TFe\nxzHGtIKO0RFcO6k3Dy3OIafoCP26xnsd6bQ1VTDyVPWU7nMXkeeBqUCyiOQDvwIiAVT1MeAeoAvw\niNM1Ro1zRVQ34HVnWgTwnKq+eyoZgsUzn+2kuq6O70/O9DqKMaYVXXtGBvOW5/L40lz+8K2RTb8h\nyDVVMP6F74ooRORVVb0s0BWr6pVNzL8euL6e6blA29+zjvKqWp5ZsYtzB3Wjb0pHr+MYY1pRl47R\nzBnXi2dX7uLO8wfQs3Os15FOS1NXSfn3W9HHzSCh6pW1+Rw8Vs3cs2z3GdMeXT8lkzqFf4TAMK5N\nFQxt4LkJQG2dMv/jPEamJTAuIzSukjDGNE9aYhwXDO/BC6v2UFbZtrs+b6pgjBSRUhE5AoxwnpeK\nyBERKW2NgG3ZB1v2kbf/KNdP6WNdmBvTjn1/ciZHKmt4afUer6OclkYLhqqGq2onVY1X1Qjn+fHX\nnVorZFv1xPJcUjvHMmuYdWFuTHs2Kr0zWb0T+cenedTWtd2TNdb7nUs+332Q1TsP8r3JmSHXY6Ux\npvm+PzmTPSXlLNq8z+sop8z+krnkieV5xMdEcMW4dK+jGGOCwPSh3UlPiuXvH+d6HeWUWcFwQcGh\nct7ZWMhV43vRMTrQIUeMMaEsPEy47oxMVu88yLo9h7yOc0qsYLjg2RW7ALh6Um+Pkxhjgsm3s9Lo\nGB3B3z9um5fYWsFoYRXVtbyweg/nDe5GWmKc13GMMUEkPiaSOePSWbihkK8OlXsdp9msYLSwt9YX\nUnK0imtQcDLwAAAUCklEQVTPyPA6ijEmCF17RgZ1qvzTORPRlljBaEGqylOf7qRf146c0beL13GM\nMUEoPSmOcwd348XVe6isqfU6TrNYwWhBn+85xIaCw1w7qbfdqGeMadDVE3tz4GgV77axcb+tYLSg\npz/dScfoCC4dE/KjzRpjTsPkfslkdInjmc/a1mkpKxgtpPhIJW9vKOTysWl2Ka0xplFhYcJ/TexN\n9q6DbP6q7fSyZAWjhbyUvYfqWrVLaY0xAfnW2HRiIsN4pg01flvBaAF1dcqLq/cwITPJxrwwxgQk\nIS6Si0f25F+fF1BaUe11nIBYwWgBn+UeYHfJMa4c38vrKMaYNuTqiRmUV9fy2pp8r6MExApGC3h+\n1W4SYiOZab3SGmOaYXhaAiPTO/Psyt2oBn8vtlYwTlPJ0Sre37SPS0enEhMZ7nUcY0wbM2dcOtuL\nyvi8DfQvZQXjNL22Np+q2jo7HWWMOSUXjehBbGQ4L2cH/+BKrhYMEZkvIkUisrGB+SIiD4lIjois\nF5ExfvOuFZHtzuNaN3OeKlXl+VW7Gd2rMwO7x3sdxxjTBsXHRHLB8B68ua6QY1XBPYSr20cYTwIz\nG5k/C+jvPOYCjwKISBLwK2ACMB74lYgE3aDYa3YdZEfxUa4cZ0cXxphTd8W4dMoqa1i4Ibjv/Ha1\nYKjqMqCkkUVmA0+rzwqgs4j0AGYAi1S1RFUPAotovPB44tW1+cRFhXPhiB5eRzHGtGHjMhLJTO4Q\n9GN+e92GkQr476F8Z1pD0/+DiMwVkWwRyS4uLnYt6Mkqqmt5a30hM4d2p4Pd2W2MOQ0iwrey0li1\ns4Tc4jKv4zTI64Jx2lR1nqpmqWpWSkpKq2138ZdFHKmo4dIx9dYxY4xplsvGpBEm8HIQ35PhdcEo\nAPwHvU5zpjU0PWi8traAbp2iOaNvstdRjDEhoFunGKYO7Mpra/OprQvOezK8LhgLgGucq6UmAodV\ntRB4D5guIolOY/d0Z1pQKDlaxZKtRcwelUp4mHVjboxpGZeMTmVfaSWr8hpr+vWOqyffReR5YCqQ\nLCL5+K58igRQ1ceAhcAFQA5wDPiuM69ERO4HVjuruk9Vg2YPvrX+K2rqlG/a6ShjTAs6f3A34qLC\neeOLAiYF4SBsrhYMVb2yifkK3NzAvPnAfDdyna5X1xYwuEcnBnXv5HUUY0wIiY0KZ8bQ7izcUMiv\nZw8lOiK4eo/w+pRUm7PrwFHW7TnEpaN7eh3FGBOCLh7Vk9KKGpZsbb2rPgNlBaOZ3t5QCMCFI6xg\nGGNa3uR+yXTpEMWCL77yOsp/sILRTG+vL2R0r86kdo71OooxJgRFhodx4YgefLBlH0eCbJwMKxjN\nsOvAUTZ9VcqFw+3ObmOMe2aPSqWypo73Nu3zOsrXWMFohuOno2ZZwTDGuGiMcxbj7fXBdVrKCkYz\nLNxQyKh0Ox1ljHGXiDBrWHc+ztkfVMO3WsEI0K4DR9lYYKejjDGtY9bw7lTXKou3FHkd5QQrGAH6\n9+koG4bVGOO+0emJdOsUzTsbC72OcoIVjAC9s2Evo9I7k5YY53UUY0w7EBYmzBjanaXbioNmYCUr\nGAEoPFzOhoLDTB/azesoxph2ZOaw7lRU17E0SG7is4IRgA+dc4jnD7aCYYxpPeMzkkjqEMU7G4Nj\nJD4rGAH4cMs+eiXF0a9rR6+jGGPakYjwMKYP6cbiL4uorKn1Oo4VjKYcq6rhkx0HOG9wN0SsK3Nj\nTOs6b3A3yiprWJ130OsoVjCasnz7fqpq6jhvcFevoxhj2qEz+nUhKiKMxV96f3mtFYwmfLhlH/Ex\nEYzLTPI6ijGmHYqLiuCMvl1Y/KX33YRYwWhEXZ2y+Msipg7sSmS47SpjjDfOGdSVnQeOkVtc5mkO\n+yvYiHX5h9hfVmWno4wxnpo20Pc3yOvTUlYwGrFs235E4Kz+KV5HMca0Y+lJcQzo1tEKRjBbvr2Y\nEakJJHaI8jqKMaadmzaoK6vySjwdI8MKRgNKK6r5fM8hJvdP9jqKMcZwzsCu1NQpH2/f71kGVwuG\niMwUka0ikiMiP61n/v8TkS+cxzYROeQ3r9Zv3gI3c9ZnxY4D1NYpU+x0lDEmCIztnUh8dATLPCwY\nEW6tWETCgYeB84F8YLWILFDVzceXUdU7/Za/FRjtt4pyVR3lVr6mLN++n7iocMb0SvQqgjHGnBAR\nHsbEvl34OMe7fqXcPMIYD+Soaq6qVgEvALMbWf5K4HkX8zTL8u3FTOrju2HGGGOCweR+yewpKWf3\ngWOebN/Nv4apwB6/1/nOtP8gIr2BTGCx3+QYEckWkRUicklDGxGRuc5y2cXFLVN5dx84xs4Dx5hi\n7RfGmCByvE11uUdHGcHy7/Mc4BVV9e9dq7eqZgFXAX8Skb71vVFV56lqlqpmpaS0THvD8R/GlAHW\nfmGMCR59kjvQIyGGT3K8acdws2AUAOl+r9OcafWZw0mno1S1wPmaCyzh6+0brvp4+356JsTQJ7lD\na23SGGOaJCJM7pfMp85FOa3NzYKxGugvIpkiEoWvKPzH1U4iMghIBD7zm5YoItHO82TgTGDzye91\ng6qyMq+EiX27WO+0xpigM7l/MoeOVbPpq8Otvm3XCoaq1gC3AO8BW4CXVHWTiNwnIhf7LToHeEFV\n/cvlYCBbRNYBHwH/4391lZu2F5VRcrSKiX26tMbmjDGmWc7o67RjeHB5rWuX1QKo6kJg4UnT7jnp\n9b31vO9TYLib2RqyMvcAABMzrWAYY4JPSnw0A7vFszKvhJunte62g6XRO2isyCuhR0IM6UmxXkcx\nxph6jc9MYs3OEmpq61p1u1Yw/KgqK3NLmJCZZO0XxpigNT4ziaNVtWwuLG3V7VrB8JO7/yj7yyqZ\nYO0XxpggNt4Z0G1VXkmrbtcKhp+Vub6dP8FG1zPGBLFunWLI6BLHSisY3lmZd4CU+Ggy7f4LY0yQ\nG5+ZxOqdJdS14v0YVjD8rM4rYby1Xxhj2oDxmV04dKya7UWtN2yrFQzH3sMVfHW4grHWO60xpg2Y\ncKId40CrbdMKhmPt7oMAjOltBcMYE/zSEmPpkRDDilZsx7CC4Vi76yDREWEM6dHJ6yjGGNMkEWFs\n70Q+33Ww1bZpBcOxdvdBhqcm2PgXxpg2Y2zvRL46XEHh4fJW2Z79dQQqa2rZWFBqp6OMMW3K8RFB\n1+461MSSLcMKBrDpq1KqausY06uz11GMMSZgQ3p2IiYy7EQbrNusYOBrvwAYbVdIGWPakMjwMEak\ndmZNK7VjWMEAPt99iNTOsXTrFON1FGOMaZaxGYlU19a1yg18rnZv3las3X2QsdZ+YYxpg348YyA/\nmTmoVbbV7gtGZU0tk/slc2a/ZK+jGGNMs7VmzxTtvmBER4Tzh2+N9DqGMcYEPWvDMMYYExArGMYY\nYwJiBcMYY0xAXC0YIjJTRLaKSI6I/LSe+deJSLGIfOE8rvebd62IbHce17qZ0xhjTNNca/QWkXDg\nYeB8IB9YLSILVHXzSYu+qKq3nPTeJOBXQBagwBrnva3Xy5YxxpivcfMIYzyQo6q5qloFvADMDvC9\nM4BFqlriFIlFwEyXchpjjAmAmwUjFdjj9zrfmXayy0RkvYi8IiLpzXwvIjJXRLJFJLu4uLglchtj\njKmH143ebwIZqjoC31HEU81dgarOU9UsVc1KSUlp8YDGGGN83LxxrwBI93ud5kw7QVX9xxZ8Avi9\n33unnvTeJU1tcM2aNftFZNcpZAVIBvaf4nvdZLmax3I1j+VqnlDM1TvQBUXVnQ6rRCQC2Aaci68A\nrAauUtVNfsv0UNVC5/mlwE9UdaLT6L0GGOMsuhYYq6qujUUoItmqmuXW+k+V5Woey9U8lqt52nsu\n144wVLVGRG4B3gPCgfmquklE7gOyVXUBcJuIXAzUACXAdc57S0TkfnxFBuA+N4uFMcaYprnal5Sq\nLgQWnjTtHr/nPwN+1sB75wPz3cxnjDEmcF43egeTeV4HaIDlah7L1TyWq3nadS7X2jCMMcaEFjvC\nMMYYE5B2VTBE5FsisklE6kQk66R5P3P6vNoqIjMaeH+miKx0lntRRKJcyPiiX99aO0XkiwaW2yki\nG5zlsls6Rz3bu1dECvyyXdDAco32H+ZCrj+IyJfOzZ+vi0jnBpZrlf0VQP9p0c7POMf5LGW4lcVv\nm+ki8pGIbHY+/7fXs8xUETns9/O9p751uZCt0Z+L+Dzk7K/1IjKmvvW0cKaBfvvhCxEpFZE7Tlqm\nVfaXiMwXkSIR2eg3LUlEFjn97C0SkXqHC3WlPz5VbTcPYDAwEN89HVl+04cA64BoIBPYAYTX8/6X\ngDnO88eAm1zO+7/APQ3M2wkkt+K+uxe4u4llwp191weIcvbpEJdzTQcinOe/A37n1f4K5PsHfgA8\n5jyfg68vNbd/dj2AMc7zeHyXu5+cayrwVmt9ngL9uQAXAO8AAkwEVrZyvnBgL9Dbi/0FnIXv9oKN\nftN+D/zUef7T+j7zQBKQ63xNdJ4nnm6ednWEoapbVHVrPbNmAy+oaqWq5gE5+PrCOkFEBDgHeMWZ\n9BRwiVtZne19G3jerW244HT6Dzslqvq+qtY4L1fgu8nTK4F8/7P5d48GrwDnOj9r16hqoaqudZ4f\nAbbQQFc7QWg28LT6rAA6i0iPVtz+ucAOVT3VG4JPi6ouw3fLgT//z1BDf4dc6Y+vXRWMRgTSd1UX\n4JDfH6cG+7dqIVOAfaq6vYH5CrwvImtEZK6LOfzd4pwWmN/AYXDAfYC55Hv4/hutT2vsr0C+/xPL\nOJ+lw/g+W63COQU2GlhZz+xJIrJORN4RkaGtFKmpn4vXn6k5NPxPmxf7C6CbOjc84zv66VbPMq7s\nt5Ab01tEPgC61zPrF6r6RmvnqU+AGa+k8aOLyapaICJdgUUi8qXz34gruYBHgfvx/YLfj+902fdO\nZ3stkev4/hKRX+C7AfTZBlbT4vurrRGRjsCrwB2qWnrS7LX4TruUOe1T/wL6t0KsoP25OG2UF1P/\nvWJe7a+vUVUVkVa71DXkCoaqnncKb2uy3yvgAL7D4QjnP8P6lmmRjOLrVuWbwNhG1lHgfC0Skdfx\nnQ45rV+0QPediPwNeKueWYHsxxbPJSLXARcB56pzAreedbT4/qpHIN//8WXynZ9zAr7PlqtEJBJf\nsXhWVV87eb5/AVHVhSLyiIgkq6qr/SYF8HNx5TMVoFnAWlXdd/IMr/aXY5843So5p+eK6lnmlPrj\na4qdkvJZAMxxrmDJxPefwir/BZw/RB8BlzuTrgXcOmI5D/hSVfPrmykiHUQk/vhzfA2/G+tbtqWc\ndN740ga2txroL76ryaLwHc4vcDnXTODHwMWqeqyBZVprfwXy/S/A99kB32dpcUNFrqU4bSR/B7ao\n6oMNLNP9eFuKiIzH97fB1UIW4M9lAXCNc7XUROCw3+kYtzV4lO/F/vLj/xlq6O/Qe8B0EUl0Th9P\nd6adHrdb+YPpge8PXT5QCewD3vOb9wt8V7hsBWb5TV8I9HSe98FXSHKAl4Fol3I+Cdx40rSewEK/\nHOucxyZ8p2bc3nfPABuA9c4HtsfJuZzXF+C7CmdHK+XKwXeu9gvn8djJuVpzf9X3/QP34StoADHO\nZyfH+Sz1aYV9NBnfqcT1fvvpAuDG458z4BZn36zDd/HAGa2Qq96fy0m5BN/InTucz1+W27mc7XbA\nVwAS/Ka1+v7CV7AKgWrnb9f38bV5fQhsBz4Akpxls4An/N77PedzlgN8tyXy2J3exhhjAmKnpIwx\nxgTECoYxxpiAWMEwxhgTECsYxhhjAmIFwxhjTECsYJiQIiKXntTT6Bfi6514lkvbu1FErnGeXyci\nPf3mPSEiQ1pgG8d7Cr6vBdY1RXy91rp6344JTXZZrQlpTv9E3wGmqWqdy9tagq9H3xbtPl1E7gXK\nVPWPLbS+DHw9rQ5rifWZ9sOOMEzIEpEBwD3A1ScXCxHJEN84Gs+KyBYReUVE4px554rI5+Ibp2G+\niEQ70//H+e98vYj80Zl2r4jcLSKX47tx6lnnqCZWRJaIM+6KiFzprG+jiPzOL0eZiDzgdGK3QkTq\n60ju5O+ro4j8w1nfehG5zG9dfxDfmBcfiMh4J0OuiFzcMnvVtGdWMExIcvpOeg64S1V3N7DYQOAR\nVR0MlAI/EJEYfHfaX6Gqw/H1t3aTiHTB11PAUFUdAfzGf0Wq+gqQDXxHVUeparlflp74xuo4BxgF\njBOR411SdwBWqOpIfH0o3RDAt/dLfF1kDHeyLPZb12JVHQoccTKe7+Q+7dNZxljBMKHqfmCTqr7Y\nyDJ7VPUT5/k/8XWhMRDIU9VtzvSn8A1icxioAP4uIt8E6u23qgHjgCWqWqy+jiufddYJUMW/O3Jc\nA2QEsL7z8HWXAYD6xjs4vq53necbgKWqWu08D2S9xjTKCoYJOSIyFbgMX38/jTm5Aa/BBj3nD/14\nfIMeXcS//zCfrmr9d0NiLafXg7T/uurw9ZmGczou5HqmNq3PCoYJKU7PnP8ArlHf6HKN6SUik5zn\nVwEf4+t8MkNE+jnTrwaWim8siQRVXQjcCYysZ31H8A2BerJVwNkikiwi4fh6QV3anO/rJIuAm4+/\nkAbGdDampVnBMKHmRqAr8OhJl9ZeUc+yW4GbRWQLvnGPH1XVCuC7wMsisgHff+qP4SsEb4nIenyF\n5Yf1rO9J4LHjjd7HJ6qvO+6f4usefx2wRk9vMK/fAIlOA/o6YNpprMuYgNlltaZdakuXltpltSZY\n2BGGMcGvDJjbUjfuAW8CrTEynAkxdoRhjDEmIHaEYYwxJiBWMIwxxgTECoYxxpiAWMEwxhgTECsY\nxhhjAmIFwxhjTED+PzCTROuSYv2mAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "z = np.linspace(zmin, zmax, 1000)\n", + "plt.plot(z, phi(z))\n", + "plt.xlabel('Z position [cm]')\n", + "plt.ylabel('Flux [n/src]')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "As you might expect, we get a rough cosine shape but with a flux depression in the middle due to the boron slab that we introduced. To get a more accurate distribution, we'd likely need to use a higher order expansion.\n", + "\n", + "One more thing we can do is confirm that integrating the distribution gives us the same value as the first moment (since $P_0(z') = 1$). This can easily be done by numerically integrating using the trapezoidal rule:" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "36.434786672754925" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.trapz(phi(z), z)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In addition to being able to tally Legendre moments, there are also functional expansion filters available for spherical harmonics (`SphericalHarmonicsFilter`) and Zernike polynomials over a unit disk (`ZernikeFilter`). A separate `LegendreFilter` class can also be used for determining Legendre scattering moments (i.e., an expansion of the scattering cosine, $\\mu$)." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.6.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/include/openmc.h b/include/openmc.h index 10c0b2de8e..1c2106432a 100644 --- a/include/openmc.h +++ b/include/openmc.h @@ -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; diff --git a/openmc/__init__.py b/openmc/__init__.py index 9b13fa6927..191528327c 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -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 * diff --git a/openmc/checkvalue.py b/openmc/checkvalue.py index 2f80ee4c0c..3b334319f6 100644 --- a/openmc/checkvalue.py +++ b/openmc/checkvalue.py @@ -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. """ diff --git a/openmc/filter.py b/openmc/filter.py index c1b9f8dc92..2bab17dbb6 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1,10 +1,9 @@ from abc import ABCMeta from collections import Iterable, OrderedDict import copy -from functools import reduce import hashlib +from itertools import product from numbers import Real, Integral -import operator from xml.etree import ElementTree as ET import numpy as np @@ -19,9 +18,12 @@ from .surface import Surface from .universe import Universe -_FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', - 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom'] +_FILTER_TYPES = ( + 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', + 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', + 'energyfunction', 'cellfrom', 'legendre', 'spatiallegendre', + 'sphericalharmonics', 'zernike' +) _CURRENT_NAMES = ( 'x-min out', 'x-min in', 'x-max out', 'x-max in', @@ -186,20 +188,15 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): def bins(self): return self._bins + @bins.setter + def bins(self, bins): + self.check_bins(bins) + self._bins = bins + @property def num_bins(self): return len(self.bins) - @bins.setter - def bins(self, bins): - # Format the bins as a 1D numpy array. - bins = np.atleast_1d(bins) - - # Check the bin values. - self.check_bins(bins) - - self._bins = bins - def check_bins(self, bins): """Make sure given bins are valid for this filter. @@ -221,8 +218,6 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): XML element containing filter data """ - - element = ET.Element('filter') element.set('id', str(self.id)) element.set('type', self.short_name.lower()) @@ -307,7 +302,7 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): Parameters ---------- - filter_bin : Integral or tuple + filter_bin : int or tuple The bin is the 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 @@ -318,13 +313,9 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): Returns ------- - filter_index : Integral + filter_index : int The index in the Tally data array for this filter bin. - See also - -------- - Filter.get_bin() - """ if filter_bin not in self.bins: @@ -332,47 +323,10 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): 'is not one of the bins'.format(filter_bin) raise ValueError(msg) - return np.where(self.bins == filter_bin)[0][0] - - def get_bin(self, bin_index): - """Returns the filter bin for some filter bin index. - - Parameters - ---------- - bin_index : Integral - The zero-based index into the filter's array of bins. The bin - index for 'material', 'surface', 'cell', 'cellborn', and 'universe' - filters corresponds to the ID in the filter's list of bins. For - 'distribcell' tallies the bin index necessarily can only be zero - since only one cell can be tracked per tally. The bin index for - 'energy' and 'energyout' filters corresponds to the energy range of - interest in the filter bins of energies. The bin index for 'mesh' - filters is the index into the flattened array of (x,y) or (x,y,z) - mesh cell bins. - - Returns - ------- - bin : 1-, 2-, or 3-tuple of Real - The bin in the Tally data array. The bin for 'material', surface', - 'cell', 'cellborn', 'universe' and 'distribcell' filters is a - 1-tuple of the ID corresponding to the appropriate filter bin. - The bin for 'energy' and 'energyout' filters is a 2-tuple of the - lower and upper energies bounding the energy interval for the filter - bin. The bin for 'mesh' tallies is a 2-tuple or 3-tuple of the x,y - or x,y,z mesh cell indices corresponding to the bin in a 2D/3D mesh. - - See also - -------- - Filter.get_bin_index() - - """ - - cv.check_type('bin_index', bin_index, Integral) - cv.check_greater_than('bin_index', bin_index, 0, equality=True) - cv.check_less_than('bin_index', bin_index, self.num_bins) - - # Return a 1-tuple of the bin. - return (self.bins[bin_index],) + if isinstance(self.bins, np.ndarray): + return np.where(self.bins == filter_bin)[0][0] + else: + return self.bins.index(filter_bin) def get_pandas_dataframe(self, data_size, stride, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. @@ -423,25 +377,24 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): class WithIDFilter(Filter): - """Abstract parent for filters of types with ids (Cell, Material, etc.).""" - - @Filter.bins.setter - def bins(self, bins): - # Format the bins as a 1D numpy array. + """Abstract parent for filters of types with IDs (Cell, Material, etc.).""" + def __init__(self, bins, filter_id=None): bins = np.atleast_1d(bins) - # Check the bin values. + # Make sure bins are either integers or appropriate objects cv.check_iterable_type('filter bins', bins, (Integral, self.expected_type)) + + # Extract ID values + bins = np.array([b if isinstance(b, Integral) else b.id + for b in bins]) + self.bins = bins + self.id = filter_id + + def check_bins(self, bins): + # Check the bin values. for edge in bins: - if isinstance(edge, Integral): - cv.check_greater_than('filter bin', edge, 0, equality=True) - - # Extract id values. - bins = np.atleast_1d([b if isinstance(b, Integral) else b.id - for b in bins]) - - self._bins = bins + cv.check_greater_than('filter bin', edge, 0, equality=True) class UniverseFilter(WithIDFilter): @@ -449,7 +402,7 @@ class UniverseFilter(WithIDFilter): Parameters ---------- - bins : openmc.Universe, Integral, or iterable thereof + bins : openmc.Universe, int, or iterable thereof The Universes to tally. Either openmc.Universe objects or their Integral ID numbers can be used. filter_id : int @@ -601,12 +554,13 @@ class MeshFilter(Filter): Attributes ---------- - bins : Integral - The Mesh ID mesh : openmc.Mesh The Mesh object that events will be tallied onto id : int Unique identifier for the filter + bins : list of tuple + A list of mesh indices for each filter bin, e.g. [(1, 1, 1), (2, 1, 1), + ...] num_bins : Integral The number of filter bins @@ -614,7 +568,13 @@ class MeshFilter(Filter): def __init__(self, mesh, filter_id=None): self.mesh = mesh - super().__init__(mesh.id, filter_id) + self.id = filter_id + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tMesh ID', self.mesh.id) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + return string @classmethod def from_hdf5(cls, group, **kwargs): @@ -643,69 +603,12 @@ class MeshFilter(Filter): def mesh(self, mesh): cv.check_type('filter mesh', mesh, openmc.Mesh) self._mesh = mesh - self.bins = mesh.id - - @property - def num_bins(self): - return reduce(operator.mul, self.mesh.dimension) - - def check_bins(self, bins): - if not len(bins) == 1: - msg = 'Unable to add bins "{0}" to a MeshFilter since ' \ - 'only a single mesh can be used per tally'.format(bins) - raise ValueError(msg) - elif not isinstance(bins[0], Integral): - msg = 'Unable to add bin "{0}" to MeshFilter since it ' \ - 'is a non-integer'.format(bins[0]) - raise ValueError(msg) - elif bins[0] < 0: - msg = 'Unable to add bin "{0}" to MeshFilter since it ' \ - 'is a negative integer'.format(bins[0]) - raise ValueError(msg) + self.bins = list(mesh.indices) def can_merge(self, other): # Mesh filters cannot have more than one bin return False - def get_bin_index(self, filter_bin): - # Filter bins for a mesh are an (x,y,z) tuple. Convert (x,y,z) to a - # single bin -- this is similar to subroutine mesh_indices_to_bin in - # openmc/src/mesh.F90. - n_dim = len(self.mesh.dimension) - if n_dim == 3: - i, j, k = filter_bin - nx, ny, nz = self.mesh.dimension - return (i - 1) + (j - 1)*nx + (k - 1)*nx*ny - elif n_dim == 2: - i, j, *_ = filter_bin - nx, ny = self.mesh.dimension - return (i - 1) + (j - 1)*nx - else: - return filter_bin[0] - 1 - - def get_bin(self, bin_index): - cv.check_type('bin_index', bin_index, Integral) - cv.check_greater_than('bin_index', bin_index, 0, equality=True) - cv.check_less_than('bin_index', bin_index, self.num_bins) - - n_dim = len(self.mesh.dimension) - if n_dim == 3: - # Construct 3-tuple of x,y,z cell indices for a 3D mesh - nx, ny, nz = self.mesh.dimension - x = (bin_index % nx) + 1 - y = (bin_index % (nx * ny)) // nx + 1 - z = bin_index // (nx * ny) + 1 - return (x, y, z) - - elif n_dim == 2: - # Construct 2-tuple of x,y cell indices for a 2D mesh - nx, ny = self.mesh.dimension - x = (bin_index % nx) + 1 - y = bin_index // nx + 1 - return (x, y) - else: - return (bin_index + 1,) - def get_pandas_dataframe(self, data_size, stride, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. @@ -783,6 +686,19 @@ class MeshFilter(Filter): return df + def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing filter data + + """ + element = super().to_xml_element() + element[0].text = str(self.mesh.id) + return element + class MeshSurfaceFilter(MeshFilter): """Filter events by surface crossings on a regular, rectangular mesh. @@ -802,36 +718,25 @@ class MeshSurfaceFilter(MeshFilter): The Mesh object that events will be tallied onto id : int Unique identifier for the filter + bins : list of tuple + + A list of mesh indices / surfaces for each filter bin, e.g. [(1, 1, + 'x-min out'), (1, 1, 'x-min in'), ...] + num_bins : Integral The number of filter bins """ - @property - def num_bins(self): - n_dim = len(self.mesh.dimension) - return 4*n_dim*reduce(operator.mul, self.mesh.dimension) + @MeshFilter.mesh.setter + def mesh(self, mesh): + cv.check_type('filter mesh', mesh, openmc.Mesh) + self._mesh = mesh - def get_bin_index(self, filter_bin): - # Split bin into mesh/surface parts - *mesh_tuple, surf = filter_bin - - # Get index for mesh alone - mesh_index = super().get_bin_index(mesh_tuple) - - # Combine surface and mesh index - n_dim = len(self.mesh.dimension) - for surf_index, name in enumerate(_CURRENT_NAMES): - if surf == name: - return 4*n_dim*mesh_index + surf_index - else: - raise ValueError("'{}' is not a valid mesh surface.".format(surf)) - - def get_bin(self, bin_index): - n_dim = len(self.mesh.dimension) - mesh_index, surf_index = divmod(bin_index, 4*n_dim) - mesh_bin = super().get_bin(mesh_index) - return mesh_bin + (_CURRENT_NAMES[surf_index],) + # Take the product of mesh indices and current names + n_dim = len(mesh.dimension) + self.bins = [mesh_tuple + (surf,) for mesh_tuple, surf in + product(mesh.indices, _CURRENT_NAMES[:4*n_dim])] def get_pandas_dataframe(self, data_size, stride, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. @@ -920,42 +825,74 @@ class RealFilter(Filter): Parameters ---------- - bins : Iterable of Real - A grid of bin values. + values : iterable of float + A list of values for which each successive pair constitutes a range of + values for a single bin filter_id : int Unique identifier for the filter Attributes ---------- - bins : Iterable of Real - A grid of bin values. + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + values for a single bin id : int Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of values indicating a + filter bin range num_bins : int The number of filter bins """ + def __init__(self, values, filter_id=None): + self.values = np.asarray(values) + self.bins = np.vstack((self.values[:-1], self.values[1:])).T + self.id = filter_id def __gt__(self, other): if type(self) is type(other): # Compare largest/smallest bin edges in filters # This logic is used when merging tallies with real filters - return self.bins[0] >= other.bins[-1] + return self.values[0] >= other.values[-1] else: return super().__gt__(other) - @property - def num_bins(self): - return len(self.bins) - 1 + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tValues', self.values) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + return string + + @Filter.bins.setter + def bins(self, bins): + Filter.bins.__set__(self, np.asarray(bins)) + + def check_bins(self, bins): + for v0, v1 in bins: + # Values should be real + cv.check_type('filter value', v0, Real) + cv.check_type('filter value', v1, Real) + + # Make sure that each tuple has values that are increasing + if v1 < v0: + raise ValueError('Values {} and {} appear to be out of order' + .format(v0, v1)) + + for pair0, pair1 in zip(bins[:-1], bins[1:]): + # Successive pairs should be ordered + if pair1[1] < pair0[1]: + raise ValueError('Values {} and {} appear to be out of order' + .format(pair1[1], pair0[1])) def can_merge(self, other): if type(self) is not type(other): return False - if self.bins[0] == other.bins[-1]: + if self.bins[0, 0] == other.bins[-1][1]: # This low edge coincides with other's high edge return True - elif self.bins[-1] == other.bins[0]: + elif self.bins[-1][1] == other.bins[0, 0]: # This high edge coincides with other's low edge return True else: @@ -968,11 +905,11 @@ class RealFilter(Filter): raise ValueError(msg) # Merge unique filter bins - merged_bins = np.concatenate((self.bins, other.bins)) - merged_bins = np.unique(merged_bins) + merged_values = np.concatenate((self.values, other.values)) + merged_values = np.unique(merged_values) # Create a new filter with these bins and a new auto-generated ID - return type(self)(sorted(merged_bins)) + return type(self)(sorted(merged_values)) def is_subset(self, other): """Determine if another filter is a subset of this filter. @@ -994,80 +931,19 @@ class RealFilter(Filter): if type(self) is not type(other): return False - elif len(self.bins) != len(other.bins): + elif self.num_bins != other.num_bins: return False else: - return np.allclose(self.bins, other.bins) + return np.allclose(self.values, other.values) def get_bin_index(self, filter_bin): - i = np.where(self.bins == filter_bin[1])[0] + i = np.where(self.bins[:, 1] == filter_bin[1])[0] if len(i) == 0: msg = 'Unable to get the bin index for Filter since "{0}" ' \ 'is not one of the bins'.format(filter_bin) raise ValueError(msg) else: - return i[0] - 1 - - def get_bin(self, bin_index): - cv.check_type('bin_index', bin_index, Integral) - cv.check_greater_than('bin_index', bin_index, 0, equality=True) - cv.check_less_than('bin_index', bin_index, self.num_bins) - - # Construct 2-tuple of lower, upper bins for real-valued filters - return (self.bins[bin_index], self.bins[bin_index + 1]) - - -class EnergyFilter(RealFilter): - """Bins tally events based on incident particle energy. - - Parameters - ---------- - bins : Iterable of Real - A grid of energy values in eV. - filter_id : int - Unique identifier for the filter - - Attributes - ---------- - bins : Iterable of Real - A grid of energy values in eV. - id : int - Unique identifier for the filter - num_bins : int - The number of filter bins - - """ - - def get_bin_index(self, filter_bin): - # Use lower energy bound to find index for RealFilters - deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1] - min_delta = np.min(deltas) - if min_delta < 1E-3: - return deltas.argmin() - 1 - else: - msg = 'Unable to get the bin index for Filter since "{0}" ' \ - 'is not one of the bins'.format(filter_bin) - raise ValueError(msg) - - def check_bins(self, bins): - for edge in bins: - if not isinstance(edge, Real): - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is a non-integer or floating point ' \ - 'value'.format(edge, type(self)) - raise ValueError(msg) - elif edge < 0.: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is a negative value'.format(edge, type(self)) - raise ValueError(msg) - - # Check that bin edges are monotonically increasing - for index in range(1, len(bins)): - if bins[index] < bins[index-1]: - msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ - 'since they are not monotonically ' \ - 'increasing'.format(bins, type(self)) - raise ValueError(msg) + return i[0] def get_pandas_dataframe(self, data_size, stride, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. @@ -1102,35 +978,103 @@ class EnergyFilter(RealFilter): # Extract the lower and upper energy bounds, then repeat and tile # them as necessary to account for other filters. - lo_bins = np.repeat(self.bins[:-1], stride) - hi_bins = np.repeat(self.bins[1:], stride) + lo_bins = np.repeat(self.bins[:, 0], stride) + hi_bins = np.repeat(self.bins[:, 1], stride) tile_factor = data_size // len(lo_bins) lo_bins = np.tile(lo_bins, tile_factor) hi_bins = np.tile(hi_bins, tile_factor) # Add the new energy columns to the DataFrame. - df.loc[:, self.short_name.lower() + ' low [eV]'] = lo_bins - df.loc[:, self.short_name.lower() + ' high [eV]'] = hi_bins + if hasattr(self, 'units'): + units = ' [{}]'.format(self.units) + else: + units = '' + + df.loc[:, self.short_name.lower() + ' low' + units] = lo_bins + df.loc[:, self.short_name.lower() + ' high' + units] = hi_bins return df + def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing filter data + + """ + element = super().to_xml_element() + element[0].text = ' '.join(str(x) for x in self.values) + return element + + +class EnergyFilter(RealFilter): + """Bins tally events based on incident particle energy. + + Parameters + ---------- + values : Iterable of Real + A list of values for which each successive pair constitutes a range of + energies in [eV] for a single bin + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + energies in [eV] for a single bin + id : int + Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of energies in [eV] + for a single filter bin + num_bins : int + The number of filter bins + + """ + units = 'eV' + + def get_bin_index(self, filter_bin): + # Use lower energy bound to find index for RealFilters + deltas = np.abs(self.bins[:, 1] - filter_bin[1]) / filter_bin[1] + min_delta = np.min(deltas) + if min_delta < 1E-3: + return deltas.argmin() + else: + msg = 'Unable to get the bin index for Filter since "{0}" ' \ + 'is not one of the bins'.format(filter_bin) + raise ValueError(msg) + + def check_bins(self, bins): + super().check_bins(bins) + for v0, v1 in bins: + cv.check_greater_than('filter value', v0, 0., equality=True) + cv.check_greater_than('filter value', v1, 0., equality=True) + class EnergyoutFilter(EnergyFilter): """Bins tally events based on outgoing particle energy. Parameters ---------- - bins : Iterable of Real - A grid of energy values in eV. + values : Iterable of Real + A list of values for which each successive pair constitutes a range of + energies in [eV] for a single bin filter_id : int Unique identifier for the filter Attributes ---------- - bins : Iterable of Real - A grid of energy values in eV. + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + energies in [eV] for a single bin id : int Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of energies in [eV] + for a single filter bin num_bins : int The number of filter bins @@ -1412,98 +1356,41 @@ class MuFilter(RealFilter): Parameters ---------- - bins : Iterable of Real or Integral - A grid of scattering angles which events will binned into. Values - represent the cosine of the scattering angle. If an Iterable is given, - the values will be used explicitly as grid points. If a single Integral - is given, the range [-1, 1] will be divided up equally into that number - of bins. + values : int or Iterable of Real + A grid of scattering angles which events will binned into. Values + represent the cosine of the scattering angle. If an iterable is given, + the values will be used explicitly as grid points. If a single int is + given, the range [-1, 1] will be divided up equally into that number of + bins. filter_id : int Unique identifier for the filter Attributes ---------- - bins : Integral - A grid of scattering angles which events will binned into. Values - represent the cosine of the scattering angle. If an Iterable is given, - the values will be used explicitly as grid points. If a single Integral - is given, the range [-1, 1] will be divided up equally into that number - of bins. + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + scattering angle cosines for a single bin id : int Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of scattering angle + cosines for a single filter bin num_bins : Integral The number of filter bins """ + def __init__(self, values, filter_id=None): + if isinstance(values, Integral): + values = np.linspace(-1., 1., values + 1) + super().__init__(values, filter_id) def check_bins(self, bins): - for edge in bins: - if not isinstance(edge, Real): - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is a non-integer or floating point ' \ - 'value'.format(edge, type(self)) - raise ValueError(msg) - elif edge < -1.: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is less than -1'.format(edge, type(self)) - raise ValueError(msg) - elif edge > 1.: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is greater than 1'.format(edge, type(self)) - raise ValueError(msg) - - # Check that bin edges are monotonically increasing - for index in range(1, len(bins)): - if bins[index] < bins[index-1]: - msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ - 'since they are not monotonically ' \ - 'increasing'.format(bins, type(self)) - raise ValueError(msg) - - def get_pandas_dataframe(self, data_size, stride, **kwargs): - """Builds a Pandas DataFrame for the Filter's bins. - - This method constructs a Pandas DataFrame object for the filter with - columns annotated by filter bin information. This is a helper method - for :meth:`Tally.get_pandas_dataframe`. - - Parameters - ---------- - data_size : int - The total number of bins in the tally corresponding to this filter - stride : int - Stride in memory for the filter - - Returns - ------- - pandas.DataFrame - A Pandas DataFrame with one column of the lower energy bound and one - column of upper energy bound for each filter bin. The number of - rows in the DataFrame is the same as the total number of bins in the - corresponding tally, with the filter bin appropriately tiled to map - to the corresponding tally bins. - - See also - -------- - Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() - - """ - # Initialize Pandas DataFrame - df = pd.DataFrame() - - # Extract the lower and upper energy bounds, then repeat and tile - # them as necessary to account for other filters. - lo_bins = np.repeat(self.bins[:-1], stride) - hi_bins = np.repeat(self.bins[1:], stride) - tile_factor = data_size // len(lo_bins) - lo_bins = np.tile(lo_bins, tile_factor) - hi_bins = np.tile(hi_bins, tile_factor) - - # Add the new energy columns to the DataFrame. - df.loc[:, self.short_name.lower() + ' low'] = lo_bins - df.loc[:, self.short_name.lower() + ' high'] = hi_bins - - return df + super().check_bins(bins) + for x in np.ravel(bins): + if not np.isclose(x, -1.): + cv.check_greater_than('filter value', x, -1., equality=True) + if not np.isclose(x, 1.): + cv.check_less_than('filter value', x, 1., equality=True) class PolarFilter(RealFilter): @@ -1511,98 +1398,44 @@ class PolarFilter(RealFilter): Parameters ---------- - bins : Iterable of Real or Integral - A grid of polar angles which events will binned into. Values represent - an angle in radians relative to the z-axis. If an Iterable is given, - the values will be used explicitly as grid points. If a single Integral - is given, the range [0, pi] will be divided up equally into that number - of bins. + values : int or Iterable of Real + A grid of polar angles which events will binned into. Values represent + an angle in radians relative to the z-axis. If an iterable is given, the + values will be used explicitly as grid points. If a single int is given, + the range [0, pi] will be divided up equally into that number of bins. filter_id : int Unique identifier for the filter Attributes ---------- - bins : Iterable of Real or Integral - A grid of polar angles which events will binned into. Values represent - an angle in radians relative to the z-axis. If an Iterable is given, - the values will be used explicitly as grid points. If a single Integral - is given, the range [0, pi] will be divided up equally into that number - of bins. + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + polar angles in [rad] for a single bin + id : int + Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of polar angles for a + single filter bin id : int Unique identifier for the filter num_bins : Integral The number of filter bins """ + units = 'rad' + + def __init__(self, values, filter_id=None): + if isinstance(values, Integral): + values = np.linspace(0., np.pi, values + 1) + super().__init__(values, filter_id) def check_bins(self, bins): - for edge in bins: - if not isinstance(edge, Real): - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is a non-integer or floating point ' \ - 'value'.format(edge, type(self)) - raise ValueError(msg) - elif edge < 0.: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is less than 0'.format(edge, type(self)) - raise ValueError(msg) - elif not np.isclose(edge, np.pi) and edge > np.pi: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is greater than pi'.format(edge, type(self)) - raise ValueError(msg) - - # Check that bin edges are monotonically increasing - for index in range(1, len(bins)): - if bins[index] < bins[index-1]: - msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ - 'since they are not monotonically ' \ - 'increasing'.format(bins, type(self)) - raise ValueError(msg) - - def get_pandas_dataframe(self, data_size, stride, **kwargs): - """Builds a Pandas DataFrame for the Filter's bins. - - This method constructs a Pandas DataFrame object for the filter with - columns annotated by filter bin information. This is a helper method for - :meth:`Tally.get_pandas_dataframe`. - - Parameters - ---------- - data_size : int - The total number of bins in the tally corresponding to this filter - stride : int - Stride in memory for the filter - - Returns - ------- - pandas.DataFrame - A Pandas DataFrame with a column corresponding to the lower polar - angle bound for each of the filter's bins. The number of rows in - the DataFrame is the same as the total number of bins in the - corresponding tally, with the filter bin appropriately tiled to map - to the corresponding tally bins. - - See also - -------- - Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() - - """ - # Initialize Pandas DataFrame - df = pd.DataFrame() - - # Extract the lower and upper angle bounds, then repeat and tile - # them as necessary to account for other filters. - lo_bins = np.repeat(self.bins[:-1], stride) - hi_bins = np.repeat(self.bins[1:], stride) - tile_factor = data_size // len(lo_bins) - lo_bins = np.tile(lo_bins, tile_factor) - hi_bins = np.tile(hi_bins, tile_factor) - - # Add the new angle columns to the DataFrame. - df.loc[:, 'polar low'] = lo_bins - df.loc[:, 'polar high'] = hi_bins - - return df + super().check_bins(bins) + for x in np.ravel(bins): + if not np.isclose(x, 0.): + cv.check_greater_than('filter value', x, 0., equality=True) + if not np.isclose(x, np.pi): + cv.check_less_than('filter value', x, np.pi, equality=True) class AzimuthalFilter(RealFilter): @@ -1610,98 +1443,43 @@ class AzimuthalFilter(RealFilter): Parameters ---------- - bins : Iterable of Real or Integral - A grid of azimuthal angles which events will binned into. Values + values : int or Iterable of Real + A grid of azimuthal angles which events will binned into. Values represent an angle in radians relative to the x-axis and perpendicular - to the z-axis. If an Iterable is given, the values will be used - explicitly as grid points. If a single Integral is given, the range + to the z-axis. If an iterable is given, the values will be used + explicitly as grid points. If a single int is given, the range [-pi, pi) will be divided up equally into that number of bins. filter_id : int Unique identifier for the filter Attributes ---------- - bins : Iterable of Real or Integral - A grid of azimuthal angles which events will binned into. Values - represent an angle in radians relative to the x-axis and perpendicular - to the z-axis. If an Iterable is given, the values will be used - explicitly as grid points. If a single Integral is given, the range - [-pi, pi) will be divided up equally into that number of bins. + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + azimuthal angles in [rad] for a single bin id : int Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of azimuthal angles + for a single filter bin num_bins : Integral The number of filter bins """ + units = 'rad' + + def __init__(self, values, filter_id=None): + if isinstance(values, Integral): + values = np.linspace(-np.pi, np.pi, values + 1) + super().__init__(values, filter_id) def check_bins(self, bins): - for edge in bins: - if not isinstance(edge, Real): - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is a non-integer or floating point ' \ - 'value'.format(edge, type(self)) - raise ValueError(msg) - elif not np.isclose(edge, -np.pi) and edge < -np.pi: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is less than -pi'.format(edge, type(self)) - raise ValueError(msg) - elif not np.isclose(edge, np.pi) and edge > np.pi: - msg = 'Unable to add bin edge "{0}" to a "{1}" ' \ - 'since it is greater than pi'.format(edge, type(self)) - raise ValueError(msg) - - # Check that bin edges are monotonically increasing - for index in range(1, len(bins)): - if bins[index] < bins[index-1]: - msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \ - 'since they are not monotonically ' \ - 'increasing'.format(bins, type(self)) - raise ValueError(msg) - - def get_pandas_dataframe(self, data_size, stride, paths=True): - """Builds a Pandas DataFrame for the Filter's bins. - - This method constructs a Pandas DataFrame object for the filter with - columns annotated by filter bin information. This is a helper method for - :meth:`Tally.get_pandas_dataframe`. - - Parameters - ---------- - data_size : int - The total number of bins in the tally corresponding to this filter - stride : int - Stride in memory for the filter - - Returns - ------- - pandas.DataFrame - A Pandas DataFrame with a column corresponding to the lower - azimuthal angle bound for each of the filter's bins. The number of - rows in the DataFrame is the same as the total number of bins in the - corresponding tally, with the filter bin appropriately tiled to map - to the corresponding tally bins. - - See also - -------- - Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() - - """ - # Initialize Pandas DataFrame - df = pd.DataFrame() - - # Extract the lower and upper angle bounds, then repeat and tile - # them as necessary to account for other filters. - lo_bins = np.repeat(self.bins[:-1], stride) - hi_bins = np.repeat(self.bins[1:], stride) - tile_factor = data_size // len(lo_bins) - lo_bins = np.tile(lo_bins, tile_factor) - hi_bins = np.tile(hi_bins, tile_factor) - - # Add the new angle columns to the DataFrame. - df.loc[:, 'azimuthal low'] = lo_bins - df.loc[:, 'azimuthal high'] = hi_bins - - return df + super().check_bins(bins) + for x in np.ravel(bins): + if not np.isclose(x, -np.pi): + cv.check_greater_than('filter value', x, -np.pi, equality=True) + if not np.isclose(x, np.pi): + cv.check_less_than('filter value', x, np.pi, equality=True) class DelayedGroupFilter(Filter): @@ -1709,7 +1487,7 @@ class DelayedGroupFilter(Filter): Parameters ---------- - bins : Integral or Iterable of Integral + bins : iterable of int The delayed neutron precursor groups. For example, ENDF/B-VII.1 uses 6 precursor groups so a tally with all groups will have bins = [1, 2, 3, 4, 5, 6]. @@ -1718,7 +1496,7 @@ class DelayedGroupFilter(Filter): Attributes ---------- - bins : Integral or Iterable of Integral + bins : iterable of int The delayed neutron precursor groups. For example, ENDF/B-VII.1 uses 6 precursor groups so a tally with all groups will have bins = [1, 2, 3, 4, 5, 6]. @@ -1728,17 +1506,10 @@ class DelayedGroupFilter(Filter): The number of filter bins """ - @Filter.bins.setter - def bins(self, bins): - # Format the bins as a 1D numpy array. - bins = np.atleast_1d(bins) - + def check_bins(self, bins): # Check the bin values. - cv.check_iterable_type('filter bins', bins, Integral) - for edge in bins: - cv.check_greater_than('filter bin', edge, 0, equality=True) - - self._bins = bins + for g in bins: + cv.check_greater_than('delayed group', g, 0) class EnergyFunctionFilter(Filter): @@ -1751,18 +1522,18 @@ class EnergyFunctionFilter(Filter): Parameters ---------- energy : Iterable of Real - A grid of energy values in eV. + A grid of energy values in [eV] y : iterable of Real - A grid of interpolant values in eV. + A grid of interpolant values in [eV] filter_id : int Unique identifier for the filter Attributes ---------- energy : Iterable of Real - A grid of energy values in eV. + A grid of energy values in [eV] y : iterable of Real - A grid of interpolant values in eV. + A grid of interpolant values in [eV] id : int Unique identifier for the filter num_bins : Integral @@ -1903,6 +1674,14 @@ class EnergyFunctionFilter(Filter): raise RuntimeError('EnergyFunctionFilters have no bins.') def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing filter data + + """ element = ET.Element('filter') element.set('id', str(self.id)) element.set('type', self.short_name.lower()) @@ -1925,10 +1704,6 @@ class EnergyFunctionFilter(Filter): # This filter only has one bin. Always return 0. return 0 - def get_bin(self, bin_index): - """This function is invalid for EnergyFunctionFilters.""" - raise RuntimeError('EnergyFunctionFilters have no get_bin() method') - def get_pandas_dataframe(self, data_size, stride, **kwargs): """Builds a Pandas DataFrame for the Filter's bins. diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py new file mode 100644 index 0000000000..872eaac9f5 --- /dev/null +++ b/openmc/filter_expansion.py @@ -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 + `_ 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 diff --git a/openmc/mesh.py b/openmc/mesh.py index cd8eb3c5ed..aed7a46d8c 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -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), ...] """ diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 4932feff1e..5bc00cbc95 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -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 diff --git a/openmc/tallies.py b/openmc/tallies.py index 6c055d1913..50398b786e 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -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: diff --git a/src/constants.F90 b/src/constants.F90 index 0e138f0ee5..b335b78c86 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -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 :: & diff --git a/src/math.F90 b/src/math.F90 index 2fa0a6ce1e..ee8cd0530b 100644 --- a/src/math.F90 +++ b/src/math.F90 @@ -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 !=============================================================================== diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index 3d6eaaef21..e484ac2315 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -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_)) diff --git a/src/tallies/tally_filter_energy.F90 b/src/tallies/tally_filter_energy.F90 index caba6755d1..d186fa62a8 100644 --- a/src/tallies/tally_filter_energy.F90 +++ b/src/tallies/tally_filter_energy.F90 @@ -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 diff --git a/src/tallies/tally_filter_header.F90 b/src/tallies/tally_filter_header.F90 index 34ad75388f..93d050b1af 100644 --- a/src/tallies/tally_filter_header.F90 +++ b/src/tallies/tally_filter_header.F90 @@ -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 !=============================================================================== diff --git a/src/tallies/tally_filter_legendre.F90 b/src/tallies/tally_filter_legendre.F90 new file mode 100644 index 0000000000..9545b86926 --- /dev/null +++ b/src/tallies/tally_filter_legendre.F90 @@ -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 diff --git a/src/tallies/tally_filter_material.F90 b/src/tallies/tally_filter_material.F90 index 778fec4b68..3b9f843f84 100644 --- a/src/tallies/tally_filter_material.F90 +++ b/src/tallies/tally_filter_material.F90 @@ -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 diff --git a/src/tallies/tally_filter_sph_harm.F90 b/src/tallies/tally_filter_sph_harm.F90 new file mode 100644 index 0000000000..5c47c2c71f --- /dev/null +++ b/src/tallies/tally_filter_sph_harm.F90 @@ -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 diff --git a/src/tallies/tally_filter_sptl_legendre.F90 b/src/tallies/tally_filter_sptl_legendre.F90 new file mode 100644 index 0000000000..d6594d05a9 --- /dev/null +++ b/src/tallies/tally_filter_sptl_legendre.F90 @@ -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 diff --git a/src/tallies/tally_filter_zernike.F90 b/src/tallies/tally_filter_zernike.F90 new file mode 100644 index 0000000000..e96a53a7f9 --- /dev/null +++ b/src/tallies/tally_filter_zernike.F90 @@ -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 diff --git a/src/tallies/tally_header.F90 b/src/tallies/tally_header.F90 index 0495ce14dc..02a578d7e0 100644 --- a/src/tallies/tally_header.F90 +++ b/src/tallies/tally_header.F90 @@ -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 diff --git a/tests/regression_tests/tally_slice_merge/test.py b/tests/regression_tests/tally_slice_merge/test.py index e52d0fde49..20ab57a077 100644 --- a/tests/regression_tests/tally_slice_merge/test.py +++ b/tests/regression_tests/tally_slice_merge/test.py @@ -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) diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py new file mode 100644 index 0000000000..488badbfaf --- /dev/null +++ b/tests/unit_tests/test_filters.py @@ -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)