From 8952865c099c2d710ede4033aaf174182bd6d2bf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 22 Jul 2016 15:35:17 -0500 Subject: [PATCH] Remove n_reaction and n_product attributes in HDF5 data. Add /nuclide/reactions/ group, change reaction group names (reaction_MT). --- docs/source/io_formats/nuclear_data.rst | 9 ++++--- openmc/data/neutron.py | 26 +++++++++---------- openmc/data/reaction.py | 12 +++++---- src/multipole.F90 | 2 +- src/nuclide_header.F90 | 34 ++++++++++++++++++------- src/physics.F90 | 2 +- src/reaction_header.F90 | 22 +++++++++++++--- 7 files changed, 70 insertions(+), 37 deletions(-) diff --git a/docs/source/io_formats/nuclear_data.rst b/docs/source/io_formats/nuclear_data.rst index e7d4f2d32..ba6a54eb1 100644 --- a/docs/source/io_formats/nuclear_data.rst +++ b/docs/source/io_formats/nuclear_data.rst @@ -12,15 +12,16 @@ Incident Neutron Data **//** :Attributes: - **Z** (*int*) -- Atomic number - - **A** (*int*) -- Mass number - - **metastable** (*int*) -- Metastable state + - **A** (*int*) -- Mass number. For a natural element, A=0 is given. + - **metastable** (*int*) -- Metastable state (0=ground, 1=first + excited, etc.) - **atomic_weight_ratio** (*double*) -- Mass in units of neutron masses - **temperature** (*double*) -- Temperature in MeV - **n_reaction** (*int*) -- Number of reactions :Datasets: - **energy** (*double[]*) -- Energy points at which cross sections are tabulated -**//reaction_/** +**//reactions/reaction_/** :Attributes: - **mt** (*int*) -- ENDF MT reaction number - **label** (*char[]*) -- Name of the reaction @@ -33,7 +34,7 @@ Incident Neutron Data :Datasets: - **xs** (*double[]*) -- Cross section values tabulated against the nuclide energy grid -**//reaction_/product_/** +**//reactions/reaction_/product_/** Reaction product data is described in :ref:`product`. diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index c77ea86a1..b6a4f03b5 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -256,14 +256,14 @@ class IncidentNeutron(object): g.attrs['metastable'] = self.metastable g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio g.attrs['temperature'] = self.temperature - g.attrs['n_reaction'] = len(self.reactions) # Write energy grid g.create_dataset('energy', data=self.energy) # Write reaction data - for i, rx in enumerate(self.reactions.values()): - rx_group = g.create_group('reaction_{}'.format(i)) + rxs_group = g.create_group('reactions') + for rx in self.reactions.values(): + rx_group = rxs_group.create_group('reaction_{:03}'.format(rx.mt)) rx.to_hdf5(rx_group) # Write total nu data if available @@ -315,18 +315,16 @@ class IncidentNeutron(object): data.energy = group['energy'].value # Read reaction data - n_reaction = group.attrs['n_reaction'] + rxs_group = group['reactions'] + for name, obj in sorted(rxs_group.items()): + if name.startswith('reaction_'): + rx = Reaction.from_hdf5(obj, data.energy) + data.reactions[rx.mt] = rx - # Write reaction data - for i in range(n_reaction): - rx_group = group['reaction_{}'.format(i)] - rx = Reaction.from_hdf5(rx_group, data.energy) - data.reactions[rx.mt] = rx - - # Read total nu data if available - if 'total_nu' in rx_group: - tgroup = rx_group['total_nu'] - rx.derived_products = [Product.from_hdf5(tgroup)] + # Read total nu data if available + if rx.mt in (18, 19, 20, 21, 38) and 'total_nu' in group: + tgroup = group['total_nu'] + rx.derived_products.append(Product.from_hdf5(tgroup)) # Read unresolved resonance probability tables if 'urr' in group: diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 8fb79f109..cadab0bf7 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -368,7 +368,6 @@ class Reaction(object): group.attrs['Q_value'] = self.q_value group.attrs['threshold_idx'] = self.threshold_idx + 1 group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 - group.attrs['n_product'] = len(self.products) if self.xs is not None: group.create_dataset('xs', data=self.xs.y) for i, p in enumerate(self.products): @@ -403,13 +402,16 @@ class Reaction(object): xs = group['xs'].value rx.xs = Tabulated1D(energy, xs) + # Determine number of products + n_product = 0 + for name in group: + if name.startswith('product_'): + n_product += 1 + # Read reaction products - n_product = group.attrs['n_product'] - products = [] for i in range(n_product): pgroup = group['product_{}'.format(i)] - products.append(Product.from_hdf5(pgroup)) - rx.products = products + rx.products.append(Product.from_hdf5(pgroup)) return rx diff --git a/src/multipole.F90 b/src/multipole.F90 index 382002536..770121b9d 100644 --- a/src/multipole.F90 +++ b/src/multipole.F90 @@ -133,7 +133,7 @@ contains accumulated_fission = .true. case default ! Search through all of our secondary reactions - do j = 1, nuc % n_reaction + do j = 1, size(nuc % reactions) if (nuc % reactions(j) % MT == MT(i)) then ! Match found diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 9cff8a5d5..9f5594138 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -1,8 +1,9 @@ module nuclide_header use, intrinsic :: ISO_FORTRAN_ENV + use, intrinsic :: ISO_C_BINDING - use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5iget_name_f + use hdf5 use h5lt, only: h5ltpath_valid_f use constants @@ -82,7 +83,6 @@ module nuclide_header type(MultipoleArray), pointer :: multipole => null() ! Reactions - integer :: n_reaction ! # of reactions type(Reaction), allocatable :: reactions(:) type(DictIntInt) :: reaction_index ! map MT values to index in reactions ! array; used at tally-time @@ -182,15 +182,20 @@ module nuclide_header integer :: i integer :: Z integer :: A - integer :: n_reaction + integer :: storage_type + integer :: max_corder + integer :: n_links integer :: hdf5_err integer(HID_T) :: urr_group, nu_group integer(HID_T) :: energy_dset + integer(HID_T) :: rxs_group integer(HID_T) :: rx_group integer(HID_T) :: total_nu integer(SIZE_T) :: name_len, name_file_len + integer(HSIZE_T) :: j integer(HSIZE_T) :: dims(1) character(MAX_WORD_LEN) :: temp + type(VectorInt) :: MTs logical :: exists ! Get name of nuclide from group @@ -206,8 +211,6 @@ module nuclide_header this % zaid = 1000*Z + A + 400*this % metastable call read_attribute(this % awr, group_id, 'atomic_weight_ratio') call read_attribute(this % kT, group_id, 'temperature') - call read_attribute(n_reaction, group_id, 'n_reaction') - this % n_reaction = n_reaction ! Read energy grid energy_dset = open_dataset(group_id, 'energy') @@ -217,13 +220,26 @@ module nuclide_header call read_dataset(this % energy, energy_dset) call close_dataset(energy_dset) + ! Get MT values based on group names + rxs_group = open_group(group_id, 'reactions') + call h5gget_info_f(rxs_group, storage_type, n_links, max_corder, hdf5_err) + do j = 0, n_links - 1 + call h5lget_name_by_idx_f(rxs_group, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, & + j, temp, hdf5_err, name_len) + if (starts_with(temp, "reaction_")) then + call MTs % push_back(int(str_to_int(temp(10:12)))) + end if + end do + ! Read reactions - allocate(this % reactions(n_reaction)) + allocate(this % reactions(MTs % size())) do i = 1, size(this % reactions) - rx_group = open_group(group_id, 'reaction_' // trim(to_str(i - 1))) + rx_group = open_group(rxs_group, 'reaction_' // trim(& + zero_padded(MTs % data(i), 3))) call this % reactions(i) % from_hdf5(rx_group) call close_group(rx_group) end do + call close_group(rxs_group) ! Read unresolved resonance probability tables if present call h5ltpath_valid_f(group_id, 'urr', .true., exists, hdf5_err) @@ -513,11 +529,11 @@ module nuclide_header write(unit_,*) ' # of grid points = ' // trim(to_str(this % n_grid)) write(unit_,*) ' Fissionable = ', this % fissionable write(unit_,*) ' # of fission reactions = ' // trim(to_str(this % n_fission)) - write(unit_,*) ' # of reactions = ' // trim(to_str(this % n_reaction)) + write(unit_,*) ' # of reactions = ' // trim(to_str(size(this % reactions))) ! Information on each reaction write(unit_,*) ' Reaction Q-value COM IE' - do i = 1, this % n_reaction + do i = 1, size(this % reactions) associate (rxn => this % reactions(i)) write(unit_,'(3X,A11,1X,F8.3,3X,L1,3X,I6)') & reaction_name(rxn % MT), rxn % Q_value, rxn % scatter_in_cm, & diff --git a/src/physics.F90 b/src/physics.F90 index 09a43aa99..a9f226395 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -346,7 +346,7 @@ contains i = i + 1 ! Check to make sure inelastic scattering reaction sampled - if (i > nuc % n_reaction) then + if (i > size(nuc % reactions)) then call write_particle_restart(p) call fatal_error("Did not sample any reaction for nuclide " & &// trim(nuc % name)) diff --git a/src/reaction_header.F90 b/src/reaction_header.F90 index 9829f6e8f..185945563 100644 --- a/src/reaction_header.F90 +++ b/src/reaction_header.F90 @@ -1,11 +1,12 @@ module reaction_header - use hdf5, only: HID_T, HSIZE_T + use hdf5 + use constants, only: MAX_WORD_LEN use hdf5_interface, only: read_attribute, open_group, close_group, & open_dataset, read_dataset, close_dataset, get_shape use product_header, only: ReactionProduct - use string, only: to_str + use string, only: to_str, starts_with implicit none @@ -34,9 +35,16 @@ contains integer :: i integer :: cm integer :: n_product + integer :: storage_type + integer :: max_corder + integer :: n_links + integer :: hdf5_err integer(HID_T) :: pgroup integer(HID_T) :: xs + integer(SIZE_T) :: name_len integer(HSIZE_T) :: dims(1) + integer(HSIZE_T) :: j + character(MAX_WORD_LEN) :: name call read_attribute(this % Q_value, group_id, 'Q_value') call read_attribute(this % MT, group_id, 'mt') @@ -51,8 +59,16 @@ contains call read_dataset(this % sigma, xs) call close_dataset(xs) + ! Determine number of products + call h5gget_info_f(group_id, storage_type, n_links, max_corder, hdf5_err) + n_product = 0 + do j = 0, n_links - 1 + call h5lget_name_by_idx_f(group_id, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, & + j, name, hdf5_err, name_len) + if (starts_with(name, "product_")) n_product = n_product + 1 + end do + ! Read products - call read_attribute(n_product, group_id, 'n_product') allocate(this % products(n_product)) do i = 1, n_product pgroup = open_group(group_id, 'product_' // trim(to_str(i - 1)))