Removed unused variables, cleaned up other code base, fixed all the bugs (seems to work now!), set default MGXS library temperature (via the API) to 294K to be consistent with the default w/in OpenMC. Looking good folks. Next up is re-implementing the user-facing option for converting Legendre to tabular (it shouldnt be in the library like it was before), implementing a script to convert from xml to h5, and switching the library format to use sparse scattering matrices instead of dense

This commit is contained in:
Adam Nelson 2016-09-04 15:40:02 -04:00
parent aeb7f0c787
commit 7e5326d705
7 changed files with 108 additions and 89 deletions

View file

@ -288,9 +288,12 @@ class Mesh(object):
openmc.YPlane(y0=self.upper_right[1],
boundary_type=bc[3])]
if twod:
zplanes = [openmc.ZPlane(z0=np.finfo(np.float).min,
# Would prefer to have the z ranges be the max supported float, but
# these values are apparently different between python and Fortran.
# Choosing a safe and sane default.
zplanes = [openmc.ZPlane(z0=-1000.,
boundary_type='reflective'),
openmc.ZPlane(z0=np.finfo(np.float).max,
openmc.ZPlane(z0=1000.,
boundary_type='reflective')]
else:
zplanes = [openmc.ZPlane(z0=self.lower_left[2],

View file

@ -907,8 +907,11 @@ class Library(object):
# the provided order or the Library's order.
xsdata.order = min(order, self.legendre_order)
# Right now only 'legendre' data and isotropic weighting is supported
self.scatter_type = 'legendre'
self.representation = 'isotropic'
if nuclide is not 'total':
xsdata.zaid = self._nuclides[nuclide][0]
xsdata.awr = self._nuclides[nuclide][1]
if subdomain is None:

View file

@ -20,13 +20,7 @@ _SCATTER_TYPE_TABULAR = 3
_SCATTER_TYPE_LEGENDRE = 4
_SCATTER_TYPE_HISTOGRAM = 5
_REPRESENTATIONS = ['isotropic', 'angle']
_REPRESENTATIONS_TEXT = {'isotropic': _REPRESENTATION_ISOTROPIC,
'angle': _REPRESENTATION_ANGLE}
_SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
_SCATTER_TYPES_TEXT = {'tabular': _SCATTER_TYPE_TABULAR,
'legendre': _SCATTER_TYPE_LEGENDRE,
'histogram': _SCATTER_TYPE_HISTOGRAM}
class XSdata(object):
"""A multi-group cross section data set providing all the
@ -43,7 +37,7 @@ class XSdata(object):
weighting). Defaults to 'isotropic'
temperatures : numpy.ndarray
Temperatures (in units of Kelvin) of the provided datasets. Defaults
to a single temperature at 300K.
to a single temperature at 294K.
Attributes
----------
@ -54,7 +48,7 @@ class XSdata(object):
of the isotope to the mass of a single neutron.
temperatures : numpy.ndarray
Temperatures (in units of Kelvin) of the provided datasets. Defaults
to a single temperature at 300K.
to a single temperature at 294K.
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure
fissionable : bool
@ -164,7 +158,7 @@ class XSdata(object):
"""
def __init__(self, name, energy_groups, temperatures=[300.],
def __init__(self, name, energy_groups, temperatures=[294.],
representation='isotropic'):
# Initialize class attributes
self.name = name
@ -377,7 +371,7 @@ class XSdata(object):
check_type('use_chi', use_chi, bool)
self._use_chi = use_chi
def set_total(self, total, temperature=300.):
def set_total(self, total, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -387,7 +381,7 @@ class XSdata(object):
Total Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -405,7 +399,7 @@ class XSdata(object):
i = self.temperatures.tolist().index(temperature)
self._total[i] = nptotal
def set_absorption(self, absorption, temperature=300.):
def set_absorption(self, absorption, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -415,7 +409,7 @@ class XSdata(object):
Absorption Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -434,7 +428,7 @@ class XSdata(object):
i = self.temperatures.tolist().index(temperature)
self._absorption[i] = npabsorption
def set_fission(self, fission, temperature=300.):
def set_fission(self, fission, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -444,7 +438,7 @@ class XSdata(object):
Fission Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -465,7 +459,7 @@ class XSdata(object):
if np.sum(npfission) > 0.0:
self._fissionable = True
def set_kappa_fission(self, kappa_fission, temperature=300.):
def set_kappa_fission(self, kappa_fission, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -475,7 +469,7 @@ class XSdata(object):
Kappa-Fission Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -498,7 +492,7 @@ class XSdata(object):
if np.sum(npkappa_fission) > 0.0:
self._fissionable = True
def set_chi(self, chi, temperature=300.):
def set_chi(self, chi, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -508,7 +502,7 @@ class XSdata(object):
Fission Spectrum
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -538,7 +532,7 @@ class XSdata(object):
if self.use_chi is not None:
self.use_chi = True
def set_scatter_matrix(self, scatter, temperature=300.):
def set_scatter_matrix(self, scatter, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -548,7 +542,7 @@ class XSdata(object):
Scattering Matrix Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -567,7 +561,7 @@ class XSdata(object):
i = self.temperatures.tolist().index(temperature)
self._scatter_matrix[i] = npscatter
def set_multiplicity_matrix(self, multiplicity, temperature=300.):
def set_multiplicity_matrix(self, multiplicity, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -577,7 +571,7 @@ class XSdata(object):
Multiplicity Matrix Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -597,7 +591,7 @@ class XSdata(object):
i = self.temperatures.tolist().index(temperature)
self._multiplicity_matrix[i] = npmultiplicity
def set_nu_fission(self, nu_fission, temperature=300.):
def set_nu_fission(self, nu_fission, temperature=294.):
"""This method sets the cross section for this XSdata object at the
provided temperature.
@ -607,7 +601,7 @@ class XSdata(object):
Nu-fission Cross Section
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
See also
--------
@ -655,7 +649,7 @@ class XSdata(object):
if np.sum(npnu_fission) > 0.0:
self._fissionable = True
def set_total_mgxs(self, total, temperature=300., nuclide='total',
def set_total_mgxs(self, total, temperature=294., nuclide='total',
xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.TotalXS or
openmc.mgxs.TransportXS to be used to set the total cross section
@ -668,7 +662,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -702,7 +696,7 @@ class XSdata(object):
msg = 'Angular-Dependent MGXS have not yet been implemented'
raise ValueError(msg)
def set_absorption_mgxs(self, absorption, temperature=300.,
def set_absorption_mgxs(self, absorption, temperature=294.,
nuclide='total', xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.AbsorptionXS
to be used to set the absorption cross section for this XSdata object.
@ -714,7 +708,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -749,7 +743,7 @@ class XSdata(object):
msg = 'Angular-Dependent MGXS have not yet been implemented'
raise ValueError(msg)
def set_fission_mgxs(self, fission, temperature=300., nuclide='total',
def set_fission_mgxs(self, fission, temperature=294., nuclide='total',
xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.FissionXS
to be used to set the fission cross section for this XSdata object.
@ -761,7 +755,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -796,7 +790,7 @@ class XSdata(object):
msg = 'Angular-Dependent MGXS have not yet been implemented'
raise ValueError(msg)
def set_nu_fission_mgxs(self, nu_fission, temperature=300.,
def set_nu_fission_mgxs(self, nu_fission, temperature=294.,
nuclide='total', xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.NuFissionXS
to be used to set the nu-fission cross section for this XSdata object.
@ -808,7 +802,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -852,7 +846,7 @@ class XSdata(object):
if np.sum(self._nu_fission) > 0.0:
self._fissionable = True
def set_kappa_fission_mgxs(self, k_fission, temperature=300.,
def set_kappa_fission_mgxs(self, k_fission, temperature=294.,
nuclide='total', xs_type='macro',
subdomain=None):
"""This method allows for an openmc.mgxs.KappaFissionXS
@ -866,7 +860,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -901,7 +895,7 @@ class XSdata(object):
msg = 'Angular-Dependent MGXS have not yet been implemented'
raise ValueError(msg)
def set_chi_mgxs(self, chi, temperature=300., nuclide='total',
def set_chi_mgxs(self, chi, temperature=294., nuclide='total',
xs_type='macro', subdomain=None):
"""This method allows for an openmc.mgxs.Chi
to be used to set chi for this XSdata object.
@ -912,7 +906,7 @@ class XSdata(object):
MGXS Object containing chi for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -954,7 +948,7 @@ class XSdata(object):
if self.use_chi is not None:
self.use_chi = True
def set_scatter_matrix_mgxs(self, scatter, temperature=300.,
def set_scatter_matrix_mgxs(self, scatter, temperature=294.,
nuclide='total', xs_type='macro',
subdomain=None):
"""This method allows for an openmc.mgxs.ScatterMatrixXS
@ -969,7 +963,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -1027,7 +1021,7 @@ class XSdata(object):
raise ValueError(msg)
def set_multiplicity_matrix_mgxs(self, nuscatter, scatter=None,
temperature=300., nuclide='total',
temperature=294., nuclide='total',
xs_type='macro', subdomain=None):
"""This method allows for either the direct use of only an
openmc.mgxs.MultiplicityMatrixXS OR
@ -1050,7 +1044,7 @@ class XSdata(object):
for the domain of interest.
temperature : float
Temperature (in units of Kelvin) of the provided dataset. Defaults
to 300K
to 294K
nuclide : str
Individual nuclide (or 'total' if obtaining material-wise data)
to gather data for. Defaults to 'total'.
@ -1119,22 +1113,26 @@ class XSdata(object):
if self.fissionable is not None:
grp.attrs['fissionable'] = self.fissionable
if self.representation is not None:
grp.attrs['representation'] = \
np.string_(_REPRESENTATIONS_TEXT[self.representation])
grp.attrs['representation'] = np.string_(self.representation)
if self.representation == 'angle':
if self.num_azimuthal is not None:
grp.attrs['num-azimuthal'] = self.num_azimuthal
if self.num_polar is not None:
grp.attrs['num-polar'] = self.num_polar
if self.scatter_type is not None:
grp.attrs['scatter-type'] = \
np.string_(_SCATTER_TYPES_TEXT[self.scatter_type])
grp.attrs['scatter-type'] = np.string_(self.scatter_type)
if self.order is not None:
grp.attrs['order'] = self.order
ktg = grp.create_group('kTs')
for temperature in self.temperatures:
temp_label = str(int(np.round(temperature))) + "K"
kT = temperature * openmc.data.K_BOLTZMANN
ktg.create_dataset(temp_label, data=kT)
# Create the temperature datasets
for i, temperature in enumerate(self.temperatures):
xsgrp = grp.create_group(str(np.round(temperature)) + "K")
xsgrp = grp.create_group(str(int(np.round(temperature))) + "K")
if self._total[i] is not None:
xsgrp.create_dataset("total", data=self._total[i],
compression=compression)

View file

@ -16,6 +16,7 @@ module input_xml
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
use mesh_header, only: RegularMesh
use mgxs_data, only: create_macro_xs, read_mgxs
use multipole, only: multipole_read
use output, only: write_message
use plot_header

View file

@ -1,6 +1,7 @@
module mgxs_data
use constants
use algorithm, only: find
use error, only: fatal_error
use global
use hdf5
@ -28,6 +29,7 @@ contains
integer :: i_nuclide ! index in nuclides array
character(20) :: name ! name of library to load
integer :: representation ! Data representation
character(MAX_LINE_LEN) :: temp_str
type(Material), pointer :: mat
type(SetChar) :: already_read
integer(HID_T) :: file_id
@ -99,9 +101,12 @@ contains
! First find out the data representation
if (check_attribute(xsdata_group, "representation")) then
call read_attribute(representation, xsdata_group, "representation")
if (representation /= MGXS_ISOTROPIC .and. &
representation /= MGXS_ANGLE) then
call read_attribute(temp_str, xsdata_group, "representation")
if (trim(temp_str) == 'isotropic') then
representation = MGXS_ISOTROPIC
else if (trim(temp_str) == 'angle') then
representation = MGXS_ANGLE
else
call fatal_error("Invalid Data Representation!")
end if
else
@ -179,7 +184,9 @@ contains
allocate(MgxsAngle :: macro_xs(i_mat) % obj)
end select
call macro_xs(i_mat) % obj % combine(kTs(i_mat), mat, nuclides_MG, &
energy_groups, max_order)
energy_groups, max_order, &
temperature_tolerance, &
temperature_method)
end do
end subroutine create_macro_xs

View file

@ -6,6 +6,7 @@ module mgxs_header
use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5iget_name_f
use h5lt, only: h5ltpath_valid_f
use algorithm, only: find, sort
use constants, only: MAX_FILE_LEN, ZERO, ONE, TWO, PI
use error, only: fatal_error
use hdf5_interface
@ -100,7 +101,8 @@ module mgxs_header
integer, intent(in) :: max_order ! Maximum requested order
end subroutine mgxs_from_hdf5_
subroutine mgxs_combine_(this, temps, mat, nuclides, groups, max_order)
subroutine mgxs_combine_(this, temps, mat, nuclides, groups, max_order, &
tolerance, method)
import Mgxs, Material, MgxsContainer, VectorReal
class(Mgxs), intent(inout) :: this ! The Mgxs to initialize
type(VectorReal), intent(in) :: temps ! Temperatures to obtain
@ -108,6 +110,8 @@ module mgxs_header
type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
integer, intent(in) :: max_order ! Maximum requested order
real(8), intent(in) :: tolerance ! Tolerance on method
integer, intent(in) :: method ! Type of temperature access
end subroutine mgxs_combine_
pure function mgxs_get_xs_(this, xstype, gin, gout, uvw, mu) result(xs_val)
@ -206,12 +210,11 @@ module mgxs_header
integer :: hdf5_err
integer(SIZE_T) :: name_len, name_file_len
integer(HID_T) :: kT_group, kT_dset
integer(HID_T) :: kT_group
character(MAX_FILE_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
logical :: exists
character(MAX_WORD_LEN) :: temp_str
real(8) :: dangle
integer :: ipol, iazi
@ -236,6 +239,7 @@ module mgxs_header
do i = 1, size(dset_names)
! Read temperature value
call read_dataset(temps_available(i), kT_group, trim(dset_names(i)))
! Convert MeV to Kelvin
temps_available(i) = temps_available(i) / K_BOLTZMANN
end do
@ -377,7 +381,7 @@ module mgxs_header
real(8), allocatable :: input_scatt(:, :, :)
real(8), allocatable :: temp_scatt(:, :, :)
real(8) :: dmu, mu, norm
integer :: order, order_dim, gin, gout, l, arr_len
integer :: order, order_dim, gin, gout, l
integer :: legendre_mu_points, imu
type(VectorInt) :: temps_to_read
integer :: t
@ -387,12 +391,13 @@ module mgxs_header
temps_to_read, order_dim)
!!!TODO: Fix
enable_leg_mu = .True.
legendre_mu_points = 33
! Load the more specific data
do t = 1, temps_to_read % size()
associate(xs => this % xs(t))
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
temp_str = trim(to_str(temps_to_read % data(t))) // "K"
xsdata_grp = open_group(xs_id, trim(temp_str))
if (this % fissionable) then
allocate(xs % nu_fission(groups))
@ -619,8 +624,8 @@ module mgxs_header
real(8), allocatable :: scatt_coeffs(:, :, :, :, :)
real(8), allocatable :: input_scatt(:, :, :, :, :)
real(8), allocatable :: temp_scatt(:, :, :, :, :)
real(8) :: dmu, mu, norm, dangle
integer :: order, order_dim, gin, gout, l, arr_len
real(8) :: dmu, mu, norm
integer :: order, order_dim, gin, gout, l
integer :: legendre_mu_points, imu, ipol, iazi
type(VectorInt) :: temps_to_read
integer :: t
@ -630,12 +635,13 @@ module mgxs_header
temps_to_read, order_dim)
!!!TODO: Fix
enable_leg_mu = .True.
legendre_mu_points = 33
! Load the more specific data
do t = 1, temps_to_read % size()
associate(xs => this % xs(t))
! Get temperature as a string
temp_str = trim(to_str(temps_to_read % data(i))) // "K"
temp_str = trim(to_str(temps_to_read % data(t))) // "K"
xsdata_grp = open_group(xs_id, trim(temp_str))
if (this % fissionable) then
if (check_dataset(xsdata_grp, "chi")) then
@ -899,16 +905,17 @@ module mgxs_header
! MGXS*_COMBINE Builds a macroscopic Mgxs object from microscopic Mgxs objects
!===============================================================================
subroutine mgxs_combine(this, temps, mat, nuclides, scatter_type, order_dim)
subroutine mgxs_combine(this, temps, mat, nuclides, max_order, &
scatter_type, order_dim)
class(Mgxs), intent(inout) :: this ! The Mgxs to initialize
type(VectorReal), intent(in) :: temps ! Temperatures to obtain
type(Material), pointer, intent(in) :: mat ! base material
type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: max_order ! Maximum requested order
integer, intent(out) :: scatter_type ! Type of scatter
integer, intent(out) :: order_dim ! Scattering data order size
integer :: t, mat_max_order, max_order, order, n_pol, n_azi
real(8) :: dangle
integer :: t, mat_max_order, order
! Fill in meta-data from material information
if (mat % name == "") then
@ -924,7 +931,7 @@ module mgxs_header
allocate(this % kTs(temps % size()))
do t = 1, temps % size()
this % kTs(t) = temps % data(i)
this % kTs(t) = temps % data(t)
end do
! Allocate the XS object for the number of temperatures
@ -1004,20 +1011,23 @@ module mgxs_header
end subroutine mgxs_combine
subroutine mgxsiso_combine(this, temps, mat, nuclides, groups, max_order)
subroutine mgxsiso_combine(this, temps, mat, nuclides, groups, max_order, &
tolerance, method)
class(MgxsIso), intent(inout) :: this ! The Mgxs to initialize
type(VectorReal), intent(in) :: temps ! Temperatures to obtain [MeV]
type(Material), pointer, intent(in) :: mat ! base material
type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
integer, intent(in) :: max_order ! Maximum requested order
real(8), intent(in) :: tolerance ! Tolerance on method
integer, intent(in) :: method ! Type of temperature access
integer :: i ! loop index over nuclides
integer :: t ! Index in to temps
integer :: gin, gout ! group indices
real(8) :: atom_density ! atom density of a nuclide
real(8) :: norm, nuscatt
integer :: mat_max_order, order, order_dim, nuc_order_dim
integer :: order_dim, nuc_order_dim
real(8), allocatable :: temp_mult(:, :), mult_num(:, :), mult_denom(:, :)
real(8), allocatable :: scatt_coeffs(:, :, :)
integer :: nuc_t
@ -1025,7 +1035,8 @@ module mgxs_header
integer :: scatter_type
! Set the meta-data
call mgxs_combine(this, temps, mat, nuclides, scatter_type, order_dim)
call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_type, &
order_dim)
! Create the Xs Data for each temperature
TEMP_LOOP: do t = 1, temps % size()
@ -1067,11 +1078,11 @@ module mgxs_header
! Determine actual temperatures to read
temp_desired = temps % data(i)
nuc_t = minloc(abs(nuc % kTs - temp_desired), dim=1)
temp_actual = nuc % kTs(i_closest)
temp_actual = nuc % kTs(nuc_t)
if (abs(temp_actual - temp_desired) >= tolerance) then
call fatal_error("MGXS library does not contain cross sections &
&for " // trim(this % name) // " at or near " // &
trim(to_str(nint(temp_desired))) // " K.")
trim(to_str(nint(temp_desired / K_BOLTZMANN))) // " K.")
end if
case (TEMPERATURE_INTERPOLATION)
@ -1174,13 +1185,16 @@ module mgxs_header
end subroutine mgxsiso_combine
subroutine mgxsang_combine(this, temps, mat, nuclides, groups, max_order)
subroutine mgxsang_combine(this, temps, mat, nuclides, groups, max_order, &
tolerance, method)
class(MgxsAngle), intent(inout) :: this ! The Mgxs to initialize
type(VectorReal), intent(in) :: temps ! Temperatures to obtain
type(Material), pointer, intent(in) :: mat ! base material
type(MgxsContainer), intent(in) :: nuclides(:) ! List of nuclides to harvest from
integer, intent(in) :: groups ! Number of E groups
integer, intent(in) :: max_order ! Maximum requested order
real(8), intent(in) :: tolerance ! Tolerance on method
integer, intent(in) :: method ! Type of temperature access
integer :: i ! loop index over nuclides
integer :: t ! temperature loop index
@ -1188,7 +1202,7 @@ module mgxs_header
real(8) :: atom_density ! atom density of a nuclide
integer :: ipol, iazi, n_pol, n_azi
real(8) :: norm, nuscatt
integer :: mat_max_order, order, order_dim, nuc_order_dim
integer :: order_dim, nuc_order_dim
real(8), allocatable :: temp_mult(:, :, :, :), mult_num(:, :, :, :)
real(8), allocatable :: mult_denom(:, :, :, :), scatt_coeffs(:, :, :, :, :)
integer :: nuc_t
@ -1196,7 +1210,8 @@ module mgxs_header
integer :: scatter_type
! Set the meta-data
call mgxs_combine(this, temps, mat, nuclides, scatter_type, order_dim)
call mgxs_combine(this, temps, mat, nuclides, max_order, scatter_type, &
order_dim)
! Get the number of each polar and azi angles and make sure all the
! NuclideAngle types have the same number of these angles
@ -1268,11 +1283,11 @@ module mgxs_header
! Determine actual temperatures to read
temp_desired = temps % data(i)
nuc_t = minloc(abs(nuc % kTs - temp_desired), dim=1)
temp_actual = nuc % kTs(i_closest)
temp_actual = nuc % kTs(nuc_t)
if (abs(temp_actual - temp_desired) >= tolerance) then
call fatal_error("MGXS library does not contain cross sections &
&for " // trim(this % name) // " at or near " // &
trim(to_str(nint(temp_desired))) // " K.")
trim(to_str(nint(temp_desired / K_BOLTZMANN))) // " K.")
end if
case (TEMPERATURE_INTERPOLATION)
@ -1325,7 +1340,7 @@ module mgxs_header
nuscatt = nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % scattxs(gin) * &
nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % energy(gin) % data(gout)
mult_num(iazi, ipol, gout, gin) = &
mult_num(iazi, ipo, gout, gin) + atom_density * &
mult_num(iazi, ipol, gout, gin) + atom_density * &
nuscatt
if (nuc % xs(nuc_t) % scatter(iazi, ipol) % obj % mult(gin) % data(gout) > ZERO) then
mult_denom(iazi, ipol, gout, gin) = &
@ -1716,18 +1731,11 @@ module mgxs_header
! sqrt(temperature), (with temperature in units of MeV)
!===============================================================================
pure subroutine mgxs_find_temperature(this, sqrtkT, tolerance)
class(Mgxs), intent(in) :: this
real(8), intent(in) :: sqrtkT ! Temperature (in units of of MeV)
real(8), intent(in) :: tolerance ! Temperature tolerance
pure subroutine mgxs_find_temperature(this, sqrtkT)
class(Mgxs), intent(inout) :: this
real(8), intent(in) :: sqrtkT ! Temperature (in units of of MeV)
real(8) :: kT
integer :: i_temp
kT = sqrtkT**2
do i_temp = 1, size(this % kTs)
if (abs(this % kTs(i_temp) - kT) < K_BOLTZMANN * tolerance) exit
end do
this % index_temp = minloc(abs(this % kTs - (sqrtkT * sqrtkT)), dim=1)
end subroutine mgxs_find_temperature

View file

@ -93,8 +93,7 @@ contains
! After every collision for the MGXS mode
if (p % material /= MATERIAL_VOID) then
! Update the temperature index
call macro_xs(p % material) % obj % find_temperature(p % sqrtkT, &
temperature_tolerance)
call macro_xs(p % material) % obj % find_temperature(p % sqrtkT)
! Get the data
call macro_xs(p % material) % obj % calculate_xs(p % g, &
p % coord(p % n_coord) % uvw, material_xs)