HexLattice cell rotation appears to break DistribcellFilter tally spatial ordering
Summary: Rotating the cell containing a HexLattice via cell.rotation as described in pull request #1232 produces correct geometry for transport, but maps the DistribcellFilter tally results to the un-rotated geometry. Has this solution been deprecated? Or is it still considered valid? I could not find this issue documented so far, and might be worth mentioning the limitations of the rotation attribute in the documentation.
Overview
When searching for how to rotate a hex lattice (Not the orientation, just the rotation) I came across a pull request #1232 that describes how to rotate a lattice by rotating it’s parent cell. Has this solution been deprecated? Or is it still considered valid?
# rotate the containing cell
lattice = openmc.HexLattice()
[...]
main_cell = openmc.Cell(fill=lattice)
main_cell.rotation = (0., 0., 60.)
If you then apply a DistribcellFilter to your fuel pins, it will return data for the un-rotated pin locations, not the current locations.
Why this happens (As far as I can tell)
DistribcellFilter assigns tally bins by cell instance number, which is determined by walking the lattice universe array in definition order (outermost ring to innermost, clockwise from “top”). This walk happens in the lattice’s own coordinate frame and is fixed at geometry build time.
cell.rotation is applied as a CSG transform after instance numbering is established. The neutron transport sees the correct rotated geometry, but instance 0 still maps to the first element in the Python ring array. After rotation, this now sits at a different physical location than expected.
The result: tally results appear to be correct in aggregate, but each bin is attributed to the wrong pin. The flux map appears rotated or scrambled relative to the true geometry.
Demonstration
import openmc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
from matplotlib.patches import RegularPolygon
# --- Materials ---
fuel = openmc.Material(name='fuel')
fuel.add_nuclide('U235', 1.0)
fuel.set_density('g/cm3', 10.0)
poison = openmc.Material(name='poison')
poison.add_nuclide('B10', 1.0)
poison.set_density('g/cm3', 2.52)
water = openmc.Material(name='water')
water.add_nuclide('H1', 2.0, 'ao')
water.add_nuclide('O16', 1.0, 'ao')
water.set_density('g/cm3', 1.0)
mats = openmc.Materials([fuel, poison, water])
mats.export_to_xml()
# --- Universes ---
r_pin = openmc.ZCylinder(r=0.4)
fuel_cell = openmc.Cell(name='fuel', fill=fuel, region=-r_pin)
poison_cell = openmc.Cell(name = 'poison', fill = poison, region = -r_pin)
fuel_water_cell = openmc.Cell(name='water in fuel pin', fill=water, region=+r_pin)
poison_water_cell = openmc.Cell(name='water in poison pin', fill=water, region=+r_pin)
fuel_pin_u = openmc.Universe(cells=[fuel_cell, fuel_water_cell])
poison_pin_u = openmc.Universe(cells=[poison_cell, poison_water_cell])
all_water_cell = openmc.Cell(name='outer water', fill=water)
outer_u = openmc.Universe(cells=[all_water_cell])
# --- Lattice ---
lat = openmc.HexLattice()
lat.center = (0., 0.)
lat.pitch = (1.25,)
lat.outer = outer_u
lat.universes = [[fuel_pin_u]*11 + [poison_pin_u], [fuel_pin_u]*6, [fuel_pin_u]]
# --- Geometry ---
boundary = openmc.ZCylinder(r=4.0, boundary_type='vacuum')
main_cell = openmc.Cell(fill=lat, region=-boundary)
# *** THE ISSUE: rotating the containing cell does not shift instance numbering ***
main_cell.rotation = (0., 0., 60.)
geom = openmc.Geometry([main_cell])
geom.export_to_xml()
# --- Settings ---
settings = openmc.Settings()
settings.batches = 50
settings.inactive = 10
settings.particles = 1000
settings.source = openmc.IndependentSource(
space=openmc.stats.Point((0, 0, 0)),
energy=openmc.stats.Discrete([2e6], [1.0])
)
settings.export_to_xml()
# --- Tally ---
# DistribcellFilter scores flux separately for each pin instance
fuel_filter = openmc.DistribcellFilter(fuel_cell)
poison_filter = openmc.DistribcellFilter(poison_cell)
tally_fuel = openmc.Tally(name='fuel_pin_flux')
tally_poison = openmc.Tally(name='poison_pin_flux')
tally_fuel.filters = [fuel_filter]
tally_poison.filters = [poison_filter]
tally_fuel.scores = ['flux']
tally_poison.scores = ['flux']
tallies = openmc.Tallies([tally_poison] + [tally_fuel])
tallies.export_to_xml()
# --- Combined Plot Figure ---
fig, ax = plt.subplots(1, 5, figsize=(24, 5))
# --- Plot 1: rotated containing cell ---
geom.plot(
basis='xy',
origin=(0.0, 0.0, 0.0),
width=(10.0, 10.0),
pixels=(1200, 1200),
color_by='material',
axes=ax[0]
)
ax[0].set_title("Cell Rotation Geometry")
# --- Run ---
openmc.run(threads=7)
# --- Extract rotated flux ---
with openmc.StatePoint('statepoint.50.h5') as sp:
t_fuel = sp.get_tally(name='fuel_pin_flux').get_pandas_dataframe()
t_poison = sp.get_tally(name='poison_pin_flux').get_pandas_dataframe()
def extract_pins(df):
df.columns = ['_'.join(filter(None, [str(c) for c in col])).strip('_')
for col in df.columns]
return pd.DataFrame({
'ix': df['level 2_lat_x'].astype(int),
'iy': df['level 2_lat_y'].astype(int),
'flux': df['mean'].astype(float),
'std_dev': df['std. dev.'].astype(float),
})
fuel_pins = extract_pins(t_fuel)
poison_pins = extract_pins(t_poison)
openmc_pins = pd.concat([fuel_pins, poison_pins], ignore_index=True)
openmc_pin = openmc_pins.groupby(['ix', 'iy'])[['flux']].sum().reset_index()
openmc_pin = openmc_pin[openmc_pin['flux'] > 0].reset_index(drop=True)
rotated_lattice_cell = openmc_pin
# --- INSTEAD ---
lat = openmc.HexLattice()
lat.center = (0., 0.)
lat.pitch = (1.25,)
lat.outer = outer_u
a_Ring = [fuel_pin_u]
b_Ring = [fuel_pin_u]*6
c_Ring = [fuel_pin_u]*11 + [poison_pin_u]
# Shift lattice 60 degrees to match rotation of the core
shiftDir = 1
shift = [2,1] * shiftDir
c_Ring = c_Ring[shift[0]:] + c_Ring[:shift[0]]
b_Ring = b_Ring[shift[1]:] + b_Ring[:shift[1]]
a_Ring = a_Ring
lat.universes = [c_Ring, b_Ring, a_Ring]
boundary = openmc.ZCylinder(r=4.0, boundary_type='vacuum')
main_cell = openmc.Cell(fill=lat, region=-boundary)
geom = openmc.Geometry([main_cell])
geom.export_to_xml()
# --- Tally ---
# DistribcellFilter scores flux separately for each pin instance
fuel_filter = openmc.DistribcellFilter(fuel_cell)
poison_filter = openmc.DistribcellFilter(poison_cell)
tally_fuel = openmc.Tally(name='fuel_pin_flux')
tally_poison = openmc.Tally(name='poison_pin_flux')
tally_fuel.filters = [fuel_filter]
tally_poison.filters = [poison_filter]
tally_fuel.scores = ['flux']
tally_poison.scores = ['flux']
tallies = openmc.Tallies([tally_poison] + [tally_fuel])
tallies.export_to_xml()
# --- Plot 2: shifted pin locations ---
geom.plot(
basis='xy',
origin=(0.0, 0.0, 0.0),
width=(10.0, 10.0),
pixels=(1200, 1200),
color_by='material',
axes=ax[1]
)
ax[1].set_title("Pin Location Rotation Geometry")
# --- Run ---
openmc.run(threads=7)
# --- Extract rotated flux ---
with openmc.StatePoint('statepoint.50.h5') as sp:
t_fuel = sp.get_tally(name='fuel_pin_flux').get_pandas_dataframe()
t_poison = sp.get_tally(name='poison_pin_flux').get_pandas_dataframe()
def extract_pins(df):
df.columns = ['_'.join(filter(None, [str(c) for c in col])).strip('_')
for col in df.columns]
return pd.DataFrame({
'ix': df['level 2_lat_x'].astype(int),
'iy': df['level 2_lat_y'].astype(int),
'flux': df['mean'].astype(float),
'std_dev': df['std. dev.'].astype(float),
})
fuel_pins = extract_pins(t_fuel)
poison_pins = extract_pins(t_poison)
openmc_pins = pd.concat([fuel_pins, poison_pins], ignore_index=True)
openmc_pin = openmc_pins.groupby(['ix', 'iy'])[['flux']].sum().reset_index()
openmc_pin = openmc_pin[openmc_pin['flux'] > 0].reset_index(drop=True)
rotated_pin_locations = openmc_pin
# --- Comparisons ---
a = rotated_lattice_cell
b = rotated_pin_locations
d = a.merge(b, on=["ix", "iy"], suffixes=("_lattice", "_pin"))
d["rel_diff"] = 100 * (d["flux_pin"] - d["flux_lattice"]) / d["flux_lattice"]
# Flat-top axial hex coordinates
d["x"] = 1.5 * d["ix"]
d["y"] = np.sqrt(3) * (d["iy"] + 0.5 * d["ix"])
def plot_hex(ax, data, values, title, cmap="viridis"):
hexes = [
RegularPolygon(
(row.x, row.y),
numVertices=6,
radius=0.55,
orientation=0,
)
for row in data.itertuples()
]
pc = PatchCollection(hexes, cmap=cmap, edgecolor="black", linewidth=0.8)
pc.set_array(data[values].to_numpy())
ax.add_collection(pc)
ax.autoscale()
ax.set_aspect("equal")
ax.set_title(title)
ax.axis("off")
plt.colorbar(pc, ax=ax, shrink=0.75)
plot_hex(ax[2], d, "flux_lattice", "Lattice Cell Rotation Flux")
plot_hex(ax[3], d, "flux_pin", "Pin Location Rotation Flux")
plot_hex(ax[4], d, "rel_diff", "Relative Difference (%)", cmap="coolwarm")
plt.tight_layout()
plt.savefig("Geometry_Comparison_All_5.png", dpi=300, bbox_inches="tight")
plt.show()
Workaround
Instead of rotating the lattice containing cell, I just shifted the ring arrays themselves before assigning them to lattice.universes. For a 60° rotation on a 6 axis symmetric lattice, shift each ring by n/6 positions (where n is the ring size):
# Instead of main_cell.rotation = (0., 0., 60.)
# Shift rings so instance ordering matches physical positions
outer_ring = outer_ring[2:] + outer_ring[:2] # 12-element ring: shift by 2
mid_ring = mid_ring[1:] + mid_ring[:1] # 6-element ring: shift by 1
# inner ring (1 element) unchanged
lat.universes = [outer_ring, mid_ring, inner_ring]
# No rotation on the containing cell
This bakes the rotation into the array definition order, so instance numbers and physical positions remain consistent with the tallies. Is this a better solution to lattice rotation? Should this be mentioned in the documentation?
Environment
- OpenMC version: 0.15.3
- Application: TRIGA reactor full-core model, hex lattice with 7 rings
- Tally type:
DistribcellFilterpin flux map