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

@@ -12,15 +12,16 @@ Incident Neutron Data
 **/<nuclide name>/**
 
 :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
 
-**/<nuclide name>/reaction_<i>/**
+**/<nuclide name>/reactions/reaction_<mt>/**
 
 :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
 
-**/<nuclide name>/reaction_<i>/product_<j>/**
+**/<nuclide name>/reactions/reaction_<mt>/product_<j>/**
 
    Reaction product data is described in :ref:`product`.
 

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:

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
 

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
 

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, &

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))

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)))