# OpenMC CSG Geometry Gotchas: Universe Fill Truncation, Overlapping Cells, and Cell vs Mesh Tallies
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.
```python
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)
```python
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:
```python
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:
```python
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).
```python
# 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).
```python
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, n_nuclides, n_scores); do not just call `mean.reshape(nz, ny, nx)`,
the axes will be physically wrong.
```python
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):
```python
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
```python
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.