mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Fixes to scatter format conversion as well as initial implementation of conversion of MG library to CE
This commit is contained in:
Adam Nelson
parent
1c367b049c
commit
7e1700227c
1 changed files with 239 additions and 14 deletions
|
|
@ -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.
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue