# OpenMC CSG geometry gotchas: universe fill truncation, overlapping cells, and cell vs mesh tallies Build correct OpenMC CSG geometry and tallies: how fill= truncates universes, why overlapping cells silently corrupt results and how to catch them, and how to set up cell, mesh, and distribcell tallies with scores, filters, and estimators for criticality or fixed-source runs. Exact reference: {"kind":"skill_version","skill_id":"skl_lDXMt_ccTcsSYYl0AF_NBA","version_id":"skv_SRIZ5x-ZTbXUYZiOzl4CUg"} Applicability: [{"constraint":">=0.11.0","technology":"openmc","version_scheme":"semver"}] # 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. ## Supporting basis and limitations 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 lower_left, 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 (show_overlaps / overlap_color 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/lower_left/upper_right) and get_reshaped_data axis ordering (x-fastest); estimator default tracklength reverting to analog when post-collision info is required from io_formats/tallies.rst; current-score constraints with MeshSurfaceFilter/SurfaceFilter from the v0.11.0 scores table. ## Change and rationale New skill: practical OpenMC CSG geometry and tally guide. Covers cell/universe/lattice roles, universe-fill truncation by the parent cell, the silently-unchecked overlapping-cells trap with plot and --geometry-debug validation, root-universe and boundary-condition requirements, and tally setup: CellFilter vs MeshFilter vs DistribcellFilter, scores, nuclides, estimator defaults, and reading mesh results back with get_reshaped_data. Includes a complete runnable 2x2 pin criticality model. OpenMC geometry errors (overlapping cells, misunderstood universe fill truncation, missing root universe) fail silently rather than raising errors, and tally misuse (CellFilter summing over lattice instances, hand-reshaped mesh data, wrong estimator expectations) produces misleading numbers. A single skill that pairs each gotcha with the concrete validation step and real API code prevents the most common silent-failure modes in one place.