Segfault when trying to write tallies.out with DistribcellFilter

[lewisgross1296]

Bug Description

Originally, I created a forum post about a segmentation fault that occurred for my model. Included is the ordinary error and a more detailed stack trace using the debugger. That model is a bit complicated, so I created a simpler example that reproduces the error.

The segfault occurs when when I use a DistribcellFilter on a cell representing the inside of an infinite ZCylinder. This cell appears in the model as part of a universe which is repeated in a HexLattice and is eventually intersected with openmc.ZPlanes to make a geometry with finite cells.

I'm pretty sure the error happens during the process of writing tallies.out, since the error vanishes when I tell my settings not to produce the tallies.out file.

Steps to Reproduce

Running this script will produce the error

import openmc
import openmc.model
import numpy as np

radius = 0.9 # cm
pin_lattice_pitch = 2.0 # cm

# zplanes to define lower and upper region
z_low = openmc.ZPlane(z0=0,boundary_type='vacuum') # z=0 cm
z_mid = openmc.ZPlane(z0=10) # z=10 cm
z_high = openmc.ZPlane(z0=20,boundary_type='vacuum') # z=20 cm
hex_boundary_region = openmc.model.hexagonal_prism(edge_length=pin_lattice_pitch,orientation='y', boundary_type='periodic')

# materials
nat_U = openmc.Material(name="nat_U")
nat_U.set_density('g/cm3',1.0)
nat_U.add_element('U',1.0)

graphite = openmc.Material(name="Graphite")
graphite.set_density('g/cm3', 1.1995)
graphite.add_element('C', 1.0)

mats = openmc.Materials([nat_U,graphite])
mats.export_to_xml()

# cylinder
cyl = openmc.ZCylinder(r=radius)
cell = openmc.Cell(region=-cyl, fill=nat_U)
outer_cell = openmc.Cell(region=+cyl, fill=graphite)
univ = openmc.Universe(cells=[cell,outer_cell])

# make outer graphite universe to give to lattice
all_graphite_cell = openmc.Cell(fill=graphite)
outer_graphite_universe = openmc.Universe(cells=(all_graphite_cell,))

# create a hexagonal lattice of compacts
active_core_lattice = openmc.HexLattice()
active_core_lattice.pitch = (pin_lattice_pitch,)
active_core_lattice.center = (0., 0.)
active_core_lattice.outer = outer_graphite_universe

center = [univ]
ring = [univ, univ, univ, univ, univ, univ]

active_core_lattice.universes = [ring, center]
lower_hex_boundary_cell = openmc.Cell(fill=active_core_lattice, region=hex_boundary_region & +z_low & -z_mid)
upper_hex_boundary_cell = openmc.Cell(fill=active_core_lattice, region=hex_boundary_region & +z_mid & -z_high)
hex_cells = [lower_hex_boundary_cell,upper_hex_boundary_cell]

geometry = openmc.Geometry(hex_cells)
geometry.export_to_xml()

tallies = openmc.Tallies()

# moderator
tally = openmc.Tally(name="tally")
tally.scores = ['flux']
filter = openmc.DistribcellFilter(cell)
tally.filters = [filter]
tallies.append(tally)
tallies.export_to_xml()

# settings
settings = openmc.Settings()
# source definition. fission source given bounding box of graphite active region
system_LL = (-pin_lattice_pitch*np.sqrt(3)/2, -pin_lattice_pitch, 0)
system_UR = (pin_lattice_pitch*np.sqrt(3)/2, pin_lattice_pitch, 20)
source_dist = openmc.stats.Box(system_LL, system_UR, only_fissionable=True)
source = openmc.Source(space=source_dist)
settings.source = source
# k-eigenvalue settings
settings.particles = 10000
settings.inactive = 25
settings.batches = 100
settings.export_to_xml()

openmc.run()

Environment

I am using OpenMC version 0.13.3 on Ubuntu 20.04.6 LTS. I installed from source in debug mode. The simulation is run using python 3.10.5.

[lewisgross1296]

Tagging @pshriwise

[pshriwise]

Thanks for the MWE @lewisgross1296! This issue looked really familiar to me and I thought we had fixed it (see #1995 and the related test), which is I was a little confused when you presented this in the forum. It appears that fix only works for rectangular lattices though. I'll think this a bit and get back to you.

[lewisgross1296]

Considering that my example could be more minimal, I've made the body of a test similar to the test you linked that causes the same error. Maybe we can chat at SW about how to turn this into a full test that goes into OpenMC. I'm not sure I fully understand exactly how to make an analogous expected_result or the cdtemp chunk, but here is a more minimal case that causes the same error:

import numpy as np
import openmc

model = openmc.Model()

# Create a single material
m = openmc.Material()
m.add_nuclide('U235', 1.0)
m.set_density('g/cm3', 10.0)
model.materials.append(m)

# Create a universe with a single infinite cell
c = openmc.Cell(fill=m)
u = openmc.Universe(cells=[c])

# create a hexagonal lattice
pin_lattice_pitch = 2.0 # cm
hex_lat = openmc.HexLattice()
hex_lat.pitch = (pin_lattice_pitch,)
hex_lat.center = (0., 0.)
hex_lat.outer = u

# define center and ring
center = [u]
ring = [u, u, u, u, u, u]
hex_lat.universes = [ring, center]

# zplanes to define lower and upper region
z_low = openmc.ZPlane(z0=0,boundary_type='vacuum') # z=0 cm
z_mid = openmc.ZPlane(z0=10) # z=10 cm
z_high = openmc.ZPlane(z0=20,boundary_type='vacuum') # z=20 cm

# hexagonal region to house hex lattice cells
hex_boundary_region = openmc.model.hexagonal_prism(edge_length=pin_lattice_pitch,orientation='y', boundary_type='periodic')
lower_hex_boundary_cell = openmc.Cell(fill=hex_lat, region=hex_boundary_region & +z_low & -z_mid)
upper_hex_boundary_cell = openmc.Cell(fill=hex_lat, region=hex_boundary_region & +z_mid & -z_high)
model.geometry = openmc.Geometry([lower_hex_boundary_cell,upper_hex_boundary_cell])

# Add necessary settings and export
model.settings.batches = 10
model.settings.inactive = 0
model.settings.particles = 100
# settings
# source definition. fission source given bounding box of graphite active region
system_LL = (-pin_lattice_pitch*np.sqrt(3)/2, -pin_lattice_pitch, 0)
system_UR = (pin_lattice_pitch*np.sqrt(3)/2, pin_lattice_pitch, 20)
source_dist = openmc.stats.Box(system_LL, system_UR, only_fissionable=True)
source = openmc.Source(space=source_dist)
model.settings.source = source
model.export_to_xml()

openmc.run()

[pshriwise]

Welp, this let me down the wonderful rabbit hole of cell instances and lattice offsets, uncovering a few different issues along the way. Thanks for your patience on this issue @lewisgross1296.

The failure point seen here is a result of getting through the following loop in distribcell_path_inner, which is called to determine the text labels of distributed cell instances written to the tallies.out file.

openmc/src/geometry_aux.cpp

Lines 526 to 545 in fb65dfd

	for (; cell_it != search_univ.cells_.crend(); ++cell_it) {
	Cell& c = *model::cells[*cell_it];
	// Material cells don't contain other cells so ignore them.
	if (c.type_ != Fill::MATERIAL) {
	int32_t temp_offset;
	if (c.type_ == Fill::UNIVERSE) {
	temp_offset = offset + c.offset_[map];
	} else {
	Lattice& lat = *model::lattices[c.fill_];
	int32_t indx = lat.universes_.size() * map + lat.begin().indx_;
	temp_offset = offset + lat.offsets_[indx];
	}
	// The desired cell is the first cell that gives an offset smaller or
	// equal to the target offset.
	if (temp_offset <= target_offset)
	break;
	}
	}

As a start, I've created a PR to report a failure if we make it through that loop (#2812).

In C++, our lattice classes hold an "offset table" that represents the instance offset for a cell used within one of the lattice universes.

It turns out that, for hexagonal lattices, this happens because the entry in the offset table being checked in this line is -1 and will always be less than the offset we're looking for when determining the instance of the path as we iterate backward through the offset table, which should result in monotonically decreasing values. However, the lattice offset table is allocated to be "square" as it reflects the allocation of the HexagonalLattice::universes_ member. Specifically, universes_ is allocated to be (2*n_rings -1)**2 * n_axial_. This results in n_rings_ - 1 empty entries at the beginning and end of the universes_ vector, making reliance on the same iterator methods (i.e. begin, end, rbegin, and rend) for both RectangularLattice and HexagonalLattice.

It turns out that the incorrect loop-exit happens before we get to the call of distribcell_path_inner that results in a segfault. This lowest call in the stack is made with an incorrect search cell, at which point the loop goes over a set of material-filled cells, triggers neither condition within the loop, and attempts to access an entry outside of the model::cells vector.

Additionally, a bug was found when more than one cell was part of the universe filling the lattice. We weren't accounting for the cell offset in the offset comparison in that loop either to make sure the correct section of the offset table is accessed when determining paths. This also showed up in how tallies are handled and will require another fix -- that one will go in soon as well.

TL;DR

There's a problem with how we iterate over hexagonal cells and index into offsets that only manifests itself when a lattice is used to fill more than one cell. I'll be submitting a fix for those shortly and this should fix the issue you're seeing with writing the tallies.out file.

PR Plan:

• Add an error message for the failure mode in distribcell_path_inner
  (Provide error message if a cell path can't be determined #2812)
• Fix distribcell tallies for models with these features
  (Fix distribcell labels for lattices used as fill in multiple cells #2813)
• Fix lattice offset iterators for hexagonal lattices
  (on deck, branch)

[lewisgross1296]

Thanks Patrick! Glad you got to the bottom of this. Tagging @gonuke for the update