diff --git a/CMakeLists.txt b/CMakeLists.txt index d3673df26c..be62dcd57b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -368,7 +368,6 @@ set(LIBOPENMC_FORTRAN_SRC src/message_passing.F90 src/mgxs_data.F90 src/mgxs_header.F90 - src/multipole.F90 src/multipole_header.F90 src/nuclide_header.F90 src/output.F90 diff --git a/docs/source/io_formats/data_wmp.rst b/docs/source/io_formats/data_wmp.rst index d645c6e8f1..6174e85ff7 100644 --- a/docs/source/io_formats/data_wmp.rst +++ b/docs/source/io_formats/data_wmp.rst @@ -28,12 +28,6 @@ Windowed Multipole Library Format ":math:`r`" and ":math:`i`" identifiers, similar to how `h5py`_ does it. - **end_E** (*double*) Highest energy the windowed multipole part of the library is valid for. - - **energy_points** (*double[]*) - Energy grid for the pointwise library in the reaction group. - - **fissionable** (*int*) - 1 if this nuclide has fission data. 0 if it does not. - - **fit_order** (*int*) - The order of the curve fit. - **formalism** (*int*) The formalism of the underlying data. Uses the `ENDF-6`_ format formalism numbers. @@ -51,18 +45,6 @@ Windowed Multipole Library Format - **l_value** (*int[]*) The index for a corresponding pole. Equivalent to the :math:`l` quantum number of the resonance the pole comes from :math:`+1`. - - **length** (*int*) - Total count of poles in `data`. - - **max_w** (*int*) - Maximum number of poles in a window. - - **MT_count** (*int*) - Number of pointwise tables in the library. - - **MT_list** (*int[]*) - A list of available MT identifiers. See `ENDF-6`_ for meaning. - - **n_grid** (*int*) - Total length of the pointwise data. - - **num_l** (*int*) - Number of possible :math:`l` quantum states for this nuclide. - **pseudo_K0RS** (*double[]*) :math:`l` dependent value of @@ -90,13 +72,6 @@ Windowed Multipole Library Format The pole to start from for each window. - **w_end** (*int[]*) The pole to end at for each window. - - **windows** (*int*) - Number of windows. - -**/nuclide/reactions/MT** - - **MT_sigma** (*double[]*) -- Cross section value for this reaction. - - **Q_value** (*double*) -- Energy released in this reaction, in eV. - - **threshold** (*int*) -- The first non-zero entry in ``MT_sigma``. .. _h5py: http://docs.h5py.org/en/latest/ .. _ENDF-6: https://www.oecd-nea.org/dbdata/data/manual-endf/endf102.pdf diff --git a/openmc/data/multipole.py b/openmc/data/multipole.py index 33078f5046..f7a78d9532 100644 --- a/openmc/data/multipole.py +++ b/openmc/data/multipole.py @@ -24,7 +24,7 @@ _RM_RF = 3 # Residue fission # Multi-level Breit Wigner indices _MLBW_RT = 1 # Residue total -_MLBW_RX = 2 # Residue compettitive +_MLBW_RX = 2 # Residue competitive _MLBW_RA = 3 # Residue absorption _MLBW_RF = 4 # Residue fission @@ -141,6 +141,12 @@ def _broaden_wmp_polynomials(E, dopp, n): class WindowedMultipole(EqualityMixin): """Resonant cross sections represented in the windowed multipole format. + Parameters + ---------- + formalism : {'MLBW', 'RM'} + The R-matrix formalism used to reconstruct resonances. Either 'MLBW' + for multi-level Breit Wigner or 'RM' for Reich-Moore. + Attributes ---------- num_l : Integral @@ -195,11 +201,9 @@ class WindowedMultipole(EqualityMixin): a/E + b/sqrt(E) + c + d sqrt(E) + ... """ - def __init__(self): - self.num_l = None - self.fit_order = None - self.fissionable = None - self.formalism = None + def __init__(self, formalism): + self._num_l = None + self.formalism = formalism self.spacing = None self.sqrtAWR = None self.start_E = None @@ -218,11 +222,15 @@ class WindowedMultipole(EqualityMixin): @property def fit_order(self): - return self._fit_order + return self.curvefit.shape[1] - 1 @property def fissionable(self): - return self._fissionable + if self.formalism == 'RM': + return self.data.shape[1] == 4 + else: + # Assume self.formalism == 'MLBW' + return self.data.shape[1] == 5 @property def formalism(self): @@ -272,35 +280,10 @@ class WindowedMultipole(EqualityMixin): def curvefit(self): return self._curvefit - @num_l.setter - def num_l(self, num_l): - if num_l is not None: - cv.check_type('num_l', num_l, Integral) - cv.check_greater_than('num_l', num_l, 1, equality=True) - cv.check_less_than('num_l', num_l, 4, equality=True) - # There is an if block in _evaluate that assumes num_l <= 4. - self._num_l = num_l - - @fit_order.setter - def fit_order(self, fit_order): - if fit_order is not None: - cv.check_type('fit_order', fit_order, Integral) - cv.check_greater_than('fit_order', fit_order, 2, equality=True) - # _broaden_wmp_polynomials assumes the curve fit has at least 3 - # terms. - self._fit_order = fit_order - - @fissionable.setter - def fissionable(self, fissionable): - if fissionable is not None: - cv.check_type('fissionable', fissionable, bool) - self._fissionable = fissionable - @formalism.setter def formalism(self, formalism): - if formalism is not None: - cv.check_type('formalism', formalism, str) - cv.check_value('formalism', formalism, ('MLBW', 'RM')) + cv.check_type('formalism', formalism, str) + cv.check_value('formalism', formalism, ('MLBW', 'RM')) self._formalism = formalism @spacing.setter @@ -337,9 +320,20 @@ class WindowedMultipole(EqualityMixin): cv.check_type('data', data, np.ndarray) if len(data.shape) != 2: raise ValueError('Multipole data arrays must be 2D') - if data.shape[1] not in (3, 4, 5): # 3 or 4 for RM, 4 or 5 for MLBW - raise ValueError('The second dimension of multipole data arrays' - ' must have a length of 3, 4 or 5') + if self.formalism == 'RM': + if data.shape[1] not in (3, 4): + raise ValueError('For the Reich-Moore formalism, ' + 'data.shape[1] must be 3 or 4. One value for the pole.' + ' One each for the total and absorption residues. ' + 'Possibly one more for a fission residue.') + else: + # Assume self.formalism == 'MLBW' + if data.shape[1] not in (4, 5): + raise ValueError('For the Multi-level Breit-Wigner ' + 'formalism, data.shape[1] must be 4 or 5. One value ' + 'for the pole. One each for the total, competitive, ' + 'and absorption residues. Possibly one more for a ' + 'fission residue.') if not np.issubdtype(data.dtype, complex): raise TypeError('Multipole data arrays must be complex dtype') self._data = data @@ -363,6 +357,12 @@ class WindowedMultipole(EqualityMixin): if not np.issubdtype(l_value.dtype, int): raise TypeError('Multipole l_value arrays must be integer' ' dtype') + + self._num_l = len(np.unique(l_value)) + + else: + self._num_l = None + self._l_value = l_value @w_start.setter @@ -428,6 +428,7 @@ class WindowedMultipole(EqualityMixin): format. """ + if isinstance(group_or_filename, h5py.Group): group = group_or_filename else: @@ -442,20 +443,12 @@ class WindowedMultipole(EqualityMixin): 'Python API expects version ' + WMP_VERSION) group = h5file['nuclide'] - out = cls() - - # Read scalar values. Note that group['max_w'] is ignored. - - length = group['length'].value - windows = group['windows'].value - out.num_l = group['num_l'].value - out.fit_order = group['fit_order'].value - out.fissionable = bool(group['fissionable'].value) + # Read scalars. if group['formalism'].value == _FORM_MLBW: - out.formalism = 'MLBW' + out = cls('MLBW') elif group['formalism'].value == _FORM_RM: - out.formalism = 'RM' + out = cls('RM') else: raise ValueError('Unrecognized/Unsupported R-matrix formalism') @@ -466,39 +459,36 @@ class WindowedMultipole(EqualityMixin): # Read arrays. - err = "WMP '{}' array shape is not consistent with the '{}' value" + err = "WMP '{}' array shape is not consistent with the '{}' array shape" out.data = group['data'].value - if out.data.shape[0] != length: - raise ValueError(err.format('data', 'length')) + + out.l_value = group['l_value'].value + if out.l_value.shape[0] != out.data.shape[0]: + raise ValueError(err.format('l_value', 'data')) out.pseudo_k0RS = group['pseudo_K0RS'].value if out.pseudo_k0RS.shape[0] != out.num_l: - raise ValueError(err.format('pseudo_k0RS', 'num_l')) - - out.l_value = group['l_value'].value - if out.l_value.shape[0] != length: - raise ValueError(err.format('l_value', 'length')) + raise ValueError(err.format('pseudo_k0RS', 'l_value')) out.w_start = group['w_start'].value - if out.w_start.shape[0] != windows: - raise ValueError(err.format('w_start', 'windows')) out.w_end = group['w_end'].value - if out.w_end.shape[0] != windows: - raise ValueError(err.format('w_end', 'windows')) + if out.w_end.shape[0] != out.w_start.shape[0]: + raise ValueError(err.format('w_end', 'w_start')) out.broaden_poly = group['broaden_poly'].value.astype(np.bool) - if out.broaden_poly.shape[0] != windows: - raise ValueError(err.format('broaden_poly', 'windows')) + if out.broaden_poly.shape[0] != out.w_start.shape[0]: + raise ValueError(err.format('broaden_poly', 'w_start')) out.curvefit = group['curvefit'].value - if out.curvefit.shape[0] != windows: - raise ValueError(err.format('curvefit', 'windows')) - if out.curvefit.shape[1] != out.fit_order + 1: - raise ValueError(err.format('curvefit', 'fit_order')) + if out.curvefit.shape[0] != out.w_start.shape[0]: + raise ValueError(err.format('curvefit', 'w_start')) - # Note that all the file 3 data (group['reactions/MT...']) are ignored. + # _broaden_wmp_polynomials assumes the curve fit has at least 3 terms. + if out.fit_order < 2: + raise ValueError("Windowed multipole is only supported for " + "curvefits with 3 or more terms.") return out @@ -661,3 +651,47 @@ class WindowedMultipole(EqualityMixin): fun = np.vectorize(lambda x: self._evaluate(x, T)) return fun(E) + + def export_to_hdf5(self, path, libver='earliest'): + """Export windowed multipole data to an HDF5 file. + + Parameters + ---------- + path : str + Path to write HDF5 file to + libver : {'earliest', 'latest'} + Compatibility mode for the HDF5 file. 'latest' will produce files + that are less backwards compatible but have performance benefits. + + """ + + # Open file and write version. + with h5py.File(path, 'w', libver=libver) as f: + f.create_dataset('version', (1, ), dtype='S10') + f['version'][:] = WMP_VERSION.encode('ASCII') + + # Make a nuclide group. + g = f.create_group('nuclide') + + # Write scalars. + if self.formalism == 'MLBW': + g.create_dataset('formalism', + data=np.array(_FORM_MLBW, dtype=np.int32)) + else: + # Assume RM. + g.create_dataset('formalism', + data=np.array(_FORM_RM, dtype=np.int32)) + g.create_dataset('spacing', data=np.array(self.spacing)) + g.create_dataset('sqrtAWR', data=np.array(self.sqrtAWR)) + g.create_dataset('start_E', data=np.array(self.start_E)) + g.create_dataset('end_E', data=np.array(self.end_E)) + + # Write arrays. + g.create_dataset('data', data=self.data) + g.create_dataset('l_value', data=self.l_value) + g.create_dataset('pseudo_K0RS', data=self.pseudo_k0RS) + g.create_dataset('w_start', data=self.w_start) + g.create_dataset('w_end', data=self.w_end) + g.create_dataset('broaden_poly', + data=self.broaden_poly.astype(np.int8)) + g.create_dataset('curvefit', data=self.curvefit) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 8c3669b2e0..918d471d3c 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -21,7 +21,6 @@ module input_xml use message_passing use mgxs_data, only: create_macro_xs, read_mgxs use mgxs_header - use multipole, only: multipole_read use nuclide_header use output, only: title, header, print_plot use plot_header @@ -4377,7 +4376,7 @@ contains allocate(nuc % multipole) ! Call the read routine - call multipole_read(filename, nuc % multipole, i_table) + call nuc % multipole % from_hdf5(filename) nuc % mp_present = .true. end associate diff --git a/src/multipole.F90 b/src/multipole.F90 deleted file mode 100644 index 0282adb7d9..0000000000 --- a/src/multipole.F90 +++ /dev/null @@ -1,89 +0,0 @@ -module multipole - - use hdf5 - - use constants - use error, only: fatal_error - use hdf5_interface - use multipole_header, only: MultipoleArray, FIT_T, FIT_A, FIT_F, & - MP_FISS, FORM_MLBW, FORM_RM - use nuclide_header, only: nuclides - - implicit none - -contains - -!=============================================================================== -! MULTIPOLE_READ Reads in a multipole HDF5 file with the original API -! specification. Subject to change as the library format matures. -!=============================================================================== - - subroutine multipole_read(filename, multipole, i_table) - character(len=*), intent(in) :: filename ! Filename of the - ! multipole library - ! to load - type(MultipoleArray), intent(out), target :: multipole ! The object to fill - integer, intent(in) :: i_table ! index in nuclides/ - ! sab_tables - - integer(HID_T) :: file_id - integer(HID_T) :: group_id - - ! Intermediate loading components - integer :: is_fissionable - character(len=10) :: version - - associate (nuc => nuclides(i_table)) - - ! Open file for reading and move into the /isotope group - file_id = file_open(filename, 'r', parallel=.true.) - group_id = open_group(file_id, "/nuclide") - - ! Check the file version number. - call read_dataset(version, file_id, "version") - if (version /= VERSION_MULTIPOLE) call fatal_error("The current multipole& - & format version is " // trim(VERSION_MULTIPOLE) // " but the file "& - // trim(filename) // " uses version " // trim(version) // ".") - - ! Load in all the array size scalars - call read_dataset(multipole % length, group_id, "length") - call read_dataset(multipole % windows, group_id, "windows") - call read_dataset(multipole % num_l, group_id, "num_l") - call read_dataset(multipole % fit_order, group_id, "fit_order") - call read_dataset(multipole % max_w, group_id, "max_w") - call read_dataset(is_fissionable, group_id, "fissionable") - if (is_fissionable == MP_FISS) then - multipole % fissionable = .true. - else - multipole % fissionable = .false. - end if - call read_dataset(multipole % formalism, group_id, "formalism") - - call read_dataset(multipole % spacing, group_id, "spacing") - call read_dataset(multipole % sqrtAWR, group_id, "sqrtAWR") - call read_dataset(multipole % start_E, group_id, "start_E") - call read_dataset(multipole % end_E, group_id, "end_E") - - ! Allocate the multipole array components - call multipole % allocate() - - ! Read in arrays - call read_dataset(multipole % data, group_id, "data") - call read_dataset(multipole % pseudo_k0RS, group_id, "pseudo_K0RS") - call read_dataset(multipole % l_value, group_id, "l_value") - call read_dataset(multipole % w_start, group_id, "w_start") - call read_dataset(multipole % w_end, group_id, "w_end") - call read_dataset(multipole % broaden_poly, group_id, "broaden_poly") - - call read_dataset(multipole % curvefit, group_id, "curvefit") - - call close_group(group_id) - - ! Close file - call file_close(file_id) - - end associate - - end subroutine multipole_read - -end module multipole diff --git a/src/multipole_header.F90 b/src/multipole_header.F90 index 143047b0b5..7fb438bce1 100644 --- a/src/multipole_header.F90 +++ b/src/multipole_header.F90 @@ -1,5 +1,12 @@ module multipole_header + use hdf5 + + use constants + use dict_header, only: DictIntInt + use error, only: fatal_error + use hdf5_interface + implicit none !======================================================================== @@ -29,9 +36,6 @@ module multipole_header FIT_A = 2, & ! Absorption FIT_F = 3 ! Fission - ! Value of 'true' when checking if nuclide is fissionable - integer, parameter :: MP_FISS = 1 - !=============================================================================== ! MULTIPOLE contains all the components needed for the windowed multipole ! temperature dependent cross section libraries for the resolved resonance @@ -42,87 +46,142 @@ module multipole_header !========================================================================= ! Isotope Properties - logical :: fissionable = .false. ! Is this isotope fissionable? - integer :: length ! Number of poles - integer, allocatable :: l_value(:) ! The l index of the pole - real(8), allocatable :: pseudo_k0RS(:) ! The value (sqrt(2*mass neutron)/reduced planck constant) - ! * AWR/(AWR + 1) * scattering radius for each l - complex(8), allocatable :: data(:,:) ! Contains all of the pole-residue data - real(8) :: sqrtAWR ! Square root of the atomic weight ratio + + logical :: fissionable ! Is this isotope fissionable? + integer, allocatable :: l_value(:) ! The l index of the pole + integer :: num_l ! Number of unique l values + real(8), allocatable :: pseudo_k0RS(:) ! The value (sqrt(2*mass neutron + ! /reduced planck constant) + ! * AWR/(AWR + 1) + ! * scattering radius for + ! each l + complex(8), allocatable :: data(:,:) ! Poles and residues + real(8) :: sqrtAWR ! Square root of the atomic + ! weight ratio + integer :: formalism ! R-matrix formalism !========================================================================= ! Windows - integer :: windows ! Number of windows - integer :: fit_order ! Order of the fit. 1 linear, 2 quadratic, etc. + integer :: fit_order ! Order of the fit. 1 linear, + ! 2 quadratic, etc. real(8) :: start_E ! Start energy for the windows real(8) :: end_E ! End energy for the windows - real(8) :: spacing ! The actual spacing in sqrt(E) space. - ! spacing = sqrt(multipole_w%endE - multipole_w%startE)/multipole_w%windows - integer, allocatable :: w_start(:) ! Contains the index of the pole at the start of the window - integer, allocatable :: w_end(:) ! Contains the index of the pole at the end of the window - real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function. (reaction type, coeff index, window index) + real(8) :: spacing ! The actual spacing in sqrt(E) + ! space. + ! spacing = sqrt(multipole_w % endE - multipole_w % startE) + ! / multipole_w % windows + integer, allocatable :: w_start(:) ! Contains the index of the pole at + ! the start of the window + integer, allocatable :: w_end(:) ! Contains the index of the pole at + ! the end of the window + real(8), allocatable :: curvefit(:,:,:) ! Contains the fitting function. + ! (reaction type, coeff index, + ! window index) integer, allocatable :: broaden_poly(:) ! if 1, broaden, if 0, don't. - !========================================================================= - ! Storage Helpers - integer :: num_l - integer :: max_w + contains - integer :: formalism + procedure :: from_hdf5 => multipole_from_hdf5 - contains - procedure :: allocate => multipole_allocate ! Allocates Multipole end type MultipoleArray contains !=============================================================================== -! MULTIPOLE_ALLOCATE allocates necessary data for Multipole. +! FROM_HDF5 loads multipole data from an HDF5 file. !=============================================================================== - subroutine multipole_allocate(multipole) - class(MultipoleArray), intent(inout) :: multipole ! Multipole object to allocate. + subroutine multipole_from_hdf5(this, filename) + class(MultipoleArray), intent(inout) :: this + character(len=*), intent(in) :: filename - ! This function assumes length, numL, fissionable, windows, fitorder, - ! and formalism are known + character(len=10) :: version + integer :: i, n_poles, n_residue_types, n_windows + integer(HSIZE_T) :: dims_1d(1), dims_2d(2), dims_3d(3) + integer(HID_T) :: file_id + integer(HID_T) :: group_id + integer(HID_T) :: dset + type(DictIntInt) :: l_val_dict - ! Allocate the pole-residue storage. - ! MLBW has one more pole than Reich-Moore, and fissionable nuclides - ! have further one more. - if (multipole % formalism == FORM_MLBW) then - if (multipole % fissionable) then - allocate(multipole % data(5, multipole % length)) - else - allocate(multipole % data(4, multipole % length)) - end if - else if (multipole % formalism == FORM_RM) then - if (multipole % fissionable) then - allocate(multipole % data(4, multipole % length)) - else - allocate(multipole % data(3, multipole % length)) - end if + ! Open file for reading and move into the /isotope group + file_id = file_open(filename, 'r', parallel=.true.) + group_id = open_group(file_id, "/nuclide") + + ! Check the file version number. + call read_dataset(version, file_id, "version") + if (version /= VERSION_MULTIPOLE) call fatal_error("The current multipole& + & format version is " // trim(VERSION_MULTIPOLE) // " but the file "& + // trim(filename) // " uses version " // trim(version) // ".") + + ! Read scalar values. + call read_dataset(this % formalism, group_id, "formalism") + call read_dataset(this % spacing, group_id, "spacing") + call read_dataset(this % sqrtAWR, group_id, "sqrtAWR") + call read_dataset(this % start_E, group_id, "start_E") + call read_dataset(this % end_E, group_id, "end_E") + + ! Read the "data" array. Use its shape to figure out the number of poles + ! and residue types in this data. + dset = open_dataset(group_id, "data") + call get_shape(dset, dims_2d) + n_residue_types = int(dims_2d(1), 4) - 1 + n_poles = int(dims_2d(2), 4) + allocate(this % data(n_residue_types+1, n_poles)) + call read_dataset(this % data, dset) + call close_dataset(dset) + + ! Check to see if this data includes fission residues. + if (this % formalism == FORM_RM) then + this % fissionable = (n_residue_types == 3) + else + ! Assume FORM_MLBW. + this % fissionable = (n_residue_types == 4) + end if + + ! Read the "l_value" array. + allocate(this % l_value(n_poles)) + call read_dataset(this % l_value, group_id, "l_value") + + ! Figure out the number of unique l values in the l_value array. + do i = 1, n_poles + if (.not. l_val_dict % has(this % l_value(i))) then + call l_val_dict % set(this % l_value(i), 0) end if + end do + this % num_l = l_val_dict % size() + call l_val_dict % clear() - ! Allocate the l value table for each pole-residue set. - allocate(multipole % l_value(multipole % length)) + ! Read the "pseudo_K0RS" array. + allocate(this % pseudo_k0RS(this % num_l)) + call read_dataset(this % pseudo_k0RS, group_id, "pseudo_K0RS") - ! Allocate the table of pseudo_k0RS values at each l. - allocate(multipole % pseudo_k0RS(multipole % num_l)) + ! Read the "w_start" array and use its shape to figure out the number of + ! windows. + dset = open_dataset(group_id, "w_start") + call get_shape(dset, dims_1d) + n_windows = int(dims_1d(1), 4) + allocate(this % w_start(n_windows)) + call read_dataset(this % w_start, dset) + call close_dataset(dset) - ! Allocate window start, window end - allocate(multipole % w_start(multipole % windows)) - allocate(multipole % w_end(multipole % windows)) + ! Read the "w_end" and "broaden_poly" arrays. + allocate(this % w_end(n_windows)) + call read_dataset(this % w_end, group_id, "w_end") + allocate(this % broaden_poly(n_windows)) + call read_dataset(this % broaden_poly, group_id, "broaden_poly") - ! Allocate broaden_poly - allocate(multipole % broaden_poly(multipole % windows)) + ! Read the "curvefit" array. + dset = open_dataset(group_id, "curvefit") + call get_shape(dset, dims_3d) + allocate(this % curvefit(dims_3d(1), dims_3d(2), dims_3d(3))) + call read_dataset(this % curvefit, dset) + call close_dataset(dset) + this % fit_order = int(dims_3d(2), 4) - 1 - ! Allocate curvefit - if(multipole % fissionable) then - allocate(multipole % curvefit(FIT_F, multipole % fit_order+1, multipole % windows)) - else - allocate(multipole % curvefit(FIT_A, multipole % fit_order+1, multipole % windows)) - end if - end subroutine + ! Close the group and file. + call close_group(group_id) + call file_close(file_id) + end subroutine multipole_from_hdf5 end module multipole_header diff --git a/tests/unit_tests/test_data_multipole.py b/tests/unit_tests/test_data_multipole.py index 4a4a9f0d21..a1cc5bc02c 100644 --- a/tests/unit_tests/test_data_multipole.py +++ b/tests/unit_tests/test_data_multipole.py @@ -17,6 +17,13 @@ def u235(): return openmc.data.WindowedMultipole.from_hdf5(filename) +@pytest.fixture(scope='module') +def u234(): + directory = os.environ['OPENMC_MULTIPOLE_LIBRARY'] + filename = os.path.join(directory, '092234.h5') + return openmc.data.WindowedMultipole.from_hdf5(filename) + + @pytest.fixture(scope='module') def fe56(): directory = os.environ['OPENMC_MULTIPOLE_LIBRARY'] @@ -24,8 +31,8 @@ def fe56(): return openmc.data.WindowedMultipole.from_hdf5(filename) -def test_evaluate(u235): - """Make sure multipole object can be called.""" +def test_evaluate_rm(u235): + """Make sure a Reich-Moore multipole object can be called.""" energies = [1e-3, 1.0, 10.0, 50.] total, absorption, fission = u235(energies, 0.0) assert total[1] == pytest.approx(90.64895383) @@ -33,6 +40,15 @@ def test_evaluate(u235): assert total[1] == pytest.approx(91.12534964) +def test_evaluate_mlbw(u234): + """Make sure a Multi-Level Breit-Wigner multipole object can be called.""" + energies = [1e-3, 1.0, 10.0, 50.] + total, absorption, fission = u234(energies, 0.0) + assert total[3] == pytest.approx(15.02827953) + total, absorption, fission = u234(energies, 300.0) + assert total[3] == pytest.approx(15.08269143) + + def test_high_l(fe56): """Test a nuclide (Fe56) with a high l-value (4).""" energies = [1e-3, 1.0, 10.0, 1e3, 1e5] @@ -40,3 +56,9 @@ def test_high_l(fe56): assert total[0] == pytest.approx(25.072619556789267) total, absorption, fission = fe56(energies, 300.0) assert total[0] == pytest.approx(27.85535792368082) + + +def test_export_to_hdf5(tmpdir, u235): + filename = str(tmpdir.join('092235.h5')) + u235.export_to_hdf5(filename) + assert os.path.exists(filename)