From 7e5326d705b90f5081eca9206388cc91fa6bafb9 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 4 Sep 2016 15:40:02 -0400 Subject: [PATCH] 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 --- openmc/mesh.py | 7 +++- openmc/mgxs/library.py | 5 ++- openmc/mgxs_library.py | 90 +++++++++++++++++++++--------------------- src/input_xml.F90 | 1 + src/mgxs_data.F90 | 15 +++++-- src/mgxs_header.F90 | 76 +++++++++++++++++++---------------- src/tracking.F90 | 3 +- 7 files changed, 108 insertions(+), 89 deletions(-) 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)