diff --git a/src/ace.F90 b/src/ace.F90 index c42e4cbecf..ed51f591f5 100644 --- a/src/ace.F90 +++ b/src/ace.F90 @@ -34,6 +34,7 @@ contains integer :: i ! index in materials array integer :: j ! index over nuclides in material + integer :: k ! index over S(a,b) tables in material integer :: i_listing ! index in xs_listings array integer :: i_nuclide ! index in nuclides integer :: i_sab ! index in sab_tables @@ -79,8 +80,9 @@ contains end if end do - if (mat % has_sab_table) then - name = mat % sab_name + do k = 1, mat % n_sab + ! Get name of S(a,b) table + name = mat % sab_names(k) if (.not. already_read % contains(name)) then i_listing = xs_listing_dict % get_key(name) @@ -93,7 +95,7 @@ contains ! Add name to dictionary call already_read % add(name) end if - end if + end do end do ! ========================================================================== @@ -102,26 +104,26 @@ contains do i = 1, n_materials mat => materials(i) - if (mat % has_sab_table) then - ! In order to know which nuclide the S(a,b) table applies to, we need - ! to search through the list of nuclides for one which has a matching - ! zaid - sab => sab_tables(mat % sab_table) + do k = 1, mat % n_sab + ! In order to know which nuclide the S(a,b) table applies to, we need to + ! search through the list of nuclides for one which has a matching zaid + sab => sab_tables(mat % i_sab_tables(k)) - do j = 1, mat % n_nuclides - nuc => nuclides(mat % nuclide(j)) - if (nuc % zaid == sab % zaid) then - mat % sab_nuclide = j + ! Loop through nuclides and find match + FIND_NUCLIDE: do j = 1, mat % n_nuclides + if (nuclides(mat % nuclide(j)) % zaid == sab % zaid) then + mat % i_sab_nuclides(k) = j + exit FIND_NUCLIDE end if - end do + end do FIND_NUCLIDE ! Check to make sure S(a,b) table matched a nuclide - if (mat % sab_nuclide == 0) then - message = "S(a,b) table " // trim(mat % sab_name) // " did not match & - &any nuclide on material " // trim(to_str(mat % id)) + if (mat % i_sab_nuclides(k) == NONE) then + message = "S(a,b) table " // trim(mat % sab_names(k)) // " did not & + &match any nuclide on material " // trim(to_str(mat % id)) call fatal_error() end if - end if + end do end do end subroutine read_xs diff --git a/src/ace_header.F90 b/src/ace_header.F90 index f5afeb171f..8506352bb1 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -195,7 +195,7 @@ module ace_header real(8) :: kappa_fission ! microscopic energy-released from fission ! Information for S(a,b) use - logical :: use_sab ! in S(a,b) energy range? + integer :: index_sab ! index in sab_tables (zero means no table) real(8) :: elastic_sab ! microscopic elastic scattering on S(a,b) table ! Information for URR probability table use diff --git a/src/cross_section.F90 b/src/cross_section.F90 index b0eac7788c..9e0b33d9f1 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -23,8 +23,9 @@ contains integer :: i ! loop index over nuclides integer :: i_nuclide ! index into nuclides array integer :: i_sab ! index into sab_tables array + integer :: j ! index in mat % i_sab_nuclides real(8) :: atom_density ! atom density of a nuclide - real(8) :: sab_threshold ! threshold for S(a,b) table + logical :: check_sab ! should we check for S(a,b) table? type(Material), pointer :: mat => null() ! current material ! Set all material macroscopic cross sections to zero @@ -43,31 +44,50 @@ contains ! Find energy index on unionized grid if (grid_method == GRID_UNION) call find_energy_index() - ! Check if there's an S(a,b) table for this material - if (mat % has_sab_table) then - sab_threshold = sab_tables(mat % sab_table) % threshold_inelastic - else - sab_threshold = ZERO - end if + ! Determine if this material has S(a,b) tables + check_sab = (mat % n_sab > 0) + + ! Initialize position in i_sab_nuclides + j = 1 ! Add contribution from each nuclide in material do i = 1, mat % n_nuclides + ! ======================================================================== + ! CHECK FOR S(A,B) TABLE + + i_sab = 0 + + ! Check if this nuclide matches one of the S(a,b) tables specified -- this + ! relies on i_sab_nuclides being in sorted order + if (check_sab .and. i == mat % i_sab_nuclides(j)) then + ! Get index in sab_tables + i_sab = mat % i_sab_tables(j) + + ! If particle energy is greater than the highest energy for the S(a,b) + ! table, don't use the S(a,b) table + if (p % E > sab_tables(i_sab) % threshold_inelastic) i_sab = 0 + + ! Increment position in i_sab_nuclides + j = j + 1 + + ! Don't check for S(a,b) tables if there are no more left + if (j > mat % n_sab) check_sab = .false. + end if + + ! ======================================================================== + ! CALCULATE MICROSCOPIC CROSS SECTION + ! Determine microscopic cross sections for this nuclide i_nuclide = mat % nuclide(i) - ! Determine whether to use S(a,b) based on energy of particle and whether - ! the nuclide matches the S(a,b) table - if (p % E < sab_threshold .and. i == mat % sab_nuclide) then - i_sab = mat % sab_table - else - i_sab = 0 - end if - ! Calculate microscopic cross section for this nuclide if (p % E /= micro_xs(i_nuclide) % last_E) then call calculate_nuclide_xs(i_nuclide, i_sab) end if + ! ======================================================================== + ! ADD TO MACROSCOPIC CROSS SECTION + ! Copy atom density of nuclide in material atom_density = mat % atom_density(i) @@ -148,7 +168,7 @@ contains micro_xs(i_nuclide) % interp_factor = f ! Initialize sab treatment to false - micro_xs(i_nuclide) % use_sab = .false. + micro_xs(i_nuclide) % index_sab = NONE micro_xs(i_nuclide) % elastic_sab = ZERO micro_xs(i_nuclide) % use_ptable = .false. @@ -231,7 +251,7 @@ contains type(SAlphaBeta), pointer :: sab => null() ! Set flag that S(a,b) treatment should be used for scattering - micro_xs(i_nuclide) % use_sab = .true. + micro_xs(i_nuclide) % index_sab = i_sab ! Get pointer to S(a,b) table sab => sab_tables(i_sab) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 9237d6bac1..a10e76d27f 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -464,9 +464,11 @@ contains dims, m % atom_density, hdf5_err) ! Write S(a,b) information if present - if (m % has_sab_table) then - call h5ltmake_dataset_string_f(temp_group, "sab_table", & - m % sab_name, hdf5_err) + if (m % n_sab > 0) then + call hdf5_write_integer(temp_group, "i_sab_nuclides", & + m % i_sab_nuclides) + call hdf5_write_integer(temp_group, "i_sab_tables", & + m % i_sab_tables) end if ! Close group for i-th material diff --git a/src/input_xml.F90 b/src/input_xml.F90 index ac3dcb2380..b03f8ecb2b 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -952,6 +952,7 @@ contains integer :: i ! loop index for materials integer :: j ! loop index for nuclides integer :: n ! number of nuclides + integer :: n_sab ! number of sab tables for a material integer :: index_list ! index in xs_listings array integer :: index_nuclide ! index in nuclides integer :: index_sab ! index in sab_tables @@ -1152,43 +1153,52 @@ contains ! ======================================================================= ! READ AND PARSE TAG FOR S(a,b) DATA - if (size(material_(i) % sab) == 1) then - ! Get pointer to S(a,b) table - sab => material_(i) % sab(1) + n_sab = size(material_(i) % sab) + if (n_sab > 0) then + ! Set number of S(a,b) tables + mat % n_sab = n_sab - ! Determine name of S(a,b) table - name = trim(sab % name) // "." // trim(sab % xs) - mat % sab_name = name + ! Allocate names and indices for nuclides and tables + allocate(mat % sab_names(n_sab)) + allocate(mat % i_sab_nuclides(n_sab)) + allocate(mat % i_sab_tables(n_sab)) - ! Check that this nuclide is listed in the cross_sections.xml file - if (.not. xs_listing_dict % has_key(name)) then - message = "Could not find S(a,b) table " // trim(name) // & - " in cross_sections.xml file!" - call fatal_error() - end if - mat % has_sab_table = .true. + ! Initialize i_sab_nuclides + mat % i_sab_nuclides = NONE - ! Find index in xs_listing and set the name and alias according to the - ! listing - index_list = xs_listing_dict % get_key(name) - name = xs_listings(index_list) % name - alias = xs_listings(index_list) % alias + do j = 1, n_sab + ! Get pointer to S(a,b) table + sab => material_(i) % sab(j) - ! If this S(a,b) table hasn't been encountered yet, we need to add its - ! name and alias to the sab_dict - if (.not. sab_dict % has_key(name)) then - index_sab = index_sab + 1 - mat % sab_table = index_sab + ! Determine name of S(a,b) table + name = trim(sab % name) // "." // trim(sab % xs) + mat % sab_names(j) = name - call sab_dict % add_key(name, index_sab) - call sab_dict % add_key(alias, index_sab) - else - mat % sab_table = sab_dict % get_key(name) - end if + ! Check that this nuclide is listed in the cross_sections.xml file + if (.not. xs_listing_dict % has_key(name)) then + message = "Could not find S(a,b) table " // trim(name) // & + " in cross_sections.xml file!" + call fatal_error() + end if - elseif (size(material_(i) % sab) > 1) then - message = "Cannot have multiple S(a,b) tables on a single material." - call fatal_error() + ! Find index in xs_listing and set the name and alias according to the + ! listing + index_list = xs_listing_dict % get_key(name) + name = xs_listings(index_list) % name + alias = xs_listings(index_list) % alias + + ! If this S(a,b) table hasn't been encountered yet, we need to add its + ! name and alias to the sab_dict + if (.not. sab_dict % has_key(name)) then + index_sab = index_sab + 1 + mat % i_sab_tables(j) = index_sab + + call sab_dict % add_key(name, index_sab) + call sab_dict % add_key(alias, index_sab) + else + mat % i_sab_tables(j) = sab_dict % get_key(name) + end if + end do end if ! Add material to dictionary diff --git a/src/material_header.F90 b/src/material_header.F90 index 3b90342669..cbfd917498 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -9,16 +9,18 @@ module material_header type Material integer :: id ! unique identifier integer :: n_nuclides ! number of nuclides - character(12), allocatable :: names(:) ! isotope names integer, allocatable :: nuclide(:) ! index in nuclides array real(8) :: density ! total atom density in atom/b-cm real(8), allocatable :: atom_density(:) ! nuclide atom density in atom/b-cm ! S(a,b) data references - logical :: has_sab_table = .false. - character(12) :: sab_name ! name of S(a,b) table - integer :: sab_table = 0 ! index in sab_tables - integer :: sab_nuclide = 0 ! index of nuclide which has S(a,b) table + integer :: n_sab = 0 ! number of S(a,b) tables + integer, allocatable :: i_sab_nuclides(:) ! index of corresponding nuclide + integer, allocatable :: i_sab_tables(:) ! index in sab_tables + + ! Temporary names read during initialization + character(12), allocatable :: names(:) ! isotope names + character(12), allocatable :: sab_names(:) ! name of S(a,b) table end type Material end module material_header diff --git a/src/output.F90 b/src/output.F90 index dfd1a6c9c0..b31d19edcd 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -624,8 +624,12 @@ contains end do ! Write information on S(a,b) table - if (mat % has_sab_table) then - write(unit_,*) ' S(a,b) table = ' // trim(mat % sab_name) + if (mat % n_sab > 0) then + write(unit_,*) ' S(a,b) tables:' + do i = 1, mat % n_sab + write(unit_,*) ' ' // trim(& + sab_tables(mat % i_sab_tables(i)) % name) + end do end if write(unit_,*) diff --git a/src/physics.F90 b/src/physics.F90 index 8c77620902..0602b6bd27 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -422,9 +422,9 @@ contains ! ======================================================================= ! ELASTIC SCATTERING - if (micro_xs(i_nuclide) % use_sab) then + if (micro_xs(i_nuclide) % index_sab /= NONE) then ! S(a,b) scattering - call sab_scatter(i_nuclide, mat % sab_table) + call sab_scatter(i_nuclide, micro_xs(i_nuclide) % index_sab) else ! get pointer to elastic scattering reaction