From 73228f29242d08306ba8b552f128dc2db3b30903 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 10 Nov 2011 11:59:58 -0500 Subject: [PATCH] Reading binary ACE libraries complete. convert_xsdir.py script updated. --- src/constants.f90 | 5 + src/cross_section.f90 | 311 +++++++++--------- src/cross_section_header.f90 | 5 +- src/input_xml.f90 | 46 ++- src/utils/convert_xsdir.py | 21 +- .../templates/cross_sections_t.xml | 5 +- 6 files changed, 200 insertions(+), 193 deletions(-) diff --git a/src/constants.f90 b/src/constants.f90 index 6ae68440dc..106a38b0a4 100644 --- a/src/constants.f90 +++ b/src/constants.f90 @@ -221,6 +221,11 @@ module constants NU_POLYNOMIAL = 1, & ! Nu values given by polynomial NU_TABULAR = 2 ! Nu values given by tabular distribution + ! Cross section filetypes + integer, parameter :: & + ASCII = 1, & ! ASCII cross section file + BINARY = 2 ! Binary cross section file + ! Maximum number of partial fission reactions integer, parameter :: PARTIAL_FISSION_MAX = 4 diff --git a/src/cross_section.f90 b/src/cross_section.f90 index 1d900d0295..c9d03053e1 100644 --- a/src/cross_section.f90 +++ b/src/cross_section.f90 @@ -21,6 +21,8 @@ module cross_section real(8), allocatable :: XSS(:) ! Cross section data integer :: XSS_index ! current index in XSS data + type(DictionaryCI), pointer :: already_read => null() + private :: NXS private :: JXS private :: XSS @@ -45,10 +47,6 @@ contains 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 @@ -93,13 +91,6 @@ contains else 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 @@ -120,8 +111,8 @@ contains 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 this S(a,b) table 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 @@ -131,13 +122,6 @@ contains 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 end if end do @@ -152,19 +136,45 @@ contains ! ========================================================================== ! READ ALL ACE CROSS SECTION TABLES - ! Get list to list of nuclides - list => dict_keys(library_dict) + call dict_create(already_read) ! Loop over all files - do while (associated(list)) - ! Read all cross sections from this file - path = list % data % key - call read_ace_tables(path) + do i = 1, n_materials + mat => materials(i) - ! Move pointer to next item - list => list % next + do j = 1, mat % n_nuclides + name = mat % names(j) + + if (.not. dict_has_key(already_read, name)) then + index = dict_get_key(xs_listing_dict, name) + index_nuclides = dict_get_key(nuclide_dict, name) + name = xs_listings(index) % name + alias = xs_listings(index) % alias + + call read_ace_table(index_nuclides, index) + + call dict_add_key(already_read, name, 0) + call dict_add_key(already_read, alias, 0) + end if + end do + + if (mat % has_sab_table) then + name = mat % sab_name + + if (.not. dict_has_key(already_read, name)) then + index = dict_get_key(xs_listing_dict, name) + index_sab = dict_get_key(sab_dict, name) + + call read_ace_table(index_sab, index) + + call dict_add_key(already_read, name, 0) + end if + end if end do + ! ========================================================================== + ! ASSIGN S(A,B) TABLES TO SPECIFIC NUCLIDES WITHIN MATERIALS + do i = 1, n_materials mat => materials(i) @@ -190,163 +200,136 @@ contains end if end do - call dict_delete(library_dict) - end subroutine read_xs !=============================================================================== -! 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. +! READ_ACE_BINARY reads a single cross section table in binary format. This +! routine reads the header data for each table and then calls appropriate +! subroutines to parse the actual data. !=============================================================================== - subroutine read_ace_tables(filename) + subroutine read_ace_table(index_table, index) - character(*) :: filename ! name of ACE library file + integer, intent(in) :: index_table + integer, intent(in) :: index - 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 + integer :: i + integer :: j, j1, j2 + integer :: record_length + integer :: location + integer :: entries + integer :: length + integer :: in = 7 + integer :: zaids(16) + integer :: filetype + real(8) :: kT + real(8) :: awrs(16) + real(8) :: awr + character(10) :: name + character(10) :: date + character(10) :: mat + character(70) :: comment + character(MAX_FILE_LEN) :: filename ! path to ACE cross section library type(Nuclide), pointer :: nuc => null() type(SAB_Table), pointer :: sab => null() + type(XsListing), pointer :: listing => null() - ! 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() + ! determine path, record length, and location of table + listing => xs_listings(index) + filename = listing % path + record_length = listing % recl + location = listing % location + entries = listing % entries + filetype = listing % filetype + + ! display message + message = "Loading ACE cross section table: " // listing % name + call write_message(6) + + if (filetype == ASCII) then + ! ======================================================================= + ! READ ACE TABLE IN ASCII FORMAT + + ! Find location of table + open(UNIT=in, FILE=filename, STATUS='old', ACTION='read') + rewind(UNIT=in) + do i = 1, location - 1 + read(UNIT=in, FMT=*) + end do + + ! Read first line of header + read(UNIT=in, FMT='(A10,2E12.0,1X,A10)') name, awr, kT, date + + ! Read more header and NXS and JXS + read(UNIT=in, FMT=100) comment, mat, & + (zaids(i), awrs(i), i=1,16), NXS, JXS +100 format(A70,A10/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/4(I7,F11.0)/& + &,8I9/8I9/8I9/8I9/8I9/8I9) + + ! determine table length + length = NXS(1) + allocate(XSS(length)) + + ! Read XSS array + read(UNIT=in, FMT='(4E20.0)') XSS + + elseif (filetype == BINARY) then + ! ======================================================================= + ! READ ACE TABLE IN BINARY FORMAT + + open(UNIT=in, FILE=filename, STATUS='old', ACTION='read', & + ACCESS='direct', RECL=record_length) + + ! Read all header information + read(UNIT=in, REC=location) name, awr, kT, date, & + comment, mat, (zaids(i), awrs(i), i=1,16), NXS, JXS + + ! determine table length + length = NXS(1) + allocate(XSS(length)) + + ! Read remaining records with XSS + do i = 1, (length + entries - 1)/entries + j1 = 1 + (i-1)*entries + j2 = min(length, j1 + entries - 1) + read(UNIT=IN, REC=location + i) (XSS(j), j=j1,j2) + end do end if - ! 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 + ! ========================================================================== + ! PARSE DATA BASED ON NXS, JXS, AND XSS ARRAYS - do - ! Reset flags for finding table to read - found_nuclide = .false. - found_sab = .false. + select case(listing % type) + case (ACE_NEUTRON) + nuc => nuclides(index_table) + nuc % name = name + nuc % awr = awr + nuc % temp = kT / K_BOLTZMANN + nuc % zaid = NXS(2) - ! 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) exit - call split_string(line, words, n) + 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) + case (ACE_THERMAL) + sab => sab_tables(index_table) + sab % name = name + sab % awr = awr + sab % temp = kT / K_BOLTZMANN + sab % zaid = zaids(1) - 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)) - 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 - - ! 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_nuclide) then - ! Set ZAID of nuclide - nuc % zaid = NXS(2) - - ! allocate storage for XSS array - allocate(XSS(n)) - - ! 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_thermal_data(sab) + end select + deallocate(XSS) if(associated(nuc)) nullify(nuc) if(associated(sab)) nullify(sab) close(UNIT=in) - end subroutine read_ace_tables + end subroutine read_ace_table !=============================================================================== ! READ_ESZ - reads through the ESZ block. This block contains the energy grid, diff --git a/src/cross_section_header.f90 b/src/cross_section_header.f90 index 4698d1c2aa..abca6504d7 100644 --- a/src/cross_section_header.f90 +++ b/src/cross_section_header.f90 @@ -158,10 +158,13 @@ module cross_section_header character(10) :: alias integer :: type integer :: zaid + integer :: filetype + integer :: location + integer :: recl + integer :: entries real(8) :: awr real(8) :: temp logical :: metastable - logical :: binary character(MAX_FILE_LEN) :: path end type XsListing diff --git a/src/input_xml.f90 b/src/input_xml.f90 index ba881eb554..15a6cf7077 100644 --- a/src/input_xml.f90 +++ b/src/input_xml.f90 @@ -937,6 +937,9 @@ contains integer :: i ! loop index integer :: n_listings ! number of listings in cross_sections.xml + integer :: filetype ! default file type + integer :: recl ! default record length + integer :: entries ! default number of entries logical :: file_exists ! does cross_sections.xml exist? character(MAX_WORD_LEN) :: directory type(XsListing), pointer :: listing => null() @@ -953,8 +956,11 @@ contains message = "Reading cross sections XML file..." call write_message(5) - ! Initialize directory_ variable + ! Initialize variables that may go unused directory_ = "" + filetype_ = "" + record_length_ = 0 + entries_ = 0 ! Parse cross_sections.xml file call read_xml_file_cross_sections_t(path_cross_sections) @@ -962,6 +968,22 @@ contains ! Copy directory information if present directory = trim(directory_) + ! determine whether binary/ascii + if (filetype_ == 'ascii') then + filetype = ASCII + elseif (filetype_ == 'binary') then + filetype = BINARY + elseif (len_trim(filetype_) == 0) then + filetype = ASCII + else + message = "Unknown filetype in cross_sections.xml: " // trim(filetype_) + call fatal_error() + end if + + ! copy default record length and entries for binary files + recl = record_length_ + entries = entries_ + ! Allocate xs_listings array if (.not. associated(ace_tables_)) then message = "No ACE table listings present in cross_sections.xml file!" @@ -981,16 +1003,19 @@ contains listing % zaid = ace_tables_(i) % zaid listing % awr = ace_tables_(i) % awr listing % temp = ace_tables_(i) % temperature + listing % location = ace_tables_(i) % location ! determine type of cross section - select case(ace_tables_(i) % type) - case ('neutron') + if (ends_with(listing % name, 'c')) then listing % type = ACE_NEUTRON - case ('thermal') + elseif (ends_with(listing % name, 't')) then listing % type = ACE_THERMAL - case ('dosimetry') - listing % type = ACE_DOSIMETRY - end select + end if + + ! set filetype, record length, and number of entries + listing % filetype = filetype + listing % recl = recl + listing % entries = entries ! determine metastable state if (ace_tables_(i) % metastable == 0) then @@ -999,13 +1024,6 @@ contains listing % metastable = .true. end if - ! determine whether binary/ascii - if (ace_tables_(i) % binary == 0) then - listing % binary = .false. - else - listing % binary = .true. - end if - ! determine path of cross section table if (starts_with(ace_tables_(i) % path, '/')) then listing % path = ace_tables_(i) % path diff --git a/src/utils/convert_xsdir.py b/src/utils/convert_xsdir.py index bac86eae38..912dd2b2f2 100755 --- a/src/utils/convert_xsdir.py +++ b/src/utils/convert_xsdir.py @@ -71,8 +71,8 @@ class Xsdir(object): table.filename = words[2] table.access = words[3] table.filetype = int(words[4]) - table.address = int(words[5]) - table.tablelength = int(words[6]) + table.location = int(words[5]) + table.length = int(words[6]) self.filetype.add(table.filetype) @@ -159,8 +159,8 @@ class XsdirTable(object): self.filename = None self.access = None self.filetype = None - self.address = None - self.tablelength = None + self.location = None + self.length = None self.recordlength = None self.entries = None self.temperature = None @@ -207,8 +207,8 @@ class XsdirTable(object): def to_xml_node(self, doc): node = doc.createElement("ace_table") node.setAttribute("name", self.name) - for attribute in ["alias", "zaid", "type", "metastable", - "awr", "temperature", "binary", "path"]: + for attribute in ["alias", "zaid", "type", "metastable", "awr", + "temperature", "path", "location"]: if hasattr(self, attribute): # Join string for alias attribute if attribute == "alias": @@ -216,24 +216,19 @@ class XsdirTable(object): continue string = " ".join(self.alias) else: - string = "{0}".format(getattr(self,attribute)) + string = str(getattr(self,attribute)) # Skip metastable and binary if 0 if attribute == "metastable" and self.metastable == 0: continue - if attribute == "binary" and self.binary == 0: - continue # Skip any attribute that is none if getattr(self, attribute) is None: continue # Create attribute node - # nodeAttr = doc.createElement(attribute) - # text = doc.createTextNode(string) - # nodeAttr.appendChild(text) - # node.appendChild(nodeAttr) node.setAttribute(attribute, string) + return node diff --git a/src/xml-fortran/templates/cross_sections_t.xml b/src/xml-fortran/templates/cross_sections_t.xml index d49d798620..07bdaddb99 100644 --- a/src/xml-fortran/templates/cross_sections_t.xml +++ b/src/xml-fortran/templates/cross_sections_t.xml @@ -13,11 +13,14 @@ - + + + +