Fixed bug in statepoint instance attribute for openmc.mgxs.Library

openmc/mgxs/library.py

@@ -1,4 +1,5 @@
 import sys
+import os
 import copy
 from numbers import Integral
 
@@ -90,6 +91,7 @@ class Library(object):
             clone._domain_type = self.domain_type
             clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
             clone._all_mgxs = self.all_mgxs
+            clone._statepoint = self._statepoint
 
             clone._all_mgxs = {}
             for domain in self.domains:
@@ -396,7 +398,7 @@ class Library(object):
         HDF5 groups from the domain type, domain id, subdomain id (for
         distribcell domains), nuclides and cross section types. Two datasets for
         the mean and standard deviation are stored for each subdomain entry in
-        the HDF5 file.
+        the HDF5 file. The number of groups is stored as a file attribute.
 
         NOTE: This requires the h5py Python package.
 
@@ -435,6 +437,19 @@ class Library(object):
                   'library since a statepoint has not yet been loaded'
             raise ValueError(msg)
 
+        import h5py
+
+        # Make directory if it does not exist
+        if not os.path.exists(directory):
+            os.makedirs(directory)
+
+        # Add an attribute for the number of energy groups to the HDF5 file
+        full_filename = os.path.join(directory, filename + '.h5')
+        full_filename = full_filename.replace(' ', '-')
+        f = h5py.File(full_filename, 'w')
+        f.attrs["# groups"] = self.num_groups
+        f.close()
+
         # Export MGXS for each domain and mgxs type to an HDF5 file
         for domain in self.domains:
             for mgxs_type in self.mgxs_types:
@@ -443,4 +458,5 @@ class Library(object):
                 if subdomains == 'avg':
                     mgxs = mgxs.get_subdomain_avg_xs()
 
-                mgxs.build_hdf5_store(filename, directory, xs_type)
+                mgxs.build_hdf5_store(filename, directory,
+                                      xs_type=xs_type, nuclides=nuclides)

openmc/mgxs/mgxs.py

@@ -903,8 +903,9 @@ class MGXS(object):
 
         print(string)
 
-    def build_hdf5_store(self, filename='mgxs', directory='mgxs', append=True,
-                         subdomains='all', nuclides='all', xs_type='macro'):
+    def build_hdf5_store(self, filename='mgxs', directory='mgxs',
+                         subdomains='all', nuclides='all',
+                         xs_type='macro', append=True):
         """Export the multi-group cross section data to an HDF5 binary file.
 
         This method constructs an HDF5 file which stores the multi-group
@@ -921,9 +922,6 @@ class MGXS(object):
             Filename for the HDF5 file. Defaults to 'mgxs'.
         directory : str
             Directory for the HDF5 file. Defaults to 'mgxs'.
-        append : boolean
-            If true, appends to an existing HDF5 file with the same filename
-            directory (if one exists). Defaults to True.
         subdomains : Iterable of Integral or 'all'
             The subdomain IDs of the cross sections to include in the report.
             Defaults to 'all'.
@@ -936,6 +934,9 @@ class MGXS(object):
         xs_type: {'macro', 'micro'}
             Store the macro or micro cross section in units of cm^-1 or barns.
             Defaults to 'macro'.
+        append : boolean
+            If true, appends to an existing HDF5 file with the same filename
+            directory (if one exists). Defaults to True.
 
         Raises
         ------
@@ -952,12 +953,7 @@ class MGXS(object):
                   'cross section has not been computed'
             raise ValueError(msg)
 
-        # Attempt to import h5py
-        try:
-            import h5py
-        except ImportError:
-            msg = 'The h5py Python package must be installed on your system'
-            raise ImportError(msg)
+        import h5py
 
         # Make directory if it does not exist
         if not os.path.exists(directory):
@@ -993,27 +989,10 @@ class MGXS(object):
 
         cv.check_value('xs_type', xs_type, ['macro', 'micro'])
 
-        '''
-        if self.by_nuclide:
-            nuclides = self.domain.get_all_nuclides()
-            densities = np.zeros(len(nuclides), dtype=np.float)
-            for i, nuclide in enumerate(nuclides):
-                densities[i] = nuclides[nuclide][1]
-        else:
-            nuclides = ['sum']
-        '''
-
         # Create an HDF5 group within the file for the domain
         domain_type_group = xs_results.require_group(self.domain_type)
         domain_group = domain_type_group.require_group(str(self.domain.id))
 
-        '''
-        if subdomains == 'all' and self.domain_type == 'distribcell':
-            subdomains = np.arange(self.num_subdomains, dtype=np.int)
-        else:
-            subdomains = [self.domain.id]
-        '''
-
         # Determine number of digits to pad subdomain group keys
         num_digits = len(str(self.num_subdomains))
 
@@ -1976,7 +1955,7 @@ class Chi(MGXS):
         nu_fission_in = self.tallies['nu-fission-in']
         nu_fission_out = self.tallies['nu-fission-out']
 
-        # Remove the coarse energy filter to keep it out of tally arithmetic
+        # Remove coarse energy filter to keep it out of tally arithmetic
         energy_filter = nu_fission_in.find_filter('energy')
         nu_fission_in.remove_filter(energy_filter)
 
@@ -2071,9 +2050,17 @@ class Chi(MGXS):
                 nu_fission_in = nu_fission_in.summation(nuclides=nuclides)
                 nu_fission_out = nu_fission_out.summation(nuclides=nuclides)
 
+                # Remove coarse energy filter to keep it out of tally arithmetic
+                energy_filter = nu_fission_in.find_filter('energy')
+                nu_fission_in.remove_filter(energy_filter)
+
                 # Compute chi and store it as the xs_tally attribute so we can
                 # use the generic get_xs(...) method
                 xs_tally = nu_fission_out / nu_fission_in
+
+                # Add the coarse energy filter back to the nu-fission tally
+                nu_fission_in.add_filter(energy_filter)
+
                 xs = xs_tally.get_values(filters=filters,
                                          filter_bins=filter_bins, value=value)