diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 05cabf4961..97bb6d96b8 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -287,7 +287,7 @@ class Library(object): def by_nuclide(self, by_nuclide): cv.check_type('by_nuclide', by_nuclide, bool) - if by_nuclide == True and self.domain_type == 'mesh': + if by_nuclide and self.domain_type == 'mesh': raise ValueError('Unable to create MGXS library by nuclide with ' 'mesh domain') @@ -297,7 +297,7 @@ class Library(object): def domain_type(self, domain_type): cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES) - if self.by_nuclide == True and domain_type == 'mesh': + if self.by_nuclide and domain_type == 'mesh': raise ValueError('Unable to create MGXS library by nuclide with ' 'mesh domain') @@ -372,7 +372,15 @@ class Library(object): @scatter_format.setter def scatter_format(self, scatter_format): - cv.check_value('scatter_format', scatter_format, openmc.mgxs.MU_TREATMENTS) + cv.check_value('scatter_format', scatter_format, + openmc.mgxs.MU_TREATMENTS) + + if scatter_format == 'histogram' and self.correction == 'P0': + msg = 'The P0 correction will be ignored since the ' \ + 'scatter format is set to histogram' + warn(msg) + self.correction = None + self._scatter_format = scatter_format @legendre_order.setter @@ -405,6 +413,13 @@ class Library(object): msg = 'The histogram bins will be ignored since the ' \ 'scatter format is set to legendre' warn(msg) + elif self.scatter_format == 'histogram': + if self.correction == 'P0': + msg = 'The P0 correction will be ignored since ' \ + 'a histogram representation of the scattering '\ + 'kernel is requested' + warn(msg, RuntimeWarning) + self.correction = None self._histogram_bins = histogram_bins @@ -471,7 +486,7 @@ class Library(object): mgxs.delayed_groups = None else: delayed_groups \ - = list(range(1,self.num_delayed_groups+1)) + = list(range(1, self.num_delayed_groups + 1)) mgxs.delayed_groups = delayed_groups # If a tally trigger was specified, add it to the MGXS @@ -1351,7 +1366,7 @@ class Library(object): 'scattering matrix is not provided.') # Total or transport can be present, but if using # self.correction=="P0", then we should use transport. - if (((self.correction is "P0") and + if (((self.correction == "P0") and ('nu-transport' not in self.mgxs_types))): error_flag = True warn('A "nu-transport" MGXS type is required since a "P0" ' diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 324cbe56c4..036feac74a 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -3775,15 +3775,16 @@ class ScatterMatrixXS(MatrixMGXS): # tally data is stored in order of increasing energies if order_groups == 'increasing': xs = xs[:, ::-1, ::-1, ...] - - # Place the histogram bins before the incoming groups - if self.scatter_format == 'histogram': - xs = np.swapaxes(xs, 3, 2) - xs = np.swapaxes(xs, 2, 1) - + # import pdb; pdb.set_trace() if squeeze: - xs = np.squeeze(xs) - xs = np.atleast_2d(xs) + # We want to squeeze out everything but the in_groups, out_groups, + # and, if needed, num_mu_bins dimension. These must not be squeezed + # so 1-group problems have the correct shape. + if self.scatter_format == 'histogram': + axes = (5, 4, 0) + else: + axes = (4, 3, 0) + xs = np.squeeze(xs, axis=axes) return xs diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index e3af59632f..5256dd102c 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -15,7 +15,7 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \ # Supported incoming particle MGXS angular treatment representations _REPRESENTATIONS = ['isotropic', 'angle'] _SCATTER_TYPES = ['tabular', 'legendre', 'histogram'] -_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", +_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]", "[G'][DG]"] @@ -134,7 +134,7 @@ class XSdata(object): Note that some cross sections can be input in more than one shape so they are listed multiple times: - [Order][G][G']: scatter_matrix + [G][G'][Order]: scatter_matrix [G]: total, absorption, fission, kappa_fission, nu_fission, prompt_nu_fission, inverse_velocity @@ -298,9 +298,9 @@ class XSdata(object): self.num_delayed_groups) self._xs_shapes["[G'][DG]"] = (self.energy_groups.num_groups, self.num_delayed_groups) - self._xs_shapes["[Order][G][G']"] \ - = (self.num_orders, self.energy_groups.num_groups, - self.energy_groups.num_groups) + self._xs_shapes["[G][G'][Order]"] \ + = (self.energy_groups.num_groups, + self.energy_groups.num_groups, self.num_orders) # If representation is by angle prepend num polar and num azim if self.representation == 'angle': @@ -718,7 +718,7 @@ class XSdata(object): """ # Get the accepted shapes for this xs - shapes = [self.xs_shapes["[Order][G][G']"]] + shapes = [self.xs_shapes["[G][G'][Order]"]] # Convert to a numpy array so we can easily get the shape for checking scatter = np.asarray(scatter) @@ -1534,11 +1534,10 @@ class XSdata(object): if self.representation == 'isotropic': if self.scatter_format == 'legendre': # Get the scattering orders in the outermost dimension - self._scatter_matrix[i] = np.zeros((self.num_orders, - self.energy_groups.num_groups, - self.energy_groups.num_groups)) + self._scatter_matrix[i] = \ + np.zeros(self.xs_shapes["[G][G'][Order]"]) for moment in range(self.num_orders): - self._scatter_matrix[i][moment, :, :] = \ + self._scatter_matrix[i][:, :, moment] = \ scatter.get_xs(nuclides=nuclide, xs_type=xs_type, moment=moment, subdomains=subdomain) else: @@ -1657,7 +1656,7 @@ class XSdata(object): if self.num_polar is not None: grp.attrs['num_polar'] = self.num_polar - grp.attrs['scatter_shape'] = np.string_("[Order][G][G']") + grp.attrs['scatter_shape'] = np.string_("[G][G'][Order]") if self.scatter_format is not None: grp.attrs['scatter_format'] = np.string_(self.scatter_format) if self.order is not None: @@ -1736,14 +1735,13 @@ class XSdata(object): # Get the sparse scattering data to print to the library G = self.energy_groups.num_groups if self.representation == 'isotropic': - g_out_bounds = np.zeros((G, 2), dtype=np.int) for g_in in range(G): if self.scatter_format == 'legendre': - nz = np.nonzero(self._scatter_matrix[i][0, g_in, :]) + nz = np.nonzero(self._scatter_matrix[i][g_in, :, 0]) elif self.scatter_format == 'histogram': nz = np.nonzero(np.sum( - self._scatter_matrix[i][:, g_in, :], axis=0)) + self._scatter_matrix[i][g_in, :, :], axis=1)) if len(nz[0]) == 0: g_out_bounds[g_in, 0] = 0 g_out_bounds[g_in, 1] = 0 @@ -1757,8 +1755,8 @@ class XSdata(object): for g_in in range(G): for g_out in range(g_out_bounds[g_in, 0], g_out_bounds[g_in, 1] + 1): - for l in range(len(matrix[:, g_in, g_out])): - flat_scatt.append(matrix[l, g_in, g_out]) + for l in range(len(matrix[g_in, g_out, :])): + flat_scatt.append(matrix[g_in, g_out, l]) # And write it. scatt_grp = xs_grp.create_group('scatter_data') @@ -1767,7 +1765,6 @@ class XSdata(object): # Repeat for multiplicity if self._multiplicity_matrix[i] is not None: - # Now create the flattened scatter matrix array matrix = self._multiplicity_matrix[i][:, :] flat_mult = [] @@ -1793,7 +1790,13 @@ class XSdata(object): for p in range(Np): for a in range(Na): for g_in in range(G): - matrix = self._scatter_matrix[i][p, a, 0, g_in, :] + if self.scatter_format == 'legendre': + matrix = \ + self._scatter_matrix[i][p, a, g_in, :, 0] + elif self.scatter_format == 'histogram': + matrix = \ + np.sum(self._scatter_matrix[i][p, a, g_in, :, :], + axis=1) nz = np.nonzero(matrix) g_out_bounds[p, a, g_in, 0] = nz[0][0] g_out_bounds[p, a, g_in, 1] = nz[0][-1] @@ -1806,8 +1809,8 @@ class XSdata(object): for g_in in range(G): for g_out in range(g_out_bounds[p, a, g_in, 0], g_out_bounds[p, a, g_in, 1] + 1): - for l in range(len(matrix[:, g_in, g_out])): - flat_scatt.append(matrix[l, g_in, g_out]) + for l in range(len(matrix[g_in, g_out, :])): + flat_scatt.append(matrix[g_in, g_out, l]) # And write it. scatt_grp = xs_grp.create_group('scatter_data') diff --git a/scripts/openmc-update-mgxs b/scripts/openmc-update-mgxs index cad48d0e31..54a553f77a 100755 --- a/scripts/openmc-update-mgxs +++ b/scripts/openmc-update-mgxs @@ -23,12 +23,10 @@ optional arguments: from __future__ import print_function import os -from shutil import move import warnings import xml.etree.ElementTree as ET import argparse -import h5py import numpy as np import openmc.mgxs_library @@ -51,7 +49,7 @@ def parse_args(): if args['output'] == '': filename = args['input'].name - extension = filenameos.path.splitext() + extension = os.path.splitext(filename) if extension == '.xml': filename = filename[:filename.rfind('.')] + '.h5' args['output'] = filename @@ -84,6 +82,8 @@ if __name__ == '__main__': temp = tree.find('group_structure').text.strip() temp = np.array(temp.split()) group_structure = temp.astype(np.float) + # Convert from MeV to eV + group_structure *= 1.e6 energy_groups = openmc.mgxs.EnergyGroups(group_structure) temp = tree.find('inverse-velocity') if temp is not None: @@ -103,7 +103,7 @@ if __name__ == '__main__': temperature = get_data(xsdata_elem, 'kT') if temperature is not None: temperature = \ - float(temperature) / openmc.data.K_BOLTZMANN + float(temperature) / openmc.data.K_BOLTZMANN * 1.E6 else: temperature = 294. temperatures = [temperature] @@ -163,7 +163,7 @@ if __name__ == '__main__': if temp is not None: temp = np.array(temp.split(), dtype=float) total = temp.astype(np.float) - total.shape = xsd[i].vector_shape + total.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_total(total, temperature) if inverse_velocity is not None: @@ -172,41 +172,55 @@ if __name__ == '__main__': temp = get_data(xsdata_elem, 'absorption') temp = np.array(temp.split()) absorption = temp.astype(np.float) - absorption.shape = xsd[i].vector_shape + absorption.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_absorption(absorption, temperature) temp = get_data(xsdata_elem, 'scatter') temp = np.array(temp.split()) scatter = temp.astype(np.float) - scatter.shape = xsd[i].pn_matrix_shape + # This is now a flattened-array of something that started with a + # shape of [Order][G][G']; we need to unflatten and then switch the + # ordering + in_shape = (order_dim, energy_groups.num_groups, + energy_groups.num_groups) + if representation == 'angle': + in_shape = (n_pol, n_azi) + in_shape + scatter.shape = in_shape + scatter = np.swapaxes(scatter, 2, 3) + scatter = np.swapaxes(scatter, 3, 4) + else: + scatter.shape = in_shape + scatter = np.swapaxes(scatter, 0, 1) + scatter = np.swapaxes(scatter, 1, 2) + xsd[i].set_scatter_matrix(scatter, temperature) temp = get_data(xsdata_elem, 'multiplicity') if temp is not None: temp = np.array(temp.split()) multiplicity = temp.astype(np.float) - multiplicity.shape = xsd[i].matrix_shape + multiplicity.shape = xsd[i].xs_shapes["[G][G']"] xsd[i].set_multiplicity_matrix(multiplicity, temperature) temp = get_data(xsdata_elem, 'fission') if temp is not None: temp = np.array(temp.split()) fission = temp.astype(np.float) - fission.shape = xsd[i].vector_shape + fission.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_fission(fission, temperature) temp = get_data(xsdata_elem, 'kappa_fission') if temp is not None: temp = np.array(temp.split()) kappa_fission = temp.astype(np.float) - kappa_fission.shape = xsd[i].vector_shape + kappa_fission.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_kappa_fission(kappa_fission, temperature) temp = get_data(xsdata_elem, 'chi') if temp is not None: temp = np.array(temp.split()) chi = temp.astype(np.float) - chi.shape = xsd[i].vector_shape + chi.shape = xsd[i].xs_shapes['[G]'] xsd[i].set_chi(chi, temperature) else: chi = None @@ -216,9 +230,9 @@ if __name__ == '__main__': temp = np.array(temp.split()) nu_fission = temp.astype(np.float) if chi is not None: - nu_fission.shape = xsd[i].vector_shape + nu_fission.shape = xsd[i].xs_shapes['[G]'] else: - nu_fission.shape = xsd[i].matrix_shape + nu_fission.shape = xsd[i].xs_shapes["[G][G']"] xsd[i].set_nu_fission(nu_fission, temperature) # Build library as we go, but first we have enough to initialize it diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index 17104d5088..4fb1f72727 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -355,7 +355,7 @@ module mgxs_header if (attribute_exists(xs_id, "scatter_shape")) then call read_attribute(temp_str, xs_id, "scatter_shape") temp_str = trim(temp_str) - if (to_lower(temp_str) /= "[order][g][g']") then + if (to_lower(temp_str) /= "[g][g'][order]") then call fatal_error("Invalid scatter_shape option!") end if end if @@ -1794,8 +1794,8 @@ module mgxs_header order_dim = order + 1 end if - ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing - ! to ScattData + ! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for + ! passing to ScattData allocate(input_scatt(energy_groups, this % n_azi, this % n_pol)) index = 1 diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index 6004cbd190..684be88b83 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -274,7 +274,7 @@ contains ! Get this by summing the un-normalized P0 coefficient in matrix ! over all outgoing groups do gin = 1, groups - this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1) + this % scattxs(gin) = sum(matrix(gin) % data(:, :)) end do allocate(energy(groups)) @@ -282,7 +282,7 @@ contains ! while also normalizing matrix itself (making CDF of f(mu=1)=1) do gin = 1, groups allocate(energy(gin) % data(gmin(gin):gmax(gin))) - do gout = 1, groups + do gout = gmin(gin), gmax(gin) norm = sum(matrix(gin) % data(:, gout)) energy(gin) % data(gout) = norm if (norm /= ZERO) then