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.
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
$$ \phi(z') = \sum\limits_{n=0}^N a_n P_n(z') $$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
$$ a_n = \frac{2n + 1}{2} \int_{-1}^1 dz' P_n(z') \phi(z').$$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.
In [6]:
%matplotlib inline
import openmc
import numpy as np
import matplotlib.pyplot as plt
# Define fuel and B4C materials
fuel = openmc.Material()
fuel.add_element('U', 1.0, enrichment=4.5)
fuel.add_nuclide('O16', 2.0)
fuel.set_density('g/cm3', 10.0)
b4c = openmc.Material()
b4c.add_element('B', 4.0)
b4c.add_nuclide('C0', 1.0)
b4c.set_density('g/cm3', 2.5)
# Define surfaces used to construct regions
zmin, zmax = -10., 10.
box = openmc.model.get_rectangular_prism(10., 10., boundary_type='reflective')
bottom = openmc.ZPlane(z0=zmin, boundary_type='vacuum')
boron_lower = openmc.ZPlane(z0=-0.5)
boron_upper = openmc.ZPlane(z0=0.5)
top = openmc.ZPlane(z0=zmax, boundary_type='vacuum')
# Create three cells and add them to geometry
fuel1 = openmc.Cell(fill=fuel, region=box & +bottom & -boron_lower)
absorber = openmc.Cell(fill=b4c, region=box & +boron_lower & -boron_upper)
fuel2 = openmc.Cell(fill=fuel, region=box & +boron_upper & -top)
geom = openmc.Geometry([fuel1, absorber, fuel2])
settings = openmc.Settings()
spatial_dist = openmc.stats.Box(*geom.bounding_box)
settings.source = openmc.Source(space=spatial_dist)
settings.batches = 10
settings.inactive = 0
settings.particles = 1000
# Create a flux tally
flux_tally = openmc.Tally()
flux_tally2 = openmc.Tally()
flux_tally.scores = ['flux']
flux_tally2.scores = ['flux']
# Create a Legendre polynomial expansion filter and add to tally
order = 8
expand_filter = openmc.SpatialLegendreFilter(order, 'z', zmin, zmax)
cell_filter = openmc.CellFilter([absorber, fuel2])
cell_filter2 = openmc.CellFilter([fuel2, fuel1])
flux_tally.filters.append(cell_filter)
flux_tally.filters.append(expand_filter)
flux_tally2.filters.append(expand_filter)
flux_tally2.filters.append(cell_filter2)
tallies = openmc.Tallies([flux_tally, flux_tally2])
model = openmc.model.Model(geometry=geom, settings=settings, tallies=tallies)
model.export_to_xml()
import openmc.capi
openmc.capi.init()
openmc.capi.run()
In [5]:
openmc.capi.finalize()
In [7]:
tallies = openmc.capi.tallies
In [8]:
tallies
Out[8]:
In [9]:
openmc.capi.cells
Out[9]:
In [15]:
results(tallies[3], 4)
In [13]:
results(tallies[3], 5)
Out[13]:
In [14]:
results(tallies[3], 6)
Out[14]:
In [18]:
results(tallies[4],4)
Out[18]:
In [19]:
results(tallies[4],5)
In [20]:
results(tallies[4],6)
Out[20]:
In [17]:
from ctypes import c_int, c_int32, POINTER
from openmc.capi.filter import SpatialLegendreFilter, ZernikeFilter, SphericalHarmonicsFilter, CellFilter
expansion_types = (SpatialLegendreFilter, ZernikeFilter, SphericalHarmonicsFilter)
def results(tally, cell_id):
filters = tally.filters
if len(filters) != 2:
raise("We expect there to be two filters, "
"one a cell filter and the other an expansion filter")
index_to_id = {}
for key,value in openmc.capi.cells.items():
index_to_id[value._index] = key
if isinstance(filters[0], CellFilter):
cells = filters[0].bins
cell_ids = [index_to_id[cell_index] for cell_index in cells]
if cell_id not in cell_ids:
raise RuntimeError("Requested cell_id not in the passed tally")
stride_integer = cell_ids.index(cell_id)
if not isinstance(filters[1], expansion_types):
raise TypeError("Expected an expansion filter "
"as the second filter")
num_bins = filters[1].order + 1
starting_point = num_bins * stride_integer
return tally.mean[starting_point:starting_point+num_bins]
elif isinstance(filters[0], expansion_types):
num_bins = filters[0].order + 1
if not isinstance(filters[1], CellFilter):
raise TypeError("Expected a cell filter as the second filter")
cells = filters[1].bins
cell_ids = [index_to_id[cell_index] for cell_index in cells]
if cell_id not in cell_ids:
raise RuntimeError("Requested cell_id not in the passed tally")
stride_integer = cell_ids.index(cell_id)
total_bins = cells.size * num_bins
return tally.mean[stride_integer:stride_integer+total_bins+cells.size:cells.size]
In [ ]:
model.run(output=True, openmc_exec='/Users/lindad/projects/Okapi/openmc/build/bin/openmc')
Now that the run is finished, we need to load the results from the statepoint file.
In [3]:
with openmc.StatePoint('statepoint.210.h5') as sp:
df = sp.tallies[flux_tally.id].get_pandas_dataframe()
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.
In [4]:
df
Out[4]:
Since the expansion coefficients are given as
$$ a_n = \frac{2n + 1}{2} \int_{-1}^1 dz' P_n(z') \phi(z')$$we just need to multiply the Legendre moments by $(2n + 1)/2$.
In [10]:
n = np.arange(order + 1)
a_n = (2*n + 1)/2 * df['mean']
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
and, for this case, $z = 10z'$, it follows that
$$ \phi(z) = \frac{\phi(z')}{10} = \sum_{n=0}^N \frac{a_n}{10} P_n(z'). $$
In [11]:
phi = np.polynomial.Legendre(a_n/10, domain=(zmin, zmax))
Let's plot it and see how our flux looks!
In [12]:
z = np.linspace(zmin, zmax, 1000)
plt.plot(z, phi(z))
plt.xlabel('Z position [cm]')
plt.ylabel('Flux [n/src]')
Out[12]:
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.
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:
In [13]:
np.trapz(phi(z), z)
Out[13]:
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$).