diff --git a/openmc/mesh.py b/openmc/mesh.py index 7d7b483f73..43ac7aef05 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -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], diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index 9f0e173f3c..a04ec8c557 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -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: diff --git a/openmc/mgxs_library.py b/openmc/mgxs_library.py index 68c0063f2b..15d80dc913 100644 --- a/openmc/mgxs_library.py +++ b/openmc/mgxs_library.py @@ -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) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index eafbc717c2..c44c6b7bda 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 7cbb4864c7..d4ad96da57 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -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 diff --git a/src/mgxs_header.F90 b/src/mgxs_header.F90 index c86685bc0d..0c0bcc94f2 100644 --- a/src/mgxs_header.F90 +++ b/src/mgxs_header.F90 @@ -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 diff --git a/src/tracking.F90 b/src/tracking.F90 index a03cde98dd..b62cfc85a0 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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)