Allow DistribcellFilter to work with apply_tally_results=True (#3667)

openmc/__init__.py

@@ -37,7 +37,7 @@ from openmc.tracks import *
 from .config import *
 
 # Import a few names from the model module
-from openmc.model import Model
+from openmc.model import Model, SearchResult
 
 from . import examples
 

openmc/data/decay.py

@@ -13,7 +13,7 @@ import openmc.checkvalue as cv
 from openmc.exceptions import DataError
 from openmc.mixin import EqualityMixin
 from openmc.stats import Discrete, Tabular, Univariate, combine_distributions
-from .data import ATOMIC_SYMBOL, ATOMIC_NUMBER
+from .data import ATOMIC_NUMBER, gnds_name
 from .function import INTERPOLATION_SCHEME
 from .endf import Evaluation, get_head_record, get_list_record, get_tab1_record
 
@@ -126,9 +126,7 @@ class FissionProductYields(EqualityMixin):
                 for j in range(n_products):
                     Z, A = divmod(int(values[4*j]), 1000)
                     isomeric_state = int(values[4*j + 1])
-                    name = ATOMIC_SYMBOL[Z] + str(A)
-                    if isomeric_state > 0:
-                        name += f'_m{isomeric_state}'
+                    name = gnds_name(Z, A, isomeric_state)
                     yield_j = ufloat(values[4*j + 2], values[4*j + 3])
                     yields[name] = yield_j
 
@@ -256,10 +254,7 @@ class DecayMode(EqualityMixin):
                         A += delta_A
                         Z += delta_Z
 
-        if self._daughter_state > 0:
-            return f'{ATOMIC_SYMBOL[Z]}{A}_m{self._daughter_state}'
-        else:
-            return f'{ATOMIC_SYMBOL[Z]}{A}'
+        return gnds_name(Z, A, self._daughter_state)
 
     @property
     def parent(self):
@@ -348,10 +343,7 @@ class Decay(EqualityMixin):
         self.nuclide['atomic_number'] = Z
         self.nuclide['mass_number'] = A
         self.nuclide['isomeric_state'] = metastable
-        if metastable > 0:
-            self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}_m{metastable}'
-        else:
-            self.nuclide['name'] = f'{ATOMIC_SYMBOL[Z]}{A}'
+        self.nuclide['name'] = gnds_name(Z, A, metastable)
         self.nuclide['mass'] = items[1]  # AWR
         self.nuclide['excited_state'] = items[2]  # State of the original nuclide
         self.nuclide['stable'] = (items[4] == 1)  # Nucleus stability flag

openmc/data/neutron.py

@@ -11,7 +11,7 @@ import h5py
 
 from . import HDF5_VERSION, HDF5_VERSION_MAJOR
 from .ace import Library, Table, get_table, get_metadata
-from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV
+from .data import ATOMIC_SYMBOL, K_BOLTZMANN, EV_PER_MEV, gnds_name
 from .endf import (
     Evaluation, SUM_RULES, get_head_record, get_tab1_record, get_evaluations)
 from .fission_energy import FissionEnergyRelease
@@ -678,11 +678,7 @@ class IncidentNeutron(EqualityMixin):
         temperature = ev.target['temperature']
 
         # Determine name
-        element = ATOMIC_SYMBOL[atomic_number]
-        if metastable > 0:
-            name = f'{element}{mass_number}_m{metastable}'
-        else:
-            name = f'{element}{mass_number}'
+        name = gnds_name(atomic_number, mass_number, metastable)
 
         # Instantiate incident neutron data
         data = cls(name, atomic_number, mass_number, metastable,

openmc/filter.py

@@ -1839,12 +1839,21 @@ class DistribcellFilter(Filter):
 
     @property
     def paths(self):
-        return self._paths
+        if self._paths is None:
+            if not hasattr(self, '_geometry'):
+                raise ValueError(
+                    "Model must be exported before the 'paths' attribute is" \
+                    "available for a DistribcellFilter.")
 
-    @paths.setter
-    def paths(self, paths):
-        cv.check_iterable_type('paths', paths, str)
-        self._paths = paths
+            # Determine paths for cell instances
+            self._geometry.determine_paths()
+
+            # Get paths for the corresponding cell
+            cell_id = self.bins[0]
+            cell = self._geometry.get_all_cells()[cell_id]
+            self._paths = cell.paths
+
+        return self._paths
 
     @Filter.bins.setter
     def bins(self, bins):

openmc/model/model.py

@@ -546,6 +546,13 @@ class Model:
                 depletion_operator.cleanup_when_done = True
                 depletion_operator.finalize()
 
+    def _link_geometry_to_filters(self):
+        """Establishes a link between distribcell filters and the geometry"""
+        for tally in self.tallies:
+            for f in tally.filters:
+                if isinstance(f, openmc.DistribcellFilter):
+                    f._geometry = self.geometry
+
     def export_to_xml(self, directory: PathLike = '.', remove_surfs: bool = False,
                       nuclides_to_ignore: Iterable[str] | None = None):
         """Export model to separate XML files.
@@ -587,6 +594,8 @@ class Model:
         if self.plots:
             self.plots.export_to_xml(d)
 
+        self._link_geometry_to_filters()
+
     def export_to_model_xml(self, path: PathLike = 'model.xml', remove_surfs: bool = False,
                             nuclides_to_ignore: Iterable[str] | None = None):
         """Export model to a single XML file.
@@ -666,6 +675,8 @@ class Model:
                 fh.write(ET.tostring(plots_element, encoding="unicode"))
             fh.write("</model>\n")
 
+        self._link_geometry_to_filters()
+
     def import_properties(self, filename: PathLike):
         """Import physical properties
 

openmc/statepoint.py

@@ -723,7 +723,7 @@ class StatePoint:
                     cell = cells[cell_id]
                     if not cell._paths:
                         summary.geometry.determine_paths()
-                    tally_filter.paths = cell.paths
+                    tally_filter._paths = cell.paths
 
         self._summary = summary
 

openmc/tallies.py

@@ -290,7 +290,7 @@ class Tally(IDManagerMixin):
 
     @property
     def num_nuclides(self):
-        return len(self._nuclides)
+        return max(len(self._nuclides), 1)
 
     @property
     def scores(self):
@@ -393,6 +393,11 @@ class Tally(IDManagerMixin):
             group = f[f'tallies/tally {self.id}']
             self._num_realizations = int(group['n_realizations'][()])
 
+            for filt in self.filters:
+                if isinstance(filt, openmc.DistribcellFilter):
+                    filter_group = f[f'tallies/filters/filter {filt.id}']
+                    filt._num_bins = int(filter_group['n_bins'][()])
+
             # Update nuclides
             nuclide_names = group['nuclides'][()]
             self._nuclides = [name.decode().strip() for name in nuclide_names]
@@ -3704,43 +3709,17 @@ class Tallies(cv.CheckedList):
             if possible. Defaults to False.
 
         """
-        if not isinstance(tally, Tally):
-            msg = f'Unable to add a non-Tally "{tally}" to the Tallies instance'
-            raise TypeError(msg)
-
         if merge:
-            merged = False
-
             # Look for a tally to merge with this one
             for i, tally2 in enumerate(self):
-
                 # If a mergeable tally is found
                 if tally2.can_merge(tally):
                     # Replace tally2 with the merged tally
                     merged_tally = tally2.merge(tally)
                     self[i] = merged_tally
-                    merged = True
-                    break
+                    return
 
-            # If no mergeable tally was found, simply add this tally
-            if not merged:
-                super().append(tally)
-
-        else:
-            super().append(tally)
-
-    def insert(self, index, item):
-        """Insert tally before index
-
-        Parameters
-        ----------
-        index : int
-            Index in list
-        item : openmc.Tally
-            Tally to insert
-
-        """
-        super().insert(index, item)
+        super().append(tally)
 
     def merge_tallies(self):
         """Merge any mergeable tallies together. Note that n-way merges are

tests/unit_tests/test_filter_distribcell.py

@@ -0,0 +1,53 @@
+import openmc
+import pandas as pd
+
+
+def test_distribcell_filter_apply_tally_results(run_in_tmpdir):
+    # Reset IDs to ensure consistent paths
+    openmc.reset_auto_ids()
+
+    mat = openmc.Material()
+    mat.add_nuclide("U235", 1.0)
+    mat.set_density("g/cm3", 1.0)
+
+    # Define 2x2 lattice with a cylinder in each universe
+    cyl = openmc.ZCylinder(r=1.0)
+    cell1 = openmc.Cell(fill=mat, region=-cyl)
+    cell2 = openmc.Cell(fill=None, region=+cyl)
+    univ = openmc.Universe(cells=[cell1, cell2])
+    lattice = openmc.RectLattice()
+    lattice.lower_left = (-3.0, -3.0)
+    lattice.pitch = (3.0, 3.0)
+    lattice.universes = [[univ, univ], [univ, univ]]
+    box = openmc.model.RectangularPrism(6., 6., boundary_type='reflective')
+    root_cell = openmc.Cell(region=-box, fill=lattice)
+    geometry = openmc.Geometry([root_cell])
+
+    # Create model and add tally with distribcell filter
+    model = openmc.Model(geometry)
+    model.settings.batches = 10
+    model.settings.particles = 1000
+    tally = openmc.Tally()
+    distribcell_filter = openmc.DistribcellFilter(cell1)
+    tally.filters = [distribcell_filter]
+    tally.scores = ['flux']
+    model.tallies = [tally]
+
+    # Run OpenMC and apply tally results
+    model.run(apply_tally_results=True)
+
+    # Check that mean and standard deviation are available on tally
+    assert tally.mean.shape == (4, 1, 1)
+    assert tally.std_dev.shape == (4, 1, 1)
+
+    # Make sure paths attribute on filter is correct
+    assert distribcell_filter.paths == [
+        'u3->c3->l2(0,0)->u1->c1',
+        'u3->c3->l2(1,0)->u1->c1',
+        'u3->c3->l2(0,1)->u1->c1',
+        'u3->c3->l2(1,1)->u1->c1',
+    ]
+
+    # Check that we can get a DataFrame from the tally
+    df = tally.get_pandas_dataframe()
+    assert isinstance(df, pd.DataFrame)