OpenMC CSG geometry gotchas: universe fill truncation, overlapping cells, and cell vs mesh tallies

TL;DR

Use it to confirm a universe fill lands where you expect. Run the executable with the -g, -geometry-debug, or --geometry-debug flag. Use fewer particles, and check the output report for how many overlap checks were performed per cell: regions that particles rarely visit are poorly checked, so adjust the source or particle count to get coverage where it matters.

When to use

You are seeing this: A cell is a region of space filled with something. Use this skill when you run into "OpenMC CSG geometry gotchas: universe fill truncation, overlapping cells, and cell vs mesh tallies".

When not to use

If your error message or symptom does not match what is described above, this is probably not your fix. Search for your exact error text instead of forcing this one to fit.

Versions

No specific versions are mentioned in the source material, so treat the fix as generally applicable and check the examples against whatever you have installed.

Use this when building OpenMC models with the Python API: cells, universes, lattices, and setting up tallies (cell, mesh, distribcell filters; scores; estimators) for a criticality or fixed-source run. It covers the mistakes that produce silently wrong answers instead of errors.

The three building blocks: cell, universe, lattice

A cell is a region of space filled with something. The fill is set with the fill attribute and may be a material, a universe, or a lattice. No fill means a void cell: particles stream through it and undergo no collisions.

fuel = openmc.Cell(name='fuel', fill=uo2_material, region=-fuel_cylinder)
gap  = openmc.Cell(name='gap', region=+fuel_cylinder & -clad_cylinder)  # void

A universe is a collection of cells that forms a repeatable unit of geometry. Create one with openmc.Universe(cells=[...]) or build it up with universe.add_cell(cell) / universe.add_cells([...]). Universes are used three ways: as the root universe assigned to openmc.Geometry, as the fill of a cell, and as the contents of lattice elements. Cells that are never placed in any universe automatically belong to the base universe (ordinary Euclidean space), so universes are optional when nothing repeats.

A lattice (openmc.RectLattice, openmc.HexLattice) is a regular array of universes. A RectLattice needs four things:

• lower_left: coordinates of the lattice's lower-left corner
• pitch: distance between centers of adjacent elements (per axis)
• universes: nested lists of universes, one per lattice element
• outer: universe filling any position outside the defined lattice (optional)

lattice = openmc.RectLattice(name='2x2 pins')
lattice.pitch = (pitch, pitch)
lattice.lower_left = (-pitch, -pitch)
lattice.universes = [[pin_univ, pin_univ],
                     [pin_univ, pin_univ]]

Row ordering for RectLattice.universes: the first row corresponds to the lattice elements with the highest y value. A 3-D lattice needs x, y, z in lower_left and pitch and a triply-nested list for universes. HexLattice instead takes center, a radial (and optional axial) pitch list, and universes ordered as rings from outermost to innermost; use openmc.HexLattice.show_indices(n) to see the element layout before filling it in.

Gotcha 1: a filled universe is truncated by its parent cell

When a cell is filled with a universe, only the part of that universe lying inside the parent cell's region exists in the geometry. The universe's cells may extend to infinity (for example a moderator cell with region +clad_surface and no outer bound), they get cut off at the parent cell or lattice element boundary. This is how a pin universe becomes one assembly element: the same universe looks different depending on the cell it fills.

Two consequences:

1. cell.rotation and cell.translation may only be set on a cell filled with a universe, never on a material-filled cell. They rotate/translate the whole filled universe into the parent cell.

2. To inspect what is actually at a point, use Universe.find((x, y, z)) or Geometry.find((x, y, z)). The returned list is the traversal path (universes, cells, lattices); the last element is the lowest-level cell at that location. Use it to confirm a universe fill lands where you expect.

Gotcha 2: overlapping cells are NOT checked in normal runs

OpenMC does not check whether cells within a universe overlap during a normal simulation. Overlaps produce incorrect results that may or may not be obvious: a particle in an overlapped region is located in one of the cells (effectively an arbitrary choice), so material assignments and tallies in that region are wrong without any error being raised.

How to catch overlaps before trusting results:

1. Plot with overlap detection. The geometry plotters (Universe.plot, Geometry.plot, or the standalone plot run) accept show_overlaps=True, which flags overlapped pixels/voxels (colored red by default via overlap_color). Plot each universe, not just the root, at a resolution high enough that thin overlap slivers are not missed, the plotter only samples the center of each pixel.

2. Geometry debug mode. Run the executable with the -g, -geometry-debug, or --geometry-debug flag. This enables overlap checks at every particle move, with considerable runtime overhead. Use fewer particles, and check the output report for how many overlap checks were performed per cell: regions that particles rarely visit are poorly checked, so adjust the source or particle count to get coverage where it matters.

Overlaps almost always come from regions built with & and | that were meant to be disjoint but are not (a cylinder radius larger than the enclosing box, a ~ complement that swallows a neighbor, lattice elements whose pitch does not match the cell region that contains them). Fix the regions, not the symptom.

Gotcha 3: every model needs a root universe

openmc.Geometry needs a root universe. You can build one explicitly:

root_cell = openmc.Cell(name='root', fill=lattice, region=-bounds)
root_universe = openmc.Universe(universe_id=0, name='root')
root_universe.add_cell(root_cell)
geometry = openmc.Geometry(root_universe)

or pass a plain list of cells and let the constructor wrap them in a root universe automatically:

geometry = openmc.Geometry([root_cell])

The root cell's region should be bounded by surfaces with boundary_type='vacuum' (or 'reflective' / 'periodic' where appropriate); with no boundary condition at all, particles that should leak instead trigger "could not be located in any cell" errors or silently vanish from tallies. A bare reflective box models an infinite lattice without any universe hierarchy at all.

Tallies: filters choose where, scores choose what

A tally is a phase-space integral of a scoring function times the flux. Filters pick the phase-space regions; scores pick the response; Tally.nuclides restricts reaction rates to specific nuclides (default is the total over all nuclides).

# Cell tally: flux and fission neutron production in the fuel cell
cell_filter = openmc.CellFilter([fuel_cell])   # or a list of cell IDs
flux_tally = openmc.Tally(name='pin flux and fission')
flux_tally.filters = [cell_filter]
flux_tally.scores = ['flux', 'fission', 'nu-fission']

# Energy-binned fission rate in the fuel
energy_filter = openmc.EnergyFilter([0.0, 0.625, 20.0e6])  # eV, ascending edges
energy_tally = openmc.Tally(name='fuel fission spectrum')
energy_tally.filters = [cell_filter, energy_filter]
energy_tally.scores = ['fission']

Mesh tally: independent of the CSG cells. Define the mesh, wrap it in a MeshFilter, and read the flat tally.mean array back with correct axis ordering after the run (see below).

mesh = openmc.RegularMesh()
mesh.dimension = (50, 50, 1)
mesh.lower_left = [-pitch, -pitch, -5.0]
mesh.upper_right = [pitch, pitch, 5.0]
mesh_filter = openmc.MeshFilter(mesh)

mesh_tally = openmc.Tally(name='mesh flux')
mesh_tally.filters = [mesh_filter]
mesh_tally.scores = ['flux']

tallies = openmc.Tallies([flux_tally, energy_tally, mesh_tally])
tallies.export_to_xml()   # required before openmc.run(); openmc.model.Model.run() does this for you

Reading a mesh result back: the flat result array iterates x-fastest, z-slowest. Use Tally.get_reshaped_data(expand_dims=True) to get an array shaped (nx, ny, nz, nnuclides, nscores); do not just call mean.reshape(nz, ny, nx), the axes will be physically wrong.

with openmc.StatePoint('statepoint.100.h5') as sp:
    t = sp.get_tally(name='mesh flux')
    flux_3d = t.get_reshaped_data(expand_dims=True)  # shape (50, 50, 1, 1, 1)

CellFilter sums over all instances; DistribcellFilter separates them

A CellFilter on a cell that appears many times (e.g. the fuel cell inside a pin universe filling a lattice) scores the sum over every instance of that cell. To score each repeated instance separately, use openmc.DistribcellFilter with the cell itself (not a list):

distrib_filter = openmc.DistribcellFilter(fuel_cell)
inst_tally = openmc.Tally(name='per-pin fission')
inst_tally.filters = [distrib_filter]
inst_tally.scores = ['fission']

If you need reaction rates broken down by nuclide, set tally.nuclides = ['U235', 'U238'] (naming follows the material convention); 'all' gives a separate bin per nuclide in the model.

Estimator: tracklength by default, but not always

Tally.estimator accepts 'tracklength', 'collision', or 'analog'. The default is tracklength, which is generally the most efficient, and OpenMC reverts to analog automatically when the tally requires it. Tracklength and collision estimators cannot be used for anything needing post-collision information, for example, a scattering tally with outgoing-energy filters must run analog. If a tally comes back with suspiciously low scores or an estimator complaint, check whether your filter/score combination forces analog and whether that is acceptable for your statistics.

The 'current' score is special: with a MeshSurfaceFilter it scores partial currents on mesh element faces and may only be combined with energy and mesh filters and no other score; with a SurfaceFilter it scores net currents on defined geometry surfaces.

Complete example: 2x2 pin criticality with cell and mesh tallies

import openmc

# --- Materials ---
fuel = openmc.Material(name='UO2 fuel')
fuel.add_nuclide('U235', 0.04)
fuel.add_nuclide('U238', 0.96)
fuel.add_nuclide('O16', 2.00)
fuel.set_density('g/cc', 10.0)

water = openmc.Material(name='water')
water.add_nuclide('H1', 2.0)
water.add_nuclide('O16', 1.0)
water.set_density('g/cc', 1.0)
materials = openmc.Materials([fuel, water])

# --- Geometry: pin universe -> 2x2 lattice -> bounded root cell ---
pitch, fuel_r, clad_r = 1.26, 0.41, 0.47
fuel_surf = openmc.ZCylinder(r=fuel_r)
clad_surf = openmc.ZCylinder(r=clad_r)

fuel_cell = openmc.Cell(name='fuel', fill=fuel, region=-fuel_surf)
gap_cell  = openmc.Cell(name='gap', region=+fuel_surf & -clad_surf)
mod_cell  = openmc.Cell(name='moderator', fill=water, region=+clad_surf)
pin_univ  = openmc.Universe(name='pin', cells=[fuel_cell, gap_cell, mod_cell])

lattice = openmc.RectLattice(name='2x2 pins')
lattice.pitch = (pitch, pitch)
lattice.lower_left = (-pitch, -pitch)
lattice.universes = [[pin_univ, pin_univ],
                     [pin_univ, pin_univ]]

bounds = openmc.model.RectangularParallelepiped(
    -pitch, pitch, -pitch, pitch, -5.0, 5.0, boundary_type='vacuum')
root_cell = openmc.Cell(name='root', fill=lattice, region=-bounds)
geometry = openmc.Geometry([root_cell])

# Sanity: verify what cell sits at a few points before running
print(geometry.find((0.0, 0.0, 0.0))[-1].name)   # expect 'fuel'

# --- Settings: k-eigenvalue run ---
settings = openmc.Settings()
settings.run_mode = 'eigenvalue'
settings.particles = 1000
settings.batches = 50
settings.inactive = 10

# --- Tallies ---
cell_filter = openmc.CellFilter([fuel_cell])
pin_tally = openmc.Tally(name='fuel flux and fission')
pin_tally.filters = [cell_filter]
pin_tally.scores = ['flux', 'fission', 'nu-fission']

mesh = openmc.RegularMesh()
mesh.dimension = (50, 50, 1)
mesh.lower_left = [-pitch, -pitch, -5.0]
mesh.upper_right = [pitch, pitch, 5.0]
mesh_tally = openmc.Tally(name='mesh flux')
mesh_tally.filters = [openmc.MeshFilter(mesh)]
mesh_tally.scores = ['flux']

tallies = openmc.Tallies([pin_tally, mesh_tally])

# --- Export and run ---
model = openmc.model.Model(geometry, materials, settings, tallies)
sp_path = model.run()   # exports XML, runs, returns statepoint path

Pre-run checklist

1. Every universe's cells are disjoint: plot with show_overlaps=True or run with --geometry-debug before production statistics.

2. The root cell is bounded by surfaces that all carry a boundary condition ('vacuum', 'reflective', or paired 'periodic'); otherwise particles leak into undefined space.

3. Lattice pitch and lower_left actually tile the region of the cell that contains the lattice, and the lattice's universes shape matches the dimensionality.

4. CellFilter vs DistribcellFilter: decide whether repeated cells should be summed or reported per instance before interpreting the results.

5. Tallies needing post-collision data (outgoing-energy scattering spectra) silently revert to the analog estimator, budget particles accordingly.

6. Mesh tally results are x-fastest in the flat array: reshape with get_reshaped_data, never by hand.

Geometry roles grounded in docs.openmc.org v0.14.0 usersguide/geometry.html: Cell.fill accepts material/universe/lattice (void when empty); universes used as Geometry root, cell fill, or lattice contents; RectLattice requires lowerleft, pitch, universes, outer; universe-filled cells truncated by parent cell boundaries and rotation/translation only valid on universe-filled cells (v0.13.1 methods/geometry.html). Overlapping cells: v0.12.0 troubleshoot section 15.2.3.1 states normal runs do not check overlaps and results may be silently wrong; plotters check pixel/voxel centers (showoverlaps / overlapcolor per pythonapi/openmc.Cell); reliable validation is -g/--geometry-debug flag with per-cell check counts. Root universe requirement and boundary types (vacuum/reflective/periodic/white) from v0.14.0 geometry guide. Tallies: filters/scores/nuclides semantics, CellFilter and DistribcellFilter usage from v0.11.0 tallies guide and the pandas-dataframes and tally-arithmetic examples; MeshFilter with RegularMesh (dimension/lowerleft/upperright) and getreshapeddata axis ordering (x-fastest); estimator default tracklength reverting to analog when post-collision info is required from ioformats/tallies.rst; current-score constraints with MeshSurfaceFilter/SurfaceFilter from the v0.11.0 scores table.