diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index b740b4dae9..6f46f48875 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -1793,7 +1793,7 @@ class XSdata(object): """ from scipy.interpolate import interp1d - from scipy.integrate import quad, simps + from scipy.integrate import simps from scipy.special import eval_legendre check_value('target_format', target_format, _SCATTER_TYPES) @@ -1827,24 +1827,34 @@ class XSdata(object): new_data[..., :order] = orig_data[..., :order] elif target_format == 'tabular': mu = np.linspace(-1, 1, xsdata.num_orders) - # Create the Legendre series and either do a point-wise - # evaluation (tabular) or bin-wise integral (histogram) - new_data[..., :] = \ - np.sum((l + 0.5) * eval_legendre(l, mu) * orig_data - for l in range(self.num_orders)) + # Evaluate the legendre on the tabular grid within mu + for imu in range(len(mu)): + new_data[..., imu] = \ + np.sum((l + 0.5) * eval_legendre(l, mu[imu]) * + orig_data[..., l] + for l in range(self.num_orders)) elif target_format == 'histogram': mu = np.linspace(-1, 1, xsdata.num_orders + 1) # This code will be written to utilize the vectorized # integration capabilities instead of having an isotropic # and angle representation path. for h_bin in range(xsdata.num_orders): - mu_fine = np.linspace(mu[i], mu[i + 1], _NMU) - new_data[..., h_bin] = \ - simps(np.sum((l + 0.5) * - eval_legendre(l, mu_fine) * - orig_data - for l in range(self.num_orders)), - mu_fine) + mu_fine = np.linspace(mu[h_bin], mu[h_bin + 1], _NMU) + table_shape = new_data.shape[:-1] + (_NMU,) + table_fine = np.zeros(table_shape) + for imu in range(len(mu_fine)): + table_fine[..., imu] = \ + np.sum((l + 0.5) * + eval_legendre(l, mu_fine[imu]) * + orig_data[..., l] + for l in range(self.num_orders)) + new_data[..., h_bin] = simps(table_fine, mu_fine) + # new_data[..., h_bin] = \ + # simps(np.sum((l + 0.5) * + # eval_legendre(l, mu_fine) * + # orig_data[..., l] + # for l in range(self.num_orders)), + # mu_fine) # Remove the very small results from numerical precision # issues (allowing conversions to be reproduced exactly) @@ -1899,7 +1909,7 @@ class XSdata(object): tab_data[..., :-1] = orig_data tab_data[..., -1] = (orig_data[..., -2] - orig_data[..., -3]) + orig_data[..., -2] - # Now get the distribution normalization factor.such + # Now get the distribution normalization factor. # This factor will be such that its average value is the # group -> group cross section (which is the sum of the # original data) @@ -1951,6 +1961,198 @@ class XSdata(object): return xsdata + def convert_to_continuous_energy(self, temperature=294.): + """Converts the XSdata object to an equivalent + openmc.data.IncidentNeutron object. + + Parameters + ---------- + temperature : float + Temperature of dataset to print; defaults to 294K + + Returns + ------- + openmc.data.IncidentNeutron + The continuous-energy IncidentNeutron data library equivalent + to the MGXS data in self + """ + + # Check if this can be performed successfully + if self.representation == 'angle': + raise ValueError("Cannot convert angle-dependent MGXS; convert to " + "an isotropic representation first") + required_types = ['absorption', 'scatter_matrix'] + for type_check in required_types: + if getattr(self, '_' + type_check)[0] is None: + raise ValueError(type_check + ' data is required') + if temperature not in self.temperatures: + raise ValueError("Invalid temperature") + + def convert_xs(group_edges, values): + cexs = openmc.data.Tabulated1D(group_edges, + np.append(values[::-1], + [values[0]]), + breakpoints=[len(group_edges)], + interpolation=[1]) + return cexs + + # Build required metadata + kTs = self.temperatures[:] * openmc.data.K_BOLTZMANN + if self.atomic_weight_ratio: + awr = self.atomic_weight_ratio + else: + awr = 1. + + # Get the temperature index + iT = np.where(self.temperatures == temperature)[0][0] + data = openmc.data.IncidentNeutron(self.name, 1, 1, 0, awr, kTs) + strT = "{}K".format(int(round(self.temperatures[iT]))) + data.energy = {strT: self.energy_groups.group_edges} + energy_midpoints = (self.energy_groups.group_edges[1:] + + self.energy_groups.group_edges[:-1]) / 2. + + # Elastic MGXS: must explicitly be 0 barns to avoid incorrect + # elastic_scatter calculation with data available to us in MG library. + el_rxn = openmc.data.Reaction(2) + el_rxn.xs[strT] = \ + convert_xs(self.energy_groups.group_edges, + np.zeros(self.energy_groups.num_groups)) + data.reactions[2] = el_rxn + + # Absorption MGXS + abs_rxn = openmc.data.Reaction(102) + abs_rxn.xs[strT] = convert_xs(self.energy_groups.group_edges, + np.subtract(self._absorption[iT], + self._fission[iT])) + data.reactions[102] = abs_rxn + + if self._nu_fission[iT] is not None and self._fission[iT] is not None: + fiss_rxn = openmc.data.Reaction(18) + fiss_rxn.xs[strT] = convert_xs(self.energy_groups.group_edges, + self._fission[iT]) + + # Get nu_fission and chi from the presence of both, or the + # nu_fission matrix + if self._nu_fission[iT].shape == self.xs_shapes["[G][G']"]: + nu_fiss = np.sum(self._nu_fission[iT], axis=1)[::-1] + chi = self._nu_fission[iT][::-1] / nu_fiss + else: + nu_fiss = self._nu_fission[iT][::-1] + chi = np.reshape(np.tile(self._chi[iT][::-1], + self.energy_groups.num_groups), + (self.energy_groups.num_groups, + self.energy_groups.num_groups)) + nu = np.divide(nu_fiss, self._fission[iT][::-1]) + + # Build a histogram distribution to represent the yield for the + # incoming groups + prod = openmc.data.Product() + prod.yield_ = \ + openmc.data.Tabulated1D(self.energy_groups.group_edges[:-1], + nu, breakpoints=[len(nu)], + interpolation=[1]) + + # Now build the outgoing energy distribution using discrete energy + # lines within each of the outgoing groups + chi_eouts = [] + for g in range(self.energy_groups.num_groups): + chi_eouts.append(openmc.stats.Discrete(energy_midpoints, + chi[g, :])) + # Ensure the distribution CDF starts with 0 + chi_eouts[-1].c = np.cumsum(chi[g, :]) - chi[g, 0] + + # Create the continuous distribution for the fission + # The angular distribution will remain None (thus isotropic) + chi_distrib = openmc.data.ContinuousTabular( + [len(chi)], [1], self.energy_groups.group_edges[:-1], + chi_eouts) + prod.distribution = \ + [openmc.data.UncorrelatedAngleEnergy(energy=chi_distrib)] + fiss_rxn.products = [prod] + fiss_rxn.center_of_mass = False + data.reactions[18] = fiss_rxn + + # Scattering Data + # First convert the data to a tabular representation + if self.scatter_format == 'tabular': + tabular = self + else: + tabular = self.convert_scatter_format('tabular', 33) + + # Calculate the isotropic scattering matrix and use that to find the + # total scattering x/s and outgoing energy distributions + isotropic_matrix = np.mean(tabular._scatter_matrix[iT], axis=-1)[::-1, + ::-1] + scatt_xs = np.sum(isotropic_matrix, axis=1) + energy = np.zeros((self.energy_groups.num_groups, + self.energy_groups.num_groups)) + for gin in range(self.energy_groups.num_groups): + energy[gin, :] = isotropic_matrix[gin, :] / scatt_xs[gin] + + # Get the anisotropic but normalized angular distribution + distrib = np.zeros((self.energy_groups.num_groups, + self.energy_groups.num_groups, + tabular.num_orders)) + for gin in range(self.energy_groups.num_groups): + for gout in range(self.energy_groups.num_groups): + distrib[gin, gout, :] = \ + np.divide(tabular._scatter_matrix[iT][gin, gout, :], + isotropic_matrix[gin, gout]) + distrib = np.nan_to_num(distrib) + + # Incorporate the scattering multiplication, if required + scatt_prod = openmc.data.Product() + if self._multiplicity_matrix[iT]: + yield_ = self._multiplicity_matrix[iT][::-1, ::-1] + scatt_prod.yield_ = \ + openmc.data.Tabulated1D(self.energy_groups.group_edges[:-1], + yield_, breakpoints=[len(yield_)], + interpolation=[1]) + + # Now build the outgoing energy distribution using discrete energy + # lines within each of the outgoing groups + scatt_eouts = [] + for g in range(self.energy_groups.num_groups): + scatt_eouts.append( + openmc.stats.Discrete(energy_midpoints, energy[g, :])) + # Ensure the distribution CDF starts with 0 + scatt_eouts[-1].c = np.cumsum(energy[g, :]) - energy[g, 0] + scatt_eouts.append(scatt_eouts[-1]) + + # Build the angular distributions associated with each outgoing energy + # group + mu = np.linspace(-1., 1., tabular.num_orders) + scatt_angles = [] + for gin in range(self.energy_groups.num_groups): + scatt_angles.append([]) + for gout in range(self.energy_groups.num_groups): + scatt_angles[gin].append( + openmc.stats.Tabular(mu, distrib[gin, gout, :])) + # Ensure the distribution CDF starts with 0 + scatt_angles[gin][-1].c = \ + np.cumsum(distrib[gin, gout, :]) - distrib[gin, gout, 0] + scatt_angles.append(scatt_angles[-1]) + + # Combine the energy and angle distributions in to a correlated + # angle/energy object + scatt_prod.distribution = \ + [openmc.data.CorrelatedAngleEnergy( + breakpoints=[len(scatt_eouts)], interpolation=[1], + energy=self.energy_groups.group_edges[:-1], + energy_out=scatt_eouts, mu=scatt_angles)] + + # Finally build the reaction with the just-calculated information + # This will be set to the (n,2n) reaction, though any scattering + # reaction would suffice + scatt_rxn = openmc.data.Reaction(16) + scatt_rxn.xs[strT] = \ + convert_xs(self.energy_groups.group_edges, scatt_xs[::-1]) + scatt_rxn.products = [scatt_prod] + scatt_rxn.center_of_mass = False + data.reactions[16] = scatt_rxn + + return data + def to_hdf5(self, file): """Write XSdata to an HDF5 file @@ -2512,6 +2714,29 @@ class MGXSLibrary(object): return library + def convert_to_continuous_energy(self, h5_filename='ce_mgxs.h5', + library_filename='ce_mgxs.xml'): + """Converts the MGXSLibrary object to an equivalent + library of openmc.data.IncidentNeutron objects + + Parameters + ---------- + h5_filename : str + HDF5 file to write with all the files + library_filename : str + cross_sections.xml file describing the HDF5 file + + """ + + library = openmc.data.DataLibrary() + data = [] + for i, xsdata in enumerate(self.xsdatas): + data.append(xsdata.convert_to_continuous_energy()) + data[-1].export_to_hdf5(h5_filename) + + library.register_file(h5_filename) + library.export_to_xml(library_filename) + def export_to_hdf5(self, filename='mgxs.h5'): """Create an hdf5 file that can be used for a simulation.