From c63ae4755b0535fc2776127c306ae16c3af53acc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Tue, 8 Nov 2011 13:12:50 -0500 Subject: [PATCH] Improved performance of cross section reading. Only one pass is made through each ACE library now. --- src/cross_section.f90 | 447 +++++++++++++++++---------------------- src/datatypes_header.f90 | 2 +- 2 files changed, 198 insertions(+), 251 deletions(-) diff --git a/src/cross_section.f90 b/src/cross_section.f90 index e5920e1e10..1d900d0295 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -3,8 +3,8 @@ module cross_section use constants use cross_section_header, only: Nuclide, Reaction, SAB_Table, XsListing use datatypes, only: dict_create, dict_add_key, dict_get_key, & - dict_has_key, dict_delete - use datatypes_header, only: DictionaryCI + dict_has_key, dict_delete, dict_keys + use datatypes_header, only: DictionaryCI, ListKeyValueCI use endf, only: reaction_name use error, only: fatal_error use fileio, only: read_line, skip_lines @@ -34,16 +34,24 @@ contains subroutine read_xs() - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - integer :: index ! index in xs_listings array - integer :: index_nuclides ! index in nuclides - integer :: index_sab ! index in sab_tables - character(10) :: key ! name of isotope, e.g. 92235.03c - type(Material), pointer :: mat => null() - type(Nuclide), pointer :: nuc => null() - type(SAB_Table), pointer :: sab => null() - type(DictionaryCI), pointer :: temp_dict => null() + integer :: i ! index in materials array + integer :: j ! index over nuclides in material + integer :: index ! index in xs_listings array + integer :: index_nuclides ! index in nuclides + integer :: index_sab ! index in sab_tables + character(10) :: name ! name of isotope, e.g. 92235.03c + character(10) :: alias ! alias of nuclide, e.g. U-235.03c + character(MAX_FILE_LEN) :: path ! path to ACE cross section library + type(Material), pointer :: mat => null() + type(Nuclide), pointer :: nuc => null() + type(SAB_Table), pointer :: sab => null() + type(DictionaryCI), pointer :: library_dict => null() + type(ListKeyValueCI), pointer :: list => null() + + call dict_create(library_dict) + + ! ========================================================================== + ! COUNT NUMBER OF TABLES AND CREATE DICTIONARIES ! First we need to determine how many continuous-energy tables and how many ! S(a,b) thermal scattering tables there are -- this loop doesn't actually @@ -58,89 +66,109 @@ contains ! materials -- if not, then increment the count for the number of ! nuclides and add to dictionary do j = 1, mat % n_nuclides - key = mat % names(j) - call lower_case(key) - if (.not. dict_has_key(nuclide_dict, key)) then - index_nuclides = index_nuclides + 1 - call dict_add_key(nuclide_dict, key, index_nuclides) + name = mat % names(j) + + ! First check that this nuclide is listed in the cross_sections.xml + ! file + if (.not. dict_has_key(xs_listing_dict, name)) then + message = "Could not find nuclide " // trim(name) // & + " in cross_sections.xml file!" + call fatal_error() + end if + + ! Find index in xs_listing and set the name and alias according to the + ! listing + index = dict_get_key(xs_listing_dict, name) + name = xs_listings(index) % name + alias = xs_listings(index) % alias + + ! If this nuclide hasn't been encountered yet, we need to add its name + ! and alias to the nuclide_dict + if (.not. dict_has_key(nuclide_dict, name)) then + index_nuclides = index_nuclides + 1 mat % nuclide(j) = index_nuclides + + call dict_add_key(nuclide_dict, name, index_nuclides) + call dict_add_key(nuclide_dict, alias, index_nuclides) else - mat % nuclide(j) = dict_get_key(nuclide_dict, key) + mat % nuclide(j) = dict_get_key(nuclide_dict, name) + end if + + ! Determine the set of ACE files which contain cross-sections used + ! in this problem + path = xs_listings(index) % path + if (.not. dict_has_key(library_dict, path)) then + call dict_add_key(library_dict, path, 0) end if end do ! Check if S(a,b) exists on other materials if (mat % has_sab_table) then - key = mat % sab_name - call lower_case(key) - if (.not. dict_has_key(sab_dict, key)) then - index_sab = index_sab + 1 - call dict_add_key(sab_dict, key, index_sab) - mat % sab_table = index_sab - else - mat % sab_table = dict_get_key(sab_dict, key) + name = mat % sab_name + + ! First check that this nuclide is listed in the cross_sections.xml + ! file + if (.not. dict_has_key(xs_listing_dict, name)) then + message = "Could not find S(a,b) table " // trim(name) // & + " in cross_sections.xml file!" + call fatal_error() + end if + + ! Find index in xs_listing and set the name and alias according to the + ! listing + index = dict_get_key(xs_listing_dict, name) + name = xs_listings(index) % name + alias = xs_listings(index) % alias + + ! If this sab hasn't been encountered yet, we need to add its name + ! and alias to the sab_dict + if (.not. dict_has_key(sab_dict, name)) then + index_sab = index_sab + 1 + mat % sab_table = index_sab + + call dict_add_key(sab_dict, name, index_sab) + call dict_add_key(sab_dict, alias, index_sab) + else + mat % sab_table = dict_get_key(sab_dict, name) + end if + + ! Determine the set of ACE files which contain cross-sections used + ! in this problem + path = xs_listings(index) % path + if (.not. dict_has_key(library_dict, path)) then + call dict_add_key(library_dict, path, 0) end if - else - mat % sab_table = 0 end if end do n_nuclides_total = index_nuclides n_sab_tables = index_sab - ! allocate arrays for ACE table storage + ! allocate arrays for ACE table storage and cross section cache allocate(nuclides(n_nuclides_total)) allocate(sab_tables(n_sab_tables)) - - ! allocate array for microscopic cross section cache allocate(micro_xs(n_nuclides_total)) - call dict_create(temp_dict) + ! ========================================================================== + ! READ ALL ACE CROSS SECTION TABLES + + ! Get list to list of nuclides + list => dict_keys(library_dict) + + ! Loop over all files + do while (associated(list)) + ! Read all cross sections from this file + path = list % data % key + call read_ace_tables(path) + + ! Move pointer to next item + list => list % next + end do - ! Now that the nuclides and sab_tables arrays have been allocated, we can - ! repeat the same loop through each table in each material and read the - ! cross sections - index_nuclides = 0 - index_sab = 0 do i = 1, n_materials mat => materials(i) - do j = 1, mat % n_nuclides - ! Get index in xs_listings array for this nuclide - index = mat % xs_listing(j) - ! Get name of nuclide - key = mat % names(j) - call lower_case(key) - - if (.not. dict_has_key(temp_dict, key)) then - index_nuclides = index_nuclides + 1 - call read_ACE_continuous(index_nuclides, index) - call dict_add_key(temp_dict, key, index_nuclides) - end if - end do - - ! Read S(a,b) table if this material has one and it hasn't been read - ! already if (mat % has_sab_table) then - ! Get name of S(a,b) table - key = mat % sab_name - - if (.not. dict_has_key(temp_dict, key)) then - ! Find the entry in xs_listings for this table - if (dict_has_key(xs_listing_dict, key)) then - index = dict_get_key(xs_listing_dict, key) - else - message = "Cannot find cross-section " // trim(key) // & - " in specified cross_sections.xml file." - call fatal_error() - end if - - ! Read the table and add entry to dictionary - index_sab = index_sab + 1 - call read_ACE_thermal(index_sab, index) - call dict_add_key(temp_dict, key, index_sab) - end if - ! 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 @@ -162,51 +190,35 @@ contains end if end do - ! delete dictionary - call dict_delete(temp_dict) + call dict_delete(library_dict) end subroutine read_xs !=============================================================================== -! READ_ACE_CONTINUOUS reads in a single ACE continuous-energy neutron cross -! section table. This routine reads the header data and then calls appropriate -! subroutines to parse the actual data. +! READ_ACE_TABLES reads in all nuclides and S(a,b) tables that are present in +! the problem from a single ACE library. This routine reads the header data for +! each table and then calls appropriate subroutines to parse the actual data. !=============================================================================== - subroutine read_ACE_continuous(index_table, index) + subroutine read_ace_tables(filename) - integer, intent(in) :: index_table ! index in nuclides array - integer, intent(in) :: index ! index in xs_listings array + character(*) :: filename ! name of ACE library file - integer :: in = 7 ! unit to read from - integer :: ioError ! error status for file access - integer :: lines ! number of lines (data - integer :: n ! number of data values - real(8) :: kT ! ACE table temperature - logical :: file_exists ! does ACE library exist? - logical :: found_xs ! did we find table in library? + integer :: index_nuclides ! index in nuclides array + integer :: index_sab ! index in sab_tables array + integer :: in = 7 ! unit to read from + integer :: ioError ! error status for file access + integer :: lines ! number of lines (data + integer :: n ! number of data values + real(8) :: kT ! ACE table temperature + logical :: file_exists ! does ACE library exist? + logical :: found_nuclide ! did we find nuclide in library? + logical :: found_sab ! did we find S(a,b) table in library? character(7) :: readable ! is ACE library readable? character(MAX_LINE_LEN) :: line ! single line to read character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line - character(MAX_FILE_LEN) :: filename ! name of ACE library file - character(10) :: tablename ! name of cross section table - type(Nuclide), pointer :: nuc => null() - - ! Check to make sure index in nuclides array and xs_listings arrays are - ! valid - if (index_table > size(nuclides)) then - message = "Index of table to read is greater than length of nuclides." - call fatal_error() - elseif (index > size(xs_listings)) then - message = "Index of xs_listing entry is greater than length of " // & - "xs_listings." - call fatal_error() - end if - - filename = xs_listings(index) % path - tablename = xs_listings(index) % name - - nuc => nuclides(index_table) + type(Nuclide), pointer :: nuc => null() + type(SAB_Table), pointer :: sab => null() ! Check if input file exists and is readable inquire(FILE=filename, EXIST=file_exists, READ=readable) @@ -219,37 +231,73 @@ contains call fatal_error() end if - ! display message - message = "Loading ACE cross section table: " // tablename - call write_message(6) - ! open file - open(file=filename, unit=in, status='old', & - & action='read', iostat=ioError) + open(FILE=filename, UNIT=in, STATUS='old', & + & ACTION='read', IOSTAT=ioError) if (ioError /= 0) then message = "Error while opening file: " // filename call fatal_error() end if - found_xs = .false. - do while (.not. found_xs) + do + ! Reset flags for finding table to read + found_nuclide = .false. + found_sab = .false. + + ! Read next line and tokenize it into strings. If we have reached the end + ! of the file, then the loop is exited call read_line(in, line, ioError) - if (ioError < 0) then - message = "Could not find ACE table " // tablename // "." - call fatal_error() - end if + if (ioError < 0) exit call split_string(line, words, n) - if (trim(words(1)) == trim(tablename)) then - found_xs = .true. + + if (dict_has_key(nuclide_dict, words(1))) then + found_nuclide = .true. + index_nuclides = dict_get_key(nuclide_dict, words(1)) + nuc => nuclides(index_nuclides) + + ! Copy header data nuc % name = words(1) nuc % awr = str_to_real(words(2)) - kT = str_to_real(words(3)) + kT = str_to_real(words(3)) nuc % temp = kT / K_BOLTZMANN + + ! display message + message = "Loading ACE cross section table: " // nuc % name + call write_message(6) + + ! Skip 5 lines + call skip_lines(in, 5, ioError) + + elseif (dict_has_key(sab_dict, words(1))) then + found_sab = .true. + index_sab = dict_get_key(sab_dict, words(1)) + sab => sab_tables(index_sab) + + ! Copy header data + sab % name = words(1) + sab % awr = str_to_real(words(2)) + kT = str_to_real(words(3)) + sab % temp = kT / K_BOLTZMANN + + ! display message + message = "Loading ACE cross section table: " // sab % name + call write_message(6) + + ! Skip 1 line + call skip_lines(in, 1, ioError) + + ! Read corresponding ZAID + call read_line(in, line, ioError) + call split_string(line, words, n) + sab % zaid = str_to_int(words(1)) + + ! Skip remaining 3 lines + call skip_lines(in, 3, ioError) + + else + call skip_lines(in, 5, ioError) end if - ! Skip 5 lines - call skip_lines(in, 5, ioError) - ! Read NXS and JXS data read(UNIT=in, FMT=*) NXS read(UNIT=in, FMT=*) JXS @@ -258,7 +306,7 @@ contains n = NXS(1) lines = (n + 3)/4 - if (found_xs) then + if (found_nuclide) then ! Set ZAID of nuclide nuc % zaid = NXS(2) @@ -267,26 +315,38 @@ contains ! Read XSS read(UNIT=in, FMT=*) XSS + + call read_esz(nuc) + call read_nu_data(nuc) + call read_reactions(nuc) + call read_angular_dist(nuc) + call read_energy_dist(nuc) + call read_unr_res(nuc) + + ! Free memory from XSS array + deallocate(XSS) + elseif (found_sab) then + ! allocate storage for XSS array + allocate(XSS(n)) + + ! Read XSS + read(UNIT=in, FMT=*) XSS + + call read_thermal_data(sab) + + ! Free memory from XSS array + deallocate(XSS) else call skip_lines(in, lines, ioError) end if - end do - call read_esz(nuc) - call read_nu_data(nuc) - call read_reactions(nuc) - call read_angular_dist(nuc) - call read_energy_dist(nuc) - call read_unr_res(nuc) - - ! Free memory from XSS array - if(allocated(XSS)) deallocate(XSS) if(associated(nuc)) nullify(nuc) + if(associated(sab)) nullify(sab) - close(unit=in) + close(UNIT=in) - end subroutine read_ACE_continuous + end subroutine read_ace_tables !=============================================================================== ! READ_ESZ - reads through the ESZ block. This block contains the energy grid, @@ -1057,119 +1117,6 @@ contains end subroutine read_unr_res -!=============================================================================== -! READ_ACE_THERMAL reads in a single ACE thermal scattering S(a,b) table. This -! routine in particular reads the header data and then calls read_thermal_data -! to parse the actual data blocks. -!=============================================================================== - - subroutine read_ACE_thermal(index_table, index) - - integer, intent(in) :: index_table ! index in sab_tables array - integer, intent(in) :: index ! index in xs_listings array - - integer :: in = 7 ! unit to read from - integer :: ioError ! error status for file access - integer :: lines ! number of lines (data - integer :: n ! number of data values - real(8) :: kT ! ACE table temperature - logical :: file_exists ! does ACE library exist? - logical :: found_xs ! did we find table in library? - character(7) :: readable ! is ACE library readable? - character(MAX_LINE_LEN) :: line ! single line to read - character(MAX_WORD_LEN) :: words(MAX_WORDS) ! words on a line - character(MAX_WORD_LEN) :: filename ! name of ACE library file - character(10) :: tablename ! name of cross section table - type(SAB_Table), pointer :: table => null() - - filename = xs_listings(index) % path - tablename = xs_listings(index) % name - - table => sab_tables(index_table) - - ! Check if input file exists and is readable - inquire(FILE=filename, EXIST=file_exists, READ=readable) - if (.not. file_exists) then - message = "ACE library '" // trim(filename) // "' does not exist!" - call fatal_error() - elseif (readable(1:3) == 'NO') then - message = "ACE library '" // trim(filename) // "' is not readable! & - &Change file permissions with chmod command." - call fatal_error() - end if - - ! display message - message = "Loading ACE cross section table: " // tablename - call write_message(6) - - ! open file - open(file=filename, unit=in, status='old', & - & action='read', iostat=ioError) - if (ioError /= 0) then - message = "Error while opening file: " // filename - call fatal_error() - end if - - found_xs = .false. - do while (.not. found_xs) - call read_line(in, line, ioError) - if (ioError < 0) then - message = "Could not find ACE table " // tablename // "." - call fatal_error() - end if - call split_string(line, words, n) - if (trim(words(1)) == trim(tablename)) then - found_xs = .true. - table%name = words(1) - table%awr = str_to_real(words(2)) - kT = str_to_real(words(3)) - table%temp = kT / K_BOLTZMANN - - ! Skip 1 line - call skip_lines(in, 1, ioError) - - ! Read corresponding ZAID - call read_line(in, line, ioError) - call split_string(line, words, n) - table % zaid = str_to_int(words(1)) - - ! Skip remaining 3 lines - call skip_lines(in, 3, ioError) - else - ! Skip 5 lines - call skip_lines(in, 5, ioError) - end if - - ! Read NXS and JXS data - read(UNIT=in, FMT=*) NXS - read(UNIT=in, FMT=*) JXS - - ! Calculate how many data points and lines in the XSS array - n = NXS(1) - lines = (n + 3)/4 - - if (found_xs) then - ! allocate storage for XSS array - allocate(XSS(n)) - - ! Read XSS - read(UNIT=in, FMT=*) XSS - else - call skip_lines(in, lines, ioError) - end if - - end do - - call read_thermal_data(table) - - ! Free memory from XSS array - if(allocated(XSS)) deallocate(XSS) - if(associated(table)) nullify(table) - - close(unit=in) - - end subroutine read_ACE_thermal - !=============================================================================== ! READ_THERMAL_DATA reads elastic and inelastic cross sections and corresponding ! secondary energy/angle distributions derived from experimental S(a,b) diff --git a/src/datatypes_header.f90 b/src/datatypes_header.f90 index 7d332cea63..491e5c42ea 100644 --- a/src/datatypes_header.f90 +++ b/src/datatypes_header.f90 @@ -17,7 +17,7 @@ module datatypes_header !=============================================================================== ! Key length for dictionary - integer, parameter :: DICT_KEY_LENGTH = 20 + integer, parameter :: DICT_KEY_LENGTH = 255 type KeyValueCI character(len=DICT_KEY_LENGTH) :: key