Discrepancy Between Legendre and Mesh Flux Tally Magnitudes

User Support
1

I am a bit confused on how the calculation of the expansion coefficients work for the spatial legendre filter. I ran the slab geometry example from the functional expansion jupyter notebook with an additional flux tally on a 1x1x40 mesh for comparison. The flux calculated from using the legendre polynomial expansion is the same shape as the flux calculated on the mesh but is about twice as large.

image
image2723×1820 210 KB

I have repeated this on a couple more geometries. The shape always matches but magnitude is always off by some (problem dependent) factor. Is there something that I am missing with regard to calculating the coefficients?

Here is the python code I ran to generate the plot from above:

import openmc
import numpy as np
import matplotlib.pyplot as plt

plt.rcParams['figure.dpi'] = 500
plt.rcParams['savefig.dpi'] = 500

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_element('C', 1.0)
b4c.set_density('g/cm3', 2.5)

zmin, zmax = -10., 10.
box = openmc.model.RectangularPrism(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.IndependentSource(space=spatial_dist)
settings.batches = 210
settings.inactive = 10
settings.particles = 1000

flux_tally = openmc.Tally()
flux_tally.scores = ['flux']

# Create a Legendre polynomial expansion filter and add to tally
order = 16
expand_filter = openmc.SpatialLegendreFilter(order, 'z', zmin, zmax)
flux_tally.filters.append(expand_filter)

flux_mesh_tally = openmc.Tally()
flux_mesh_tally.scores = ['flux']
mesh = openmc.RegularMesh()
mesh.dimension = (1,1,40)
mesh.lower_left = geom.bounding_box[0]
mesh.upper_right = geom.bounding_box[1]

meshfilter = openmc.MeshFilter(mesh)

flux_mesh_tally.filters.append(meshfilter)

tallies = openmc.Tallies([flux_tally,flux_mesh_tally])
model = openmc.model.Model(geometry=geom, settings=settings, tallies=tallies)

sp_file = model.run(output=False)

#%%
with openmc.StatePoint(sp_file) as sp:
    df1 = sp.tallies[flux_tally.id].get_pandas_dataframe()
    df2 = sp.tallies[flux_mesh_tally.id].get_pandas_dataframe()
    
n = np.arange(order + 1)
a_n = (2*n + 1)/2 * df1['mean']


# divided by 10 for change of variables
phi10 = np.polynomial.Legendre(a_n/10, domain=(zmin, zmax))
# divided by 20 to match mesh tally
phi20 = np.polynomial.Legendre(a_n/20, domain=(zmin, zmax))


z = np.linspace(zmin, zmax, 1000)

fig,ax = plt.subplots()
ax.plot(z,phi10(z),label='expansion tally',alpha=0.75)
ax.plot(z,phi20(z),label='expansion tally halved',alpha=0.75)

z2 = np.linspace(zmin, zmax,41)
z2 = z2[:-1] + np.diff(z2)/2

ax.plot(z2,df2['mean'],label='mesh tally',alpha=0.75)
ax.set_xlabel('z(cm)')
ax.set_ylabel('flux')
plt.legend()