mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Merge pull request #983 from smharper/wmp_cleanup
Clean up the multipole file format
This commit is contained in:
commit
507c155560
7 changed files with 246 additions and 247 deletions
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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<i>**
|
||||
- **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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue